OPERA Surface Displacement Modeling

• OPERA Surface Displacement is a computational pipeline that models geodetically observable ground movements from subsurface fault slip using state-of-the-art forward solvers and adjoint techniques.
• It integrates analytical models (e.g., Okada) with numerical methods like BEM, FEM, and DG to account for topography and heterogeneous elasticity.
• The approach employs iterative, gradient-based optimization to update fault geometry and slip distributions, achieving millimeter- to centimeter-scale modeling precision.

OPERA Surface Displacement refers to the forward and inverse modeling of geodetically observable ground displacements caused by subsurface elastic dislocations, calculated within the “OPERA” pipeline—a modular system designed for high-precision, physics-based geophysical inference. In this context, “OPERA” denotes not a specific algorithm but a pipeline paradigm in geophysical surface-displacement inversion, integrating state-of-the-art forward solvers, adjoint methods, and gradient-based optimization capable of resolving millimeter- to centimeter-scale surface motion. Accurate OPERA surface displacement modeling is essential for interpreting InSAR and GNSS data, reconstructing active fault geometry, and inferring slip distributions in aseismic and co-seismic deformation scenarios.

1. Physical and Mathematical Formulation

OPERA surface displacement modeling adopts an inhomogeneous, linear elasticity framework for the Earth's crust, embedding seismic faults as open surfaces $\Gamma$ with a prescribed slip vector field $s$. The displacement field $u$ satisfies: $$\begin{cases} -\nabla\cdot\left(\mathbb{C}(x):\nabla^s u\right) = 0 & \text{in } \Omega \setminus \Gamma, \ [\![ u ]\!]_\Gamma = s, \quad [\![ \mathbb{C}:\nabla^s u ]\!]_\Gamma = 0, & \text{across } \Gamma, \ \mathbb{C}(x):\nabla^s u \cdot \nu = 0 & \text{on } S_{\rm surf}, \ u \to 0 & |x| \to \infty \text{ (or } u = 0 \text{ on a distant boundary)} \end{cases}$$ with $\mathbb{C}(x)$ the fourth-order stiffness tensor and $\nabla^s u$ the symmetric gradient. The weak formulation seeks $u$ with $[\![u]\!]_\Gamma = s$ and $u=0$ on the Dirichlet boundary, such that

$$a_\Gamma(u,v) := \int_{\Omega\setminus\Gamma} \mathbb{C}\nabla^s u : \nabla^s v\, dx = 0$$

for all admissible test functions $v$.

These equations describe displacement in response to fault slip, forming the computational core for both forward simulation and, via adjoint methods, the construction of inverse algorithms for slip and geometry recovery (Aspri et al., 2024).

2. Surface Displacement Misfit and Shape Sensitivity

The inversion process in OPERA minimizes the least-squares surface–displacement misfit: $$J(\Gamma, s) = \frac12 \int_{\mathcal{S}} |u(\Gamma, s)(x) - u^{\rm obs}(x)|^2\, dS(x)$$ over fault geometries $\Gamma$ and slip distributions $s$, where $u^{\rm obs}$ is geodetic data on a subset $\mathcal{S}$ of the Earth's surface (Aspri et al., 2024).

The key to efficient optimization is the shape derivative of $J$, derived using material and adjoint field techniques. For small geometric and slip perturbations $(U, h)$: $$D_{U,h}J(\Gamma, s) = \int_{\Omega\setminus\Gamma} \mathbb{M}[\nabla^s u] : \nabla^s p\, dx + \int_\Gamma (\mathbb{C}:\nabla^s p\, n)\cdot h\, dS$$ where $p$ is the adjoint field (solution of an elastic problem with a surface misfit right-hand side), and $\mathbb{M}$ is a polarization tensor involving both the elasticity tensor and perturbation field. For $2$D problems with constant elasticity, explicit boundary expressions are available (Aspri et al., 2024). This machinery exposes sensitivities of surface displacements to variations in fault geometry and slip, supporting gradient-based updates within the inversion pipeline.

3. Analytical and Numerical Approximations

Classical surface displacement calculations use what is commonly termed the "flat-slab" (Okada 1985) model: an isotropic, homogeneous half-space with a planar free surface. In this setting, surface displacements $u(x,y,0) = (u_x, u_y, u_z)$ are available in closed form for rectangular faults parameterized by slip magnitude, fault strike, dip, and depth. For example, the vertical displacement due to pure dip-slip is

$$u_z^{(2)}(\xi, \eta) = -\frac{D_2}{2\pi(1-\nu)} \left[ \frac{\xi\eta}{R(R+\eta)} + (1-2\nu) \arctan \frac{\xi\eta}{qR} \right]_{(\xi, \eta) = \text{corners}}$$

with parameters as defined in the model (Lara et al., 2019). Mogi-type formulas are used for point sources.

Yet, flat-slab models neglect variations in topography, lateral heterogeneities in $\mathbb{C}(x)$, and non-planarity of the free surface. To address these deficiencies, OPERA modules incorporate numerical methods:

• Boundary-integral equations (BEM) are solved directly on triangulated representations of real-Earth topography (TINs from digital elevation models).
• Finite element or discontinuous Galerkin methods discretize elastostatic equations on conforming meshes containing the fault (Aspri et al., 2024).

The choice of formulation impacts achievable modeling accuracy, computational cost, and the ability to incorporate complex fault and slip geometries.

4. Topographic Corrections and Accuracy Assessment

Comparisons between flat-slab and topography-resolving forward models reveal systematic accuracy limitations in the former when applied to regions of significant relief or strong horizontal slip (Lara et al., 2019). For several large earthquakes, the maximum difference between flat and topography-corrected surface displacements reaches up to 11.4% of the peak displacement, with lateral (horizontal) displacement errors generally exceeding vertical errors for events dominated by vertical motion, and the reverse for horizontal-slip events. The empirical relationship

$$\epsilon_{\rm local}(x) \simeq \alpha \cdot |\nabla h(x)|$$

with $\alpha \approx 3$–6% per 100 m of surface relief, describes the correlation of local errors to topographic gradient.

To achieve the millimeter- to centimeter-scale modeling error targets required by OPERA, topography must be fully resolved, either via:

• Numerical boundary-integral solutions on high-resolution surface meshes derived from DEMs.
• Finite-element solutions using surface-conforming meshes.

A pragmatic workflow is to localize slip with rapid Okada-based inversions and subsequently refine surface displacement predictions by topography-corrected forward simulation, eliminating the residual systematic bias (Lara et al., 2019).

Model Type	Error Magnitude	Topography Inclusion	Usage in OPERA
Okada (flat-slab)	Up to 10–19%	No	Rapid inversion
Topography-corrected	<1 cm, regionally	Yes	High-precision modeling

5. Iterative Inverse Algorithms in OPERA

OPERA implements surface displacement inversion as a PDE-constrained optimization process. Using the shape derivative formalism, the geometry (parameterized via mesh vertices or spline control points) and, optionally, slip can be iteratively updated. At each iteration (Aspri et al., 2024):

1. A forward elastostatic problem is solved (typically by DG-FEM) for current $(\Gamma^k,s^k)$.
2. The misfit $J^k$ is assembled from predicted and observed surface displacements.
3. The adjoint elastic problem is solved, with forcing determined by the misfit.
4. For each geometric parameter (e.g., vertex position), the gradient $dJ/d\tau$ is computed via the shape derivative.
5. Descent directions and step sizes are determined (line search or trust-region).
6. Geometry and slip parameters are updated.

Convergence is monitored via decreases in $J$ and maximal parameter displacement, typically stalling at the noise floor of the data. In 2D test settings, convergence to sub-micron accuracy can be achieved in 100–200 iterations, assuming sufficiently rich data (Aspri et al., 2024).

6. Numerical Implementation and Test Cases

A typical OPERA-style implementation for inversion leverages a fine triangular mesh for the forward solution, a coarser one for shape updates, and enforces slip as a non-homogeneous jump on mesh edges overlying the fault. The piecewise polynomial approximation space and interior-penalty terms are key components in the DG discretization (Aspri et al., 2024). Penalization coefficients are tuned with respect to the stiffness tensor and mesh size.

Numerical tests demonstrate:

• Robust recovery for faults with constant slip.
• Modestly degraded but still successful inversion for variable slip distributions.
• Sensitivity of reconstruction stability to observation geometry—at least two-sided measurement scenarios are needed for reliable inversion; one-sided data yields unstable solutions.
• Monotonic decrease in $J$ until convergence to the data noise floor, with typical mesh vertex updates dropping below $10^{-7}$ after 150–200 iterations.

7. Integration of Forward, Adjoint, and Gradient Modules

The modular nature of OPERA allows wrappers around forward, adjoint, and gradient solvers, enabling seamless construction of data-fitting and parameter-updating loops. The advantage of the distributed shape-derivative approach is the compatibility with:

• Heterogeneous, spatially variable elasticity $\mathbb{C}(x)$,
• Arbitrary (including nonplanar) measurement configurations,
• Arbitrary fault geometry initializations.

All essential equations—forward, adjoint, and shape derivative—are encapsulated as plugins or modules, supporting rapid deployment in large-scale geodetic inversion workflows. This structure fosters systematic uncertainty quantification and model refinement as demanded by contemporary seismogeodetic applications (Aspri et al., 2024).

• A. Aspri, E. Beretta, E. Lee, A. Mazzucato, "A shape derivative algorithm for reconstructing elastic dislocations in geophysics" (Aspri et al., 2024)
• R. Armijo et al., "Topography effect on the seismogenic deformation of the earth's surface" (Lara et al., 2019)