Mineral Inversion

Summary

Mineral inversion finds the volumes of all minerals and fluids at once by fitting the logs with a forward model, with the volumes constrained to add to one. It is a method of the Porosity step. It takes Clay volume from Clay Volume and Kerogen volume from TOC Analysis as inputs and does not compute them. Use it where the mineralogy varies, the lithology is complex, or porosity depends on a matrix that a single density cannot describe.

Inputs and outputs

Equations

Forward model. Each log responds to the sum of the contributions of the components in the rock. For log \(i\) and component \(j\), with \(A_{ij}\) the response of the pure component (its endmember) and \(V_j\) its volume per unit bulk volume:

The volumes must add to one (closure) and cannot be negative:

Logs that mix linearly in volume are used directly: bulk density, neutron porosity (approximately), sonic slowness (as in the time-average relation) and the volumetric cross section \(U\). The photoelectric factor PE does not mix linearly with volume and must be converted to \(U\) first.

Constraints enter as extra equations. The clay volume \(\Vcl\) from Clay Volume is a row that has 1 under the clay component and 0 elsewhere, with data value \(\Vcl\). The kerogen volume \(\Vk\) from TOC Analysis is a row with 1 under kerogen, with data value \(\Vk\). Closure is a row of ones with data value 1 and a very small uncertainty. These rows are weighted exactly like log rows.

Estimate. The volumes are the solution of a weighted, constrained least-squares problem, with \(\sigma_i\) the uncertainty assigned to row \(i\):

Without the sign constraint, the solution is \(\hat{\mathbf V} = (\mathbf A^T \mathbf W \mathbf A)^{-1}\mathbf A^T \mathbf W\,\mathbf d\) with \(\mathbf W = \mathrm{diag}(1/\sigma_i^2)\).

Outputs. Porosity is the sum of the fluid volumes. The grain density is the volume-weighted density of the solids:

The calculator below is the smallest useful case: quartz, calcite and a fluid, solved from the bulk density and \(U\), plus closure. With three unknowns and three equations the solution is exact. Cramer's rule gives, with \(D\) the determinant of the system:

The three volumes add to one by construction. No limit is applied, so a negative value shows that the measurement is outside what the model can describe. If \(D\) is zero the minerals cannot be separated and the calculator returns no value.

Single-value calculator

Behavior

The plot holds \(U\) at 5.72 barns/cm³ and sweeps the bulk density. Porosity falls as density rises, by about 0.06 per 0.1 g/cm³ (0.281 at 2.2 and 0.035 at 2.6 g/cm³), while quartz rises from 0.479 to 0.846 and calcite falls from 0.240 to 0.119. Density mostly sets the porosity. \(U\) does the opposite job: at a bulk density of 2.332, raising \(U\) from 4 to 8 to 10 barns/cm³ moves porosity only from 0.193 to 0.209 to 0.218, but moves calcite from 0.005 to 0.458 to 0.684 and quartz from 0.802 to 0.333 to 0.098. In this system the density log carries the porosity and the photoelectric log carries the mineral split, which is why quartz-calcite separation fails without a lithology-sensitive log. The reason is in the conditioning, covered on the uncertainty page.

Parameter guidance

Choose the components. The mineral list is the main decision. Include only components that the data can support and the geology allows. The number of unknowns must not exceed the number of independent equations: logs plus constraint rows plus closure. Every extra component with a similar response to another one makes the problem worse conditioned. A typical shale-gas set is quartz, calcite, dolomite, clay, kerogen, pyrite and water.

Supply the constraints from the upstream steps. The clay volume (Clay volume, also called VWCL) comes from Clay Volume and the kerogen volume (Kerogen volume) from TOC Analysis. The inversion uses them as data with an uncertainty. It does not recompute them, and if either is wrong the inversion inherits the error. Pyrite is sometimes tied to kerogen volume by a fixed ratio. That is a regional empirical rule and should be validated.

Endmembers are on the endmember page, and the uncertainties that weight each equation are on the uncertainty page. Give the closure row a very small uncertainty (for example 0.001) so that it acts as a hard constraint, and the constraint rows an uncertainty that reflects the real accuracy of the upstream result (for example 0.03 for clay volume and 0.01 for kerogen volume).

Choose a porosity output. Porosity can be taken as the sum of the fluid volumes (the strictest, depends on the whole solution), from the density log with the inverted grain density (the usual choice, because the density log is the most reliable porosity measurement), or averaged with the neutron porosity. Compare them.

Worked example

Five components (quartz, calcite, clay, kerogen, water) from four logs (density, neutron, sonic and U), a clay volume and a kerogen volume, plus closure. That is seven equations for five unknowns. The data are made from a known rock and then given log errors and a clay-volume error. The solver enumerates the sets of non-zero components and keeps the best non-negative solution. The output shows porosity (water) within 0.02 of the truth while the quartz-calcite split is off by more than 0.1, a result explained on the uncertainty page:

import itertools
import numpy as np

# Endmember responses. Illustrative mid-range values, not a recommendation.
#             density  neutron  sonic   U
#             g/cm3    v/v      us/ft   barns/cm3
END = {
    "quartz":  [2.65, -0.02,  55.5,  4.8],
    "calcite": [2.71,  0.00,  47.5, 13.8],
    "clay":    [2.55,  0.35, 110.0,  9.0],
    "kerogen": [1.26,  0.60, 160.0,  0.25],
    "water":   [1.00,  1.00, 189.0,  0.4],
}
LOGS = ["rhob", "nphi", "dt", "u"]
SIGMA = {"rhob": 0.03, "nphi": 0.03, "dt": 5.0, "u": 1.0, "vcl": 0.03, "vker": 0.01, "unity": 0.001}


def system(minerals, curves=LOGS, constraints=("vcl", "vker")):
    """Rows: the log curves, then the constraint rows, then the closure row."""
    rows, sig, labels = [], [], []
    for c in curves:
        rows.append([END[m][LOGS.index(c)] for m in minerals])
        sig.append(SIGMA[c]); labels.append(c)
    for c, name in (("vcl", "clay"), ("vker", "kerogen")):
        if c in constraints:
            rows.append([1.0 if m == name else 0.0 for m in minerals])
            sig.append(SIGMA[c]); labels.append(c)
    rows.append([1.0] * len(minerals)); sig.append(SIGMA["unity"]); labels.append("unity")
    return np.array(rows), np.array(sig), labels


def solve(A, sig, d):
    """Weighted least squares with V >= 0: try every active set, keep the best feasible one."""
    Aw, dw = A / sig[:, None], d / sig
    n = A.shape[1]
    best, best_cost = None, np.inf
    for k in range(1, n + 1):
        for subset in itertools.combinations(range(n), k):
            x = np.linalg.lstsq(Aw[:, subset], dw, rcond=None)[0]
            if (x < -1e-12).any():
                continue
            v = np.zeros(n)
            v[list(subset)] = x
            cost = np.sum((Aw @ v - dw) ** 2)
            if cost < best_cost:
                best, best_cost = v, cost
    return best, best_cost


truth = {"quartz": 0.45, "calcite": 0.20, "clay": 0.15, "kerogen": 0.05, "water": 0.15}
minerals = list(truth)
A, sig, labels = system(minerals)
v_true = np.array([truth[m] for m in minerals])
d = A @ v_true                       # noise-free data: four logs, VWCL, kerogen volume, closure
noise = np.array([0.02, -0.02, 4.0, -0.8, 0.02, 0.0, 0.0])   # errors on the logs and on the clay volume
d_obs = d + noise
v, cost = solve(A, sig, d_obs)
print("component   true  inverted")
for m, t, x in zip(minerals, v_true, v):
    print(f"{m:10s} {t:5.2f} {x:8.3f}")
print(f"sum of volumes = {v.sum():.4f}")
print()
pred = A @ v
print("row      measured  predicted  residual/sigma")
for lab, m_, p_, s_ in zip(labels, d_obs, pred, sig):
    print(f"{lab:7s} {m_:9.3f} {p_:10.3f} {(m_ - p_) / s_:10.2f}")
print(f"weighted cost = {cost:.3f}")

Output

component   true  inverted
quartz      0.45    0.564
calcite     0.20    0.078
clay        0.15    0.174
kerogen     0.05    0.050
water       0.15    0.134
sum of volumes = 1.0000

row      measured  predicted  residual/sigma
rhob        2.350      2.346       0.14
nphi        0.203      0.214      -0.35
dt         91.325     87.531       0.76
u           5.543      5.413       0.13
vcl         0.170      0.174      -0.12
vker        0.050      0.050      -0.01
unity       1.000      1.000      -0.02
weighted cost = 0.751

Assumptions and limitations

• Every log responds linearly to the component volumes. Density and U do. Neutron and sonic are approximately linear, but the neutron porosity depends on lithology scale and hydrogen index, and the sonic on compaction and texture.
• The component list is complete. A mineral that is present and not in the model is absorbed by the others, often without a large misfit (see the model-error QC page).
• The endmembers are right. Clay and kerogen endmembers in particular vary with type and maturity.
• Log errors are independent and have the assigned standard deviations. Washouts and borehole effects are not independent random errors.
• The clay volume and the kerogen volume passed in are accurate to their assigned uncertainty.
• Pore fluids are treated as components with a fixed response, so the invaded-zone fluid, not the virgin fluid, is what the logs see.

QC checks

• Volumes sum to one, are non-negative and lie in a physical range. A large number of components at exactly zero shows the active set changing from sample to sample, which makes curves noisy.
• Predicted logs overlay the measured logs within the assigned uncertainty. See the model-error QC page.
• The inverted clay and kerogen volumes stay close to the values supplied. A large difference means the logs disagree with the upstream result.
• Porosity compares with core porosity and with the neutron-density porosity in a clean zone. Grain density compares with core grain density.
• Volumes are stable when an endmember is changed by a small amount. If they are not, the system is poorly conditioned.

Going Deeper

Multi-mineral log inversion, solving for volumes by least squares with several logs at once, came into routine use in the 1980s as computers became available to service companies and operators, replacing sequences of crossplots. The unconventional-resource era added the need to include kerogen and pyrite as components, and pushed the field toward constraining the problem with outside information (clay volume from a gamma ray, kerogen volume from TOC) instead of leaving everything to the logs. That makes inversion a downstream user of those steps and is the reason it sits within Porosity and not in front of it. The simple solver here (an exhaustive search over sets of non-zero components) is chosen for clarity and works for a handful of components; production software uses active-set or interior-point methods.

References

1. Quirein, J., Kimminau, S., LaVigne, J., Singer, J. and Wendel, F., 1986. A coherent framework for developing and applying multiple formation evaluation models. Transactions of the SPWLA 27th Annual Logging Symposium, Paper DD.
2. Doveton, J.H., 1994. Geologic Log Analysis Using Computer Methods. AAPG Computer Applications in Geology No. 2, American Association of Petroleum Geologists, Tulsa, OK.
3. Lawson, C.L. and Hanson, R.J., 1974. Solving Least Squares Problems. Prentice-Hall, Englewood Cliffs, NJ (reprinted by SIAM, 1995).

Python reference implementation

Related methods

• Inversion Endmember Parameters
• Inversion Uncertainty and Confidence
• Inversion Model-Error QC
• Density Porosity
• Gas and Organic-Matter Corrections