Invert geology

In previous posts we have only considered the rheological properties of lithologies as inversion variables - that is, the geometry of the hull that encloses these lithologies are fixed. But our knowledge of subsurface geology significantly diminishes with depth. It is intuitive, then, that the variation of these hulls be considered within our inverse framework. Inversion strategies that operate in the forward mode, like Monte Carlo techniques, can be flexible in the way geological variation is coded. However, using gradient-based approaches limit us to parameterisations that are continuously differentiable and in which the gradient vectors are smooth.

Problem statement: We need a way to relate the spatial configuration of nodes to lithology that is differentiable.

An obvious way forward, without much additional work, is to treat each node in the entire mesh as separate inversion variables. The priors for these mesh variables will require regularisation using a prior covariance matrix to reduce spurious values from arising. A Gaussian covariance matrix is common in this situation, where the value of $\sigma$ controls the spatial extent to which adjacent nodes are correlated. However, the result will be a smooth field that is not geologically feasible: lithologies are separated by sharp unconformities such as faults or intrusions; they rarely transition gradually from one rock type to another.

We propose a solution where the vertical thickness of each layer is posed as an inversion variable. This is described by a ratio of lithologies that comprise each column in the domain. The number of inversion variables thus becomes the product of the number of lithologies and the number of columns in the domain. We demonstrate an example where we relate the geometry of heat-producing bodies to the vertically integrated heat flux.

Vertically integrated heat flux

We know that in 1D surface heat flow, $q_s$, is the sum of heat flux from the base $q_m$ and integrated heat sources, $H$:

$$ q_s = \int_z H , \mathrm{d}z + q_m $$

which is approximately true for 2D and 3D simulations, where a smaller fraction of heat flux is attributed to a lateral component. Consider a domain partitioned into two lithologies, $A$ and $B$, that are stacked vertically. We can expand the above expression to:

$$ \begin{align} q_s =& q_A + q_B + q_m \\ =& H_A \Delta z_A + H_B \Delta z_B + q_m \end{align} $$

where $\Delta z_A$ and $\Delta z_B$ are the thicknesses of lithology $A$ and $B$, respectively. Here is a python snippet of the forward model:

def forward_model(x):
    r = x
    
    # normalise ratios
    rN = r/sum(r)
    
    nc = rN*N
    Z = nc*hz
    qs = sum(H*Z)
    c = sum((qs - qobs)**2/sigma_qobs**2)
    return c

Next we derive the tangent linear model. The derivatives of $q_s$ with respect to $\Delta z_A$ and $\Delta z_B$ are:

$$ \begin{align} \frac{\partial q_s}{\partial \Delta z_A} &= H_A \\ \frac{\partial q_s}{\partial \Delta z_B} &= H_B \end{align} $$

def tangent_linear(x, dx):
    r = x
    dr = dx
    
    # normalise ratios
    rN = r/sum(r)
    drNdr = 1.0/sum(r)
    drNdrsum = -r/sum(r)**2
    drN = drNdr*dr + drNdrsum*sum(dr)
    
    nc = rN*N
    dncdrN = N
    dnc = dncdrN*drN
    
    Z = nc*hz
    dZdnc = hz
    dZ = dZdnc*dnc
    
    qs = sum(H*Z)
    dqsdZ = H
    dqs = sum(dqsdZ*dZ)
    
    c = sum((qs - qobs)**2/sigma_qobs**2)
    dcdqs = (2.0*qs - 2.0*qobs)/sigma_qobs**2
    dc = dcdqs*dqs
    return c, dc

The adjoint model simply backward-propagates the derivatives backwards through the forward model:

def adjoint_model(x, df=1.0):
    r = x
    
    # normalise ratios
    rN = r/sum(r)
    
    nc = rN*N
    Z = nc*hz
    qs = sum(H*Z)
    c = sum((qs - qobs)**2/sigma_qobs**2)
    
    # ADJOINT PART
    dcdqs = (2.0*qs - 2.0*qobs)/sigma_qobs**2
    dqs = dcdqs*df
    
    dqsdZ = H
    dZ = dqsdZ*dqs
    
    dZdnc = hz
    dnc = dZdnc*dZ
    
    dncdrN = N
    drN = dncdrN*dnc
    
    drNdr = 1.0/sum(r)
    drNdrsum = -r/sum(r)**2
    dr = drNdr*drN
    drsum = drNdrsum*drN
    
    drsumdr = 1.0/N
    dr += sum(drsum)
    
    return c, dr

There you have it! This is an adjoint to the geometry of heat-producing bodies constrained by integrated vertical heat flux. Download the Python code which tests the validity of the adjoint with the tangent linear and finite-differences of the forward model.