Smooth surface reconstruction of earthquake faults from distributed potency beachballs

Source geometry inference consists of the following three model equations.

The first model equation is the representation theorem for linear continua. Displacements $u_{i}$ of data components $i$ are expressed by linear combinations of potency areal densities $D(\boldsymbol{\xi})$ at coordinates $\boldsymbol{\xi}$ on faults $\Gamma$ (Kikuchi & Kanamori, 1991):

where $G_{i}(\boldsymbol{\xi})$ denotes Green’s function relating $u_{i}$ and $D(\boldsymbol{\xi})$. The same formulation of geometry inference applies to time-dependent problems, and we drop time dependence until our teleseismic applications.

The second model equation is the transformation of double couple sources to fault slips. When the potency density $D$ is double-coupled, it is expressed by the symmetrized tensor product of a slip vector ${\bf s}$ and a fault normal ${\bf n}$:

where the superscript ${\rm T}$ denotes the transpose. Equation (2) tells that an $n$-vector is either of the nodal planes of potency in the beachball diagram standard in seismic moment representations. It is equivalent to one of the normalized sums of the maximum and minimum eigenvectors of $D$, accounting for the sign reversals of these eigenvectors. In this paper, we suppose that the normal $n^{\circ}$ for a point source solution (e.g., the first-motion solution) is known and assume that the $n$-vector for $D$ is its closest nodal plane to $n^{\circ}$. That is, the tilt of the $n$-vector from the reference is assumed to be within 45 degrees. Increasing the number of reference point sources may satisfy this assumption. Note that seismic moment is ill-defined in a strict sense across a stiffness discontinuity (e.g., a plate boundary), but potency is well-defined even across it.

The source inversion via eq. (1) and the nodal-plane extraction by eq. (2) have established methods of analysis. The potency density tensor inversion of eq. (1) is similar to the slip inversion (Yabuki & Matsu’ura, 1992; Yagi & Fukahata, 2011). The nodal-plane extraction of eq. (2) is the same as that for moment beachballs.

What remains is the correspondence between $n$-vectors and shapes of surfaces, and it is given by the definition of the $n$-vector: the normal vector ${\bf n}(\boldsymbol{\xi})$ at a coordinate $\boldsymbol{\xi}$ on a surface $\Gamma$ is orthogonal to the infinitesimal shift $d\boldsymbol{\xi}$ along $\Gamma$ at $\boldsymbol{\xi}$ as

In this paper, a ‘smooth surface’ denotes a surface whose shape is defined by eq. (3), known as a ‘regular surface’ in differential geometry. The orthogonality to the tangential plane defines the normal vector of the curved surface, and a relative location between two arbitrary points on a smooth surface is given by the line integral $\int d\boldsymbol{\xi}$ that satisfies eq. (3). Since the $n$-vector field is equivalent to the shape of the smooth surface excepting its rigid-body translation, it entirely records the macroscopic smooth shape of a fault. Hereafter, one location $\boldsymbol{\xi}^{\circ}$ on (or at least adjacent to) the fault (e.g., the rupture initiation point) is assumed to be known, and its estimation errors are neglected.

Let us suppose that the potency density $D$ and the associated $n$-vector field are evaluated given the well-known equations (1) and (2). Our focus is on how to translate them into $\Gamma$ through eq. (3). However, the domain of definition for the observed potency density $D$ depends on the fault shape $\Gamma$ (eq. 1). To perform the inference of $D$ prior to reconstructing $\Gamma$, the representation theorem (eq. 1) must be projected onto some fixed surface, termed the reference surface $\Gamma^{\circ}$ (Fig. 1).

We now introduce $\Gamma^{\circ}$ as a flat plane perpendicular to ${\bf n}^{\circ}$ passing through $\boldsymbol{\xi}^{\circ}$ (Shimizu et al., 2020, also see Fig. 1); if necessary, multiple flat planes are incorporated. Coordinates $\boldsymbol{\xi}$ on $\Gamma$ are parametrized as a function of the location $(x,y)$ taken on $\Gamma^{\circ}$ and the elevation $z$ of $\Gamma$ from $\Gamma^{\circ}$. The positive direction of $z$ is directed to ${\bf n}^{\circ}$. Then, basis function expansions are applied to the potency density by using basis functions that depend only on $(x,y)$. The $q$-component $D_{q}$ of the potency density tensor is expressed by bases $X_{j}(x,y)$ ($j=1,...,M$), and the coefficients are stored in a vector ${\bf a}$ as $a_{j+(q-1)M}$:

Substituting it into eq. (1), we have

where $H_{ij}(\Gamma)=\int_{\Gamma}d\Sigma G_{i}(\boldsymbol{\xi})X_{j}(x,y)$. Note the domain of the surface integral of $H_{ij}(\Gamma)$ is still the unknown $\Gamma$. In the PDTI, $H_{ij}(\Gamma)$ is regularly approximated by a known surface $\Gamma^{\prime}$ (e.g., $\Gamma^{\circ}$) or an estimate $\hat{\Gamma}$ [i.e., $H(\Gamma)\simeq H(\hat{\Gamma})$, see §6.1 for details]. In the above manner, the object to simultaneously estimate (fault potency and fault geometry) is decomposed into model parameters ${\bf a}$ and $\Gamma$. They are formally independent in the above formulation, and their dependencies are described by an additional constraint, which is the theme of this study.

The governing equation of the problem is now the definitional identity of the $n$-vector (eq. 7). It indicates the orthogonality between the $n$-vector and a curved surface, which in three dimensions means more than that a 90-degree rotation (cotangent) of the $n$-vector gives the downward surface slope; introducing

we can rewrite eq. (7) as

Namely, the $n$-vector indicates the direction of the slope by its horizontal projection ${\bf n}_{\parallel}$ ($|{\bf n}_{\parallel}|=1$) as well as the downward slope by its horizontal-vertical ratio $n_{\perp}/n_{z}$ (Fig. 3).

The remaining analysis is straightforward when the $n$-vector is parametrized in a spherical coordinate. Now the polar angle $\theta\in[0,\pi]$ is set as the angle between ${\bf n}^{\circ}$ and ${\bf n}$, and the azimuthal angle $\phi\in[0,2\pi)$ is measured from the $x$ axis. For this setting,

Note $n_{x}=\sin\theta\cos\phi$ and $n_{y}=\sin\theta\sin\phi$. Plugging eq. (21) into eq. (20),

Hereafter, the basis vectors for the position $(\theta,\phi)$ on a unit sphere pointing toward ${\bf n}$ are denoted as ${\bf n}$ (along the $n$-vector), ${\bf e}_{\theta}$ (parallel to the $\theta$ axis), and ${\bf e}_{\phi}$ (parallel to the $\phi$ axis).

The shape of $\Gamma$ is equivalent to both the $z$ field and the $n$ field, but formally, eq. (7) can be read as a partial differential equation of $z$ with a random external force $-(n_{\perp}/n_{z}){\bf n}_{\parallel}$ imposed. The probability $P(\Gamma|{\bf a};\kappa)$ (eq. 14) describes fluctuations of the $n$-vector from the expectation $n_{*}$. We derive an analytical solution of this stochastic partial differential equation. For this derivation, we assume non-integrable components (eq. 12) of $n_{*}$ are small; the full-order derivation becomes complicated (Appendix B).

First, we rewrite the Cartesian forms of $P(\Gamma|{\bf a};\kappa)$ (eq. 14) and the partial differential equation of the elevation (eq. 7) in the spherical coordinate. In a spherical coordinate system at a unit spherical position ($\theta_{*}$, $\phi_{*}$) of ${\bf n}_{*}$, ${\bf n}$ is expressed as

with associated coefficients ($\chi_{n}^{*}$, $\chi_{\theta}^{*}$, $\chi_{\phi}^{*}$). The first order of the Taylor series with respect to angle fluctuations $\theta-\theta_{*}$ and $\phi-\phi_{*}$ around $n_{*}$ yields

where $\mathcal{O}(\delta\Omega^{k})$ denotes the $k$-th order angle fluctuations from $n_{*}$. Note $\partial{\bf n}(\theta,\phi)/\partial\theta={\bf e}_{\theta}$ and $\partial{\bf n}(\theta,\phi)/\partial\phi={\bf e}_{\phi}\sin\theta$. The normalization condition ${\bf n}\cdot{\bf n}=1$ of ${\bf n}$ and the above expressions of $\chi_{\theta}^{*}$ and $\chi_{\phi}^{*}$ result in

Therefore, for the von Mises-Fisher distribution of eq. (14), we have

and

where $\mathcal{Z}^{\prime}$ represents normalization. Only the second-order fluctuations from $n_{*}$ matter in the dislocation limit $\kappa\to\infty$ if non-integrable components (eq. 12) of $n_{*}$ are small.

We next evaluate the partial differential equation of the elevation (eq. 7), which is expressed as eq. (22) in the spherical coordinate. The Taylor series of eq. (22) in angle fluctuations from $n_{*}$ leads to

where $z_{*,\alpha}:=-n_{\alpha}/n_{z}$. Using eq. (21) again and taking the inner and cross products of $\nabla z$ and ${\bf n}_{*\parallel}$, we find

or explicitly, recalling ${\bf n}_{\parallel}\cdot\nabla z=-n_{\perp}/n_{z}$ and ${\bf n}_{\parallel}\times\nabla z=0$ (eq. 20),

Plugging eqs. (31) and (32) into $P(\Gamma|{\bf a})$ (eq. 27 with $\kappa\to\infty$), we arrive at

Here $\sin\theta=n_{\perp}$ (eq. 21) is also used. The area element $d\Sigma(\boldsymbol{\xi})$ along the fault is larger than the area element $dxdy$ along the reference plane because of the tilt of the fault normal by $\theta$:

Note $|n_{z}|=n_{z}$ for $\theta$ below 45 degrees now considered. In the dislocation limit, with small non-integrable components (eq. 12) of $n_{*}$, the higher order fluctuations from $n_{*}$ are dropped, and the most probable estimate $\Gamma_{*}$ for given $n_{*}$ is expressed as:

where ${\rm argmin}$ denotes the arguments of the minima. To avoid zero divisions at $n_{*z}=0$, $n_{*z}^{2}$ among the weight $|n_{*z}^{3}|$ in the first term was placed inside the brackets. The first term from the $\theta$ fluctuation constrains the ${\bf n}_{\parallel}$-parallel component of the fault slope $\nabla z$, and the second term from the $\phi$ fluctuation constrains the ${\bf n}_{\parallel}$-perpendicular component of $\nabla z$. Their weights $|n_{*z}^{3}|$ and $|n_{*z}|$ derive from the independence of $P(\Gamma|{\bf a})$ from the reference plane orientation, and consequently, the constraints become weaker for steeper slopes (i.e., smaller $n_{*z}$ for larger $\theta(<\pi/4)$: $\partial n_{*z}/\partial\theta_{*}=-\sin\theta_{*}<0$).

Thus, equations to be solved (eqs. 7 and 14) are reduced to a least square form of eq. (35). The rest of analysis is the ordinary linear inversion, which we omit to detail. We have solved it using a basis function expansion in this paper, instead of solving it in the continuous space, where the $x$ and $y$ derivatives of the $n$-vector appear, complicating implementations outside differentiable programming. The basis functions employed are the quadratic B-spline functions as in the PDTI, with intervals 1.5 times wider than that of the PDTI. The surface integral of the loss function is collocated at the $x$-$y$ positions of the spline knots adopted in the PDTI. We have removed the rigid-body translation of the fault by specifying a point with zero elevation. This translation of a fault is unconstrained in eq. (35) and does not affect fault elevations regardless of the way of removal.

The resultant estimation remains the least square using eqs. (31) and (32) even when generalizing $P(\Gamma|{\bf a})$ to account for the expected elevation field $\bar{z}$ such that

Our PDTI application includes such terms for computational acceleration (§6.2). The scope of eq. (36) covers a generalization using the probabilities of source locations, such as surface reconstruction from the locations and mechanisms of distributed point sources. It is also applicable to incorporating the constraints on fault edge locations, such as accounting for fault surface traces and jointing multiple fault surfaces. Note $P(\Gamma|{\bf a})$ other than eq. (14) is reduced to the quadratic form of eq. (36) when the higher orders are dropped. The explicit representation of those higher orders dropped in the dislocation limit is given by the non-perturbative expression of the $n$-vector as a function of $\nabla z$ (eq. 37 in the next section).

The first test treats the uniaxial twist shown in Fig. 4. Its elevation is symmetrical about the short side $y$ (see Fig. 4 caption). Except for the center line $y=0$, where $n_{x}=0$ identically, both $n_{x}$ and $n_{y}$ are non-zero, and the $n$-vector varies three-dimensionally. The two-dimensional vizualizaiton of the shape is a birds-eye view hereafter.

The estimation results with a noise standard deviation of 0.05 are indicated in Fig. 5. Despite that the input $n$-vectors are noisy, the estimated elevation is consistent with the original shape. Concavities and convexities at the four corners of the original shape are reproduced while the contours shift. As mentioned in the previous section, the noise component ($\Delta n=n_{*}-n_{0}$) almost surely does not satisfy the slope integrability condition (eq. 11). Therefore, when $n_{*}$ is disturbed from the $n$-vector of a surface, the most probable surface $\Gamma_{*}$ (eq. 35) for given $n_{*}$ could selectively remove the noise in $n_{*}$.

The second test treats a bowl-shaped surface shown in Fig. 6, which is similar to the result of our actual data analysis concerning the 2013 Balochistan earthquake (§6.3). The elevation $z$ is expressed by a function that is parabolic along the long side $x$ and linear along the short side $y$ (see Fig. 6 caption). When $x$ is interpreted as the horizontal axis and $y$ as the vertical, the surface horizontally bends well at the top and increases vertical tilts at the horizontal end (larger $|z_{,x}|$ for larger $y$ and larger $z_{,y}$ for larger $|x|$, roughly speaking).

The results for a standard deviation of 0.05 are indicated in Fig. 7 as in the first test. The reconstructed shape seems more stable than the reconstructed uniaxial twist. This difference suggests that when the $n$-vector noise is at the same level, short-wavelength features such as topography are relatively susceptible to the $n$-vector noise. Considering that wider-regional elevations are controlled by longer-wavelength slopes, we ascribe this tendency to the spatial averaging that low-passes the slope error in a long-wavelength scale. The higher surface reproducibility at longer wavelengths might be reinforced in the application to the PDTI, to which the spatial potency smoothing is commonly applied.

Figure 8 summarizes the errors and residuals of the above two tests with variable noise levels. The estimation errors $1-\hat{n}\cdot n_{0}$ and data residuals $1-\hat{n}\cdot n_{*}$ for normal vectors (Fig. 8a, b) are both linear with the noise amplitude $1-n_{0}\cdot n^{*}$. The data residuals are visibly smaller than the noise amplitude, indicating significant variance reduction. The estimation error $1-\hat{n}\cdot n_{0}$ is smaller than the data residual $1-\hat{n}\cdot n_{*}$, suggesting that the output denoises the given $n$-vector rather than explains it. The spatial average of elevation errors $|z-\bar{z}|$ ($\Gamma_{*}$ in Fig. 8c, d) is more scattered than that of normal vectors (Fig. 8a, b), and its average is roughly proportional to the square root of the noise amplitude of the normal vectors $1-n_{0}\cdot n^{*}$ that approximates squared angle fluctuations. The elevation residuals are roughly the same between the two shapes (Fig. 8c, d), and these same-order absolute residuals indicate a significant difference in relative residuals, supporting higher reproducibility of a bowl-shaped curve compared to uniaxially twisting roughness observed in the previous subsections.

For comparison, we have also evaluated quasi-two-dimensional surface construction from Shimizu et al. (2021) (quasi 2D in Fig. 8c, d). It evaluates the elevation as $z(x,y)\simeq z(\bar{x},\bar{y})+\int_{\bar{x}}^{x}dx^{\prime}\bar{z}_{,x}(x^{\prime})+\bar{z}_{,y}(x)(y-\bar{y})$ by using the mean slope $\bar{z}_{,\alpha}=-\int d_{y}(n^{*}_{\alpha}/n^{*}_{z})/\int dy$ integrated over the short side $y$. We assume $z(\bar{x},\bar{y})$ at the reference plane center $(\bar{x},\bar{y})$ is known in this quasi-two-dimensional reconstruction. The stochastic surface reconstruction $\Gamma_{*}$ works well as quasi-two-dimensional reconstruction for small noise, and for large noise, rather $\Gamma_{*}$ has smaller estimation errors. Even for a uniaxial twist, for which quasi-two-dimensional approximation is errorless, the accuracy of $\Gamma_{*}$ is higher, though the difference is less than accuracy variations (Fig. 8c). Meanwhile, for a bowl shape similar to the fault shape inferred in our application (§6.3), the accuracy of $\Gamma_{*}$ is a few times greater than the quasi-two-dimensional benchmark for the maximum noise amplitude in this test (Fig. 8d).

The fault we use here has a narrow width, so quasi-two-dimensional interpolation from Shimizu et al. (2021) provides a good approximation if no noise exists. Even for such high-aspect-ratio source faults, nonetheless, the results of this study indicate that a probabilistic approach is beneficial for surface reconstruction from noisy $n$-vector fields. Rather than reducing noise through the law of large numbers associated with averaging along the short sides, denoising based on the a priori knowledge that $n$-vector components that do not satisfy the slope integrability condition (eq. 11) are fictitious is likely to better reproduce the original shape.

For the application, we combine the prior constraint on fault geometry $P(\Gamma|{\bf a})$ with the posterior probability of the potency density in the orthodox PDTI. We have considered static problems so far, and hereafter we also consider dynamic problems to invert teleseismic data. The result obtained thus far is applicable as is even when the representation theorem becomes time-dependent and the basis functions are spanned spatiotemporally.

The formulation of the orthodox PDTI, both for static and dynamic problems, consists of the potency likelihood $P({\bf d}|{\bf a},\Gamma;\boldsymbol{\sigma}^{2})$, which is given by data ${\bf d}$ through the representation theorem, and the prior $P({\bf a};\boldsymbol{\rho}^{2})$ for the potency density subject to the basis function expansions ${\bf a}$. The scales of variances for the likelihood $\boldsymbol{\sigma}^{2}$ and weight factors for prior loss functions $\boldsymbol{\rho}^{2}$ are simultaneously estimated. Usually, $\boldsymbol{\sigma}^{2}$ represents the scales of observation errors and Green’s function errors (Yagi & Fukahata, 2011), and $\boldsymbol{\rho}^{2}$ represents the smoothing intensities with respect to space and time (Shimizu et al., 2020). The problem to be solved here consists of these $P({\bf d}|{\bf a},\Gamma;\boldsymbol{\sigma}^{2})$ and $P({\bf a};\boldsymbol{\rho}^{2})$ and the prior for the surface $P(\Gamma|{\bf a})$, and its formal solution is the joint posterior of ${\bf a}$ and $\Gamma$:

where

The analysis for the dislocation limit is simplified when the joint posterior $P({\bf a},\Gamma|{\bf d};\boldsymbol{\sigma}^{2},\boldsymbol{\rho}^{2})$ is decomposed into the marginal posterior of the potency $P({\bf a}|{\bf d};{\boldsymbol{\sigma}}^{2},{\boldsymbol{\rho}}^{2})=\int d\Gamma P({\bf a},\Gamma|{\bf d};{\boldsymbol{\sigma}}^{2},{\boldsymbol{\rho}}^{2})$ and the conditional posterior of the surface $P(\Gamma|{\bf a},{\bf d};{\boldsymbol{\sigma}}^{2},{\boldsymbol{\rho}}^{2})$:

From eq. (16),

That is, the evaluation of $P({\bf a},\Gamma|{\bf d};\boldsymbol{\sigma}^{2},\boldsymbol{\rho}^{2})$ is reduced to the alternate recursion of the most probable estimation of $\Gamma$ given ${\bf a}$ returning $\Gamma_{*}$ (eq. 43) and the inference of ${\bf a}$ given $\Gamma_{*}$ (eq. 42).

We evaluate the hyperparameter optimals $\hat{\boldsymbol{\sigma}}^{2}$ and $\hat{\boldsymbol{\rho}}^{2}$ by Akaike’s Bayesian Information Criterion (ABIC; Akaike, 1980):

The associated optimals of the potency density and fault surface are evaluated by the posterior maximum:

Neglecting the third- and higher-order moments from the probability peaks as in Shimizu et al. (2021), eqs. (40) and (42) are simplified as

and

Besides, from eqs. (41) and (43),

Thus, our optimum-value estimation is the alternate recursion of solving for $\hat{\Gamma}$ given $\hat{\bf a}$ and solving for $(\hat{\bf a},\hat{\boldsymbol{\sigma}}^{2},\hat{\boldsymbol{\rho}}^{2})$ given $\hat{\Gamma}$. The analysis with a fixed $\Gamma$ is exactly the orthodox source inversion with fixed geometry (Yagi & Fukahata, 2011; Shimizu et al., 2020), and the analysis of $\Gamma$ given ${\bf a}$ is what this study has examined. Once $\hat{\Gamma}$ and $(\hat{\boldsymbol{\sigma}}^{2},\hat{\boldsymbol{\rho}}^{2})$ are obtained, the uncertainty of ${\bf a}$ is evaluable through eq. (47) by using $P({\bf a}|{\bf d};\hat{\boldsymbol{\sigma}}^{2},\hat{\boldsymbol{\rho}}^{2})$. Hyperparameter uncertainty is similarly evaluable using $P({\bf d};{\boldsymbol{\sigma}}^{2},{\boldsymbol{\rho}}^{2})$ for given $\hat{\Gamma}$, though it is a sharply peaked well-behaved distribution and follows the law of large numbers exceptionally rapidly when ${\boldsymbol{\sigma}}^{2}$ and ${\boldsymbol{\rho}}^{2}$ are exponential loss weights (Sato et al., 2022).

The above procedure does not include explicit $\Gamma$ uncertainty evaluations. In the present problem setting using a delta-functional $P(\Gamma|{\bf a})$, the estimation error of $\Gamma$ solely derives from error propagation from the potency density estimation, and its evaluation requires another joint posterior decomposition $P({\bf a},\Gamma|{\bf d};{\boldsymbol{\sigma}}^{2},{\boldsymbol{\rho}}^{2})=P(\Gamma|{\bf d};{\boldsymbol{\sigma}}^{2},{\boldsymbol{\rho}}^{2})P({\bf a}|\Gamma,{\bf d};{\boldsymbol{\sigma}}^{2},{\boldsymbol{\rho}}^{2})$ and the evaluation of $P(\Gamma|{\bf d};{\boldsymbol{\sigma}}^{2},{\boldsymbol{\rho}}^{2})=\int d{\bf a}P({\bf a},\Gamma|{\bf d};{\boldsymbol{\sigma}}^{2},{\boldsymbol{\rho}}^{2})$. Due to the slope integrability (eq. 11) on a smooth surface, the posterior of the fault surface $P(\Gamma|{\bf d};{\boldsymbol{\sigma}}^{2},{\boldsymbol{\rho}}^{2})$ is deviated from the posterior of the $n$-vector $n_{*}$ field expected from the potency density (Appendix A). The evaluation of $P({\bf a}|\Gamma,{\bf d};{\boldsymbol{\sigma}}^{2},{\boldsymbol{\rho}}^{2})$ corresponds to the conventional slip inversion, where the potency density ${\bf a}$ is forced to be a displacement discontinuity field across given $\Gamma$.

The computations of $\hat{\bf a}$, $\hat{\Gamma}$, and $(\hat{\boldsymbol{\sigma}}^{2},\hat{\boldsymbol{\rho}}^{2})$ (eqs. 44 and 45) as per eqs. (46)–(48) are formally the same as those of Shimizu et al. (2021) that alternately evaluate $({\bf a},{\boldsymbol{\sigma}}^{2},{\boldsymbol{\rho}}^{2})$ for given $\Gamma$ and $\Gamma$ for given ${\bf a}$. Therefore, simply updating their fault reconstruction stage enables their routine to reconstruct three-dimensional fault geometry.

Figure 9 outlines our computational routine following Shimizu et al. (2021). (I: Fault model initialization) A structured grid is set along the reference plane, and the spatial knots for potency density interpolations are placed at grid centers; non-planar initializations following stage V below are optionally applied. (II: PDTI with Green’s function computations) The PDTI is performed, with Green’s function approximated by point-source representations; when the point sources are allocated to grid centers, they are weighted by the areas of gridded subsurfaces. (III: Nodal plane extraction) The $n$-vector estimates of point sources are computed as nodal planes of estimated potency beachballs; the $n$-vector field estimate is approximately computed by their interpolation. (IV: Surface reconstruction) The elevation field is reconstructed from the $n$-vector field estimate. (V: PDTI with Green’s function updates) The aforementioned structured grid is vertically projected onto the reconstructed fault surface, and the vertical positions of sources and the areas of gridded subsurfaces are updated accordingly; the PDTI updates the potency beachballs and the loop returns to the stage III.

In summary, once the PDTI is performed, potency beachballs set an $n$-vector field, which in turn reconstructs the fault surface, the continuous information of which is subdivided as a digital elevation model to update Green’s function for the next PDTI. The digital elevation model is now computed by the Visualization Toolkit (Schroeder et al., 1998), and we call it by using a Python helper PyVista (Sullivan & Kaszynski, 2019). We execute the following two operations via PyVista: (A) the projection of the structured grid generated by numpy.meshgrid onto the reconstructed fault surface (pyvista.StructuredGrid) and (B) the computation of source locations and subdivided areas on vertically projected structured grids (compute_cell_sizes, cell_centers).

We solve this recursion by a gradient descent, where the $\Gamma_{*}$ determination (eq. 35) is modified by a penalty term on the residual elevation from a given shape $\bar{\Gamma}$ with elevation $\bar{z}$ as

and each surface update accounts for the pre-update state:

The elevation $\hat{z}^{(n)}$ of the $n$-th iteration $\hat{\Gamma}^{(n)}$ is updated to $\hat{z}^{(n+1)}=\hat{z}^{(n)}+\lambda\partial[\kappa^{-1}\ln P(\Gamma|\hat{\bf a}^{(n)})]/\partial z$, and $\lambda$ is the learning rate of this gradient descent. When iterations converge, $\hat{z}^{(n)}=\hat{z}^{(n+1)}$ is met, and the surface optimal $\hat{\Gamma}=\hat{\Gamma}^{(n+1)}$ is the solution that satisfies $\partial\ln P(\Gamma|\hat{\bf a})/\partial z=0$ for $\hat{\bf a}=\hat{\bf a}^{(n)}$.