# The mechanics of static non-planar faults in infinitesimal strain theory

Pierre Romanet, Tatsuhiko Saito, Eiichi Fukuyama

► To cite this version:

Pierre Romanet, Tatsuhiko Saito, Eiichi Fukuyama. The mechanics of static non-planar faults in infinitesimal strain theory. Geophysical Journal International, 2024, 239 (3), pp.1664-1693. 10.1093/gji/ggae337. hal-04844954

HAL Id: hal-04844954

https://hal.science/hal-04844954v1

Submitted on 20 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The mechanics of static non-planar faults in infinitesimal strain theory

Pierre Romanet<sup>1,2</sup> Tatsuhiko Saito<sup>1,3</sup> and Eiichi Fukuyama<sup>1,3</sup>

Accepted 2024 September 12. Received 2024 September 3; in original form 2023 December 22

## SUMMARY

Fault geometry is a key factor in controlling the mechanics of faulting. However, there is currently limited theoretical knowledge regarding the effect of non-planar fault geometry on earthquake mechanics. Here, we address this gap by introducing an expansion of the relation between fault traction and slip, up to second order, relative to the deviation from a planar fault geometry. This expansion enables the separation of the effects of non-planarities from those of planar faults. This expansion is realized in the boundary integral equation, assuming a small fault slope. It provides an interpretation for the effect of complex fault geometry on fault traction, for any fault geometry and any slip distribution. Hence, the results are also independent of the friction that applies on the fault. The findings confirm that fault geometry has a strong influence on in-plane faulting (mode II) by altering the normal traction on the fault and making it more resistant to slipping for any fault geometry. On the contrary, for out-of-plane faulting (mode III), fault geometry has a much smaller influence. Additionally, we analyse some singularities that arise for specific fault geometries often used in earthquake simulations and provide guidelines for their elimination. To conclude this study, we discuss the limits of the infinitesimal strain theory when non-planar faults are considered.

**Key words:** Earthquake dynamics; Mechanics, theory, and modelling; Dynamics and mechanics of faulting.

## 1 INTRODUCTION

Understanding the mechanical effect of fault geometry is an important question for seismology, as it can influence various aspects of earthquake mechanics, such as the areas of the fault that will nucleate, the size of the rupture, or the slip distribution (Aki 1979; King & Nabelek 1985; Schwartz & Sibson 1989; Klinger et al. 2006; Wesnousky 2006, 2008; Milliner et al. 2015). Previous theoretical efforts to understand the effect of fault non-planarities have been primarily focusing on particular fault geometries such as sinusoidale geometry (Saucier et al. 1992; Chester & Chester 2000), or rough faults (Dieterich & Smith 2009; Dunham et al. 2011; Fang & Dunham 2013; Sagy & Lyakhovsky 2019; Morad et al. 2022; Maurer 2024), limiting the generalizability of the results.

In this paper, we extract the effect of non-planarities from the effect of planar fault by performing an expansion of the boundary integral equation. It makes the results presented in this paper very general as they do not depend on one specific geometry, or one specific slip distribution. It also makes the results presented here independent of the friction law that applies on the fault.

Our emphasis will be on the two modes of shearing fault: in-plane faulting and out-of-plane faulting. However, the results presented in this paper can also be applied to opening faults, and the associated findings can be found in Appendix G. We will use a slight abuse of language by sometimes calling mode I, mode II and mode III what should respectively be called in-plane opening, in-plane shearing and out-of-plane shearing. Mode I corresponds to pure opening, Mode II corresponds to a pure shear fault, where the direction of rupture is parallel to the direction of slip. Mode III corresponds to a pure shear fault where the direction of rupture is perpendicular to the direction of slip (see Fig. 1). The primary objective of this paper is to provide a comprehensive explanation of the main effects of non-planar fault geometry on the traction that applies on the fault.

<footnote n="1"><sup>1</sup>Department of Earth Sciences, La Sapienza University of Rome, Piazzale Aldo Moro 5, 00185 Roma, Italy. E-mail: romanet@geoazur.unice.fr</footnote>

<footnote n="2"><sup>2</sup>Université Côte d'Azur, CNRS, IRD, Observatoire de la Côte d'Azur, Géoazur, Sophia-Antipolis, 06560 Valbonne, France</footnote>

<footnote n="3"><sup>3</sup>Earthquake and Tsunami research division, NIED, 3-1 Tennodai, Tsukuba-shi, Ibaraki, 305-0006, Japan</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Description of the different modes of slip for a non-planar fault. Please note that we are using a slight abuse of language by calling them mode I, mode II and mode III.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Definition of the parameters used in this study. The flat fault approximation is when the slope between any two points of the fault is considered small. This figure is modified from Romanet et al. (2020).</figcaption></figure>

## 2 TRACTIONS ON A NON-PLANAR FAULT

This paper derives exact semi-analytical solutions for the elastic strains and stresses within a homogeneous, infinite, isotropic, linear, static, 2D medium caused by slip on a finite, non-planar fault (see Appendix G, H and I). It can be viewed as the static, 2D extension of the previous work by Romanet et al. (2020), where the semi-analytical solution was derived for a fully dynamic, 3D medium using the Einstein notation. The 0th order solutions for in-plane (mode II), in space domain, were previously derived in Romanet et al. (2020). Here, we performed a comprehensive derivation for all mode of slip, in both space and spectral domain, and up to the 1st order, and we provide a physical interpretation of each term.

A complete mechanical model would need to solve the equilibrium for shear traction on the fault $\tau_f = \tau_{el} + \tau_{load}$, where $\tau_f$ is the frictional resistance of the fault, $\tau_{load}$ is the shear loading on the fault, and $\tau_{el}$ is the elastic shear traction response due to slip on the fault. Unless stated otherwise, we are not solving the equilibrium of shear traction on the fault in this paper, but only assuming the slip distribution, and examining how the elastic shear $\tau_{el}$ and normal $\sigma_{el}$ traction respond to slip on a non-planar fault.

The geometry of the fault is described by a function $h(y_1)$ representing the height of the fault at given position $y = (y_1, y_2) = (y_1, h(y_1))$ (Fig. 2). The derivative of the height with respect to the coordinate $y_1$ (the fault slope) is denoted as $m(y_1) = \frac{d}{dy_1}h(y_1)$. When the fault slope is small ($m(y_1) \ll 1$), the second order derivative of height with respect to the coordinate $y_1$ can be linked to the curvature along the fault $\kappa^t$ by:

$$ \begin{aligned} \kappa^t(y_1) &= \frac{\frac{d}{dy_1}m(y_1)}{(1 + m^2(y_1))^2} \\ &\approx \frac{d}{dy_1}m(y_1) = \frac{d^2}{dy_1^2}h(y_1), \text{ if } m(y_1) \ll 1. \end{aligned} \quad (1) $$

The absolute value of the fault curvature $|\kappa^t|$ can be interpreted geometrically as the inverse radius of the tangential circle to the fault (see Fig. 2). When the fault is locally flat, the curvature is zero ($\kappa^t(y_1) = 0$).
<figure><img src="image_4.png" /><figcaption>**Figure 3.** A non-planar fault can be made equivalent to a planar fault with additional on-fault stresses correction. This is not working for off-fault stresses.</figcaption></figure>

This paper demonstrates that the semi-analytical solution for the on-fault stresses can be significantly simplified when assuming a small fault slope. **For the on-fault stresses**, it is possible to identify the exact mechanical corrections that renders a planar fault mechanically equivalent to a non-planar fault by including additional terms that account for the non-planar geometry (see Fig. 3). In order not to distract the reader from the main points of this paper, most of the derivation is provided in Appendix B (see also Romanet et al. (2020) for a more general derivation). In the following, we just provide a summary of the principal steps involved in the derivation:

(i) We start from the representation theorem (Aki & Richards 2002, chapter 2.5) which allows to calculate the displacement $u_k$ anywhere in the medium by evaluating an integral over the fault involving the Green's function $G_{kp}$ (Tada & Yamashita 1997) and the displacement discontinuity $\Delta u_i = u_i^+ - u_i^-$ across the fault. $u_i^+$ and $u_i^-$ represent the displacement on each side of the fault, as shown in Fig. 2. The subscripts refer to specific component of the vector, for example $u_k$ is the $k^{th}$ component of the slip vector **u**. $c_{ijpq}$ is the Hooke tensor, and **n** is the normal vector to the fault. Finally, **x** is the location at which the slip is evaluated and **y** is the variable over which the integration along the fault is performed:

$$u_k(\mathbf{x}) = - \int_{\text{fault}} c_{ijpq} \Delta u_i(\mathbf{y}) n_j(\mathbf{y}) \frac{\partial}{\partial x_q} G_{kp}(\mathbf{x}, \mathbf{y}) d\xi(\mathbf{y}). \quad (2)$$

(ii) By using the strain definition $\epsilon_{cd} = \frac{1}{2} (\frac{\partial}{\partial x_c} u_d + \frac{\partial}{\partial x_d} u_c)$ and the Hooke's law ($\sigma_{ab} = c_{abcd} \epsilon_{cd}$), we can obtain the stresses $\sigma_{ab}$ at any point within the medium:

$$\sigma_{ab}(\mathbf{x}) = -c_{abcd} \int_{\text{fault}} c_{ijpq} \Delta u_i(\mathbf{y}) n_j(\mathbf{y}) \frac{\partial}{\partial x_q} \frac{\partial}{\partial x_d} G_{cp}(\mathbf{x}, \mathbf{y}) d\xi(\mathbf{y}). \quad (3)$$

It leads to an integral linking the stresses to the Green's function and the displacement discontinuity along the fault. Unfortunately the resulting integral is hypersingular for the on-fault stresses. This hypersingularity arises because of the second-order derivative of the Green's functions. As a result, conventional numerical integration techniques cannot be used (Koller et al. 1992; Tada & Yamashita 1997).

(iii) We regularize this integral using the tangential differential operator (Bonnet 1999; Sato et al. 2020) and project the slip vector onto the fault (Romanet et al. 2020). This introduces the curvature term and the gradient term into the equation:

$$\begin{aligned} \sigma_{ab}(\mathbf{x}) = & \int_{\text{fault}} \underbrace{c_{abcd} c_{ijpq} \frac{\partial}{\partial x_q} G_{cp}(\mathbf{x}, \mathbf{y}) [n_d(\mathbf{y}) t_j(\mathbf{y}) - n_j(\mathbf{y}) t_d(\mathbf{y})] t_i(\mathbf{y})}_{\text{Kernel}} \frac{\partial}{\partial y^t} \Delta u(\mathbf{y}) d\xi(\mathbf{y}) \\ & + \underbrace{\text{Gradient term}}_{\text{Kernel}} + \underbrace{\text{Curvature} \times \text{Slip}}_{\text{Curvature term}} \\ & + \int_{\text{fault}} \underbrace{c_{abcd} c_{ijpq} \frac{\partial}{\partial x_q} G_{cp}(\mathbf{x}, \mathbf{y}) [n_d(\mathbf{y}) t_j(\mathbf{y}) - n_j(\mathbf{y}) t_d(\mathbf{y})] n_i(\mathbf{y}) \kappa^t(\mathbf{y}) \Delta u(\mathbf{y}) d\xi(\mathbf{y}).}_{\text{Curvature term}} \end{aligned} \quad (4)$$

This equation is now Cauchy integrable, and **no longer** hypersingular. In the gradient term, $\frac{\partial}{\partial y^t} \Delta u(\mathbf{y})$ represents the derivative along the fault direction (see Appendix B1). In the curvature term, the local curvature of the fault $\kappa^t$ can be seen. The $t$ upper-script emphasizes that the curvature is the one in the direction of the tangential vector to the fault **t**.

(iv) **For on-fault shear ($\tau_{el} = t_i \sigma_{ij} n_j$) and normal ($\sigma_{el} = n_i \sigma_{ij} n_j$) tractions only**, we can develop the integrand that consists of Green's function and Hooke's tensor while making the approximation that fault geometry slope is small (see Fig. 2, Romanet et al. 2020; Romanet & Ozawa 2022). This approach allows to derive a simplified expansion for the **elastic tractions on the fault** due to a slip distribution:

$$\begin{aligned} \underbrace{\tau_{el}(\mathbf{x})}_{\text{Full solution}} &= \underbrace{\tau_{el}^0(x_1)}_{\text{0th order}} + \underbrace{\tau_{el}^1(x_1)}_{\text{1st order}} + \dots \\ \underbrace{\sigma_{el}(\mathbf{x})}_{\text{Full solution}} &= \underbrace{\sigma_{el}^0(x_1)}_{\text{0th order}} + \underbrace{\sigma_{el}^1(x_1)}_{\text{1st order}} + \dots \end{aligned} \quad (5)$$
The zeroth-order terms in the previous equation represent terms that are independent of the slope between any two points of the fault. The first-order terms are terms whose integrands are proportional to the fault slope $\propto \frac{x_2-y_2}{x_1-y_1}$, and the second-order terms are terms whose integrands are proportional to the square of the fault slope $\propto (\frac{x_2-y_2}{x_1-y_1})^2$ and so on. Note that assuming that the slope is small between any two points of the fault is equivalent to assume that the local slope $m(x_1) = \frac{d}{dy_1}h(y_1)$ is small, hence $m(x_1)$ is also a first-order term, and $m^2(x_1)$ a second-order term. Finally, please also note that assuming the fault slope to be small does not mean that the curvature of the fault is small. In this process, we do not make any assumption about the slip distribution, or the curvature along the fault. This is different from the previous studies using boundary perturbation theory (Chester & Chester 2000; Dunham et al. 2011; Fang & Dunham 2013), that were assuming the slope as well as the curvature to be small, and that required an expansion of the slip. This decomposition process may appear complex, but it offers two significant advantages when it comes to decomposing shear and normal traction:

(a) **Interpretability:** Breaking down the tractions into terms based on fault slope provides a clearer understanding of the physical processes at play. Each order of terms corresponds to a specific level of importance regarding the effect of fault geometry, allowing for easier interpretation and analysis.

(b) **Approximation:** The small slope approximation in the decomposition process enables us to simplify the mathematical expressions and calculations involved. This approximation is often valid in many practical scenarios, because faults are more or less linear structures. As a result, more computationally efficient numerical methods can be employed while still capturing essential aspects of fault behaviour (Romanet & Ozawa 2022).

The numerical calculation of the stress integrals (eq. 4) and the on-fault shear and normal tractions (eq. 5) is either done using a piecewise linear discretization of the fault geometry, together with a piecewise constant discretization of the slip (see Appendix C), or the spectral representation of the integrals (see Appendix G, H and I).

## 2.1 In-plane shear (mode II)

By employing the small slope approximation, the elastic shear and normal tractions due to slip on a non-planar fault can be expanded as:

$$ \underbrace{\tau_{\text{el}}}_{\text{Elastic shear traction of a non-planar fault}} = \underbrace{\tau_{\text{el}}^{0}}_{\text{Planar fault response}} + \underbrace{\tau_{\text{el}}^{1}}_{\text{Shear traction drag}} + \underbrace{\dots}_{\text{Higher order terms}}, \quad (6) $$

$$ \underbrace{\sigma_{\text{el}}}_{\substack{\text{Elastic normal traction} \\ \text{Only if the fault is non-planar}}} = \underbrace{\sigma_{\text{el}}^{0}}_{\substack{\text{Normal traction perturbation} \\ \quad \mathop{\text{\textit{c}}^{\text{i}}\text{t}}\Delta u^{\text{i}}}} + \underbrace{\sigma_{\text{el}}^{1}}_{\substack{\text{1st order} \\ \text{Normal traction perturbations}}} + \underbrace{\dots}_{\text{Higher order terms}}, \quad (7) $$

with the expressions:

$$ \tau_{\text{el}}^{0}(x_1) = -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u^{\text{i}}(y_1) \right] dy_1, \quad (8) $$

$$ \sigma_{\text{el}}^{0}(x_1) = -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{1}{x_1 - y_1} \kappa^{\text{i}}(y_1) \Delta u^{\text{i}}(y_1) \right] dy_1, \quad (8) $$

$$ \tau_{\text{el}}^{1}(x_1) = \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{m(y_1)}{x_1 - y_1} - \frac{x_2 - y_2}{(x_1 - y_1)^2} \right] \kappa^{\text{i}}(y_1) \Delta u^{\text{i}}(y_1) dy_1, \quad (9) $$

$$ \sigma_{\text{el}}^{1}(x_1) = \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{2m(x_1)}{x_1 - y_1} - \frac{m(y_1)}{x_1 - y_1} - \frac{x_2 - y_2}{(x_1 - y_1)^2} \right] \frac{d}{dy_1} \Delta u^{\text{i}}(y_1) dy_1, \quad (9) $$

where $\mu$ is the shear modulus and $\nu$ is the Poisson's ratio. Contrary to the general expression given in eq. (4), the integrals in this decomposition are over the line ($dy_1$) and not over the fault ($d\xi(y)$). An example of this decomposition of shear and normal traction into zeroth-order term and first-order term is illustrated in Fig. 4 for a rough fault and in Fig. 5 for a seamount fault geometry. The examples consider a given fault geometry with an assumed slip distribution $\Delta u(y_1) = (1 - 4\frac{y_1^2}{L^2})^{3/2}$, where $L$ is the length of the fault. We chose this slip distribution because it leads to an analytical expression for the shear traction, and it does not create any singularity at tips of the fault (see Appendix E). It can be observed that the zeroth- and first-order terms are already capturing a significant portion of the physics of the shear and normal traction for both geometries. The full solution, calculated without any approximation on the geometry, can always be decomposed as the sum of the gradient term, and the curvature term (Romanet et al. (2020), see also Appendix G2, H1 and I1). In the subsequent section, a physical interpretation of the zeroth- and first-order terms for both shear and normal traction is provided. Additionally, analytical results and possible other approximations will be presented.

### 2.1.1 The zeroth-order contribution to shear traction $\tau_{\text{el}}^0$: the planar fault response

The zeroth-order response of the shear traction, denoted as $\tau_{\text{el}}^0$, is exactly the same as if the fault was flat (Segall 2010, section 4.7). In other words, It means that the main contribution for the shear traction of a non-planar fault is the flat fault response. In Fig. 6, the shear traction
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Assuming a rough fault geometry (as shown by the x-axis) and the slip distribution $\Delta u(y_1) = (1 - 4\frac{y_1^2}{L^2})^{3/2}$ (red curve), $L$ being the length of the fault, the shear traction (upper panels) and the normal traction (lower panels) term can be decomposed into zeroth-order term and first-order term. The 0th order and the 1st order (thick blue lines), are respectively calculated using the expressions (8) and (9) (more precisely the spectral version of these equations as given by H3 and H5 in appendix H). The full solution (without any approximation on the fault geometry) and associated gradient and curvature terms (thin black lines) are calculated using eq. (H1).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Assuming a seamount fault geometry (as shown by the x-axis) and the slip distribution $\Delta u(y_1) = (1 - 4\frac{y_1^2}{L^2})^{3/2}$ (red curve), $L$ being the length of the fault, the shear traction (upper panels) and the normal traction (lower panels) term can be decomposed into zeroth-order term and first-order term. The 0th order and the 1st order (thick blue lines), are respectively calculated using the expressions (8) and (9) (more precisely the spectral version of these equations as given by H3 and H5 in appendix H). The full solution (without any approximation on the fault geometry) and associated gradient and curvature terms, in thin black lines, are calculated using eq. (H1).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Assuming a fault geometry (a. a rough fault, and b. a seamount fault geometry) and the slip distribution $\Delta u(y_1) = (1 - 4\frac{y_1^2}{L^2})^{3/2}$, L being the length of the fault, the shear traction can be calculated with the full solution as given by eq. (H1), and the zeroth-order solution as calculated by eq. (H3) in the spectral domain. The zeroth-order solution for the shear traction is exactly the one for a planar fault. It is compared with the analytical solution for the equivalent flat fault (black dash line, see Appendix E). The slight difference between the analytical solution and the zeroth-order solution are coming from the fact that the calculation was performed in spectral domain, hence there is a periodic replication of the fault that slightly increases the shear traction of the zeroth-order solution.</figcaption></figure>

response for a rough fault (Fig. 6a.) and for a seamount fault geometry (Fig. 6b.) that is subject to a slip distribution $\Delta u(y_1) = (1 - 4\frac{y_1^2}{L^2})^{3/2}$ in meter, L being the length of the fault, can be seen. This represents the shear traction for a right lateral fault. Inside the fault, the shear traction exhibits a globally negative value, resulting in decreasing the shear traction on the fault. This overall shear traction reduction is usually called the shear stress drop in seismology. On the other hand, outside the fault, the presence of the fault increases the shear traction.

### 2.1.2 The zeroth-order contribution to normal traction $\sigma_{el}^0$: the main source of normal traction variations

One of the main effects of non-planar geometry on an in-plane shear fault is the introduction of normal traction variations. This effect is well-documented in the literature (Nielsen & Knopoff 1998; Dunham et al. 2011; Romanet et al. 2020; Cattania & Segall 2021). Fig. 7 provides a closer look at the normal traction depicted in Fig. 4 for a rough fault and in Fig. 5 for a seamount fault geometry. It can be seen that the local maximums and minimums of normal traction correspond to the areas where the fault is locally flat ($\kappa^t = 0$). This relationship can be analytically demonstrated with the zeroth-order elastic normal traction (which is the main contribution to the normal traction) given by:

$$ \sigma_{el}^{0}(x_1) = -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \frac{1}{x_1 - y_1} \kappa^{t}(y_1) \Delta u^{t}(y_1) dy_1. \quad (10) $$

To identify the areas of minimums and maximums of normal traction, we can differentiate the previous expression and find where the derivative equals zero:

$$ \frac{d}{dx_1} \sigma_{el}^{0}(x_1) = \underbrace{\frac{\mu}{2\pi(1-\nu)}}_{>0} \int_{-\infty}^{+\infty} \underbrace{\frac{1}{(x_1 - y_1)^2}}_{>0} \kappa^{t}(y_1) \Delta u^{t}(y_1) dy_1. \quad (11) $$

Since the slip is always positive, the only way that the previous integral to be 0 is if the curvature reverses sign, or if the fault is only planar. Moreover, due to the weight $1/(x_1 - y_1)^2$, the maximums and minimums of normal traction occur very close to the areas where the curvature changes sign. Although it may not precisely coincide with the location of the curvature reversing sign, it generally occurs very close. Similar reasoning can be applied for the extremes of the shear traction for pure opening faults (see Appendix G).

Lastly, it is worth noting an interesting approximation for a rough fault that is not applicable to general fault geometries. If the slip distribution in Fourier domain has mainly low wavelength, so that it exists $k_c$ such that $\Delta u^t(k > k_c) = 0$, and that the curvature along the fault in Fourier domain has mainly high wavelength such that $\kappa^t(k < k_c) = 0$, then the Hilbert transform of the curvature that multiplies the slip can be simplified using Bedrosian's theorem (Bedrosian 1963):

$$ \sigma_{el}^{0}(x_1) \approx -\frac{\mu}{2\pi(1-\nu)} \Delta u^{t}(x_1) \int_{-\infty}^{+\infty} \frac{\kappa^{t}(y_1)}{x_1 - y_1} dy_1. \quad (12) $$

Although this approximation is not entirely accurate as the slip and curvature do have overlapping bandwidth, it can be considered a reasonable assumption for rough faults, as shown in Fig. 8. This approximation allows easy interpretation for the elastic normal traction change, because
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) Zoom on the normal traction for a rough fault in Fig. 4. (b) Zoom on the normal traction for a seamount fault geometry in Fig. 5. It can be seen that minimums and maximums of the normal traction correspond to area where the fault is locally flat (where the curvature is zero $κ' = 0$). The calculation was performed using the full elastic normal traction (eq. H1).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Assuming the fault geometry and the slip distribution $Δu(y_1) = (1 - 4\frac{y_1^2}{L^2})^{3/2}$, the full solution for normal traction (black curved) can be calculated using eq. (H1). The solution using the Bedrosian's approximation eq. (12) (blue dash curve) is very close to the full solution and allows to give some interpretations.</figcaption></figure>
Shear traction drag for different fault geometries

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Calculation of the first-order effect of the shear traction as expressed by eq. (13) (continuous red line), and comparison with the theoretical shear traction drag as calculated by eq. (16) (black crosses) for different fault geometries. The fault geometry is shown as the x-axis. (a) A finite rough fault, (b) a finite fault with a seamount geometry, (c) a finite fault with arctan geometry, and (d) a finite fault whose geometry is the sum of two sinusoidal functions. For the sum of two sinusoidal functions, the shear traction drag is not shown because the equation is too lengthy. The calculation is done for the prescribed slip distribution $\Delta u(y_1) = (1 - 4\frac{y_1^2}{L^2})^{3/2}$ in meter, where L is the length of the fault. The constant slip required for the calculation of the shear traction drag is taken as the average slip.</figcaption></figure>

it means that the normal traction depends only on the local slip ($\Delta u'(x_1)$) and not on the global slip inside an integral.

2.1.3 The first-order contribution to shear traction $\tau_{el}^1$: the shear traction drag

Fig. 9 illustrates the first-order term for shear traction:

$$ \tau_{el}^{1}(x_1) = \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{m(y_1)}{x_1 - y_1} - \frac{x_2 - y_2}{(x_1 - y_1)^2} \right] \kappa'(y_1) \Delta u'(y_1) dy_1, \quad (13) $$

for different fault geometries when assuming the slip distribution $\Delta u(y_1) = (1 - 4\frac{y_1^2}{L^2})^{3/2}$. This figure shows that additional shear resistance is coming from the first-order effect of the shear traction for any fault geometry (Fig. 9). It can be seen that as long as the fault is non-planar, this term is resisting to movement. It leads to the very intuitive result that it is harder to slide a fault if it is non-planar.

This additional shear resistance is very similar to the key theoretical result on rough fault obtained by Dieterich & Smith (2009) and Fang & Dunham (2013). It states that a rough fault is harder to slip than a flat fault. In other words, it says that a rough fault has an additional shear resistance when compared to a flat fault. A quantification of this term was obtained in Fang & Dunham (2013) and referred as the “roughness drag”. It is possible to derive a similar expression from the 1st order solution for shear traction. Assuming constant slip over the fault and that the fault slope is small ($m \ll 1$), the previous exact equation for the shear traction at 1st order (eq. 13) simplifies to (by writing $h(y_1) = y_2$ and $h(x_1) = x_2$):

$$ \tau_{el}^{1}(x_1) = \frac{\mu \Delta u'}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{m(y_1)}{x_1 - y_1} - \frac{x_2 - y_2}{(x_1 - y_1)^2} \right] \frac{d^2}{dy_1^2} m(y_1) dy_1. \quad (14) $$
To obtain the average $\langle \tau_{el}^1 \rangle$ of the previous equation, it is possible to Fourier transform it, and then apply it at the null wavelength ($k = 0$):

$$ \begin{aligned} \langle \tau_{el}^1 \rangle &= \lim_{T \to +\infty} \left[ \frac{1}{T} \int_{-T/2}^{+T/2} \tau_{el}^1(y_1) dy_1 \right] = \lim_{T \to +\infty} \left[ \frac{1}{T} \int_{-\infty}^{+\infty} \tau_{el}^1(y_1) \chi_T(y_1) e^0 dy_1 \right] \\ &= \lim_{T \to +\infty} \left[ \frac{1}{T} F[\tau_{el}^1 \chi_T](k=0) \right] \\ &= -\frac{\mu \Delta u^t}{4\pi(1-\nu)} \lim_{T \to +\infty} \left[ \frac{1}{T} \int_{-\infty}^{+\infty} H_T(-\tau) H_T(\tau) |\tau|^3 d\tau \right] \\ &= -\frac{\mu \Delta u^t}{4\pi(1-\nu)} \lim_{T \to +\infty} \left[ \frac{1}{T} \int_{-\infty}^{+\infty} H_T^*(\tau) H_T(\tau) |\tau|^3 d\tau \right] \\ &= -\frac{\mu \Delta u^t}{2\pi(1-\nu)} \int_0^{+\infty} P_h(\tau) |\tau|^3 d\tau, \end{aligned} \quad (15) $$

where $\langle \tau_{el}^1 \rangle$ denotes the average of $\tau_{el}^1$, $H_T$ is the Fourier transform of the height of the fault on the interval $[-T/2, T/2]$, and 0 elsewhere ($F[h_T] = H_T$), * denotes the complex conjugate. Finally, $\chi_T(y_1) = [H(y_1 + T/2) - H(y_1 - T/2)]$ is the rectangular function and $P_h = \lim_{T \to \infty} \frac{1}{T} H_T^* H_T$ is the power spectra density of the height of the fault. The above expression allows to calculate a theoretical estimate of the shear traction drag ($\tau_{drag}^1$) for any fault geometry:

$$ \tau_{drag}^1 = -\frac{\mu \Delta u^t}{2\pi(1-\nu)} \int_0^{+\infty} P_h(\tau) |\tau|^3 d\tau. \quad (16) $$

In this paper, we choose to name $\tau_{drag}^1$ the “shear traction drag”, which is the average of the first-order response of the elastic shear traction ($\tau_{el}^1$) when the slip is constant and the fault is “sufficiently” long. The “shear traction drag” ($\tau_{drag}^1$) shares a lot of similarities with the “roughness drag” of Fang & Dunham (2013), but contrary to the latter expression, it is really a traction along the fault (the expression in Fang & Dunham (2013) was integrated along the x-axis), and it is independent of the friction that applies on the fault.

For a rough self-similar fault, the power spectra density is given by $P_h(k) = (2\pi)^3 \alpha^2 |k|^{-3}$, which leads to the same expression as the roughness drag calculated in Fang & Dunham (2013):

$$ \tau_{drag}^1 = -8\pi^3 \alpha^2 \frac{\mu}{1-\nu} \frac{\Delta u^t}{\lambda_{min}}, \quad (17) $$

where $\alpha$ represents the amplitude to wavelength ratio, and $\lambda_{min}$ is the minimum roughness wavelength.

Here are some properties of the first-order term for shear traction:

(i) The first-order term for shear traction is valid for any geometry and any slip distribution. There is more resistance to slip for a non-planar fault compared to a planar fault (for in-plane shear -mode II-).

(ii) It is only associated with in-plane faulting (mode II), there is no additional shear resistance for out-of-plane shearing (mode III).

(iii) Since the stresses in the Fourier domain in three dimensional can be understood as a combination of mode II and mode III (Geubelle & Rice 1995), the analytical result mentioned above can be used to calculate the first-order term for shear traction in three dimensions.

(iv) Since the curvature term for shear traction as given by eq. (4), and the first-order term (eq. 9) are the same at order 1, the two terms can both be used to calculate the shear traction resistance.

(v) An interpretation is that for equivalent slip distribution, the shear stress drop will be higher for a non-planar fault than for a planar fault. This interpretation may seem to contradict the fact that the fault is more difficult to slip, but this is not the case. To get an equivalent slip, the loading traction on a non-planar fault will be required to be higher than on a planar fault, hence it is more difficult to slip on a non-planar fault.

(vi) The shear traction drag can be used to obtain the same result as the roughness drag expression obtained by Fang & Dunham (2013) on a rough fault: $\tau_{drag}^1 = -8\pi^3 \alpha^2 \frac{\mu}{1-\nu} \frac{\Delta u^t}{\lambda_{min}}$. The previous result was derived for a constant slip on an infinite frictionless fault,

To further validate our derivations, we conducted numerical tests comparing the elastic shear traction at first-order $\tau_{el}^1$ and the shear traction drag $\tau_{drag}^1$ for an infinite non-planar fault with constant slip. To achieve the infinite fault, we calculated the shear traction drag in the spectral domain (eq. H5), hence it is infinite with periodic replication. The results can be seen in Fig. 10(a), for a rough fault and in Fig. 10(b) for a sinusoidal fault. The mean value of the first-order term for shear traction ($\langle \tau_{el}^1 \rangle$, in dashed red line) aligns perfectly with the shear traction drag ($\tau_{drag}^1$, black crosses) for both fault geometries. This numerical comparison further supports the validity and accuracy of our derived expressions.

### 2.1.4 The first-order contribution to normal traction $\sigma_{el}^1$: slightly reduce maximum and minimums of normal traction

The main effect of the first order on normal traction $\sigma_{el}^1$ is to slightly oppose the effect of the 0th order, so that the full normal traction at maximum and minimums are slightly reduced compare to the zeroth-order effect. This effect can be seen on Fig. 11, where the maximums and minimums of normal traction at zeroth-order $\sigma_{el}^0$ are slightly reduced by the first order of normal traction $\sigma_{el}^1$ for both a rough fault (Fig. 11a.) and a seamount fault geometry (Fig. 11b.).
<figure><img src="image_12.png" /><figcaption>**Figure 10.** (a) The first-order solution for shear traction $\tau_{el}^1$ -red line- (from eq. 13) and the shear traction drag from Fang & Dunham (2013) -black crosses- (from eq. 17) for an infinite rough fault with constant slip. (b) The first-order solution for shear traction $\tau_{el}^1$ -red line- (from eq. 13) and the shear traction drag -black crosses- (eq. 16) for an infinite sinoidical fault with constant slip. The first-order solution for shear traction $\tau_{el}^1$ for both fault geometries was calculated in the Fourier domain (eq. H5, in appendix H) to make it infinite with periodic replication. It can be seen that the average of shear traction at first-order $\langle \tau_{el}^1 \rangle$ and the shear traction drag $\tau_{drag}^1$ are perfectly overlapping.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** The first-order effect on normal traction for (a) a rough fault and (b) a seamount fault geometry. If can seen that the main effect is to reduce the minimums and maximums of the normal traction. The full solution was calculated with the eq. (H1) in Appendix H, and the zeroth and first order were respectively calculated by eqs. (H3) and (H5), in Appendix H.</figcaption></figure>

### 2.1.5 Scaling of slip versus geometry

First and foremost, please note that in this section, we are introducing a new approach by solving the mechanical problem using the equilibrium of shear traction on the fault, without assuming a specific slip distribution along the fault. If the fault follows the Coulomb friction law ($f$ being the friction coefficient) and that the loading can be considered constant at 0th order along the fault, we can write the equilibrium of shear traction on the fault at 0th order as:

$$ \underbrace{-f(\sigma_{el}^0 + \sigma_{load}^0)}_{\text{Friction}} = \underbrace{\tau_{el}^0}_{\text{Shear traction due to elasticity}} + \underbrace{\tau_{load}^0}_{\text{Loading}}. \qquad (18) $$
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Scaling of geometry and slip distribution for a static, in-plane shear (mode II) fault subject to a constant loading and that follows Coulomb friction. The black dash line is the planar fault solution. The solution for a non-planar fault geometry is shown as a blue line: (a) a rough fault, (b) a seamount fault geometry, (c) an arctan fault geometry and (d) a sum of two sinusoidal functions geometry. The slip distribution is inverted using the equilibrium of shear traction on the fault. The full elastic traction was used to obtain the slip (eq. H1).</figcaption></figure>

By utilizing the equilibrium equation mentioned above, it becomes possible to solve for the distribution of slip along the fault.

This approach was applied for the different non-planar fault geometries in Fig. 12. The figure illustrates a clear scaling relationship between the geometry of the fault and the slip distribution: locally flat areas exhibit extrema of slip while highly curved areas show a strong slip gradient. It can be seen that this scaling is not working near the edges of the fault. Remarkably, the scaling appears to be linear and preserves the wavelength of the fault geometry. Indeed, in Fig. 12(d), the fault geometry is a sum of two sinusoids with different amplitude and wavelength. The amplitude and wavelength of the fault geometry can be retrieved in the resulting slip distribution.

This scaling has been analytically derived in Romanet et al. (2020), and is presented again in this paper for consistency. This scaling arises from the relationship between the shear traction and the normal traction through the friction law. For an in-plane shear fault, as shown in the previous section, the shear traction at 0th, $\tau_{el}^0$, is controlled by the gradient term along the fault while the normal traction at 0th order, $\sigma_{el}^0$, is controlled by the curvature term:

$$ 
\begin{aligned}
\tau_{el}^{0}(x_1) &= -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u'(y_1) \right] dy_1, \\
\sigma_{el}^{0}(x_1) &= -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{1}{x_1 - y_1} \kappa'(y_1) \Delta u'(y_1) \right] dy_1.
\end{aligned}
\qquad (19) $$

In a global stress state as defined as:

$$ \vec{\sigma} = \begin{bmatrix} \sigma_{11} & \sigma_{12} \\ \sigma_{12} & \sigma_{22} \end{bmatrix}, \qquad (20) $$

the shear and normal traction loadings projected onto the non-planar fault are respectively $\tau_{load} = \mathbf{t} \cdot \vec{\sigma} \cdot \mathbf{n}$ and $\sigma_{load} = \mathbf{n} \cdot \vec{\sigma} \cdot \mathbf{n}$. At 0th order, it can be shown that the shear and normal traction loads simply reduce to $\tau_{load}^0 = \sigma_{12}$ and $\sigma_{load}^0 = \sigma_{22}$, which are constant and independent of the fault geometry. The Coulomb friction on the fault links the curvature term (the fault geometry) that controls normal traction and the gradient term (the slip distribution gradient) that controls shear traction, and because at 0th order, the loading can be considered constant, it makes possible to invert analytically the two Hilbert's transforms, by using Chebyshev polynomials (Segall 2010, section 4.1), to obtain a
linear ordinary differential equation:

$$ \underbrace{\frac{d}{dx_1} \Delta u'(x_1) + f \kappa'(x_1) \Delta u'(x_1)}_{\text{Gradient of slip}} - \underbrace{\frac{2(1-\nu)}{\mu} \frac{x_1 - L/2}{L \sqrt{1 - \frac{4(x_1 - L/2)^2}{L^2}}} (f \sigma^{\text{load}} + \tau^{\text{load}})}_{\text{Effect of geometry}} \underbrace{\overbrace{\text{Effect of loading}}^{\text{Effect of geometry}}} \quad (21) $$

In the centre of fault, the effect of loading cancels so that we obtained the scaling:

$$ \frac{d\Delta u'}{\Delta u'} = -f dm, \quad (22) $$

where the property $\kappa' dx = dm$ has been used. This scaling states that the relative variation of slip $\frac{d\Delta u'}{\Delta u'}$ is equal to minus the friction coefficient $f$ that multiplies the variation of the slope of the fault $dm$.

This scaling, eq. (22), exists only for in-plane shear fault (mode II) as shown in Fig. 12, but not for out-of-plane shear fault (mode III) as depicted in Fig. 14. Although his scaling has not been yet confirmed in observations (Bruhat et al. 2020), it was verified for fully dynamic simulations using rate and state friction (Romanet & Ozawa 2022).

This scaling is an important theoretical finding as it provides a link between three crucial parameters of the fault, i.e. the fault geometry, the slip distribution and the friction coefficient.

## 2.2 Out-of-plane shear (mode III)

Similarly as for the in-plane (mode II) case, the out-of-plane (mode III) shear and normal traction can be decomposed into zeroth- and first-order responses:

$$ \underbrace{\tau_{\text{el}}}_{\substack{\text{Elastic shear traction} \\ \text{of non-planar fault}}} = \underbrace{\tau_{\text{el}}^{0}}_{\substack{\text{Planar fault response}}} + \underbrace{\tau_{\text{el}}^{1}}_{\substack{=0}} + \underbrace{\dots}_{\substack{\text{Higher order terms}}} \quad (23) $$

$$ \underbrace{\sigma_{\text{el}}}_{\substack{\text{Elastic normal traction} \\ =0}} = \underbrace{\sigma_{\text{el}}^{0}}_{\substack{=0}} + \underbrace{\sigma_{\text{el}}^{1}}_{\substack{=0}} + \underbrace{\dots}_{\substack{\text{Higher order terms} \\ =0}} \quad (24) $$

However, in this case, only the shear traction at 0th order is non-zeros:

$$ \tau_{\text{el}}^{0}(x_1) = \frac{\mu}{2\pi} \int_{-\infty}^{+\infty} \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u^s(y_1) dy_1, \quad (25) $$

$$ \sigma_{\text{el}}^{0}(x_1) = 0, \quad (26) $$

$$ \tau_{\text{el}}^{1}(x_1) = 0, \quad (26) $$

$$ \sigma_{\text{el}}^{1}(x_1) = 0. \quad (26) $$

For an out-of-plane shear fault, the contribution of non-planar fault geometry to the stresses is only a second-order effect. This means that fault geometry has a small effect for an out-of-plane shear fault (mode III), only on the shear traction (see Fig. 13). There is no normal traction variation for an out-of-plane shear fault (mode III). Fig. 14 illustrates the slip distribution for non-planar fault geometries with Coulomb friction and constant loading, using the same parameters as in the in-plane shear case (mode II) shown in Fig. 12. In the case of out-of-plane shear (mode III), there is no scaling of the slip distribution with respect to fault geometry. The slip distribution for the planar fault (dashed black line) perfectly overlaps with the slip distributions for the non-planar faults (continuous blue lines).

# 3 DISCUSSION

## 3.1 Physical singularities that appear on non-planar fault

### 3.1.1 Non-zeros slip at kinks of the fault

Many literature sources make the mistake of modelling kinks while assuming infinitesimal strain theory and no opening (Tada & Yamashita 1996; Aochi et al. 2000a; Duan & Oglesby 2005; Ely et al. 2009; Lozos et al. 2011; Fukuyama & Hok 2015; Sathiakumar & Barbot 2021). A kink on a fault with non-zeros slip and no opening create a $1/r$ singularity, where $r$ is the distance from the kink. This is a non-physical singularity that results to infinite strain energy near the kink. We can demonstrate this result in the small slope approximation, for the 0th order on normal traction:

$$ \sigma_{\text{el}}^{0}(x_1) = -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{1}{x_1 - y_1} \Delta u'(y_1) \frac{d}{dy_1} m(y_1) \right] dy_1. \quad (27) $$
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Total elastic traction for an out-of-plane fault (mode III). Assuming the fault geometry (as shown by the x-axis) and the slip distribution $\Delta u(y_1) = (1 - 4\frac{y_1^2}{L^2})^{3/2}$ in meter (red curve), the full solution, the gradient term and the curvature term (black lines) can be calculated using eq. (11) (more precisely using the spectral expressions eqs. 13 and 15). The 0th order and the 1st order (blue lines), are respectively calculated using the expressions (25) and (26).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Scaling of geometry and slip distribution for a static, out-of-plane shear (mode III) fault subject to a constant loading and that follows Coulomb friction. The black dash line is the planar fault solution. The slip distribution does not show any scaling with fault geometry. Fault geometry: (a) a rough fault, (b) a seamount fault geometry, (c) an arctan fault geometry and (d) a sum of two sinusoidal functions geometry.</figcaption></figure>

If we assume that the curvature is $\kappa'(y_1) = \frac{d}{dy_1}m(y_1)$. For a kink at the position $y_1 = 0$, the slope is discontinuous so that we can assume that $m(y_1) = AH(y_1)$, where H is the Heaviside function and A the fault slope after the kink. In this case, the previous equation becomes:

$$ 
\begin{aligned}
\sigma_{el}^{0}(x_1) &= -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{A\Delta u(y_1)}{x_1 - y_1} \delta(y_1) \right] dy_1 \\
&= -\frac{\mu}{2\pi(1-\nu)} \frac{A\Delta u'(0)}{x_1}.
\end{aligned} 
\quad (28) $$
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Figure showing the grid dependence of a kinked faults. (a) The geometry of the fault and the momentum equation. (b) The slip distribution versus position along the fault. The slip distribution is changing with gridsize. Please note the only physically acceptable slip at the kink is zero. (c) The shear traction distribution versus the position along the fault. The shear traction is closely following the normal traction because they are linked by the friction law. It is also diverging with the refinement of the mesh. (d) The normal traction distribution versus position along the fault. The normal traction is diverging with the refinement of the mesh.</figcaption></figure>

This reveals a $1/x_1$ singularity of normal traction on the fault. One straightforward way to eliminate this singularity is by imposing the slip to be null at the kink $\Delta u'(0) = 0$, which is evident from the previous equation. Another approach would be to allow for fault opening or adding a third fault and imposing the closure of slip at the kink (Andrews 1989). In the boundary element method, the maximum normal traction concentration due to the discretization is typically proportional to the gridsize $\Delta s$. Thus, the maximum normal traction on the fault is:

$$ \sigma_{el}^{0}(x_1 = 0) \propto -\frac{\mu}{2\pi(1-\nu)} \frac{A\Delta u'(0)}{\Delta s}. \qquad (29) $$

The problem is that the usual gridsize for the modelling of fault is of the order $\Delta s \sim 1000$ m, as a result, the stress concentration at the kink, which is normally infinite, is of the order of $\sigma_{el}^{0}(x_1) \simeq 1$ MPa for typical shear modulus $\mu$ (30 GPa), P-wave velocity (5000 m s$^{-1}$) and S-wave velocity (3000 m s$^{-1}$) and a slip at the kink of $\Delta u'(0) = 1$ m. This is why the issue of diverging traction at kinks has remained mostly unknown in numerical work so far. An example of this effect is shown in Fig. 15, where the slip, normal traction and shear traction was calculated for a fault following Coulomb friction with constant loading. The maximum values of normal traction and shear traction are largely underestimate (~1 MPa) for the gridsize $\Delta s = 500$ m. It is worth noting that adding plasticity or viscosity is not an entirely satisfying solution, because in most of seismological research, the infinitesimal strain theory is used together with plasticity/viscosity and the plastic/viscous effects are driven by the linear elastic interaction. If the elastic interactions are underestimated, the plastic/viscous effect will also be significantly underestimated. It means all the quantitative work done using kinks in linear elasticity (infinitesimal strain theory) is grid-dependent and does not converge numerically.

Finally, it is still possible to observe that modelling a kink is still possible in out-of-plane faulting (mode III), because there is no curvature term.

### 3.1.2 Discontinuous curvature along the fault

A discontinuity in the curvature of the fault also creates a singularity in stresses. However, this singularity can be compensated by a singularity in the slip derivative (hence the slip will still be continuous). Using the scaling of slip versus curvature—that is valid only when the fault follows Coulomb friction—(eq. 22):

$$ \frac{d}{dx_1} \Delta u'(x_1) = -f \kappa'(x_1) \Delta u'(x_1). \qquad (30) $$
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Figure showing the grid dependence of a bended fault with discontinuous curvature. (a) The geometry of the fault and the momentum equation. (b) The slip distribution versus position along the fault. The slip distribution is changing with gridsize. (c) The shear traction distribution versus the position along the fault. The shear traction is closely following the normal traction because they are linked by the friction law. It is also diverging with the refinement of the mesh. (d) The normal traction distribution versus position along the fault. The normal traction is diverging with the refinement of the mesh.</figcaption></figure>

So that if the curvature has a discontinuity of Amplitude $B$ at $x_1 = 0$, $\kappa'(x_1) = BH(x_1)$, the discontinuity in slip gradient will be: $\frac{d}{dx_1} \Delta u'(x_1) = -fBH(x_1)\Delta u'(x_1)$. This effect is well observed on Fig. 16(b), where there is a discontinuity in slip gradient due to the discontinuity of the curvature along the fault. The discontinuity in curvature creates a singularity of the stresses that is $\propto \log(x_1)$ when $x_1 \to 0$. Because this singularity is less strong than $1/\sqrt{(x_1)}$, this discontinuity in curvature does not lead to infinite strain energy, hence it is a physically acceptable singularity. Let's show the previous result using the small slope approximation and assuming constant slip distribution. For a bended fault like in Fig. 16(a), the curvature can be written $\kappa'(y_1) = B(H(y_1 - a) - H(y_1 + b))$, where $B$ is the curvature in the bended portion of the fault. It leads to the normal traction distribution (Hilbert transform of a characteristic function):

$$ \sigma_{el}^{0}(x_1) = -\frac{\mu}{2\pi(1-\nu)} \Delta u' \int_{-\infty}^{+\infty} \left[ \frac{1}{x_1 - y_1} B(H(y_1 - a) - H(y_1 + b)) \right] dy_1 \\ = -\frac{\mu}{2\pi(1-\nu)} B \Delta u' \log \left| \frac{x_1 - a}{x_1 - b} \right| \quad (31) $$

Where we found a log singularity when $x_1 \to a$ or $x_1 \to b$. The last comment is that, again, there is no such a problem in out-of-plane shear (mode III), because the stresses and strains are independent of the curvature term (see eq. I1 in appendix I).

## 3.2 Limits of infinitesimal strain theory

The infinitesimal strain theory is commonly used in seismology and is necessary to linearize the strain tensor with the displacement field. This theory assumes that both displacement and strains are small. However, when considering non-planar fault geometry, the strains and stresses will keep increasing with on-going slip breaking the small strain approximation. At a certain point, the small strain approximation becomes invalid as depicted in Fig. 17. In such cases, finite elasticity should be considered to prevent the strains and stresses from growing without bounds (Romanet et al. 2020; Tal 2023). Various approaches have been proposed in the literature to address the issue of unbounded increase or decrease in normal traction with ongoing slip (Duan & Oglesby 2005; Dunham et al. 2011; Heimisson 2020; Cattania & Segall 2021), for example the inclusion of viscosity or plasticity while retaining the small strain approximation of infinitesimal strain theory. However, this approach may lead to a significant overestimation of the effect of plasticity and viscosity in fault mechanics.

While we acknowledge the likelihood of nonlinear anelastic phenomena occurring (such as damage, plasticity or viscosity), we disagree with the interpretation that they are a necessary condition for preventing the stresses and strains from growing indefinitely with on-going slip on non-planar faults (Dieterich & Smith 2009; Shi & Day 2013), as this does not consider the effect of finite elasticity (Wise & Tal 2024).
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Figure showing the limits of infinitesimal strain theory (linear elasticity) in the modelling of non-planar fault geometries.</figcaption></figure>

<table><caption>**Table 1.** Summary of the results. Please note that there is an effect of geometry for out-of-plane shear, but this effect only appears at 2nd order.</caption><thead><tr><th></th><th>Opening fault (mode I)</th><th>In-plane shear fault (mode II)</th><th>Out-of-plane shear fault (mode III)</th></tr></thead><tbody><tr><td>Shear traction τ<sub>el</sub><sup>0</sup> (0<sup>th</sup> order)</td><td>Mainly controlled by the fault geometry ∝ κ<sup>l</sup> Δu<sup>h</sup></td><td>Same as the planar fault response</td><td>Same as the planar fault response</td></tr><tr><td>Shear traction τ<sub>el</sub><sup>1</sup> (1<sup>st</sup> order)</td><td>Nearly no effect (depends on the gradient of slip)</td><td>Shear traction drag, the fault is resisting slip</td><td>= 0 No effect of fault geometry</td></tr><tr><td>Normal traction σ<sub>el</sub><sup>0</sup> (0<sup>th</sup> order)</td><td>Same as the planar fault response</td><td>Mainly controlled by the fault geometry ∝ κ<sup>l</sup> Δu<sup>l</sup></td><td>= 0 No effect of fault geometry</td></tr><tr><td>Normal traction σ<sub>el</sub><sup>1</sup> (1<sup>st</sup> order)</td><td>Normal traction drag, the fault is resisting opening</td><td>Reduce the minimums and maximums of normal traction changes</td><td>= 0 No effect of fault geometry</td></tr></tbody></table>

It is important to recognize that small strain elasticity is violated when dealing with very small asperities, and even regularization techniques may not provide accurate results. One possible direction is to move towards finite elasticity and thoroughly test the limits of linear elasticity in current earthquake simulations on non-planar faults.

## 4 CONCLUSION

This work advances the limited theoretical knowledge regarding the effect of non-planar fault geometries on earthquake mechanics. By expanding the relation of fault traction and slip, up to first order relative to the deviation from a planar fault, this study allows for the interpretation of complex fault geometry and its impact on fault traction (see Table 1).

The results of this study confirm that fault geometry plays a significant role in in-plane faulting (mode II) by modifying the normal traction on the fault and increasing its resistance to slipping. We provide a general expression independent of fault geometry and fault slip for the shear traction drag (Fang & Dunham 2013), making it a general result. We also provide some useful simplification for the effect of rough fault on normal traction, which can be simplified as an Hilbert transform of the curvature by using Bedrosian's theorem.

Conversely, for out-of-plane faulting (mode III), the influence of fault geometry is negligible. There is no effect of fault geometry up to the second order. We also showed that, in this case, there is no scaling between the slip distribution and the geometry.

The paper also examines singularities that arise in specific fault geometries commonly used in earthquake simulations and provides guidelines for their elimination.

Ultimately, this study highlights the limitations of the small strain approximation when considering non-planar faults, emphasizing the need to consider finite elasticity for more accurate modelling of non-planar faults.

## ACKNOWLEDGMENTS

We would like to thank Robert Viesca for providing the non-singular analytical solution for stresses and slip. PR would also like to acknowledge the strong support by Jean-Paul Ampuero, Frédéric Cappa and Marco M. Scuderi. This research was mainly funded by the National Research
Institute for Earth Science and Disaster Resilience (NIED) project on the generation mechanism of large earthquakes and partly from the European Research Council (ERC) Starting Grants 101040600 (HYQUAKE). We would like to thank three anonymous reviewers, and Eric Dunham, whose reviews contributed a lot to improve this paper.

## SUPPORTING INFORMATION

Supplementary data are available at *GJIRAS* online.

mathematica.zip

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## DATA AVAILABILITY

The mathematica code used to develop the equation is available as supplementary material. No data was used in this research.

## REFERENCES

* Aki, K., 1979. Characterization of barriers on an earthquake fault, *J. Geophys. Res.*, **84**(B11), 6140–6148.

* Aki, K. & Richards, P.G., 2002. *Quantitative Seismology*, University Science Books.

* Andrews, D., 1989. Mechanics of fault junctions, *J. Geophys. Res.* **94**(B7), 9389–9397.

* Aochi, H., Fukuyama, E. & Matsu'ura, M., 2000a. Spontaneous rupture propagation on a non-planar fault in 3d elastic medium, *Pure Appl. Geophys.*, **157**, 2003–2027.

* Bedrosian, E., 1963. A product theorem for hilbert transforms, *Proc. IEEE*, **51**(5), 868–869.

* Bonnet, M., 1999. *Boundary Integral Equation Methods for Solids and Fluids*, Vol. 34, pp. 301–302, Springer.

* Bruhat, L., Klinger, Y., Vallage, A. & Dunham, E.M., 2020. Influence of fault roughness on surface displacement: from numerical simulations to coseismic slip distributions, *Geophys. J. Int.*, **220**(3), 1857–1877.

* Cattania, C. & Segall, P., 2021. Precursory slow slip and foreshocks on rough faults, *J. Geophys. Res.* **126**(4), e2020JB020430.

* Chester, F.M. & Chester, J.S., 2000. Stress and deformation along wavy frictional faults, *J. Geophys. Res.*, **105**(B10), 23 421–23 430.

* Dieterich, J.H. & Smith, D.E., 2009. Nonplanar faults: Mechanics of slip and off-fault damage, in *Mechanics, Structure and Evolution of Fault Zones*, pp. 1799–1815, eds Ben-Zion, Y. & Sammis, C., Springer. https://doi.org/10.1007/s00024-009-0517-y.

* Duan, B. & Oglesby, D.D., 2005. Multicycle dynamics of nonplanar strike-slip faults, *J. Geophys. Res.*, **110**, https://doi.org/10.1029/2004JB003298.

* Dunham, E.M., Belanger, D., Cong, L. & Kozdon, J.E., 2011. Earthquake ruptures with strongly rate-weakening friction and off-fault plasticity, part 2: Nonplanar faults, *Bull. Seism. Soc. Am.*, **101**(5), 2308–2322.

* Ely, G.P., Day, S.M. & Minster, J.B., 2009. A support-operator method for 3-d rupture dynamics, *Geophys. J. Int.*, **177**(3), 1140–1150.

* Fang, Z. & Dunham, E.M., 2013. Additional shear resistance from fault roughness and stress levels on geometrically complex faults, *J. Geophys. Res.: Solid Earth*, **118**(7), 3642–3654.

* Fukuyama, E. & Hok, S., 2015. Dynamic overshoot near trench caused by large asperity break at depth, *Pure Appl. Geophys.*, **172**(8), 2157–2165.

* Geubelle, P.H. & Rice, J.R., 1995. A spectral method for three-dimensional elastodynamic fracture problems, *J. Mech. Phys. Solids*, **43**(11), 1791–1824.

* Heimisson, E.R., 2020. Crack to pulse transition and magnitude statistics during earthquake cycles on a self-similar rough fault, *Earth Planet. Sci. Lett.*, **537**, 116–202.

* King, G.C.P. & Nabelek, J., 1985. The role of bends in faults in the initiation and termination of earthquake rupture, *Science*, **228**, 984–987.

* Klinger, Y., Michel, R. & King, G.C.P., 2006. Evidence for an earthquake barrier model from mw ~ 7.8 kokoxili (tibet) earthquake slip-distribution, *Earth Planet. Sci. Lett.*, **242**, 354–364.

* Koller, M.G., Bonnet, M. & Madariaga, R., 1992. Modelling of dynamical crack propagation using time-domain boundary integral equations, *Wave Motion*, **16**(4), 339–366.

* Lozos, J.C., Oglesby, D.D., Duan, B. & Wesnousky, S.G., 2011. The effects of double fault bends on rupture propagation: A geometrical parameter study, *Bull. Seism. Soc. Am.*, **101**(1), 385–398.

* Mason, J.C. & Handscomb, D.C., 2002. *Chebyshev Polynomials*, Chapman and Hall/CRC.

* Maurer, J., 2024. Statistical distribution of static stress resolved onto geometrically-rough faults, *Seismica*, **3**(2), https://doi.org/10.26443/seismica.v3i2.1206.

* Milliner, C.W., Dolan, J.F., Hollingsworth, J., Leprince, S., Ayoub, F. & Sammis, C.G., 2015. Quantifying near-field and off-fault deformation patterns of the 1992 mw 7.3 l anders earthquake, *Geochem. Geophys. Geosyst.*, **16**(5), 1577–1598.

* Morad, D., Lyakhovsky, V., Hatzor, Y.H. & Sagy, A., 2022. Stress heterogeneity and the onset of faulting along geometrically irregular faults, *Geophys. Res. Lett.*, **49**(17), e2021GL097,591.

* Nielsen, S.B. & Knopoff, L., 1998. The equivalent strength of geometrical barriers to earthquakes, *J. Geophys. Res.*, **103**(B5), 9953–9965.

* Romanet, P. & Ozawa, S., 2022. Fully dynamic earthquake cycle simulations on a nonplanar fault using the spectral boundary integral element method (sbiem), *Bull. Seism. Soc. Am.*, **112**(1) 78–97.

* Romanet, P., Sato, D.S. & Ando, R., 2020. Curvature, a mechanical link between the geometrical complexities of a fault: application to bends, kinks and rough faults, *Geophys. J. Int.*, **223**, 211–232.

* Sagy, A. & Lyakhovsky, V., 2019. Stress patterns and failure around rough interlocked fault surface, *J. Geophys. Res.*, **124**(7), 7138–7154.

* Sathiakumar, S. & Barbot, S., 2021. The stop-start control of seismicity by fault bends along the main himalayan thrust, *Commun. Earth Environ.*, **2**(1), 1–11.

* Sato, D.S., Romanet, P. & Ando, R., 2020. Paradox of modelling curved faults revisited with general non-hypersingular stress green's functions, *Geophys. J. Int.*, **223**(1), 197–210.

* Saucier, F., Humphreys, E. & Weldon, R., 1992. Stress near geometrically complex strike-slip faults: Application to the san andreas fault at cajon pass, southern california, *J. Geophys. Res.*, **97**(B4), 5081–5094.

* Schwartz, D.P. & Sibson, R.H., 1989. Fault segmentation and controls of rupture initiation and termination. *Open-File Rep, U.S. Geol. Surv.*, **89315**, 445.

* Segall, P., 2010. *Earthquake and Volcano Deformation*, Princeton University Press.

* Shi, Z. & Day, S.M., 2013. Rupture dynamics and ground motion from 3-d rough-fault simulations, *J. Geophys. Res.: Solid Earth*, **118**(3), 1122–1141.

* Tada, T. & Yamashita, T., 1996. The paradox of smooth and abrupt bends in two-dimensional in-plane shear-crack mechanics, *Geophys. J. Int.*, **127**(3), 795–800.
Tada, T. & Yamashita, T., 1997. Non-hypersingular boundary integral equations for two dimensional non-planar crack analysis, Geophys. J. Int., 130(2), 269–282.

Tal, Y., 2023. Modeling earthquake cycles and wear on rough faults with the mortar finite element method, Geophys. J. Int., 234(1), 190–209.

Wesnousky, S.G., 2006. Predicting the endpoints of earthquake ruptures, Nature, 444(7117), 358–360.

Wesnousky, S. G., 2008. Displacement and geometrical characteristics of earthquake surface ruptures: Issues and implications for seismic-hazard analysis and the process of earthquake rupture, Bull. Seism. Soc. Am., 98(4), 1609–1632.

Wise, L. & Tal, Y., 2024. The evolution of stresses and shear resistance on rough faults at large slip, Tectonophysics, 880, 230335, https://doi.org/10.1016/j.tecto.2024.230335.

APPENDIX A: TABLE OF PARAMETERS

Table A1 is listing the parameters used in this article.

<table><caption>Table A1. Table of parameters.</caption><thead><tr><th>Symbol</th><th>Description</th></tr></thead><tbody><tr><td>τ<sub>f</sub></td><td>Frictional shear traction resistance on the fault</td></tr><tr><td>τ<sub>el</sub></td><td>Total elastic shear traction along the fault</td></tr><tr><td>τ<sub>load</sub></td><td>Shear traction loading on the fault</td></tr><tr><td>h</td><td>Height of the fault</td></tr><tr><td>m</td><td>Slope of the fault</td></tr><tr><td>x</td><td>Point of evaluation</td></tr><tr><td>x<sub>d</sub></td><td>d<sup>th</sup> component of the point x</td></tr><tr><td>y</td><td>Point of evaluation, usually over which an integral is performed</td></tr><tr><td>y<sub>d</sub></td><td>d<sup>th</sup> component of the point y</td></tr><tr><td>κ<sup>t</sup></td><td>Curvature of the fault</td></tr><tr><td>u<sub>k</sub></td><td>k<sup>th</sup> component of the displacement vector</td></tr><tr><td>c<sub>ijpq</sub></td><td>ijpq component of the Hooke tensor</td></tr><tr><td>Δu<sub>i</sub></td><td>i<sup>th</sup> component of displacement discontinuity across the fault</td></tr><tr><td>u<sub>i</sub><sup>+</sup></td><td>i<sup>th</sup> displacement on the positive side of the fault</td></tr><tr><td>u<sub>i</sub><sup>-</sup></td><td>i<sup>th</sup> displacement on the negative side of the fault</td></tr><tr><td>Δu<sub>i</sub><sup>n</sup></td><td>i<sup>th</sup> component of the slip vector for opening (mode I)</td></tr><tr><td>Δu<sub>i</sub><sup>t</sup></td><td>i<sup>th</sup> component of the slip vector for in-plane shear (mode II)</td></tr><tr><td>Δu<sub>i</sub><sup>s</sup></td><td>i<sup>th</sup> component of the slip vector for out-of-plane shear (mode III)</td></tr><tr><td>n<sub>j</sub></td><td>j<sup>th</sup> of the vector normal to the fault</td></tr><tr><td>G<sub>kp</sub></td><td>kp component of the Green's functions</td></tr><tr><td>ξ</td><td>Arc length along the fault</td></tr><tr><td>ε<sub>cd</sub></td><td>cd component of the strain tensor</td></tr><tr><td>σ<sub>ab</sub></td><td>ab component of the stress tensor</td></tr><tr><td>t<sub>j</sub></td><td>j<sup>th</sup> component of the tangential vector to the fault</td></tr><tr><td>∂<sub>j</sub>/∂y<sub>j</sub></td><td>Derivative along the fault (∂<sub>j</sub>/∂x<sub>1</sub> = t<sub>1</sub>∂<sub>1</sub>/∂x<sub>1</sub> + t<sub>2</sub>∂<sub>2</sub>/∂x<sub>2</sub>)</td></tr><tr><td>τ<sub>el</sub><sup>0</sup></td><td>Elastic shear traction at 0<sup>th</sup> order along the fault</td></tr><tr><td>τ<sub>el</sub><sup>1</sup></td><td>Elastic shear traction at 1<sup>th</sup> order along the fault</td></tr><tr><td>σ<sub>el</sub><sup>0</sup></td><td>Total elastic normal traction along the fault</td></tr><tr><td>σ<sub>el</sub><sup>1</sup></td><td>Elastic normal traction at 0<sup>th</sup> order along the fault</td></tr><tr><td>σ<sub>el</sub><sup>2</sup></td><td>Elastic normal traction at 1<sup>th</sup> order along the fault</td></tr><tr><td>L</td><td>Length of the fault, used for the prescribed slip definition</td></tr><tr><td>μ</td><td>Shear modulus</td></tr><tr><td>ν</td><td>Poisson's ratio</td></tr><tr><td>τ<sub>drag</sub><sup>1</sup></td><td>Shear traction drag</td></tr><tr><td>l</td><td>Lengthscale for the definition of some fault geometries</td></tr><tr><td>H</td><td>Fourier transform of the height of the fault (H = F[h])</td></tr><tr><td>α</td><td>Amplitude to wavelength ratio for a self-similar fault</td></tr><tr><td>λ<sub>min</sub></td><td>Minimum roughness wavelength</td></tr><tr><td>H</td><td>Heaviside function</td></tr><tr><td>χ<sub>T</sub></td><td>Rectangular function</td></tr><tr><td>P<sub>h</sub></td><td>Power spectra density of the height of the fault</td></tr><tr><td>σ<sub>load</sub><sup>0</sup></td><td>Constant normal traction loading</td></tr><tr><td>τ<sub>load</sub><sup>0</sup></td><td>Constant shear traction loading</td></tr><tr><td>f</td><td>Coulomb friction parameter</td></tr><tr><td>σ&#773;</td><td>Constant stress state tensor</td></tr><tr><td>σ<sub>11</sub>, σ<sub>22</sub>, and σ<sub>12</sub></td><td>Components of the constant stress state tensor</td></tr><tr><td>δ</td><td>Dirac function</td></tr><tr><td>λ</td><td>Lamé's first parameter</td></tr><tr><td>θ</td><td>Angle between the axis 1 and the fault (see Fig. A2)</td></tr></tbody></table>
<figure><img src="image_20.png" /><figcaption>**Figure A1.** Figure showing the discretization of the fault.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure A2.** Figure showing the simplification that happens in the small slope approximation for the normal and tangential vectors.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure A3.** Assuming the fault geometry and the opening distribution $\Delta u(y_1) = (1 - 4\frac{y_1^2}{L^2})^{3/2}$ in meter (red curve), the full solution, the gradient term and the curvature term can be calculated using eq. (G6). The 0th order and the 1st order, are respectively calculated using the expressions (G3) and (G4).</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure A4.** Example of solving a free boundary condition on non-planar crack. (a) Geometry and boundary condition. (b) Slip and opening on the crack. Please notice that shear and opening are components of a slip-opening vector and that the modulus of this vector seems to give the planar crack opening response. (c) Shear traction on the crack. Please notice that the sum of each contribution is null (free boundary condition). (d) Normal traction on the crack. Please notice that the sum of each contribution is null (free boundary condition).</figcaption></figure>

## APPENDIX B: DERIVATION

In the following, the convention for stress is chosen as tension positive and compression negative.

### B1 Formulas

Here, we just recall some formulas that will be used in the following section for the derivation. The momentum balance equation:

$$c_{ijpq} \frac{\partial}{\partial x_j} \frac{\partial}{\partial x_q} G_{pn}(\mathbf{x}, \mathbf{y}) = 0. \qquad (B1)$$

Symmetry of the Green's function:

$$\frac{\partial}{\partial x_q} G_{ij}(\mathbf{x}, \mathbf{y}) = - \frac{\partial}{\partial y_q} G_{ij}(\mathbf{x}, \mathbf{y}). \qquad (B2)$$

Definition of the tangential differential operator on a function f (Bonnet 1999):

$$D_{ij}[f(\mathbf{y})] = n_i(\mathbf{y}) \frac{\partial}{\partial y_j} f(\mathbf{y}) - n_j(\mathbf{y}) \frac{\partial}{\partial y_i} f(\mathbf{y}). \qquad (B3)$$

Definition of the derivative along the fault, if t is the tangential vector to the fault:

$$\frac{\partial}{\partial y^t} f = t_1 \frac{\partial}{\partial y_1} f + t_2 \frac{\partial}{\partial y_2} f. \qquad (B4)$$

The expression for the Green's function can be found in Tada and Yamashita (1997):

$$G_{ij}(\mathbf{x}, \mathbf{y}) = \frac{1}{4\pi\mu} \left[ \frac{\gamma_i \gamma_j}{2(1-\nu)} - \frac{3-4\nu}{2(1-\nu)} \delta_{ij} \log(r) \right], \quad \text{for } i, j \in \{1, 2\}, \qquad (B5)$$

$$G_{33}(\mathbf{x}, \mathbf{y}) = \frac{-1}{2\pi\mu} \log(r).$$

where $\gamma_1 = \frac{x_1 - y_1}{r}$, $\gamma_2 = \frac{x_2 - y_2}{r}$, and $r = \sqrt{(x_1 - y_1)^2 + (x_2 - y_2)^2}$ is the distance between the points x and y.
## B2 Derivation

We start from the representation theorem (Tada and Yamashita 1997):

$$u_k(\mathbf{x}) = - \int_{\text{fault}} c_{ijpq} n_j(\mathbf{y}) \Delta u_i(\mathbf{y}) \frac{\partial}{\partial x_q} G_{kp}(\mathbf{x}, \mathbf{y}) d\xi(\mathbf{y}). \quad (B6)$$

From the previous equation, by deriving under the integral it can be obtained:

$$\frac{\partial}{\partial x_l} u_k(\mathbf{x}) = - \int_{\text{fault}} c_{ijpq} n_j(\mathbf{y}) \Delta u_i(\mathbf{y}) \frac{\partial}{\partial x_q} \left[ \frac{\partial}{\partial x_l} G_{kp}(\mathbf{x}, \mathbf{y}) \right] d\xi(\mathbf{y}). \quad (B7)$$

The reason to do so is that both the strains ($\epsilon_{kl} = \frac{1}{2}(\frac{\partial}{\partial x_k}u_l + \frac{\partial}{\partial x_l}u_k)$) and the stresses ($\sigma_{ab} = c_{abkl}\epsilon_{kl}$) can be obtained by linear combination of the previous equation. To regularize this integral we need to work on the integrand of the previous integral. Using the symmetry of second derivatives, we can replace the derivatives with respect to $x$ with derivatives with respect to $y$:

$$\frac{\partial}{\partial x_l} u_k(\mathbf{x}) = - \int_{\text{fault}} c_{ijpq} n_j(\mathbf{y}) \Delta u_i(\mathbf{y}) \frac{\partial}{\partial y_q} \left[ \frac{\partial}{\partial y_l} G_{kp}(\mathbf{x}, \mathbf{y}) \right] d\xi(\mathbf{y}). \quad (B8)$$

Then, we used the definition of the tangential differential operator together with the momentum balance equation (eq. B1):

$$\begin{align}
c_{ijpq} n_j(\mathbf{y}) \frac{\partial}{\partial y_l} \frac{\partial}{\partial y_q} G_{kp} &= c_{ijpq} \left[ D_{jl} \frac{\partial}{\partial y_q} G_{kp} + n_l(\mathbf{y}) \frac{\partial}{\partial y_j} \frac{\partial}{\partial y_q} G_{kp} \right] \tag{B9} \\
&= c_{ijpq} D_{jl} \frac{\partial}{\partial y_q} G_{kp}. \nonumber
\end{align}$$

Please note that this derivation is due to Daisuke Sato (Sato et al. 2020) for the 2D case. The equivalent 3D case derivation was done by Marc Bonnet (Bonnet 1999), and can also be found in Romanet et al. (2020).

We then developed the differential tangential operator in the local coordinate system of the fault, and used the fact that the derivative perpendicular to the fault direction are null:

$$\begin{align}
D_{jk}[f(\mathbf{y})] &= \left[ n_j(\mathbf{y}) \frac{\partial}{\partial y_k} - n_k(\mathbf{y}) \frac{\partial}{\partial y_j} \right] f(\mathbf{y}) \nonumber \\
&= \left[ n_j(\mathbf{y}) \left( t_k(\mathbf{y}) \frac{\partial}{\partial y^t} + n_k(\mathbf{y}) \frac{\partial}{\partial y^n} \right) - n_k(\mathbf{y}) \left( t_j(\mathbf{y}) \frac{\partial}{\partial y^t} + n_j(\mathbf{y}) \frac{\partial}{\partial y^n} \right) \right] f(\mathbf{y}) \tag{B10} \\
&= \left[ n_j(\mathbf{y}) t_k(\mathbf{y}) - n_k(\mathbf{y}) t_j(\mathbf{y}) \right] \frac{\partial}{\partial y^t} f(\mathbf{y}). \nonumber
\end{align}$$

It can be checked that the previous equation is 0 for any pair in (1,2), because $t_1 = n_2$ and $t_2 = -n_1$.

Then, we replace the integrand in eq. (B8), using eq. (B9) and eq. (B10):

$$\frac{\partial}{\partial x_l} u_k(\mathbf{x}) = - \int_{\text{fault}} c_{ijpq} \Delta u_i [n_j(\mathbf{y})t_l(\mathbf{y}) - n_l(\mathbf{y})t_j(\mathbf{y})] \frac{\partial}{\partial y^t} \frac{\partial}{\partial y_q} G_{kp}(\mathbf{x}, \mathbf{y}) d\xi(\mathbf{y}). \quad (B11)$$

Finally, it is possible to perform an integration by parts to regularize the hypersingular integral given previously (eq. B11):

$$\frac{\partial}{\partial x_l} u_k(\mathbf{x}) = \int_{\text{fault}} c_{ijpq} [n_j(\mathbf{y})t_l(\mathbf{y}) - n_l(\mathbf{y})t_j(\mathbf{y})] \frac{\partial}{\partial y^t} \Delta u_i \frac{\partial}{\partial y_q} G_{kp}(\mathbf{x}, \mathbf{y}) d\xi(\mathbf{y}). \quad (B12)$$

For numerical calculation and the discretization of the fault slip, it is best to use eq. (B11). However for facilitate interpretation, the slip vector can be projected on the local basis as:

$$\begin{align}
\Delta u_i &= \Delta u^n n_i, && \text{for mode I,} \nonumber \\
\Delta u_i &= \Delta u^t t_i, && \text{for mode II,} \tag{B13} \\
\Delta u_i &= \Delta u^s s_i, && \text{for mode III,} \nonumber
\end{align}$$

where $\Delta u^n$, $\Delta u^t$, and $\Delta u^s$ are respectively the opening, the in-plane slip and the out-of-plane slip. The derivative of slip with respect to the direction of the fault $\frac{\partial}{\partial y^t} \Delta u_i$ can then be replaced by:

$$\begin{align}
\frac{\partial}{\partial y^t} [\Delta u^n n_i] &= n_i \frac{\partial}{\partial y^t} \Delta u^n - t_i \Delta u^n \kappa^t, && \text{for mode I,} \nonumber \\
\frac{\partial}{\partial y^t} [\Delta u^t t_i] &= t_i \frac{\partial}{\partial y^t} \Delta u^t + n_i \Delta u^t \kappa^t, && \text{for mode II,} \tag{B14} \\
\frac{\partial}{\partial y^t} [\Delta u^s s_i] &= s_i \frac{\partial}{\partial y^t} \Delta u^s = \delta_{3i} \frac{\partial}{\partial y^t} \Delta u^s, && \text{for mode III,} \nonumber
\end{align}$$

where $\delta_{3i}$ is the Kronecker delta.
## APPENDIX C: DISCRETIZATION

The discretization is using the same strategy as in Romanet et al. (2020) (see Fig. A1). It consists of discretizing the slip and tangential vector to be constant over a straight line. It is convenient for this part to use the curvilinear abcisse $\xi$ instead of the position along the fault $y(\xi)$. Also, to simplify the demonstration, we will just keep the scalar $t$ instead of one component of the tangential vector $t_i$. Discretized tangential vector and slip can be written:

$$ 
\begin{aligned}
t(\xi) &= \sum_j t_j [H(\xi - \xi_j) - H(\xi - \xi_{j+1})], \\
\Delta u(\xi) &= \sum_j \Delta u_j [H(\xi - \xi_j) - H(\xi - \xi_{j+1})].
\end{aligned} 
\qquad (C1) $$

Then the integration with any kernel $K(x, y)$, will give:

$$ 
\begin{aligned}
\int_{\text{fault}} K(\mathbf{x}, \mathbf{y}) \frac{\partial}{\partial y^t} [t \Delta u] \, d\xi(\mathbf{y}) &= \\
\int_{\text{fault}} K(\mathbf{x}, \mathbf{y}) t \frac{\partial}{\partial y^t} [\Delta u] \, d\xi(\mathbf{y}) + \\
\int_{\text{fault}} K(\mathbf{x}, \mathbf{y}) \Delta u \frac{\partial}{\partial y^t} [t] \, d\xi(\mathbf{y}) &= \\
\int_{\text{fault}} K(\mathbf{x}, \mathbf{y}) t \frac{\partial}{\partial y^t} \left[ \sum_j \Delta u_j [H(\xi - \xi_j) - H(\xi - \xi_{j+1})] \right] \, d\xi(\mathbf{y}) + \\
\int_{\text{fault}} K(\mathbf{x}, \mathbf{y}) \Delta u \frac{\partial}{\partial y^t} \left[ \sum_j t_j [H(\xi - \xi_j) - H(\xi - \xi_{j+1})] \right] \, d\xi(\mathbf{y}) &= \\
\sum_j [K[\mathbf{x}, \mathbf{y}(\xi_j)] t(\xi_j) - K[\mathbf{x}, \mathbf{y}(\xi_{j+1})] t(\xi_{j+1})] \Delta u_j + \\
\sum_j [K[\mathbf{x}, \mathbf{y}(\xi_j)] \Delta u(\xi_j) - K[\mathbf{x}, \mathbf{y}(\xi_{j+1})] \Delta u(\xi_{j+1})] t_j.
\end{aligned} 
\qquad (C2) $$

We can then replace $t(\xi_j)$ and $\Delta u(\xi_j)$ by $(H(0) = 1/2)$:

$$ 
\begin{aligned}
t(\xi_j) &= (1 - 0.5)t_{j-1} + 0.5t_j = \frac{t_{j-1} + t_j}{2}, \\
\Delta u(\xi_j) &= \frac{\Delta u_{j-1} + \Delta u_j}{2},
\end{aligned} 
\qquad (C3) $$

which gives:

$$ 
\begin{aligned}
\int_{\text{fault}} K(\mathbf{x}, \mathbf{y}) \frac{\partial}{\partial y^t} [t \Delta u] \, d\xi(\mathbf{y}) &= \\
&\underbrace{\sum_j \left[ K[\mathbf{x}, \mathbf{y}(\xi_j)] \frac{t_j + t_{j-1}}{2} - K[\mathbf{x}, \mathbf{y}(\xi_{j+1})] \frac{t_{j+1} + t_j}{2} \right]}_{\text{Gradient term}} \Delta u_j + \\
&\underbrace{\sum_j \left[ K[\mathbf{x}, \mathbf{y}(\xi_j)] \frac{\Delta u_j + \Delta u_{j-1}}{2} - K[\mathbf{x}, \mathbf{y}(\xi_{j+1})] \frac{\Delta u_{j+1} + \Delta u_j}{2} \right]}_{\text{Curvature term}} t_j.
\end{aligned} 
\qquad (C4) $$

This way of discretizing the integral allows for the separation of the curvature and gradient terms.

## APPENDIX D: THE SMALL SLOPE APPROXIMATION

### D1 It transforms the integrals along the fault to integral along a straight line

Starting from a regularized boundary element method, and applying the small slope approximation, it is possible to obtain the zeroth and first-order solution of the stress on the fault due to a slip distribution for the different mode of slip. Please notice that all the integrals that were along the fault become along the $y_1$ axis only. This is due to the fact that if $y_2 = h(y_1)$ and the slope derivative is small $\frac{d}{dy_1}h(y_1) \ll 0$,
and changing the variable of the integral:

$$ 
\begin{aligned}
\int_{\text{fault}} f(y_1)d\xi(y_1) &= \int_{-\infty}^{+\infty} f(y_1)\sqrt{1+h^2(y_1)}dy_1 \\
&\approx \int_{-\infty}^{+\infty} f(y_1)(1+h^2(y_1)/2)dy_1 \qquad (\text{D1}) \\
&\approx \int_{-\infty}^{+\infty} f(y_1)dy_1.
\end{aligned}
$$

## D2 It provides an expression for the normal and tangential vector to the fault

The normal and tangential vectors can be simplified in the small slope approximation (see Fig. A2), indeed:

$$ m(y_1) = \tan[\theta(y_1)] \approx \theta(y_1), \qquad (\text{D2}) $$

so that in the small slope approximation:

$$ 
\begin{aligned}
\mathbf{n}(y_1) &= \begin{pmatrix} -\sin(\theta(y_1)) \\ \cos(\theta(y_1)) \end{pmatrix} \approx \begin{pmatrix} -m(y_1) \\ 1 \end{pmatrix}, \\
\mathbf{t}(y_1) &= \begin{pmatrix} \cos(\theta(y_1)) \\ \sin(\theta(y_1)) \end{pmatrix} \approx \begin{pmatrix} 1 \\ m(y_1) \end{pmatrix}.
\end{aligned}
\qquad (\text{D3}) $$

# APPENDIX E: A PARTICULAR SOLUTION FOR THE SEMI-ANALYTICAL FORMULATION

The slip/opening distribution for the result section is chosen as following:

$$ \Delta u(y_1) = (1 - y_1^2)^{3/2}. \qquad (\text{E1}) $$

The reason for this particular form of slip (or opening) distribution compared to the classic distribution in linear fracture mechanics ($\Delta u(y_1) = \sqrt{1 - y_1^2}$) is because it brings finite stress concentrations outside the fault (hence the numerical result are converging) and that there is an analytical solution for the stresses in the case of planar fault:

$$ \frac{1}{2\pi} \int_{-1}^{1} \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u(y_1) dy_1 = \frac{3}{4} \begin{cases} 1 - 2x_1^2 - 2x_1\sqrt{x_1^2 - 1}, & \text{if } x_1 < -1 \\ 1 - 2x_1^2, & \text{if } -1 \le x_1 \le 1, \text{ inside the fault} \\ 1 - 2x_1^2 + 2x_1\sqrt{x_1^2 - 1}, & \text{if } x_1 > 1 \end{cases} \qquad (\text{E2}) $$

The previous solution was used to check our numerical scheme against the analytical formulation. Obtaining this solution is not straightforward. For the case $-1 < x_1 < 1$, the solution can be obtained by using Chebychev polynomial (Mason and Handscomb 2002, section 9.5.1). The solution outside $|x_1| > 1$ can be obtained by doing a first change of variable $y_1 \to \sin(\theta)$ and checking that $\frac{\sin(\theta)\cos^2(\theta)}{x_1 - \sin(\theta)}$ is symmetric at $\pi/2$. That allows us to make the integration over the whole circle and take half the value. If we start from the following integral:

$$ \frac{1}{2\pi} \int_{-1}^{1} \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u(y_1) dy_1 = -\frac{3}{2\pi} \int_{-1}^{1} \frac{y_1 \sqrt{1 - y_1^2}}{x_1 - y_1} dy_1. \qquad (\text{E3}) $$

And then do the first change of variable:

$$ 
\begin{aligned}
\frac{1}{2\pi} \int_{-1}^{1} \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u(y_1) dy_1 &= -\frac{3}{2\pi} \int_{-\pi/2}^{\pi/2} \frac{\sin(\theta) \cos(\theta) |\cos(\theta)|}{x_1 - \sin(\theta)} d\theta \\
&= -\frac{3}{2\pi} \int_{-\pi/2}^{\pi/2} \frac{\sin(\theta) \cos^2(\theta)}{x_1 - \sin(\theta)} d\theta \qquad (\text{E4}) \\
&= -\frac{3}{4\pi} \int_{-\pi/2}^{3\pi/2} \frac{\sin(\theta) \cos^2(\theta)}{x_1 - \sin(\theta)} d\theta.
\end{aligned}
$$

A second change of variable is using $z = e^{i\theta}$, where $\sin(\theta) = \frac{z-1/z}{2i}$ and $\cos(\theta) = \frac{z+1/z}{2}$.

$$ -\frac{3}{4\pi} \int_{-\pi/2}^{3\pi/2} \frac{\sin(\theta) \cos^2(\theta)}{x_1 - \sin(\theta)} d\theta = -\frac{3i}{16\pi} \int_C \frac{(z^2 - 1)(1 + z^2)^2}{z^3[z - (ix + \sqrt{1-x^2})][z - (ix - \sqrt{1-x^2})]} dz. \qquad (\text{E5}) $$
For $x < -1$, there are two poles at $z = 0$ and $z = i(x + \sqrt{x^2 - 1})$, so that applying the residue theorem gives:

$$ \begin{aligned} & \frac{1}{2\pi} \int_{-1}^{1} \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u(y_1) dy_1 = 2\pi i \operatorname{Res}(f, z=0) + 2\pi i \operatorname{Res}(f, z=i(x_1 + \sqrt{x_1^2-1})) \\ & = \frac{3}{4}(1-2x_1^2-2x_1\sqrt{x_1^2-1}), \end{aligned} \quad (E6) $$

where $\operatorname{Res}(f,z)$ is the residue of $f$ in $z$.

For $x > 1$, there are also two poles at $z = 0$ and $z = i(x_1 - \sqrt{x_1^2 - 1})$, so that the residue theorem yields:

$$ \begin{aligned} & \frac{1}{2\pi} \int_{-1}^{1} \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u(y_1) dy_1 = 2\pi i \operatorname{Res}(f, z=0) + 2\pi i \operatorname{Res}(f, z=i(x_1 - \sqrt{x_1^2-1})) \\ & = \frac{3}{4}(1-2x_1^2+2x_1\sqrt{x_1^2-1}). \end{aligned} \quad (E7) $$

## APPENDIX F: DEFINITION OF THE FOURIER TRANSFORM

In the following part, the zeroth and first-order solution for both normal and shear traction are provided. Their expression depends on the choice of the definition of the Fourier transform. We define the Fourier transform as following:

$$ F[f](k) = \int_{-\infty}^{+\infty} f(x_1)e^{-ikx_1}dx_1, \quad (F1) $$

$$ F^{-1}[f](x_1) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} f(x_1)e^{ikx_1}dk. \quad (F2) $$

where $k$ is the wavenumber. In particular, we will use the Fourier transform of the two functions $u(x_1) = 1/x_1$ and $v(x_1) = 1/x_1^2$:

$$ \begin{aligned} U(k) &= F[1/x_1](k) = -i\pi \operatorname{sign}(k), \\ V(k) &= F[1/x_1^2](k) = -\pi |k|. \end{aligned} \quad (F3) $$

## APPENDIX G: IN-PLANE OPENING (MODE I)

The mode I correspond to opening. One common assumption for an opening fault is the traction free condition at the surface.

By making the small slope approximation, it can be shown that the main contribution for the shear and normal tractions are:

$$ \begin{array}{l} \underbrace{\tau_{\text{el}}}_{\substack{\text{Elastic shear traction} \\ \text{Only if the fault is non-planar}}}= \underbrace{\tau_{\text{el}}^{0}}_{\substack{\text{Shear traction perturbation} \\ \propto \alpha^t \Delta u^t}} + \underbrace{\tau_{\text{el}}^{1}}_{\substack{\text{1st order} \\ \text{Shear traction perturbation}}} + \underbrace{\dots}_{\substack{\text{Higher order terms}}} , \end{array} \quad (G1) $$

$$ \begin{array}{l} \underbrace{\sigma_{\text{el}}}_{\substack{\text{Elastic normal} \\ \text{Of non-planar fault}}} = \underbrace{\sigma_{\text{el}}^{0}}_{\substack{\text{Planar fault response}}} + \underbrace{\sigma_{\text{el}}^{1}}_{\substack{\text{Normal traction drag}}} + \underbrace{\dots}_{\substack{\text{Higher order terms}}} , \end{array} \quad (G2) $$

with:

$$ \tau_{\text{el}}^{0}(x_1) = \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \frac{1}{x_1 - y_1} \kappa^{t}(y_1) \Delta u^{n}(y_1) dy_1, \quad (G3) $$

$$ \sigma_{\text{el}}^{0}(x_1) = -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u^{n}(y_1) dy_1, \quad (G3) $$

$$ \tau_{\text{el}}^{1}(x_1) = \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{m(y_1)}{x_1 - y_1} - \frac{x_2 - y_2}{(x_1 - y_1)^2} \right] \frac{d}{dy_1} \Delta u^{n}(y_1) dy_1, \quad (G4) $$

$$ \sigma_{\text{el}}^{1}(x_1) = \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{-2m(x_1)}{x_1 - y_1} + \frac{m(y_1)}{x_1 - y_1} + \frac{x_2 - y_2}{(x_1 - y_1)^2} \right] \kappa^{t}(y_1) \Delta u^{n}(y_1) dy_1. \quad (G4) $$

On Fig. A3, an example of an in-plane opening fault is given for a rough geometry. An interesting feature is that the on-fault shear traction mainly (0th order) depends upon the geometry, because the zeroth-order shear traction depends on the local curvature (a geometrical parameter) that multiplies the opening. On the contrary, the normal traction depends mainly on the derivative of opening along the fault. That means that it is more dependent on the opening distribution than the fault geometry.

### G1 The special case of traction free condition on a fully opened non-planar crack

By looking at the zeroth-order solution (eq. G3), it is interesting to see that the result of opening on a non-planar fault leads to the apparition of shear traction on this fault. If we assume traction free condition, the only solution to remove the shear traction that appears is that there is
also some in-plane slip along the fault (mode II). One result from that is that pure opening (mode I) on a non-planar fault, with the traction free condition cannot exist without in-plane shearing (mode II). Hence, the traction free boundary condition can be written as:

$$ 
\begin{aligned}
\tau_{\text{load}} + \tau_{\text{el}}(\Delta u^n) + \tau_{\text{el}}(\Delta u^t) &= 0 \\
\sigma_{\text{load}} + \sigma_{\text{el}}(\Delta u^n) + \sigma_{\text{el}}(\Delta u^t) &= 0
\end{aligned} 
\qquad (G5) $$

On Fig. A4, we are showing an example of an opening non-planar crack. It can be seen that there is indeed shear slip as well as pure opening. The sum of all contribution, the elastic traction due to shear slip, the elastic traction due to opening, and finally the traction due to loading cancel on the fault.

## G2 Full solution

In the following, $\gamma_1 = \frac{x_1-y_1}{r}$, $\gamma_2 = \frac{x_2-y_2}{r}$ and * represent the convolution operation. These equations are developed by using the definition of the strain $\epsilon_{kl} = \frac{\partial}{\partial x_k} u_l + \frac{\partial}{\partial x_l} u_k$ and the stress ($\sigma_{ab} = c_{abkl} \epsilon_{kl}$) tensors, as well as eq. (B12), and the according slip mode eq. (B14). We also recall that the Fourier transform of $u(x) = 1/x$ and $v(x) = 1/x^2$, are respectively noted $U = F[u]$ and $V = F[v]$. Finally, $\lambda$ is the Lamé's first parameter.

$$ 
\begin{aligned}
\epsilon_{11}(\mathbf{x}) = & \frac{1}{2\pi(\lambda + 2\mu)} \int_{\text{fault}} \left[ n_1 \frac{\gamma_2}{r} (2\lambda\gamma_1^2 + \mu(3\gamma_1^2 + \gamma_2^2)) + n_2 \frac{\gamma_1}{r} (2\lambda\gamma_2^2 + \mu(\gamma_2^2 - \gamma_1^2)) \right] \frac{d}{dy'} \Delta u^n(\mathbf{y}) \, d\xi(\mathbf{y}) \\
& + \frac{1}{2\pi(\lambda + 2\mu)} \int_{\text{fault}} \left[ -n_2 \frac{\gamma_2}{r} (2\lambda\gamma_1^2 + \mu(3\gamma_1^2 + \gamma_2^2)) + n_1 \frac{\gamma_1}{r} (2\lambda\gamma_2^2 + \mu(\gamma_2^2 - \gamma_1^2)) \right] \kappa'(\mathbf{y}) \Delta u^n(\mathbf{y}) \, d\xi(\mathbf{y})
\end{aligned}
$$

$$ 
\begin{aligned}
\epsilon_{22}(\mathbf{x}) = & -\frac{1}{2\pi(\lambda + 2\mu)} \int_{\text{fault}} \left[ n_2 \frac{\gamma_1}{r} (2\lambda\gamma_2^2 + \mu(\gamma_1^2 + 3\gamma_2^2)) + n_1 \frac{\gamma_2}{r} (2\lambda\gamma_1^2 + \mu(\gamma_1^2 - \gamma_2^2)) \right] \frac{d}{dy'} \Delta u^n(\mathbf{y}) \, d\xi(\mathbf{y}) \\
& - \frac{1}{2\pi(\lambda + 2\mu)} \int_{\text{fault}} \left[ n_1 \frac{\gamma_1}{r} (2\lambda\gamma_2^2 + \mu(\gamma_1^2 + 3\gamma_2^2)) - n_2 \frac{\gamma_2}{r} (2\lambda\gamma_1^2 + \mu(\gamma_1^2 - \gamma_2^2)) \right] \kappa'(\mathbf{y}) \Delta u^n(\mathbf{y}) \, d\xi(\mathbf{y})
\end{aligned}
$$

$$ 
\begin{aligned}
\epsilon_{12}(\mathbf{x}) = & -\frac{\lambda + \mu}{2\pi(\lambda + 2\mu)} \int_{\text{fault}} \left[ (n_1 \frac{\gamma_1}{r} + n_2 \frac{\gamma_2}{r}) (\gamma_1^2 - \gamma_2^2) \right] \frac{d}{dy'} \Delta u^n(\mathbf{y}) \, d\xi(\mathbf{y}) \\
& -\frac{\lambda + \mu}{2\pi(\lambda + 2\mu)} \int_{\text{fault}} \left[ (-n_2 \frac{\gamma_1}{r} + n_1 \frac{\gamma_2}{r}) (\gamma_1^2 - \gamma_2^2) \right] \kappa'(\mathbf{y}) \Delta u^n(\mathbf{y}) \, d\xi(\mathbf{y})
\end{aligned}
$$

$$ 
\begin{aligned}
\sigma_{11}(\mathbf{x}) = & \frac{\mu}{2\pi(1-\nu)} \int_{\text{fault}} \left[ n_1(\mathbf{y}) (3\gamma_1^2 + \gamma_2^2) \frac{\gamma_2}{r} - n_2(\mathbf{y}) (\gamma_1^2 - \gamma_2^2) \frac{\gamma_1}{r} \right] \frac{d}{dy'} \Delta u^n(\mathbf{y}) \, d\xi(\mathbf{y}) \\
& - \frac{\mu}{2\pi(1-\nu)} \int_{\text{fault}} \left[ n_2(\mathbf{y}) (3\gamma_1^2 + \gamma_2^2) \frac{\gamma_2}{r} + n_1(\mathbf{y}) (\gamma_1^2 - \gamma_2^2) \frac{\gamma_1}{r} \right] \kappa'(\mathbf{y}) \Delta u^n(\mathbf{y}) \, d\xi(\mathbf{y})
\end{aligned} 
\qquad (G6) $$

$$ 
\begin{aligned}
\sigma_{22}(\mathbf{x}) = & -\frac{\mu}{2\pi(1-\nu)} \int_{\text{fault}} \left[ n_2(\mathbf{y}) (\gamma_1^2 + 3\gamma_2^2) \frac{\gamma_1}{r} + n_1(\mathbf{y}) (\gamma_1^2 - \gamma_2^2) \frac{\gamma_2}{r} \right] \frac{d}{dy'} \Delta u^n(\mathbf{y}) \, d\xi(\mathbf{y}) \\
& - \frac{\mu}{2\pi(1-\nu)} \int_{\text{fault}} \left[ n_1(\mathbf{y}) (\gamma_1^2 + 3\gamma_2^2) \frac{\gamma_1}{r} - n_2(\mathbf{y}) (\gamma_1^2 - \gamma_2^2) \frac{\gamma_2}{r} \right] \kappa'(\mathbf{y}) \Delta u^n(\mathbf{y}) \, d\xi(\mathbf{y})
\end{aligned}
$$

$$ 
\begin{aligned}
\sigma_{12}(\mathbf{x}) = & -\frac{\mu}{2\pi(1-\nu)} \int_{\text{fault}} \left[ (n_1\gamma_1 + n_2\gamma_2) \left(\frac{\gamma_1^2}{r} - \frac{\gamma_2^2}{r}\right) \right] \frac{d}{dy'} \Delta u^n(\mathbf{y}) \, d\xi(\mathbf{y}) \\
& -\frac{\mu}{2\pi(1-\nu)} \int_{\text{fault}} \left[ (n_1\gamma_2 - n_2\gamma_1) \left(\frac{\gamma_1^2}{r} - \frac{\gamma_2^2}{r}\right) \right] \kappa'(\mathbf{y}) \Delta u^n(\mathbf{y}) \, d\xi(\mathbf{y})
\end{aligned}
$$

$$ 
\begin{aligned}
\tau_{\text{el}}(\mathbf{x}) = & n_1(\mathbf{x}) n_2(\mathbf{x}) (\sigma_{11}(\mathbf{x}) - \sigma_{22}(\mathbf{x})) + (n_2^2(\mathbf{x}) - n_1^2(\mathbf{x})) \sigma_{12}(\mathbf{x}) \\
& \sigma_{\text{el}}(\mathbf{x}) = n_1^2(\mathbf{x}) \sigma_{11}(\mathbf{x}) + 2n_1(\mathbf{x}) n_2(\mathbf{x}) \sigma_{12}(\mathbf{x}) + n_2^2(\mathbf{x}) \sigma_{22}(\mathbf{x})
\end{aligned}
$$

## G3 Zeroth order

The zeroth-order solutions are obtained from the full solution by making the small slope approximation and keeping only the terms that are independent from the fault slope.
### G3.1 Space domain

$$
\begin{align*}
\sigma_{11}^{0}(x_1) &= -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u''(y_1) dy_1 \\
\sigma_{22}^{0}(x_1) &= -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u''(y_1) dy_1 \\
\sigma_{12}^{0}(x_1) &= \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \frac{1}{x_1 - y_1} \kappa'(y_1) \Delta u''(y_1) dy_1 \\
\tau_{el}^{0}(x_1) &= \sigma_{12}^{0} \\
&= \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \frac{1}{x_1 - y_1} \kappa'(y_1) \Delta u''(y_1) dy_1 \\
\sigma_{el}^{0}(x_1) &= \sigma_{22}^{0} \\
&= -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u''(y_1) dy_1
\end{align*}
\tag{G7}
$$

### G3.2 Spectral domain

$$
\begin{align*}
\sigma_{11}^{0}(k) &= -\frac{\mu}{2\pi(1-\nu)} \pi |k| F[\Delta u^n] \\
\sigma_{22}^{0}(k) &= -\frac{\mu}{2\pi(1-\nu)} \pi |k| F[\Delta u^n] \\
\sigma_{12}^{0}(k) &= -\frac{\mu}{2\pi(1-\nu)} i \pi \operatorname{sign}(k) F[\kappa'(y_1) \Delta u^n(y_1)] \tag{G8} \\
\tau_{el}^{0}(k) &= -\frac{\mu}{2\pi(1-\nu)} i \pi \operatorname{sign}(k) F[\kappa'(y_1) \Delta u^n(y_1)] \\
\sigma_{el}^{0}(k) &= -\frac{\mu}{2\pi(1-\nu)} \pi |k| F[\Delta u^n]
\end{align*}
$$

## G4 First order

The first-order solutions are obtained from the full solution by making the small slope approximation and keeping only the terms that are proportional to the fault slope.

### G4.1 Space domain

$$
\begin{align*}
\sigma_{11}^{1}(x_1) &= -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ 3 \frac{x_2 - y_2}{(x_1 - y_1)^2} - \frac{m(y_1)}{x_1 - y_1} \right] \kappa'(y_1) \Delta u''(y_1) dy_1 \\
\sigma_{22}^{1}(x_1) &= -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ -\frac{m(y_1)}{x_1 - y_1} - \frac{x_2 - y_2}{(x_1 - y_1)^2} \right] \kappa'(y_1) \Delta u''(y_1) dy_1 \\
\sigma_{12}^{1}(x_1) &= -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ -\frac{m(y_1)}{x_1 - y_1} + \frac{x_2 - y_2}{(x_1 - y_1)^2} \right] \frac{d}{dy_1} \Delta u''(y_1) dy_1 \\
\tau_{el}^{1}(x_1) &= -m(\sigma_{11}^{0} - \sigma_{22}^{0}) + \sigma_{12}^{1} \\
&= \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{m(y_1)}{x_1 - y_1} - \frac{x_2 - y_2}{(x_1 - y_1)^2} \right] \frac{d}{dy_1} \Delta u''(y_1) dy_1 \\
\sigma_{el}^{1}(x_1) &= -2m\sigma_{12}^{0} + \sigma_{22}^{1} \\
&= \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{-2m(x_1)}{x_1 - y_1} + \frac{m(y_1)}{x_1 - y_1} + \frac{x_2 - y_2}{(x_1 - y_1)^2} \right] \kappa'(y_1) \Delta u''(y_1) dy_1
\end{align*}
\tag{G9}
$$
## G4.2 Spectral domain

$$
\begin{align*}
\sigma_{11}^{1}(k) &= -\frac{\mu}{2\pi(1-\nu)} \left( \frac{3}{2\pi} F[h] * (VF[\kappa' \Delta u^n] - 3VF[h\kappa' \Delta u^n] - UF[m\kappa' \Delta u^n]) \right) \\
\sigma_{22}^{1}(k) &= -\frac{\mu}{2\pi(1-\nu)} \left( -UF[m\kappa' \Delta u^n] - \frac{1}{2\pi} F[h] * (VF[\kappa' \Delta u^n] + VF[h\kappa' \Delta u^n]) \right) \\
\sigma_{12}^{1}(k) &= -\frac{\mu}{2\pi(1-\nu)} \left( -UF \left[ m \frac{d}{dy_1} \Delta u^n \right] + \frac{1}{2\pi} F[h] * \left( VF \left[ \frac{d}{dy_1} \Delta u^n \right] \right) - VF \left[ h \frac{d}{dy_1} \Delta u^n \right] \right) \\
\tau_{el}^{1}(k) &= \frac{\mu}{2\pi(1-\nu)} \left( UF \left[ m \frac{d}{dy_1} \Delta u^n \right] - \frac{1}{2\pi} F[h] * \left( VF \left[ \frac{d}{dy_1} \Delta u^n \right] \right) + VF \left[ h \frac{d}{dy_1} \Delta u^n \right] \right) \\
\sigma_{el}^{1}(k) &= \frac{\mu}{2\pi(1-\nu)} \left( -\frac{1}{\pi} F[m] * (VF[\kappa' \Delta u^n]) + UF[m\kappa' \Delta u^n] \right. \\
&\qquad \left. + \frac{1}{2\pi} F[h] * (VF[\kappa' \Delta u^n]) - VF[h\kappa' \Delta u^n] \right)
\end{align*}
\tag{G10}
$$

# APPENDIX H: IN-PLANE SHEAR (MODE II)

## H1 Full solution

In the following, $\gamma_1 = \frac{x_1 - y_1}{r}$, $\gamma_2 = \frac{x_2 - y_2}{r}$ and $*$ represent the convolution operation. These equations are developed by using the definition of the strain $\epsilon_{kl} = \frac{1}{2} (\frac{\partial}{\partial x_k} u_l + \frac{\partial}{\partial x_l} u_k)$ and the stress $(\sigma_{ab} = c_{abkl} \epsilon_{kl})$ tensors, as well as eq. (B12), and the according slip mode eq. (B14). We also recall that the Fourier transform of $u(x) = 1/x$ and $v(x) = 1/x^2$, are respectively noted $U = F[u]$ and $V = F[v]$.

$$
\begin{align}
\epsilon_{11}(\mathbf{x}) &= \frac{1}{2\pi(\lambda + 2\mu)} \int_{\text{fault}} \left[ n_2 \frac{\gamma_2}{r} (2\lambda\gamma_1^2 + \mu(3\gamma_1^2 + \gamma_2^2)) - n_1 \frac{\gamma_1}{r} (2\lambda\gamma_2^2 + \mu(\gamma_2^2 - \gamma_1^2)) \right] \frac{d}{dy'} \Delta u'(\mathbf{y}) \, d\xi(\mathbf{y}) \nonumber \\
&\quad + \frac{1}{2\pi(\lambda + 2\mu)} \int_{\text{fault}} \left[ n_1 \frac{\gamma_2}{r} (2\lambda\gamma_1^2 + \mu(3\gamma_1^2 + \gamma_2^2)) + n_2 \frac{\gamma_1}{r} (2\lambda\gamma_2^2 + \mu(\gamma_2^2 - \gamma_1^2)) \right] \kappa'(\mathbf{y}) \Delta u'(\mathbf{y}) \, d\xi(\mathbf{y}) \nonumber \\
\epsilon_{22}(\mathbf{x}) &= \frac{1}{2\pi(\lambda + 2\mu)} \int_{\text{fault}} \left[ n_1 \frac{\gamma_1}{r} (2\lambda\gamma_2^2 + \mu(\gamma_1^2 + 3\gamma_2^2)) - n_2 \frac{\gamma_2}{r} (2\lambda\gamma_1^2 + \mu(\gamma_1^2 - \gamma_2^2)) \right] \frac{d}{dy'} \Delta u'(\mathbf{y}) \, d\xi(\mathbf{y}) \nonumber \\
&\quad - \frac{1}{2\pi(\lambda + 2\mu)} \int_{\text{fault}} \left[ n_2 \frac{\gamma_1}{r} (2\lambda\gamma_2^2 + \mu(\gamma_1^2 + 3\gamma_2^2)) + n_1 \frac{\gamma_2}{r} (2\lambda\gamma_1^2 + \mu(\gamma_1^2 - \gamma_2^2)) \right] \kappa'(\mathbf{y}) \Delta u'(\mathbf{y}) \, d\xi(\mathbf{y}) \nonumber \\
\epsilon_{12}(\mathbf{x}) &= -\frac{\lambda + \mu}{2\pi(\lambda + 2\mu)} \int_{\text{fault}} \left[ \left(n_2 \frac{\gamma_1}{r} - n_1 \frac{\gamma_2}{r}\right) (\gamma_1^2 - \gamma_2^2) \right] \frac{d}{dy'} \Delta u'(\mathbf{y}) \, d\xi(\mathbf{y}) \nonumber \\
&\quad - \frac{\lambda + \mu}{2\pi(\lambda + 2\mu)} \int_{\text{fault}} \left[ \left(n_1 \frac{\gamma_1}{r} + n_2 \frac{\gamma_2}{r}\right) (\gamma_1^2 - \gamma_2^2) \right] \kappa'(\mathbf{y}) \Delta u'(\mathbf{y}) \, d\xi(\mathbf{y}) \nonumber \\
\sigma_{11}(\mathbf{x}) &= \frac{\mu}{2\pi(1-\nu)} \int_{\text{fault}} \left[ n_2 (3\gamma_1^2 + \gamma_2^2) \frac{\gamma_2}{r} + n_1 (\gamma_1^2 - \gamma_2^2) \frac{\gamma_1}{r} \right] \frac{d}{dy'} \Delta u'(\mathbf{y}) \, d\xi(\mathbf{y}) \nonumber \\
&\quad + \frac{\mu}{2\pi(1-\nu)} \int_{\text{fault}} \left[ n_1 (3\gamma_1^2 + \gamma_2^2) \frac{\gamma_2}{r} - n_2 (\gamma_1^2 - \gamma_2^2) \frac{\gamma_1}{r} \right] \kappa'(\mathbf{y}) \Delta u'(\mathbf{y}) \, d\xi(\mathbf{y}) \tag{H1}
\end{align}
$$

$$
\sigma_{22}(\mathbf{x}) = \frac{\mu}{2\pi(1-\nu)} \int_{\text{fault}} \left[ n_1 \left( \gamma_1^2 + 3\gamma_2^2 \right) \frac{\gamma_1}{r} - n_2 \left( \gamma_1^2 - \gamma_2^2 \right) \frac{\gamma_2}{r} \right] \frac{d}{dy'} \Delta u'(\mathbf{y}) \, d\xi(\mathbf{y}) \\
- \frac{\mu}{2\pi(1-\nu)} \int_{\text{fault}} \left[ n_2 \left( \gamma_1^2 + 3\gamma_2^2 \right) \frac{\gamma_1}{r} + n_1 \left( \gamma_1^2 - \gamma_2^2 \right) \frac{\gamma_2}{r} \right] \kappa'(\mathbf{y}) \Delta u'(\mathbf{y}) \, d\xi(\mathbf{y})
$$

$$
\sigma_{12}(\mathbf{x}) = \frac{\mu}{2\pi(1-\nu)} \int_{\text{fault}} \left[ \left(n_1 \frac{\gamma_2}{r} - n_2 \frac{\gamma_1}{r}\right) \left(\gamma_1^2 - \gamma_2^2\right) \right] \frac{d}{dy'} \Delta u'(\mathbf{y}) \, d\xi(\mathbf{y}) \\
- \frac{\mu}{2\pi(1-\nu)} \int_{\text{fault}} \left[ \left(n_1 \frac{\gamma_1}{r} + n_2 \frac{\gamma_2}{r}\right) \left(\gamma_1^2 - \gamma_2^2\right) \right] \kappa'(\mathbf{y}) \Delta u'(\mathbf{y}) \, d\xi(\mathbf{y})
$$

$$
\tau_{el}(\mathbf{x}) = n_1(\mathbf{x}) n_2(\mathbf{x}) (\sigma_{11}(\mathbf{x}) - \sigma_{22}(\mathbf{x})) + (n_2^2(\mathbf{x}) - n_1^2(\mathbf{x})) \sigma_{12}(\mathbf{x})
$$

$$
\sigma_{el}(\mathbf{x}) = n_1^2(\mathbf{x}) \sigma_{11}(\mathbf{x}) + 2n_1(\mathbf{x}) n_2(\mathbf{x}) \sigma_{12}(\mathbf{x}) + n_2^2(\mathbf{x}) \sigma_{22}(\mathbf{x})
$$

## H2 Zeroth order

The zeroth-order solutions are obtained from the full solution by making the small slope approximation and keeping only the terms that are independent from the fault slope.
## H2.1 Space domain

$$
\begin{align*}
\sigma_{11}^{0}(x_1) &= -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{1}{x_1 - y_1} \kappa'(y_1) \Delta u'(y_1) \right] dy_1 \\
\sigma_{22}^{0}(x_1) &= -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{1}{x_1 - y_1} \kappa'(y_1) \Delta u'(y_1) \right] dy_1 \\
\sigma_{12}^{0}(x_1) &= -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u'(y_1) \right] dy_1 \\
\tau_{el}^{0}(x_1) &= \sigma_{12}^{0} \\
&= -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u'(y_1) \right] dy_1 \\
\sigma_{el}^{0}(x_1) &= \sigma_{22}^{0} \\
&= -\frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{1}{x_1 - y_1} \kappa'(y_1) \Delta u'(y_1) \right] dy_1
\end{align*}
\tag{H2}
$$

## H2.2 Spectral domain

$$
\begin{align*}
\sigma_{11}^{0}(k) &= \frac{\mu}{2\pi(1-\nu)} i\pi \operatorname{sign}(k) F[\kappa' \Delta u'] \\
\sigma_{22}^{0}(k) &= \frac{\mu}{2\pi(1-\nu)} i\pi \operatorname{sign}(k) F[\kappa' \Delta u'] \\
\sigma_{12}^{0}(k) &= -\frac{\mu}{2\pi(1-\nu)} \pi |k| F[\Delta u'] \\
\tau_{el}^{0}(k) &= -\frac{\mu}{2\pi(1-\nu)} \pi |k| F[\Delta u'] \\
\sigma_{el}^{0}(k) &= \frac{\mu}{2\pi(1-\nu)} i\pi \operatorname{sign}(k) F[\kappa' \Delta u']
\end{align*}
\tag{H3}
$$

## H3 First order

The first-order solutions are obtained from the full solution by making the small slope approximation and keeping only the terms that are proportional to the fault slope.

### H3.1 Space domain

$$
\begin{align*}
\sigma_{11}^{1}(x_1) &= \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ 3\frac{x_2-y_2}{(x_1-y_1)^2} - \frac{m(y_1)}{x_1-y_1} \right] \frac{d}{dy_1} \Delta u'(y_1) dy_1 \\
\sigma_{22}^{1}(x_1) &= \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ -\frac{m(y_1)}{x_1-y_1} - \frac{x_2-y_2}{(x_1-y_1)^2} \right] \frac{d}{dy_1} \Delta u'(y_1) dy_1 \\
\sigma_{12}^{1}(x_1) &= \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{m(y_1)}{x_1-y_1} - \frac{x_2-y_2}{(x_1-y_1)^2} \right] \kappa'(y_1) \Delta u'(y_1) dy_1 \\
\tau_{el}^{1}(x_1) &= -m(\sigma_{11}^{0} - \sigma_{22}^{0}) + \sigma_{12}^{1} \\
&= \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ \frac{m(y_1)}{x_1-y_1} - \frac{x_2-y_2}{(x_1-y_1)^2} \right] \kappa'(y_1) \Delta u'(y_1) dy_1 \\
\sigma_{el}^{1}(x_1) &= -2m\sigma_{12}^{0} + \sigma_{22}^{1} \\
&= \frac{\mu}{2\pi(1-\nu)} \int_{-\infty}^{+\infty} \left[ 2m(x_1) \frac{m(y_1)}{x_1-y_1} - \frac{m(y_1)}{x_1-y_1} - \frac{x_2-y_2}{(x_1-y_1)^2} \right] \frac{d}{dy_1} \Delta u'(y_1) dy_1
\end{align*}
\tag{H4}
$$
## H3.2 Spectral domain

$$
\begin{align*}
\sigma_{11}^{1}(k) &= \frac{\mu}{2\pi(1-\nu)} \left( \frac{3}{2\pi} F[h] * \left( VF \left[ \frac{d}{dy_1} \Delta u' \right] \right) - 3VF \left[ h \frac{d}{dy_1} \Delta u' \right] - UF \left[ m \frac{d}{dy_1} \Delta u' \right] \right) \\
\sigma_{22}^{1}(k) &= \frac{\mu}{2\pi(1-\nu)} \left( -UF \left[ m \frac{d}{dy_1} \Delta u' \right] - \frac{1}{2\pi} F[h] * \left( VF \left[ \frac{d}{dy_1} \Delta u' \right] \right) + VF \left[ h \frac{d}{dy_1} \Delta u' \right] \right) \\
\sigma_{12}^{1}(k) &= \frac{\mu}{2\pi(1-\nu)} \left( UF \left[ m\kappa' \Delta u' \right] - \frac{1}{2\pi} F[h] * \left( VF \left[ \kappa' \Delta u' \right] \right) + VF \left[ h\kappa' \Delta u' \right] \right) \\
\tau_{el}^{1}(k) &= \frac{\mu}{2\pi(1-\nu)} \left( UF \left[ m\kappa' \Delta u' \right] - \frac{1}{2\pi} F[h] * \left( VF \left[ \kappa' \Delta u' \right] \right) + VF \left[ h\kappa' \Delta u' \right] \right) \\
\sigma_{el}^{1}(k) &= \frac{\mu}{2\pi(1-\nu)} \left( \frac{1}{\pi} F[m] * \left( UF \left[ \frac{d}{dy_1} \Delta u' \right] \right) - UF \left[ m \frac{d}{dy_1} \Delta u' \right] \right. \\
&\qquad \left. - \frac{1}{2\pi} F[h] * \left( VF \left[ \frac{d}{dy_1} \Delta u' \right] \right) + VF \left[ h \frac{d}{dy_1} \Delta u' \right] \right)
\end{align*}
\tag{H5}
$$

# APPENDIX I: OUT-OF-PLANE (MODE III)

## I1 Full solution

In the following, $\gamma_1 = \frac{x_1-y_1}{r}$, $\gamma_2 = \frac{x_2-y_2}{r}$, and * represent the convolution operation. These equations are developed by using the definition of the strain $\epsilon_{kl} = \frac{1}{2} (\frac{\partial}{\partial x_k} u_l + \frac{\partial}{\partial x_l} u_k)$ and the stress ($\sigma_{ab} = c_{abkl}\epsilon_{kl}$) tensors, as well as eq. (B12), and the according slip mode eq. (B14). We also recall that the Fourier transform of $u(x) = 1/x$ and $v(x) = 1/x^2$, are respectively noted $U = F[u]$ and $V = F[v]$.

$$
\begin{align}
\epsilon_{13}(\mathbf{x}) &= \frac{1}{4\pi} \int_{\text{fault}} \frac{\gamma_2}{r} \frac{d}{dy'} \Delta u^s(\mathbf{y}) d\xi(\mathbf{y}) \nonumber \\
\epsilon_{23}(\mathbf{x}) &= -\frac{1}{4\pi} \int_{\text{fault}} \frac{\gamma_1}{r} \frac{d}{dy'} \Delta u^s(\mathbf{y}) d\xi(\mathbf{y}) \nonumber \\
\sigma_{13}(\mathbf{x}) &= \frac{\mu}{2\pi} \int_{\text{fault}} \frac{\gamma_2}{r} \frac{d}{dy'} \Delta u^s(\mathbf{y}) d\xi(\mathbf{y}) \nonumber \\
\sigma_{23}(\mathbf{x}) &= -\frac{\mu}{2\pi} \int_{\text{fault}} \frac{\gamma_1}{r} \frac{d}{dy'} \Delta u^s(\mathbf{y}) d\xi(\mathbf{y}) \nonumber \\
\tau_{el}(\mathbf{x}) &= \frac{\mu}{2\pi} \int_{\text{fault}} \frac{\gamma_1}{r} \frac{d}{dy'} \Delta u^s(\mathbf{y}) d\xi(\mathbf{y}) \nonumber \\
\sigma_{el}(\mathbf{x}) &= 0 \tag{11}
\end{align}
$$

## I2 Zeroth order

The zeroth-order solutions are obtained from the full solution by making the small slope approximation and keeping only the terms that are independent from the fault slope.

### I2.1 Space domain

$$
\begin{align}
\sigma_{13}^{0}(x_1) &= 0 \nonumber \\
\sigma_{23}^{0}(x_1) &= -\frac{\mu}{2\pi} \int_{-\infty}^{+\infty} \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u^s(y_1) dy_1 \nonumber \\
\tau_{el}^{0}(x_1) &= -\frac{\mu}{2\pi} \int_{-\infty}^{+\infty} \frac{1}{x_1 - y_1} \frac{d}{dy_1} \Delta u^s(y_1) dy_1 \nonumber \\
\sigma_{el}^{0}(x_1) &= 0 \tag{12}
\end{align}
$$
## 12.2 Spectral domain

$$ 
\begin{aligned}
\sigma_{13}^{0}(k) &= 0 \\
\sigma_{23}^{0}(k) &= -\frac{\mu}{2\pi}\pi|k|F[\Delta u^s] \\
\tau_{el}^{0}(k) &= -\frac{\mu}{2\pi}\pi|k|F[\Delta u^s] \\
\sigma_{el}^{0}(k) &= 0
\end{aligned}
\qquad (13) $$

## 13 First order

The first-order solutions are obtained from the full solution by making the small slope approximation and keeping only the terms that are proportional to the fault slope.

### 13.1 Space domain

$$ 
\begin{aligned}
\sigma_{13}^{1}(x_1) &= \frac{\mu}{2\pi} \int_{-\infty}^{+\infty} \frac{x_2 - y_2}{(x_1 - y_1)^2} \frac{d}{dy_1} \Delta u^s(y_1) dy_1 \\
\sigma_{23}^{1}(x_1) &= 0 \\
\tau_{el}^{1}(x_1) &= 0 \\
\sigma_{el}^{1}(x_1) &= 0
\end{aligned}
\qquad (14) $$

### 13.2 Spectral domain

$$ 
\begin{aligned}
\sigma_{13}^{1}(k) &= \frac{\mu}{2\pi} \left( \frac{1}{2\pi} F[h] * \left( VF \left[ \frac{d}{dy_1} \Delta u^s \right] \right) - VF \left[ h \frac{d}{dy_1} \Delta u^s \right] \right) \\
\sigma_{23}^{1}(k) &= 0 \\
\tau_{el}^{1}(k) &= 0 \\
\sigma_{el}^{1}(k) &= 0
\end{aligned}
\qquad (15) $$