MCMC

The Probabilistic Model

In addition to the point estimates afforded by optimization, celeri also provides uncertainty quantification via Markov Chain Monte Carlo (MCMC), which returns samples from a posterior probability distribution over the parameters of interest, governed by a probabilistic model. In this model, parameters of interest $\theta$ are given a joint prior distribution, $p(\theta)$, which is then combined with the observed data to become a posterior distribution. The model is composed of two parts: the prior, denoted $p(\theta)$, and a likelihood $p(v | \theta)$, corresponding to the likelihood of observing data $v$ given parameter values $\theta$. Multiplying these two components together constitutes the (unnormalized) posterior density $p(v|\theta)p(\theta)$, which describes the full data generating process in the Bayesian paradigm: first, parameter values are drawn from the prior, and subsequently, observations are drawn from the data distribution conditioned on those parameter values.

The geophysical forward model in celeri maps the geodynamic parameters (state vector) $\theta$ to predicted ground velocities $v$ at GPS stations. The predicted velocity field at the $N$ stations is

$$ v_\text{total}(\theta) \in \mathbb{R}^{2N}, $$

where each station contributes two observed components (east and north).

The state vector $\theta$ is assembled from multiple physical sources (strain rates, block rotations, Mogi sources, elastic slip, etc.), each contributing additively to the total expected velocity:

$$ v_\text{total}(\theta) = v_\text{strain} + v_\text{rot} + v_\text{rot,okada} + v_\text{mogi} + v_\text{elastic}. $$

The velocity components are built from linear maps $O_\text{strain}$, $O_\text{rot}$, $O_\text{rot,okada}$, $O_\text{mogi}$, and $O_\text{kinematic,fault}$ that encode elastic properties. These linear maps don't depend on $\theta$ and are precomputed by celeri.

$v_\text{total}(\theta)$ decomposes into a sum

$$ v_\text{total}(\theta) = O_\text{strain} \, \theta_\text{strain} + O_\text{rot} \, \theta_\text{rot} + \cdots $$

At $N$ GPS sites we observe noisy (east, north) velocities

$$ y_i \in \mathbb{R}^2, \quad i=1,\dots,N. $$

We model this noise with a shared variance parameter $\sigma^2$, such that the full likelihood is

$$ y_i \sim \mathcal{N}\left(v_\text{total, $i$}(\theta), \sigma^2 I_2\right). $$

The MCMC sampler generates samples of values for $\theta$ that plausibly fit the observed noisy velocities $y_i$.

Geophysical Model Components

Block strain rates

Each deformable block has a homogeneous 2D strain tensor with 3 parameters. We use a standard normal prior on the strain parameters after rescaling the operator columns. The contribution to station velocities is linear:

$$ v_\text{strain} = O_\text{strain} \theta_\text{strain}. $$

Block rotations

Each block rotates about the Earth’s center (3 parameters), again with a standard normal prior after rescaling. The contribution has two parts:

• Station velocities from rigid rotation: $v_\text{rot} = O_\text{rot} \, \theta_\text{rot}$.
• Slip along faults induced by rotation, mapped through Okada dislocation operators $v_\text{rot,okada} = O_\text{rot,okada} \, \theta_\text{rot}$.

Mogi sources

Point sources of pressure at fixed locations, with a scalar intensity parameter each. Again, we use standard normal prior after rescaling. The contribution to the velocity is again linear:

$$ v_\text{mogi} = O_\text{mogi} \, \theta_\text{mogi}. $$

Elastic slip on faults (TDEs)

On meshed fault segments, we want to model how much of the kinematically expected slip rate (from rigid block rotations) is actually accommodated elastically at the fault.

Coupling field as a Gaussian Process (GP)

A Gaussian Process is a stochastic process from which any finite collection of random variables drawn is distributed as a multivariate normal. Read more here.

Define the coupling ratio field

$$ c(s) = \frac{v_\text{e,fault}(s)}{v_\text{k,fault}(s)}, \quad s \in \text{fault surface}, $$

where $v_\text{k,fault}(s)$ is the kinematic slip rate implied by block rotations and $v_\text{e,fault}(s)$ is the actual elastic slip rate at the fault. The kinematic velocities are computed as a linear function of the rotation parameters: $v_\text{k,fault} = O_\text{k,fault} \theta_\text{rot}$.

$c(s)$ is given a Gaussian process prior, encoding smoothness along the fault mesh:

$$ c(s) \sim \mathcal{GP}(0, K(s, s')). $$

This GP is discrete, taking one value for each triangle of the mesh. For $K$ we use a Matérn covariance kernel with $\nu=5/2$ and constant variance and length scale hyperparameters. The Matérn-5/2 covariance between two triangle elements at distance $d$ is

$$ C_{5/2}(d)=\sigma^2\left(1+\frac{\sqrt{5}d}{\rho}+\frac{5d^2}{3\rho^2}\right)\exp\left(-\frac{\sqrt{5}d}{\rho}\right) $$

for lengthscale $\rho$ and variance $\sigma^2$. The lengthscale can be set by the user as either an absolute value (in km) or in units of mesh diameter, defined as the longest pairwise distance between triangle centroids on a mesh. The variance is user-defined. The full covariance matrix $K$ is defined such that $K_{ij}=C_{5/2}(d(i,j))$ where $d$ is the Euclidean distance between the centroids of triangles $i$ and $j$.

We use separate GPs for each fault. At each fault, we compute one GP for the strike-slip component and one for the dip-slip component, and assume independence between them.

Eigen-expansion (Karhunen–Loève truncation)

The GP covariance matrix $K \in \mathbb{R}^{M \times M}$ is computed on the mesh nodes ($M$ = number of fault DOFs). But with $\mathcal{O}(10^2)$ meshes each with $M=\mathcal{O}(10^3)$ triangular elements leads to an infeasible $\mathcal{O}(10^5)$ parameters to estimate, using huge matrices with $\mathcal{O}(10^6)$ entries. To avoid this very high-dimensional inference, we use low-rank approximations of the covariance kernels $K$ of each mesh. To do this, we perform an eigen-decomposition

$$ K = U \Lambda U^\top, \quad \Lambda = \text{diag}(\lambda_1, \dots, \lambda_M), $$

and truncate after $n$ leading eigenmodes (largest $\lambda_i$). Thus,

$$ c(s) \approx \sum_{i=1}^n \alpha_i u_i(s), $$

with coefficients $\alpha_i \sim \mathcal{N}(0, \lambda_i)$, where $u_i$ are the eigenfunctions defined on the mesh. The value of $n$ can be set separately for strike-slip and dip-slip by the user on a mesh-by-mesh basis. The coefficients $\alpha_i$ are the parameters in the model. This reduces the high-dimensional GP to a low-dimensional expansion (also still a GP) that captures the bulk of the smooth variation.

Resulting elastic slip and station velocities

The elastic slip rate field on the fault surface is then

$$ v_\text{e,fault}(s) = c(s) v_\text{k,fault}(s). $$

This field is then mapped through a linear operator to yield predicted velocities at GPS stations:

$$ v_\text{e} = \sum_{\text{faults}}\sum_{\text{ds, ss}} O_\text{tde} v_\text{e,fault}. $$

Alternative direct elastic formulation

Depending on the model configuration, we also support modeling the elastic slip field as a GP, instead of a coupling ratio. In this case, the elastic slip field itself is given a GP prior and eigen-expansion:

$$ v_\text{e,fault}(s) \approx \sum_{i=1}^n \beta_i , u_i(s), \quad \beta_i \sim \mathcal{N}(0, \lambda_i). $$

Which formulation is used is configured in the model setup, and can differ by fault or dip/strike component if desired.

Constraints and Regularization

Bounded transform for constraints

Physically, coupling ratios (or elastic slips) should fall within plausible limits, e.g.

$$ 0 \leq c(s) \leq 1, $$

or some more general $(\ell, u)$ bounds.

Since the MCMC sampler requires an unconstrained and differentiable probability density, we use smooth transforms of $c(s)$ to enforce these bounds. The way this works is that instead of sampling values of the parameter $c(s)$, we sample from a transformed parameter, and apply the inverse transformation when doing prediction. We use a logistic (sigmoid) if both lower and upper bounds are finite, and a softplus if only a single bound is given. So for instance,

$$ c(s) = \text{expit}(c_\text{unc}(s)) (u - \ell) + \ell, $$

where $c_\text{unc}$ is the parameter we infer over, before being transformed to the estimated coupling ratio.

Regularization

In addition to the "hard" constraints imposed on the coupling field, we'd also like to be able to impose "soft" regularization on some quantities of interest, as in the SQP formulation. In a probabilistic model, this can be done by expanding the likelihood function $p(x|\theta)$ to $p(x|\theta) \cdot f(x_c)$, where $f$ is a probability density centered around plausible values of $x_c$. In other words, we expose a latent parameter to the output of the forward model, treating it as a part of the data $x$ to discourage values outside an acceptable range. We use this method to regularize certain un-meshed fault segment slip rates where the rate is known through other observational methods to be within a certain mm/yr range. In this case, these segment slip rates truly are data with a likelihood function. But instead of a normal distribution around point observation, we impose a smooth, boxcar-like function over the acceptable interval. This interval can have both $(l,u)$ bounds or one-sided inequality constraint. For a lower bound $l$, the likelihood is a censored normal distribution where $l$ is the upper bound. This creates the effect that sampled segment slip rate values $s\leq{l}$ have likelihood

$$ 1 - \Phi(\frac{l-s}{\sigma}), $$

and vice versa for an upper bound, degrading rapidly (as governed by $\sigma$) as the slip rate moves further from the interval.

To all other slip rates without observation-based credible intervals, we impose a student's t likelihood with $\nu=5$ to regularize towards 0.

Mesh boundary conditions

We also impose soft zero-slip boundary conditions on the meshes via a pseudo-observations of 0 slip on triangular elements on mesh boundaries. That is, as in the framework above, we add terms $\phi(x|\mu=0, \sigma^2)$ to the likelihood function, where $x$ is the boundary element value and $\phi$ is the Normal PDF with a user-chosen $\sigma^2$.

Likelihood and Posterior

The complete likelihood is

$$ \begin{gathered} \sigma \sim \text{HalfNormal}(\text{scale}=2) \\ y \sim \mathcal{N}\left(v_\text{total}, \sigma^2 I \right). \end{gathered} $$

Since all inequality constraints (e.g. on slip rates or coupling ratios) are enforced by bounding transforms rather than hard truncation, the posterior density is smooth and differentiable, so we can use gradient-based MCMC.

The posterior is sampled with HMC/NUTS (via PyMC + nutpie). Since there are highly correlated parameters (e.g. slip on adjacent fault segments), we use the diagonal plus low-rank mass matrix adaptation scheme from nutpie.

Output

The posterior predictive mean ${\hat v_\text{total}(\theta)}$ is projected back into celeri’s state vector format for downstream analysis. Individual draws can be accessed as estimation.mcmc_draw(chain_idx, draw_idx).