---
title: 'OkadaTorch: Differentiable PyTorch Okada Model'
url: https://www.emergentmind.com/topics/okadatorch
type: topic
---

# OkadaTorch: Differentiable PyTorch Okada Model

Searching arXiv for the specified paper and closely related references.
OkadaTorch is a PyTorch implementation of the Okada model in which the entire computation is differentiable, enabling gradients with respect to inputs to be computed using automatic differentiation. It is presented as a direct translation of the original Okada model into PyTorch together with a wrapper interface for efficiently computing gradients and Hessians with respect to either observation station coordinates or fault parameters [2507.17126]. Within the stated scope of a 3D, homogeneous, isotropic elastic half-space with a traction-free surface, it reproduces the functionality of the standard subroutines SPOINT, SRECTF, DC3D0, and DC3D for point and finite rectangular dislocation sources at or below the surface, while making the full computational graph available to PyTorch autograd [2507.17126].

## 1. Position within elastic dislocation modeling

The Okada model is a widely used analytical solution for displacements and strains caused by a point or rectangular dislocation source in a 3D elastic half-space [2507.17126]. OkadaTorch specifically translates the classic Okada analytical solutions for elastic deformation due to shear and tensile dislocations into PyTorch tensor operations, preserving the original subroutine structure while exposing differentiability with respect to model inputs [2507.17126].

The physical setting is the standard half-space formulation: a finite rectangular dislocation embedded in an elastic half-space, with the Earth modeled as homogeneous, isotropic, linearly elastic, and bounded by a traction-free surface [2507.17126]. The source is parameterized by a rectangular fault with length $L$ along strike and width $W$ down dip, orientation given by strike $\theta$, dip $\delta$, and rake $\lambda$, and depth $d$ defined at the top edge or fault center depending on configuration [2507.17126]. Material dependence is expressed through Poisson’s ratio $\nu$, with Okada’s parameter
$$
\alpha = \frac{1 - 2\nu}{1 - \nu},
$$
while shear modulus $\mu$ determines stress scaling but is not required for displacements in the standard formulation [2507.17126].

A central motivation is methodological rather than geophysical novelty. The implementation was developed to support first- and second-order optimization, sensitivity analysis, Bayesian inference with HMC/NUTS, and end-to-end scientific machine learning models without hand-derived symbolic derivatives [2507.17126]. This suggests that its main contribution lies in computational differentiation and workflow integration rather than in altering the underlying elastic theory.

## 2. Fault geometry, coordinates, and kinematic parameterization

OkadaTorch adopts a right-handed global coordinate system with $x=\text{east}$, $y=\text{north}$, and $z=\text{up}$; at the surface $z=0$, and below the surface $z<0$ in its global convention [2507.17126]. Observation stations are specified by coordinates $(x,y,z)$, and the source geometry is referenced to a fault location $(x_0,y_0,d)$ with configurable origin semantics [2507.17126].

The horizontal transformation from global to strike-aligned coordinates is
$$
x' = (x-x_0)\sin\theta + (y-y_0)\cos\theta,
$$
$$
y' = (x-x_0)\cos\theta - (y-y_0)\sin\theta,
$$
with unit vectors
$$
\hat{s} = [\sin\theta,\cos\theta], \qquad \hat{n} = [\cos\theta,-\sin\theta],
$$
where $\hat{s}$ is along strike and $\hat{n}$ is to the right of strike [2507.17126]. Introducing $z_\downarrow=-z$ so that depth is positive downward, the dip-plane transformation is
$$
p = y'\cos\delta + z_\downarrow\sin\delta,
$$
$$
q = y'\sin\delta - z_\downarrow\cos\delta,
$$
with $x=x'$ in Okada’s local coordinates [2507.17126].

For a finite rectangle, the along-strike and down-dip coordinates are denoted $\xi\in[\xi_1,\xi_2]$ and $\eta\in[\eta_1,\eta_2]$ [2507.17126]. In common conventions, $\xi_1=0$, $\xi_2=L$, $\eta_1=0$, and $\eta_2=W$, but OkadaTorch allows different origins via `fault_origin` [2507.17126]. At the wrapper level, `fault_origin` selects whether the fault location refers to the top-left corner (`"topleft"`, default) or the rectangle center (`"center"`) [2507.17126].

Slip is decomposed into three components: strike-slip (`DISL1`), dip-slip (`DISL2`), and tensile/opening (`DISL3`) [2507.17126]. The rake $\lambda$ determines this decomposition, and displacements scale with slip [2507.17126]. In the wrapper interface, the parameter dictionary contains scalar tensors for `"x_fault"`, `"y_fault"`, `"depth"`, `"length"`, `"width"`, `"strike"`, `"dip"`, `"rake"`, and `"slip"` [2507.17126].

## 3. Mathematical structure of the forward model

OkadaTorch follows the analytical structure of Okada’s displacement and strain fields as superpositions of Green’s functions for the three slip modes [2507.17126]. In fault-local coordinates, a finite rectangular source is evaluated by summing signed contributions from its four corners. The representative displacement expression given for surface points is
$$
u_i(x,y,z=0)=\sum_{\xi\in\{\xi_1,\xi_2\}}\sum_{\eta\in\{\eta_1,\eta_2\}} \operatorname{sgn}(\xi)\operatorname{sgn}(\eta)\,[u_i^{ss}+u_i^{ds}+u_i^{ts}](x-\xi,y-\eta,q;\delta,\alpha),
$$
where the superscripts denote strike-slip, dip-slip, and tensile components, respectively [2507.17126].

The constituent kernels are built from elementary terms involving
$$
R=\sqrt{x^2+y^2+q^2},
$$
together with $\log(R+\eta)$ and $\arctan(xy/(qR))$, as well as dip-dependent and $\alpha$-dependent coefficients [2507.17126]. A representative strike-slip kernel is written as
$$
u_i^{ss}=\frac{\mathrm{DISL1}}{2\pi}\left[c_1^i \log(R+\eta)+c_2^i \arctan\!\left(\frac{xy}{qR}\right)+c_3^i F_1(x,y,q,\eta,R)+c_4^i F_2(x,y,q,\eta,R)\right],
$$
with analogous expressions for dip-slip and tensile components [2507.17126]. The full closed forms are attributed to Okada (1985) for surface points and Okada (1992) for subsurface points [2507.17126].

After evaluating local components, displacements are rotated back to global coordinates according to
$$
\begin{bmatrix}
u_x\\
u_y
\end{bmatrix}
=
\begin{bmatrix}
\sin\theta & \cos\theta\\
\cos\theta & -\sin\theta
\end{bmatrix}
\begin{bmatrix}
u_{\text{along strike}}\\
u_{\text{across strike}}
\end{bmatrix},
$$
with $u_z$ handled consistently with the sign convention $z$ up in OkadaTorch [2507.17126].

The strain tensor is defined as
$$
\varepsilon_{ij}=\tfrac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right),
$$
and the translated subroutines compute all spatial derivatives $\partial u_i/\partial x_j$; symmetrization can then be applied if strictly symmetric strain components are required [2507.17126].

## 4. Differentiable programming design

A defining feature of OkadaTorch is that all primitive operations—such as addition, multiplication, square roots, logarithms, and arctangents—are implemented as PyTorch tensor operations, so autograd tracks the entire computational graph [2507.17126]. No Python loops over stations are required, and vectorized tensors enable efficient execution on CPU or GPU [2507.17126].

If `is_degree=True`, which is the default, angular parameters are converted through differentiable scale factors $(\pi/180)$, so the rotation matrices remain differentiable functions of strike and dip [2507.17126]. The wrapper exposes derivative computation in two forms. `OkadaWrapper.gradient` uses `torch.func.jacfwd` to compute derivatives of outputs with respect to a single argument such as `"x"`, `"depth"`, or `"strike"` [2507.17126]. `OkadaWrapper.hessian` computes second derivatives with respect to two arguments by nesting `jacfwd`, allowing both diagonal and mixed entries [2507.17126]. For scalar losses, ordinary reverse-mode autodiff via `loss.backward()` computes gradients with respect to any parameter whose tensor has `requires_grad=True`, and double backward is supported for second derivatives [2507.17126].

The paper also notes the internal use of `torch.func`—formerly `functorch`—in PyTorch 2.x for efficient vectorization through `vmap` during derivative evaluation [2507.17126]. The broader PyTorch framework is described in the cited systems literature [1912.01703]. A plausible implication is that the implementation is intended to align geophysical forward modeling with standard differentiable programming idioms already established in machine learning and scientific computing.

The treatment of non-differentiable behavior is conservative. The core algorithms avoid branching on the differentiated quantities, while functions such as `atan2` and `log` are described as differentiable almost everywhere, with the usual caveats near singular points [2507.17126].

## 5. Interfaces, outputs, and computational workflow

OkadaTorch provides two interface layers: direct translations of the original subroutines and a higher-level wrapper [2507.17126]. The direct layer includes `SPOINT` for surface point sources, `SRECTF` for surface rectangular sources, `DC3D0` for subsurface point sources, and `DC3D` for subsurface rectangular sources [2507.17126]. Compared with the original Fortran, the translated functions add the keyword arguments `compute_strain` and `is_degree` [2507.17126].

For `DC3D`, the documented call signature is
```python
out, IRET = DC3D(ALPHA, X, Y, Z, DEPTH, DIP, AL1, AL2, AW1, AW2, DISL1, DISL2, DISL3, compute_strain=False, is_degree=True)
```
where `ALPHA = (1 − 2ν)/(1 − ν)`, `out` is `[ux, uy, uz]` or displacements plus strains if requested, each with the shape of `X, Y, Z`, and `IRET` is the status integer inherited from the original code [2507.17126]. The values `AL1`, `AL2`, `AW1`, and `AW2` correspond to the fault-plane extent parameters of Okada’s `DC3D` interface and appear explicitly in the core functions [2507.17126].

The high-level `OkadaWrapper` provides:
- `compute(coords, params, compute_strain=True, is_degree=True, fault_origin="topleft", nu=0.25)`
- `gradient(coords, params, arg=..., compute_strain=True, is_degree=True, fault_origin="topleft", nu=0.25)`
- `hessian(coords, params, arg1=..., arg2=..., compute_strain=True, is_degree=True, fault_origin="topleft", nu=0.25)` [2507.17126]

The `coords` argument is a dictionary containing tensors `"x"`, `"y"`, and optionally `"z"`, all of identical shape and interpreted as batched stations [2507.17126]. The `params` argument is a dictionary of scalar torch tensors, each optionally carrying `requires_grad=True` [2507.17126]. Units are required to be consistent for coordinates and geometry, for example kilometers throughout, while slip is specified in meters; outputs are then displacements in meters and strains in displacement per unit length, such as m/km [2507.17126].

If `compute_strain=True`, the wrapper returns a list of 12 tensors:
$$
[u_x,u_y,u_z,\partial u_x/\partial x,\partial u_x/\partial y,\partial u_x/\partial z,\partial u_y/\partial x,\ldots,\partial u_z/\partial z]
$$
from which strain components such as
$$
\varepsilon_{xx}=\partial u_x/\partial x,\quad
\varepsilon_{yy}=\partial u_y/\partial y,\quad
\varepsilon_{zz}=\partial u_z/\partial z,\quad
\varepsilon_{xy}=0.5(\partial u_x/\partial y+\partial u_y/\partial x)
$$
can be assembled [2507.17126]. These derivatives are themselves differentiable with respect to parameters, so quantities such as $\partial \varepsilon_{ij}/\partial p$ or second derivatives can be obtained through the same wrapper methods [2507.17126].

## 6. Performance, validation, and numerical behavior

The implementation is vectorized over stations: the core subroutines accept tensor fields `X`, `Y`, and `Z`, eliminating Python loops and permitting efficient handling of large grids on CPU or GPU [2507.17126]. Its reported computational complexity is roughly $O(N_{\text{stations}}\times N_{\text{faults}})$, with the wrapper applying the four-corner summation internally for each rectangular fault [2507.17126]. GPU acceleration follows directly from the use of PyTorch tensor operations; moving inputs and parameters to CUDA tensors is sufficient to execute on GPU [2507.17126]. The paper also notes a memory-compute trade-off: setting `compute_strain=False` avoids computing and storing intermediate variables needed only for derivatives, reducing compute time and memory when only displacements are needed [2507.17126].

Validation is described primarily in qualitative and workflow terms. The core is stated to be a direct translation of the standard Fortran `DC3D` family, and the wrapper is said to demonstrate physically consistent displacement and strain fields at the surface [2507.17126]. The paper shows forward fields for a representative rectangular fault, including $u_x$, $u_y$, $u_z$, and all nine $\partial u_i/\partial x_j$ fields; first derivatives with respect to each parameter, such as $\partial u_z/\partial p$, illustrating expected sensitivities including linearity in slip; and a synthetic inversion example in which Adam optimization recovers parameters close to ground truth with an interpretable loss history [2507.17126].

At the same time, the reported evidence has stated limits. Unit tests and formal accuracy metrics versus legacy codes are not detailed in the paper [2507.17126]. This suggests that the validation emphasis is on direct translation fidelity, qualitative field behavior, and end-to-end inversion demonstrations rather than on a published benchmark suite.

The numerical cautions closely mirror those of the classical Okada solution. Because the kernels involve logarithms, arctangents, and $1/R$-like terms, evaluations very near a fault edge or at zero distance can be ill-conditioned [2507.17126]. The implementation preserves the original subroutine logic, including the `IRET` status code, and is described as robust for standard use, but users are advised to avoid placing stations exactly on fault edges and to exercise care with very shallow sources and observation points directly above the upper edge [2507.17126]. For gradient-based inversion, regularization or parameter bounds such as minimum depth or minimum width can help when gradients become large near singular configurations [2507.17126].

## 7. Applications, scope, and limitations

The differentiable design enables several use cases explicitly identified in the paper. These include gradient-based optimization of fault parameters from GNSS or InSAR displacements using `torch.optim`, gradient-informed Bayesian inference such as HMC/NUTS, sensitivity analysis and Fisher/Hessian-based uncertainty quantification including Laplace approximations around optima, and integration with scientific machine learning pipelines in which neural networks output slip distributions and a differentiable Okada layer computes the resulting elastic fields [2507.17126]. The paper further notes possible coupling with differentiable tsunami or crustal flow models [2507.17126].

A concise summary of these use cases is given below.

| Use case | Mechanism in OkadaTorch | Stated purpose |
|---|---|---|
| Gradient-based inversion | `loss.backward()` or wrapper derivatives | Fault parameter optimization |
| Bayesian inference | Exact model gradients | HMC/NUTS workflows |
| Uncertainty analysis | Gradients and Hessians | Fisher/Hessian and Laplace approximations |
| SciML integration | End-to-end backpropagation through physics layer | Neural-network-driven slip or coupled models |

The scope is nonetheless restricted by the assumptions of the underlying Okada model. OkadaTorch assumes a homogeneous, isotropic, linear elastic half-space with a traction-free surface [2507.17126]. It does not include layering, topography, viscoelastic or poroelastic effects, rate/state friction, or inelastic processes [2507.17126]. For layered media or time-dependent rheologies, different forward models are required [2507.17126].

The paper also notes that classical identifiability issues remain: trade-offs such as those between depth and slip, or between geometry and rake, are not removed by differentiability, and gradient-based methods can converge to local minima without good initialization and appropriate regularization or priors [2507.17126]. Future extensions mentioned include multi-layered media, topography, smoothing parameterizations for slip distributions, and richer priors for Bayesian workflows [2507.17126].

Practically, the code is available in a public repository together with examples, documentation, and example notebooks for reproducing forward and derivative fields and the synthetic inversion [2507.17126]. The repository is intended for use with PyTorch, with PyTorch 2.x recommended for `torch.func.jacfwd` [2507.17126]. Publication was permitted by NIED, and the programs are stated to differ from those on the NIED site [2507.17126].

Source: https://www.emergentmind.com/topics/okadatorch