Old dogs and new tricks volume 3 - simultaneous log inversion in Excel, plotted in Python

Conventional ("Deterministic") log analysis

Conventional log analysis uses a sequential process to build up a petrophysical model of the formation. For example, we can use the GR to estimate shale volume, then use the density log with the shale volume to compute porosity. After that, we use the shale volume, porosity and the resistivity log to compute water saturation. In many clastic reservoirs this will give acceptable results, as in the North Sea gas well example I showed in a previous LI post.

Simultaneous inversion ("Stochastic") log analysis

An alternative method is to write down the response equations for each of the logs as a matrix, and then invert the matrix to evaluate shale volume, porosity and water saturation from all the logs simultaneously. The response equations look like this:

GR = GRshale*Vshale + GRsand*Vsand + GRwater*Vwater + GRgas*Vgas

where Vshale, Vsand, Vwater and Vgas are the volume fractions of the rock and fluid components of the formation. Similar linear equations can be formed for the density and neutron logs. Other logs, like resistivity and sonic, will be non-linear. The matrix can be inverted to compute the volume fractions, but some advanced mathematics is needed. The inversion requires minimizing an error function with the following form:

error = (GR - GRestimate)^2 * GRweight / GRrange

The error is summed for all the logs. Here, GRestimate is the GR value that is obtained by entering the assumed volume fractions into the response equation. GRweight is a weighting factor that can be used to control the effect of the GR error on the overall error, and GRrange is the difference between the maximum and minimum GR values. The error is zero when the estimate is the same as the actual log, but most of the time there won't be a solution that gives zero error for all the logs.

This approach, which is the basis for commercial programmes like Elan, MinSolv and Multimin, has been around since the early 80's. It has many advantages over the sequential method, but is sometimes viewed with suspicion, as a "black box". The simultaneous method can be applied to any reservoir and on any log data set, but it is particularly valuable in complex mineralogies, like carbonates and shale reservoirs.

Simultaneous inversion for Old Dogs

How can we use the simultaneous inversion method if we don't have access to high-end petrophysical software? I'm sure that it can be done in Python, but I found that I can use the Excel Solver function to do it instead. I first set up a table of matrix components and log responses, looking like this:

I then took one depth interval and the logs, and assigned a starting value for all the formation components:

The first four columns are the logs - neutron, GR, deep resistivity and density. The other columns are the formation rock and fluid components - illite, quartz, water and gas. I also have clay-bound water, VBWA, which is a fixed fraction of the illite volume. Since the density, neutron and GR measure the invaded zone, while the deep resistivity sees the virgin zone, I need to have water and gas fractions for both. VUGA and VUWA are the gas and water in the virgin zone, VXGA and VXWA are for the invaded zone.

At this point we can control the process by imposing some constraints. The volume fractions must sum to 1, and the porosity is the same for both the virgin and invaded zones, so we have:

VILL + VQUA + VBWA + VUGA + VUWA = 1

VUGA + VUWA = VXGA + VXWA

Also, we can't have more gas in the invaded zone than in the virgin zone, so we can use:

VUGA >= VXGA

Now we can reconstruct the logs from the components. That's easy for the GR, density and neutron as they have linear response equations. It's more complicated for the resistivity as the response equation is not linear. In this example I used the Dual Water equation, which uses Vcl as well as porosity and resistivity, so I have to first compute starting estimates of Vcl and porosity from my components:

Vcl = VILL + VBWA

phit = VBWA + VUGA + VUWA

I also need Sw, which is:

Swt = (VBWA + VUWA)/(VBWA + VUGA + VUWA)

Having reconstructed the logs, we can compute the error for each log as above and sum them. I also add a term that forces the components to sum to 1, and the porosity to match between the invaded and virgin zones. With this initial guess for the volumes, I get these reconstructed values:

You can see that the reconstructed logs (NEU_R etc) aren't very close to the input logs. So now I run the Solver with this dialog box:

In this box, the objective cell is the summed reconstruction errors, the variable cells are the rock and fluid components, and the constraints prevent the solver finding negative values. After solving, the reconstructed logs look like this:

Now the reconstructed logs are quite close to the inputs, and the error is very small. The formation components can be combined to give log analysis outputs:

Now we have clay volume, porosity and saturation. At this level the formation is reasonably clean with fair porosity, and is gas-bearing. There is significant invasion.

The final step is to write a Visual Basic macro that moves the data around and runs the Solver. Mine looks like this:

Note that I turn off screen updating at the start. That's because this is a very slow process! 400 depth intervals take about 13 minutes, way too long to be practical, but this is intended as a demonstration. So how does it look? Here's a display showing the input and reconstructed logs, with the volumetric composition:

The black curves are the input logs and the red are reconstructed. I get a pretty good reconstruction in the sand, but not so good in the shale. The volumetric composition, porosity and saturations are very similar to the conventional log analysis. Another way to QC the results is with cross plots:

The plots show the input log on the x-axis and the reconstructed log on the y-axis. You can see that the density reconstructs well, the GR and neutron need work, and the resistivity may have some depth issues! These plots show me where I need to modify my parameters to get a more satisfactory result.

What's next?

I would like to try this method in carbonates and unconventional reservoirs, but I can't find any data. If you have an example, please send it to me and I'll see what comes out.

Andy Brickell

Petrogeek Consulting LLC