modeled as an open surface S (an open curve in the two-dimensional case) across which 18 the elastic displacement jumps, while the normal stress, the traction, remains continuous. The jump in the displacement, a vector field g on S called the slip, measures how much 20 the rock has slipped on each side of the fault surface.

The forward or direct problem consists in determining the elastic displacement outside which we will call the acquisition manifold (see Fig. 1). The surface displacement can 26 be measured from satellite interferometric data or data from global positioning systems 27 (GPS). (For a general introduction to the geophysical model we refer to [11,25] and 28 references therein.)

While the study of elastic dislocations is by now classic, starting with the seminal work 30 of Volterra and Somigliana at the turn of the last century (see [34] and the historical 31 perspective in [12]), and elastic dislocations have been used to model creeping tectonic 32 faults during the interseismic period, especially at subduction zones in the Earth's crust 33 (among the vast literature, we mention [15,30,35]), until recently most of the available 34 results in the literature dealt with an idealized setting, in which the region of interest in 35 the upper layer of the Earth's crust is modeled as an infinite half space (or half plane in 36 the two-dimensional case) with constant elasticity parameters, the fault has the simple 37 geometry of a rectangle (or segment in the two-dimensional case) with either horizontal or 38 oblique orientation, and the slip is constant. In this specialized case then Okada obtained 39 an explicit formula for the elastic displacement [22], using the elastic Green's function in 40 a half space derived by Mindlin [20].

In recent work, it has been observed that the forward problem can be recast as a singu-42 lar source problem or as a nonhomogeneous transmission problem that can be addressed 43 using variational techniques or potential theory depending on the regularity of the coef-44 ficients, of the data, and of the domain, including the fault. This has led to establishing The minimization of the misfit functional follows a steepest descent algorithm, which we implement only in the two-dimensional setting. The numerical experiments tend to be expensive computationally, because two boundary-value problems for the system of elasticity, the forward and the adjoint problem, must be solved at each step of the descent. Therefore, as a proof of concept, we present the numerical results in a simplified setting, that of a square as the domain and where the fault S is a line segment. Furthermore, we mostly focus on the reconstruction of the fault surface assuming the slip is known. Indeed, the inverse problem is linear in the slip vector and hence stable, while the inverse problem for the fault surface is nonlinear. We then borrow the approach developed in [9] for reconstructing polygonal partitions in EIT to do the descent.

In this approach, the vertices of the dislocation line are assumed to be nodes of a coarse domain partition and are moved individually using the shape derivative calculated along linear elements on the partition. The optimal position of each vertex is computed separately and used to update the dislocation segment. The advantage of this method is that it is local to each node of the partition; hence, it can be parallelized to increase efficiency. The drawback is its accuracy. To calculate the shape derivative, we employ a DG method to compute both the forward and the adjoint problem on a finer conforming mesh. The steepest descent algorithm is stopped when a given tolerance is reached. As with most iterative methods, we have a good reconstruction only if the initial guess is sufficiently close to the true data for the problem (the fault surface and the slip vector).

Moreover, the limited accuracy of our approach, the inherent sensitivity of the shape derivative to small changes, and the underlying ill-posedness of the inverse problem make the reconstruction algorithm numerically less stable. Even though we rigorously derive the shape derivative only when the slip vanishes at the boundary, the discretized numerical algorithm can be applied in more generality (see also Remark 2.4). In our numerical experiments, as expected, we find that if the slip vector is constant, hence singularities appear in the solution at the tips of the dislocation segment, the reconstruction is more accurate, since the singularities help in locating the fault boundary. On the other hand, the error is significantly larger if the slip vector is allowed to vanish at the endpoints (see [7] where a reconstruction algorithm for linear cracks in dimension two is presented and tested). Consequently, in this latter situation we find that we may have to make more than one measurement or measure on a sufficiently large portion of the boundary, while from the uniqueness result for the inverse problem one would expect that one measurement on a small open patch at the boundary should be enough to reconstruct.

We do not address the stability of the reconstruction, nor the convergence of the method analytically, issues that we reserve to address in future work. In the case of planar faults and homogeneous elastic parameters, Lipschitz stability for the inverse fault problem was proved in [27], albeit leading to a nonquantitative estimate, which therefore cannot be directly employed to establish convergence. Quantitative Lipschitz stability estimates for a linear crack were obtained in dimension 2 in [8].

We close the Introduction with a brief overview of the paper. In Sect. 2, we discuss the set up for both the forward and the inverse problem and recall the well-posedness results for the forward problem, as well as the uniqueness results for the inverse problem in [4,6]. Section 3 is devoted to the rigorous derivation of the shape derivative. Lastly, in Sect. 4, we discuss the numerical tests and their outcomes.

We start by discussing the set-up for the forward problem and the assumptions we make on the fault surface and on the slip vector. We follow the notation and conventions in [4,6]. We only discuss the three-dimensional case with obvious modification in the twodimensional case, which is the setting for our numerical experiments.

We start by introducing some notation. We employ the blackboard font, e.g., A, for tensors, the boldface font, e.g., A, for matrices, bold italics font, e.g., a, for vectors, and italics font, e.g., a, for scalars. We denote matrix-vector multiplication simply by A b, while

• and : represent the inner product between vectors and matrices, respectively. Lastly, we use Einstein notation of summation over repeated indices, e.g.,

we similarly define the product of a tensor with a vector A a.

As in [4], we are concerned with faults buried in the Earth's crust that undergo a creeping aseismic movement. Such a situation is typical of areas of microseismicity, but can lead to earthquakes of sizeable magnitude [14]. In this context, faults tend to be shallow, and it is reasonable to work in a patch of the Earth's crust, modeled as a bounded domain with Lipschitz boundary. In the numerical implementation, is of polygonal type.

The boundary of , ∂ , is subdivided into two parts, one that represent the buried part and the other, ∂ \ , that lies on the surface of the Earth. We assume that is a closed subset of ∂ .

We will also employ a less standard space on S, the space H 1 2 00 (S), which can be defined as the weighted space:

where δ(x) =dist(x, ∂S) for x ∈ S and the distance is computed using geodesic length on S. The space H Finally, if f defined a.e. on is regular enough to have a nontangential limit on S, we let [f ] S denote the jump on S, defined as:

where f ± is the limit taken nontangentially in ± , respectively.

183 In the regime of creeping faults, where notable slip occurs without earthquakes, a quasi-184 static approximation is justified. Therefore, we model the Earth's crust as a linearly elas-185 tic, but inhomogeneous, material with time-independent elastic parameters. The elastic 186 response is encoded in a fourth-order tensor C, the elasticity tensor, which we assume 187 to be totally symmetric (that is, the material is hyperelastic) and uniformly Lipschitz 188 continuous in : 189 C ∈ W 1,∞ ( ). (2.3) 190 In fact, we assume that C is isotropic for the inverse problem. This regularity assumption 191 is not optimal in the context of geophysical applications, since the Earth's crust is typically 192 modeled as a layered medium, where material properties can vary abruptly from one layer 193 to the next. However, Lipschitz regularity is needed to compute the shape derivative, and 194 for the solvability of the forward problem in some cases as well. We discuss this point 195 further later in this section. 196 We also impose the standard (uniform) strong convexity condition, which reads 197 C(x)A : A ≥ γ |A| 2 , a.e in , (2.4) 198 for all symmetric matrices A for some constant γ > 0. This condition guarantees that the 199 energy functional is strictly positive. 200 The forward or direct fault problem consists in finding the elastic displacement u that 201 solves the homogeneous elastostatics system subject to suitable boundary conditions on 202 ∂ and suitable transmission conditions on S. On the part of ∂ that lies on the surface 203 of the Earth ∂ \ , it is appropriate to assume that there is no elastic load; hence, ∂ \ 204 is traction free. The traction is the normal component of the stress σ = C ∇u, with ∇ the 205 symmetric part of the gradient, and represents the elastic force acting on surfaces. 206

On the buried part , there are several boundary conditions that can be imposed.

For simplicity, we choose to impose homogeneous Dirichlet conditions, that is, is 208 displacement free. This choice is physically motivated by our assumption that the fault is 209 at a positive distance from the boundary, so the boundary is undisturbed by the slip on the 210 fault. For numerical purposes, absorbing boundary conditions could also be considered.

On the fault surface S, we assume continuity of the traction, while there is a jump in the 212 displacement, encoding the slippage of the rock on the two sides of the fault, given by a 213 vector field g on S, the so-called Burger's vector. The slip needs not be tangential to S and 214 can be oriented arbitrarily.

The displacement u must, hence, satisfy the following mixed-boundary-value-interface 216 problem:

[u] S = g, (C ∇u)n S = 0.

(2.5)

The strong convexity condition (2.4) implies strong ellipticity for the elastostatics operator As a matter of fact, even when and S are smooth, if g does not vanish at the boundary of S, i.e., if g ∈ H 1 2 (S), but g / ∈ H 1 2 00 (S), then singularities develop at ∂S that preclude any variational formulation. Instead, owing to the Lipschitz regularity of the elastic coefficients, provided the boundary ∂ is smooth enough, a very weak solution to (2.5) can be constructed via a duality argument, viewing the mixed-boundary-value interface problem as a singular source problem in the whole of with source

where δ S is the unit measure concentrated on S, so that f S ∈ H -3/2-, > 0. To use duality, one needs to have sufficiently regular solutions of the adjoint problem, which is a problem of the same form in our case. This is possible, since regularity holds for the system of elastostatics when is sufficiently smooth. In particular, if div(C ∇v) = F in and F ∈ L 2 ( ), then v ∈ H 2 ( ). Furthermore, it can be shown using a layer potential representation that the solution is more regular than obtained via duality, namely, it belongs to H 1-( \ S) for any > 0. Hence, the homogeneous Dirichlet condition on holds in trace sense. Somewhat informally, by a very weak solution of problem (2.5) we then mean that u satisfies

for all v ∈ D( ), where as customary D( ) is the space of smooth functions that are compactly supported in , and , (X ,X) denotes the duality pairing between a Banach space X and its dual X . In reality, one needs to work with the dual space of a suitable subspace of H 1 ( ).

For a precise definition of very weak solution and more details on the dual formulation of problem (2.5), we refer to [6]. There the domain is the lower half-space and weights must be imposed to ensure integrability at infinity. However, essentially the same arguments apply in bounded domains, without the need of weights, again if the boundary is sufficiently regular (e.g., of class C 2,α , α > 0). We recall the following well-posedness results for very weak solutions (see [6,Theorem 3.11]).

There exists a very weak solution

2 (S) and any > 0.

By contrast, when g ∈ H 1 2 00 (S), extending the slip by zero to yields a standard nonhomogeneous interface problem that admits a variational solution in H 1 ( \ S). The variational formulation only requires the elastic parameters to be bounded uniformly.

We now recall this variational or weak formulation, which will be the starting point for computing the shape derivative. Below, g ∈ H

1 2 ( ) denotes the extension of g ∈ H 1 2 00 (S)

by zero. We recall the following standard lemma (see e.g., [4], Lemma 3.3, and references therein).

The following identification holds:

Then, u ∈ H 1 ( \ S) is a weak or variational solution of the mixed-boundary-valueinterface problem (2.5) if, for all v ∈ H 1 ( ):

and [u] = g. We note that by the symmetry properties of C, the bilinear form a defined above can be written equivalently as:

The following well-posedness result holds (see [4, Theorems 3.5 and 3.6]).

Theorem 2.3 There exists a unique weak solution u ∈ H 1 ( \ S) to problem (2.5), which depends continuously on the slip vector g ∈ H 1 2 00 (S).

Very briefly, the proof consists in solving two problems, a mixed Dirichlet-Neumann problem in + and a purely Neumann problem in -, with an unknown vector potential φ on , which represents the common value of the traction on the two sides of , and subsequently solving for φ from knowledge of the slip field g using a combination of Neumann-to-Dirichlet operators for the domains ± . We refer the reader to [4] for more details and the complete proof.

We consider both the case g ∈ H 1/2 00 (S) and the general case g ∈ H 1/2 (S) for the numerical experiments in Sect. 4.

We present a detailed derivation of the shape derivative in the case g ∈ H 1 2 00 (S), that is, when the solution u ∈ H 1 ( \ S), since the variational setting is the best suited for numerical applications. Moreover, assuming that g vanishes at ∂S is more physically justified (only a patch on a fault is typically active). However, the boundary shape derivative can also be obtained using a layer potential representation of the elastic displacement u, which we established in [6], following the approach used in [8]. This derivation applies to any slip field g ∈ H Theorem 2.5 (Theorem 5.1 in [6]) Let S 1 , S 2 be two piecewise linear fault surfaces that in addition are graphs with respect to the same coordinate system. Let g i be bounded tangential fields in H 1/2 (S i ) with supp g i = S i , for i = 1, 2. Let u i , for i = 1, 2, be the unique very weak solution of (2.5) corresponding to g = g i and S = S i . If

then S 1 = S 2 and g 1 = g 2 .

The conditions on S i and g i in the theorem above can be weakened somewhat. Theorem 2.6 (Theorem 4.1 in [4]) Let S 1 , S 2 be two fault surfaces that in addition are 301 graphs with respect to the same coordinate system. Let g i ∈ H 1 2 00 (S i ), for i = 1, 2, with 302 Supp g i = S i , for i = 1, 2, and let u i , for i = 1, 2, be the unique weak solution of (2.5) 303 corresponding to g = g i and S = S i . If u 1 = u 2 , then S 1 = S 2 and g 1 = g 2 .

Briefly, the two theorems above follow both from an argument by contradiction, using 305 unique continuation for the elasticity system (this is why we have assumed that C is 306 Lipschitz continuous and isotropic, but see [5] for certain anisotropic cases). As a matter . The jump conditions allow only to conclude that w is an infinitesimal rigid motion, and 313 then are the geometric conditions imposed on S i and the conditions on g i , i = 1, 2, that 314 imply w can only be the trivial motion. Therefore, u 1 ≡ u 2 on the entire , a contradiction 315 in view of nonhomogeneous jumps g i on S i . Again, we refer to [4,6] for complete proofs.

We observe that a single measurement on is sufficient for the unique determination 317 of both the fault and the slip. The assumption that surfaces are graphs with respect to 318 a fixed, but arbitrary, coordinate system is not too restrictive in the geophysical context,

where nearby faults tend to follow the same orientation. For numerical purposes, we are 320 interested in polyhedral surfaces and, hence, both uniqueness theorems apply. In the rest of the paper we present an iterative reconstruction algorithm based on the 323 minimization of a suitable misfit functional. For ease of notation, we use to denote ≤ C,

where C is a constant that may depend on the data, i.e., the Lipschitz norm of C and the 325 strong convexity constant γ , on , S, and g, which are fixed in computing the derivative 326 (Fig. 2). where u solves (2.5) with fault surface S and slip vector g and σ (x) denotes surface area.

We remark that g depends also on S, so we have slightly abused notation. However, since the inverse problem is linear in the slip vector if S is fixed, the key step consists of reconstructing the fault surface S. Given the uniqueness results for the inverse problem and the well-posedness of the forward problem, J has a unique global minimum at S = S 0 and g = g 0 .

To find this minimum, we propose a steepest descent algorithm based on a shape derivative, defined as the derivative of J with respect to infinitesimal movements of the fault S and under infinitesimal changes of the slip vector g. Due to the ill-posedness of the inverse problem, we regularize the algorithm by minimizing only over piecewise polygonal or polyhedral faults. We do not study the convergence of this algorithm theoretically, in part because we do not know the convexity properties of J . Furthermore, there are measurement errors and we compute the solution of the forward problem numerically.

We discuss the influence of these two types of errors in the section on the numerical tests.

In this section, we tackle the computation of the shape derivative in terms of a material derivative, which gives the change of the solution u of (2.5) under infinitesimal changes of the fault S and the slip g. We encode the movement of the fault S in terms of a bi-Lipschitz diffeomorphism φ τ : → , depending on a small parameter τ ∈ [0, 1) with the following properties:

(1) φ 0 = I, and φ τ ∂ = I, where I is the identity map;

(2) φ τ is linear in τ .

Since S is at positive distance from ∂ , there exists δ > 0 such that dist(φ τ (S), ∂ ) > δ for all τ ∈ [0, 1). Therefore, we can assume without loss of generality that:

where U : → is a Lipschitz map with support in a neighborhood V of . (We recall that is at positive distance from ∂ , so V always exists.) Then:

Similarly, we encode the infinitesimal change of the slip g in terms of a vector field g τ , given by:

where h ∈ H 1 2 00 (S), and again denote with g τ its extension by zero to . The actual change in the slip is given by the composition of g τ with φ -1 τ .

For notational convenience, we set:

so that S 0 ≡ S and 0 ≡ . We recall that we assume that the elasticity tensor C is known in both the forward and the inverse problem; it is defined on the whole of and globally Lipschitz continuous.

We first obtain a distributed derivative, which can be computed in the framework of weak solutions of (2.5), assuming the slip g ∈ H 1 2 00 (S). We hence work with τ and the extensions g τ , h of the slip vectors g τ , h by zero to τ .

For ease of notation, we write J (τ ) := J ( τ , g τ ). We can now define the shape deriva-373 tive as the directional derivative:

375

Next, we let u τ ∈ H 1 ( \ S τ ) be the (weak) solution of the perturbed problem, which 376 exists and is unique by virtue of Theorem 2.6.:

)

or equivalently:

for every v ∈ H 1 ( ), where y = φ τ (x).

To compute the material derivative, we transform the perturbed problem (3.6) into a 383 problem on the fixed domain \ S. We then let

385

By the change of variable formula, which holds in Lipschitz domains, ȗτ ∈ H 1 ( \ S)

satisfies the variational problem:

It will be convenient to switch to a weak formulation with homogeneous transmission 397 conditions, so as to have the same function space for solutions and test functions. To this 398 end, we introduce a lift of the slip vector on using a Newtonian potential. Given any 399 vector field a ∈ H 1 2 ( ), we define its lift ℓ a ∈ H 1 ( \ ) as any solution to the problem: lift in H 1 ( \ ) can be obtained using a double-layer potential representation for the solution, which is well-known for the Laplacian (alternatively, one can apply the same techniques used to prove Theorem 2.3). We let:

Then u τ ∈ H 1 ( \ S) and satisfies:

since [ ȗτ ] = g τ = g + τ h. Furthermore:

We now define the material derivative u ∈ H 1 ( ) as:

where the convergence is strongly in H 1 ( ), if the limit exists.

Before tackling the existence of the material derivative, we discuss some useful properties of C and φ τ . We recall that φ τ and C are uniformly Lipschitz on . Since φ τ = I + τ U, there exists τ small enough, say τ ∈ [0, τ 0 ], for some τ 0 1, such that φ -1 τ ≈ Iτ U and lim

, and a straightforward calculation shows that lim

also in L ∞ ( ). Additionally:

and similarly for C • φ -1 τ . For ease of notation, we let

We begin with showing that a weak limit exists. We observe that, since the lift of the jump on satisfies a linear problem, ℓ τ h = τ ℓ h . Consequently, l h = ℓ h .

Lemma 3.1 Let u be the unique weak solution of (2.5) with g ∈ H 1 2 00 (S) and let u τ be given in (3.10). Then u τ -u τ converges weakly in H 1 ( ) and strongly in L 2 ( ) as τ → 0 + . Remark 3.2 Although convergence will be along subsequences (not relabeled), the whole sequence will converge as we will show uniqueness of the limit later on.

Proof By Poincaré's inequality, it is enough to establish a uniform bound on

We recall that [Dφ τ ] -1 ≈ Iτ DU for τ ∈ [0, τ 0 ], so that

Next, using the strong convexity assumption (2.4) and the symmetry of C, we have that for a.e. x ∈ ,

It then follows from Korn's inequality that 1

Next we note that, since u τ -u τ ∈ H 1 ( ), it can be used as a test function in (2.7) for u and in (3.7) for ȗτ . Hence, for τ > 0

We then write using the fact that we can integrate u τ -u equivalently on or \ or

We tackle the terms A and B separately, exploiting (3.16). For the first term, we have:

For the second term, we recall the definition of u τ (3.10) to obtain:

For the first term on the right, using the second equation in (3.16):

Hence, for τ ∈ [0, τ 0 ] we have the bound:

473

For the second term on the right, we have similarly

) 475 using also that by linearity / τ h /τ = / h . From (3.15), (3.17),(3.18),(3.19), we then obtain 476 the desired uniform bound: 477

)

where we have employed also the stability estimates for problems (2.5), (3.6), and (3.9).

We hence have that, upon passing to a subsequence, ȗτu τ converges weakly to 481 an element v ∈ H 1 ( ). By Rellich Compactness Theorem, the subsequence converges 482 strongly in L 2 ( ) (in fact, in H s ( ) for s < 1), and by weak continuity of the derivative:

485 Remark 3.3 For the validity of (3.20), it is not necessary to have C globally Lipschitz in .

It is sufficient that C is uniformly Lipschitz in a neighborhood of , which also contains 487 τ for τ sufficiently small.

We next establish uniqueness of the weak limit and strong convergence, by character-489 izing the material derivative as the solution of a certain variational problem. Remark 3.5 In general, the 4th-order tensor B does not enjoy the same symmetry properties as C, which reflects the fact that the system of elastostatics is not covariant under changes of coordinates.

Proof We begin by showing that any weak limit of u τ -u τ satisfies (3.21). This fact will also ensure its uniqueness.

We start by writing a τ ( u τ -u τ , w) for a generic test function w ∈ H 1 ( ). We proceed as in the computation of terms A and B in the proof of Lemma (3.1), again exploiting the weak formulation (3.7) for ȗτ :

For notational convenience, we write the above identity abstractly as:

Our goal is to pass to the limit τ → 0 + in each term. In all three terms, we will use the assumptions on C and φ τ , in particular (3.13) and (3.14). We also recall that the composition of two Lipschitz maps is still Lipschitz with constant that can be bounded by the product of the two Lipschitz constants, which implies in particular that C τ → C in W 1,∞ ( ) as τ → 0. Then,

a.e. x ∈ and the limit is in L ∞ ( ), because both U and C are Lipschitz.

We first investigate the behavior as τ → 0 + of the integral G τ on the lefthand side. By Lemma 3.1 the left-hand side consists of the pairing between an L 2 weakly convergent (sub)sequence ∇u τ -∇u τ and an L 2 strongly convergent sequence

. This last expression indicates the action of a 4th-order tensor onto a 2nd-order tensor yielding a 2nd-order tensor (up to the scalar factor J (φ τ )), and strong convergence follows from Hölder's inequality. Then, by weak-strong convergence:

Next, we tackle the first integral E τ on the right-hand side. We split into three parts, taking advantage of the fact that u solves (2.7):

τ dx

Using once again the assumptions on C and φ τ , the regularity of u and w, and Hölder's inequality, we can pass to the limit τ → 0 + in each term above, as follows:

where we have also employed (3.13)- (3.14) and defined B by (3.22). Convergence in the last term F τ follows straightforwardly from the convergence of φ τ to I and of C τ to C in W 1,∞ ( ): We can now show that the convergence is actually strong in H 1 ( ). We first observe that, since u ∈ H 1 ( ), we can equivalently integrate on the right-hand side of (3.21) over or \ , and we can use u as test function w. Hence, we can write (3.21) in compact form as:

where a is defined below (2.7). Next, similarly to (3.15), by Poincaré's inequality and the 564 strong convexity of C, we have the estimate:

with a τ as in (3.7).

We will show that the right-hand side of the expression above goes to zero as τ → 0 + .

To this end, we split it into three parts:

The convergence of C τ to C and of φ τ to I in W 1,∞ ( ) again implies that

.28) 576 From (3.28), we immediately have that 577 a τ ( u, u) -→ τ →0 + a ( u, u), (3.29) 578

Also, the first part of the proof using u as test function w readily yields:

We cannot directly pass to the limit

, because we only 581 have weak convergence of the arguments at this point. However, we recall the second 582 identity in (3.16), which allows us to rewrite this term as:

.34) 607 Finally, by employing (3.29), (3.30), and (3.34) in (3.27), we can conclude that:

The proof is now complete.

We recall that uτ = ȗτℓ τ h and that ℓ τ h τ = ℓ h is constant in τ . Then, from Theorem where w ∈ H 1 ( ) is the unique weak solution of the adjoint problem:

where we used the definition of the material derivative (3.35).

Let now w solve the adjoint problem (3.38) (this is a standard variational problem that admits a unique weak solution by the Lax-Milgram Theorem, for instance). Then, we can equivalently write the expression above on the right as:

We recall that we can integrate over \ or \ S, because neither u nor w jump across \ S, and that splits into an outer region + that touches ∂ and an inner core -.

From (3.38), although u cannot be used as a test function, integrating first in + and then in -gives that

recalling that by convention n is the outer normal to -. But we know that the jump

[u] = h (so, in fact, u only jumps across S), so that we finally have:

where we also used that w can be taken as test function for (3.36).

We close this section by deriving a boundary shape derivative from the distributed shape derivative in case C is constant, i.e., the rock is homogeneous in the region of interest around the fault. While this formula, which contains only integrals over the fault S, may lead to a less stable numerical implementation, it allows a more transparent geometric interpretation. We work with strong solutions of (2.5) and (3.38), assuming that ∂ , S and g, h are regular enough. By a strong solution u of (2.5), we mean that u ∈ H 2 ( \S) ∩ H 1 ( \S) and that the problem (2.5) is attained at least pointwise a.e.

Similarly, a strong solution w ∈ H 2 ( ) ∩ H 1 ( ) solves (3.38) at least pointwise a.e. . From standard regularity theory, weak solutions are strong solutions if ∂ and S are at least of class C 1,α , α > 0, and g, h ∈ H s 0 (S), s ≥ 3/2. We only compute the derivative in 2 space dimensions to simplify the resulting expression.

Our starting point is formula (3.39). We show that we can write M[∇u] : ∇w as the divergence of a certain vector b. Since we assume that C is constant, B ≡ 0. Furthermore, since both u and w satisfy the homogeneous system of elasticity a.e. in \ , it is enough to show that:

Therefore, using (3.49) and (3.50) in the second and third terms on the right-hand side of the last equality in (3.48), we obtain

Inserting this expression into (3.48) finally gives formula (3.47).

In this section, we implement a reconstruction algorithm based on a gradient descent method that uses the shape derivative (3.37) of the functional defined in (3.1). While the algorithm is applicable in both two and three space dimensions, our numerical experiments are done in the two-dimensional case, and hence we focus on describing the method in 2D. The two-dimensional problem can be seen as looking at vertical slices of the full 3D problem and neglecting any stress across different slices. This is a simplification, but our numerical experiments are a proof of concept that the algorithm can recover the fault and the slip under suitable assumptions. We reserve to address a more realistic set-up for the numerical implementation in future work.

The reconstruction procedure is based on the gradient descent algorithm proposed in [9] to reconstruct a piecewise-constant conductivity that jumps across a polygonal partition in a composite material from current-to-voltage boundary measurements. There are significant differences between the case considered in [9] and here. First of all, we deal with a system of PDEs and not a scalar PDE. Furthermore, we have nonhomogeneous jumps across the fault. In particular the solution is not H 1 regular across the fault. Lastly, we measure only on part of boundary and we measure only once (twice in some unstable cases). These additional challenges make the numerical implementation more sensitive, especially to how the jump in the displacement u across the fault is treated.

We concentrate primarily on recovering the fault geometry S, assuming the slip field g is known. Indeed, once S is determined, the problem of reconstructing g is linear. We will consider both situations addressed in our prior works [4,6], that is, the case when g vanishes at the endpoints of the line S and is in H 1 2 00 (S), and the case that g is constant on S, even though the theoretical framework for the distributed shape derivative exploits a variational formulation and hence applies to the first case.

To solve both the forward and adjoint problems (see (2.5) and (3.38)), we employ the discontinuous Galerkin (DG) method, which is nonconforming and lends itself to tackling problems with inhomogeneous jumps. Other numerical methods have been used for shape optimization. In the context of elasticity, among the many results we mention the recent works [13,18], where the boundary immerse method and immersed interface finite elements are employed.

We next discuss the DG method in the context of our problem. Fig. 4 An example of the meshes employed for generating measurements and for solving the optimization problem

and homogeneous Dirichlet conditions on the fourth side identified with the set = {(x, y) y = 1, -1 ≤ x ≤ 1} (See Fig. 3).

We consider a family of shape-regular partitions T h where 0 < h 1. These partitions consist of nonoverlapping triangles K such that = ∪ K∈T h K. Specifically, h = max K∈T h h K , where h K = diam(K). The mesh conforms to the fault S as follows. We construct a trapezoidal region inside with one side lying on the displacement-free part of ∂ and one side that agrees with the fault segment S. We then align the computational grid with S (see Fig. 4a).

We define F I as the set of all interior sides. An interior side γ in F I is characterized by the existence of two adjacent elements K + and K -in T h such that γ = ∂K + ∩ ∂K -. For clarity in presenting the DG formulation, we distinguish between the sides that belongs to S from the others, which belong to F I .

Similarly, we label the sides of the elements that touch the boundary ∂ according to Then the DG approximate solution u h ∈ V r h of the forward problem (2.5) must satisfy 885 (4.5) for all v ∈ V r h , see [1-3,23,24]. We use a similar penalized formulation to compute 886 the DG solution w h ∈ V r h of the adjoint problem (3.38). The numerical implementation 887 then follows the standard approach for finite elements methods. 888 In Fig. 5, we plot an example of the numerical solution for the forward problem (2.5), 889 obtained by using linear elements, that is, polynomials of degree 1 on each element K , 890 and setting the mesh size h = 0.01. We consider the case of a horizontal dislocation with 891 vertices at (-0.4, 0) and (0.4, 0) and a slip field vanishing at the endpoint of the dislocation 892 segment given by 893 g(s) = 10 -100 39 12 s 12 + 1 χ (-39 100 , 39 100 ) (s), 20 -100 39 12 s 12 + 1 χ (-39 100 , 39 100 ) (s) , 894 (4.6) 895 4.2 Reconstruction algorithm 896

For the reconstruction algorithm, we adopt the gradient descent approach introduced in 897 [9]. Below, we briefly recall the method as it applies to our problem.

To highlight the dependence of the shape derivative on the map U ∈ W 1,∞ ( ), identified Consequently, for notational convenience we simply write 908 dJ dτ τ =0 = DJ (S), U .

We also use ((, )) to denote the L 2 -inner product.

We also continue to specialize to the case that the dislocation line S is a straight segment, 911 so it is uniquely determined by its two vertices P l ∈ R 2 , l = 1, 2.

A key step in the descent algorithm is the computation of the descent direction at each 913 step. Let X ⊂ W 1,∞ ( , R 2 ) be a suitable subspace. Then, at iteration level k ∈ N, we solve for every δθ ∈ X, where θ k ∈ X denotes the descent direction and S k is the dislocation 917 line computed at the k-th iteration. The subspace X is chosen so that S k is a segment 918 contained in for each k. For solving (4.7), we discretize the problem and employ the 919 idea already contained in [9]. That is, given that the dislocation line is uniquely defined by 920 its vertices, we compute the descent direction for each vertex individually before updating 921 the partition at each iteration k.

922 Specifically, we recall that we have two meshes. A coarse mesh used for the measure-923 ments and a fine mesh used for the DG approximation. The vertex of the dislocation 924 segment is assumed to be on both meshes. At iteration level k, for each vertex P k l , l = 1, 2, 925 of S k , we calculate the gradient of the functional J with respect to the position of the 926 vertex P k l as follows. We select two maps

with the properties 927 that:

(1) U k l,1 , U k l,2 have support strictly contained within .

(2) U k l,1 and U k l,2 are piecewise linear along the sides of the (coarse) partition and satisfy 930 U k l,1 (P k j ) = (δ jl , 0), U k l,2 (P k j ) = (0, δ jl ),

where δ jl denotes the Kronecker delta and P k j denotes any vertex of the mesh.

Then, the gradient of J with respect to the position of the vertex P k l is given by:

934 and hence the steepest descent at the k-th iteration is represented by the vector

936

For the DG approximation, we define the maps U k l,1 and U k l,2 as follows:

where ϕ k l is the hat function associated with the node P k l of the coarse mesh (see, for Solve the forward problem (2.5)

Compute gradient of the solution.

6:

Extract boundary data for adjoint problem 7:

Solve the adjoint problem (3.38) 8:

Compute gradient of the adjoint solution. 9:

Calculate shape derivatives using U k 1 and U k 2 .

10:

Update the vertices of the dislocation, using (4.8), i.e., P k+1 l = P k lαθ k l .

11:

if k > nIter and max(θ k l ) < tol then where ε is a uniform random variable drawn from the interval (-a, a), with a > 0 deter-955 mined based on the desired noise level (unless specified otherwise, we set a = 0.0007). In 956 the numerical implementation, the noise is added pointwise to each node of the coarse 957 mesh used for the measurements. To quantify the noise level, we compute the relative L 2 As expected, the reconstruction with one measurement is poor, and two measurements 986 are needed in this case.

We include the results for the following numerical tests:

988 Test 1: Horizontal dislocation with constant slip field g (Fig. 6). Test 2: Horizontal dislocation with slip field g vanishing at the endpoints (Fig. 7) 990 Fig. 7 Reconstruction of a horizontal dislocation with vertices at (-0.4, 0) and (0.4, 0), using only one boundary measurement with 0.86% noise level. The initial guess for the reconstruction algorithm is represented by the black dashed line. The slip field g is chosen as in (4.9) Test 3: Oblique dislocation with constant slip field g (Fig. 8) 991