Forward modelling#

A fundamental aspect of any quantitative science discipline is to make observations. These observations are generated by some physical (or chemical, biological, …) process. A forward model is a mathematical description of this process, which is often expressed as a differential equation; while the observations that we make can look extremely complex (e.g. planetary objects orbiting a star), the dynamics underlying this behaviour are often much simpler to describe (Newton’s laws of motion). We can then generate synthetic observations by solving the dynamics forward in time, given some set of parameters and initial conditions.

Quantitatively, the dynamics of a system with multiple variables or dimensions can be written as:

(1)#\[\frac{\mathrm{d} \vec{X}}{\mathrm{d} t} = f\left( t, \vec{X}, \dots \right)\]

The synthetic observations are then generated by integrating Eq. (1) over time \(t \in \left[ t_0, t \right)\):

(2)#\[\vec{X}(t) = \vec{X}(t=t_0) + \int_{t_0}^{t} f\left(t', \vec{X}, \dots \right) \mathrm{d} t'\]

There exists a plethora of strategies to perform the integration as indicated above in a stable and accurate manner: Runge-Kutta, Rosenbrock, LSODA, … DiaBayes uses the Diffrax library[1], which implements several of these solvers in JAX. A good general-purpose solver is the Tsit5 algorithm, which is the default for DiaBayes.

Friction models#

A typical forward model describing a rock friction experiment comprises 3 components: a friction law that describes the slip rate of the fault in response to a given state of stress (or “friction”), a state evolution law that controls the time- and history-dependent internal state of the fault, and a model that describes how the externally imposed stress is transmitted to the fault (usually through elasticity). To make this more concrete, let’s consider a classical spring-block model governed by rate-and-state friction[2]:

(3)#\[\begin{split}\frac{\mathrm{d} \vec{X}}{\mathrm{d} t} = \begin{cases} \dot{\mu} = k \left(v_{lp} - v(\mu, \theta) \right) \quad \text{(spring-block stress transfer)} \\ \dot{\theta} = 1 - \frac{v \theta}{D_c} \quad \text{(Ageing law state evolution)} \end{cases}\end{split}\]

with

\[v(\mu, \theta) = v_0 \exp \left( \frac{1}{a} \left[ \mu - \mu_0 - b \ln \left( \frac{v_0 \theta}{D_c} \right) \right] \right)\]

In the above expressions, \(a\), \(b\), and \(D_c\) are frictional parameters that are often of interest, and can be inverted for from measured friction curves (\(\mu(t)\)). The friction parameters \(\mu_0\) and \(v_0\) are typically taken as constants as they can be inferred directly from the data. The spring-block stiffness \(k\) and the load-point velocity \(v_{lp}\) are determined by the experimental set-up. The forward model is completed by a description of the evolution of the state \(\theta\) (here the Ageing law was chosen). Integrating Eq. (3) over time yields \(\vec{X}(t) = \left[ \mu(t), \theta(t) \right]^\intercal\).

The formulations above are merely a classical example of a friction forward model. There are several other models for the friction law, state evolution, and the stress transfer, that can be used (and in some cases combined in different ways). For more background on the various models that have been implemented, see the sections on rate-and-state models, the Chen-Niemeijer-Spiers model, and spring-block models.

The remainder of this section is dedicated to the practical implementation of forward models in DiaBayes.

Implementation in DiaBayes#

DiaBayes provides a generalised interface that takes the three (compatible) forward components as an input, and connects them together into a form equivalent to Eq. (1). This interface is accessed through the diabayes.forward_models.Forward class. The forward model expressed by Eq. (3) is constructed as follows:

from diabayes.forward_models import Forward, rsf, ageing_law, springblock
forward = Forward(
    friction_model=rsf,
    state_evolution={"theta": ageing_law},
    block_type=springblock
)

The Forward class implements a __call__ method that, when given the appropriate variables (\(\vec{X}\)), parameters (\(\left\{ a, b, D_c \right\}\)), friction constants (\(\left\{ \mu_0, v_0 \right\}\)), and spring-block constants (\(\left\{ k, v_{lp} \right\}\)), yields \(\mathrm{d} \vec{X} / \mathrm{d} t\). The forward class instance is passed to an ODE solver, which will successively call forward(...) to obtain the time derivatives for integration.

To enable the assembly of a forward model from the three components, each component needs to follow a specific call signature. These are:

def friction_model(
    variables: Variables, parameters: 
    _Params, friction_constants: _Constants
    ) -> Float:
    velocity: Float = ...
    return velocity

def state_evolution(
    velocity: Float, variables: Variables, 
    parameters: _Params, friction_constants: _Constants
    ) -> Float:
    state_change: Float = ...
    return state_change

def block_type(
    time: Float, velocity: Float, velocity_derivs: Variables, 
    variables: Variables, state_rate: Float[Array, "..."], 
    constants: _BlockConstants
    ) -> Float:
    friction_change: Float = ...
    return friction_change

The type annotations Variables, _Params, _Constants, and _BlockConstants are defined in diabayes.typedefs, and they represent special containers holding the expected items that can be accessed by name. The contents of the containers are specific to a given model (as each model has its own variables and parameters). Note that each function above returns a scalar; vectorised computations are handled at a higher level by vmaping the forward model.

Any new implementations of physics (see Adding custom physics) will have to obey these call signatures to ensure a consistent computational framework.

The forward model is passed to the ODE solver upon instantiating the ODESolver class:

from diabayes.solver import ODESolver
solver = ODESolver(forward_model=forward)

The ODESolver class implements several wrappers around Forward.__call__ to make it compatible with SciPy’s solve_ivp routine for single forward calculations, or the (vmapped) diffrax.diffeqsolve routine for parallelised forward passes during inversion. The inversion methods implemented by this class are discussed in the Inverse modelling section.

References