Discontinuities without discontinuity:The Weakly-enforced Slip Method

# Discontinuities without discontinuity: The Weakly-enforced Slip Method

G.J. van Zwieten111email: g.j.v.zwieten@tue.nl, E.H. van Brummelen, K.G. van der Zee,
M.A. Gutiérrez, R.F. Hanssen

Eindhoven University of Technology, Department of Mechanical Engineering, P.O. Box 513, 5600 MB Eindhoven, The Netherlands Eindhoven University of Technology, Department of Mathematics & Computer Science, P.O. Box 513, 5600 MB Eindhoven, The Netherlands Delft University of Technology, Department of Mechanical, Maritime and Materials Engineering, P.O. Box 5058, 2600 GB Delft, The Netherlands Delft University of Technology, Department of Civil Engineering and Geosciences, P.O. Box 5048, 2600 GA Delft, The Netherlands
###### Abstract

Tectonic faults are commonly modelled as Volterra or Somigliana dislocations in an elastic medium. Various solution methods exist for this problem. However, the methods used in practice are often limiting, motivated by reasons of computational efficiency rather than geophysical accuracy. A typical geophysical application involves inverse problems for which many different fault configurations need to be examined, each adding to the computational load. In practice, this precludes conventional finite-element methods, which suffer a large computational overhead on account of geometric changes. This paper presents a new non-conforming finite-element method based on weak imposition of the displacement discontinuity. The weak imposition of the discontinuity enables the application of approximation spaces that are independent of the dislocation geometry, thus enabling optimal reuse of computational components. Such reuse of computational components renders finite-element modeling a viable option for inverse problems in geophysical applications. A detailed analysis of the approximation properties of the new formulation is provided. The analysis is supported by numerical experiments in 2D and 3D.

Keywords: Volterra dislocation, Finite Element Method, weak imposition, linear elasticity, tectonophysics.

## 1 Introduction

The world is perpetually reminded of the fact that seismic hazard is still beyond reach of prediction — as it was most recently by the disaster that struck Japan. The difficulty is not just to predict the exact moment of failure, which, as argued by some [9], might never reach a level of practicality. It is also the nature of the risk, and the extent to which stress is accumulating, that turns out to be surprisingly difficult to constrain. The 2011 Tohoku-Oki earthquake demonstrated a great lack of understanding of ongoing tectonics [12]. Arguably, a better understanding could have reduced the secondary effects if such information had led to more apt measures and regulations.

The main reason for this poor state of information can be traced to the absence of direct measurements. The primary quantities of interest, being the magnitude and orientation of the stress tensor in the earth’s crust, can be obtained only through tedious, expensive, point-wise measurements. A viable broad scale method to directly measure the global stress field does not exist. For this reason information is obtained mostly from secondary observables, earthquakes themselves being an important source. Earthquakes represent significant, near instantaneous changes in the global stress field. By accurately determining the location of the segment of the fault that collapsed, one can progressively update the stress field and evolve it in time. This way the tectonic evolution is monitored, and hazardous areas can be identified as regions where stress accumulates. For successful tracking of stress development, however, it is essential to understand the tectonic mechanism behind every earthquake. This includes the location and geometry of the section of the fault that collapsed, and the direction and magnitude of fault slip. It is increasingly popular to base such analyses on local co-seismic surface displacements. This type of information has become available since the nineties with the advent of space borne interferometric SAR measurements of the earth’s surface, and with the widespread availability of GPS measurements [19]. Analysis of this data has in recent years seen rapid adoption and is now routinely performed for all major earthquakes.

A mechanical model is required to connect observations to physics. Most commonly (if not exclusively) used is an elastic dislocation model, based on the assumption that on short time scales, nonlinear (plastic) effects are negligible. The model embeds a displacement discontinuity of given location and magnitude in an elastic medium, causing the entire medium to deform under the locked-in stress. Many different solution methods have been developed for this particular problem, based on analytical solutions or numerical approximations; see for instance [21] for an overview of the most prominent methods. However, methods founded on analytical solutions generally dictate severe model simplifications, such as elastic homogeneity or generic geometries, which restricts their validity. The computational complexity of methods based on numerical approximation, on the other hand, is typically prohibitive in practical applications. Because in practice the surface displacements are given, and the dislocation parameters are the unknowns, the computational setting is always that of an inverse problem. A typical inversion requires several thousands of evaluations of the forward model, and therefore computational efficiency is a key requirement. Moreover, the forward problems in the inversion process are essentially identical, except for the fault geometry. Reuse of computational components, such as approximate factors of the system matrix, is imperative for efficiency of the inversion. Current numerical methods for seismic problems do not offer such reuse options.

Finite-element methods provide a class of numerical techniques that are particularly versatile in terms of modeling capabilities in geophysics. Finite-element methods allow for elastic heterogeneity, anisotropy, and topography; all things that can not well be accounted for with currently used analytical and semi-analytical methods. In geophysical practice, finite-element methods are however often rejected for reasons of computational cost. The high computational cost can be retraced to the condition that the geometry of the fault coincides with element edges, which is a requirement engendered by the strong enforcement of the dislocation; see [14]. Consequently, the mesh geometry depends on the fault, which in turn implies that mesh-dependent components such as the stiffness matrix and approximate factorizations of that matrix cannot be reused for different fault geometries and must be recomputed whenever the geometry of the fault changes. The recomputation of these components in each step of a nonlinear inversion process leads to a prohibitive overall computational complexity.

To overcome the complications of standard finite-element techniques in nonlinear inversion processes in tectonophysics, this paper introduces the Weakly-enforced Slip Method (WSM), a new numerical method in which displacement discontinuities are weakly imposed. The WSM formulation is similar to Nitsche’s variational principle for enforcing Dirichlet boundary conditions [15]. The weak imposition of the discontinuity in WSM decouples the finite element mesh from the geometry of the fault, which renders the stiffness matrix and derived objects such as approximate factors independent of the fault and enables reuse of these objects. Therefore, even though the computational work required for a single realisation of the fault geometry is comparable to that of standard FEM, reuse of components makes WSM significantly more efficient when many different fault geometries are considered. This makes finite-element computations based on WSM a viable option for nonlinear inverse problems.

A characteristic feature of WSM is that it employs standard continuous finite-element approximation spaces, as opposed to the conventional FEM split-node approach [14] which introduces actual discontinuities in the approximation space. Instead, WSM approximations feature a ‘smeared out’ jump with sharply localized gradients. We will establish that the error in the WSM approximation converges only as in the -norm as the mesh width  tends to zero and that the error diverges as  in the energy norm. In addition, however, we will show that the WSM approximation displays optimal local convergence in the energy norm, i.e., optimal convergence rates are obtained on any subdomain excluding a neighborhood of the dislocation. The numerical experiments convey that WSM also displays optimal local convergence in the -norm.

The remainder of this manuscript is organized as follows. Section 2 presents strong and weak formulations of Volterra’s dislocation problem, and derives the corresponding lift-based finite-element formulation. Section 3 introduces the Weakly-enforced Slip Method based on two formal derivations, viz., by collapsing the support of the lift onto the fault and by application of Nitsche’s variational principle. In Section 4, we examine the approximation properties of the WSM formulation. Section 5 verifies and illustrates the approximation properties on the basis of numerical experiments for several 2D and 3D test cases. In addition, to illustrate the generality of WSM, in the numerical experiments we consider several test cases that violate the conditions underlying the error estimates in Section 4, such as discontinuous slip distributions and rupturing dislocations. Section 6 presents concluding remarks.

## 2 Problem formulation

In this section we define Volterra’s dislocation problem. We will postulate the strong formulation in 2.1, followed by a derivation of the weak formulation in 2.2. The latter will serve as a basis for the derivation of the Finite Element approximations, for which we lay foundations in Section 2.3.

### 2.1 The strong form

We start by defining the geometric setup. We consider an open bounded domain () with Lipschitz boundary . An -dimensional Lipschitz manifold , referred as the fault, divides the domain in two disjoint open subdomains and , such that . We equip with a unit normal vector directed into the subdomain . The fault supports a slip distribution , corresponding to a dislocation. The fault is referred to as a non-rupturing fault if the dislocation is compact in , and a rupturing fault otherwise. Let us note that in tectonics, rupturing faults correspond to intersections of the slip plane with the surface of the earth. Figure 1 illustrates this setup for . It is to be noted that need not be tangential to the fault. If has a non-vanishing normal component then the fault is opening, such as may be caused by an intruding material.

The displacement field generated by the dislocation is represented by . For convenience, we restrict our considerations to linear elasticity. We denote by the map the strain tensor corresponding to the displacement field , according to

 ϵ(u):=12[∇u+(∇u)T], (1)

under the assumption that is differentiable on . The constitutive behavior corresponds to Hooke’s law:

 σ(u):=A:ϵ(u), (2)

i.e., , where we adhere to the convention on summation on repeated indices. The tensors and are referred to as the stress tensor and the elasticity tensor, respectively. The elasticity tensor is subject to the usual symmetries . Moreover, we assume that it is bounded and satisfies a strong positivity condition, i.e., there exist positive constants and  such that:

 c–Aeijeij≤Aijkleijekl≤¯¯cAeijeij (3)

for all tensors . The elasticity tensor is in principle allowed to vary over the domain , subject to the above conditions. Auxiliary smoothness conditions on  will be introduced later.

To facilitate the formulation, we denote by the traction on a boundary corresponding to the displacement field . Moreover, we introduce the jump operator  and average operator according to:

 [[v]] :Γ∋x↦v+(x)−v−(x), (4a) {v} :Γ∋x↦12[v+(x)+v−(x)], (4b)

where and represent the traces of from within  and , respectively. Given a partition of the boundary into and such that , we define the Volterra dislocation problem as follows:

Strong formulation: given a body force , a displacement and a traction , and slip , find displacement field such that (5a) (5b) (5c) (5d) (5e)

Equation (5a) is the usual equilibrium condition, which applies everywhere in except on the manifold . Equations (5b) and (5c) respectively express that the displacements at the boundaries of and differ by the slip vector , and that the tractions at the boundaries of and are in static equilibrium, i.e., equal and opposite. It is to be remarked that this condition corresponds to a standard linear approximation in the small-slip limit, as the traction equilibrium at the fault occurs in fact in the deformed configuration. The boundary conditions  (5d) and (5e) correspond to Dirichlet and Neumann boundary conditions, representing a prescribed displacement and a prescribed traction, respectively.

To facilitate the presentation, we note that the solution to the Volterra dislocation problem (5) can be separated into a discontinuous part with homogeneous data and a continuous part with inhomogeneous data:

 (6)

The sum satisfies (5). Therefore, the inhomogeneous data in (5) can be treated separately in a standard continuous elasticity problem, and without loss of generality we can restrict our consideration to homogeneous data. We retain to identify the right member of (5a).

For rupturing faults, some compatibility conditions arise with respect to the boundary conditions. In particular, an intersection of the dislocation with the boundary of the domain is only admissible at the Neumann boundary . Otherwise, an inadmissible incompatibility between the jump condition  and the homogeneous Dirichlet condition  ensues.

### 2.2 The weak form

To derive the weak formulation of (5), we note that for any piecewise smooth function from into , we have the identities:

 ∫Ω∖Γv⋅divσ(u) =∮∂Ω−v⋅σn(u)−∫Ω−σ(u):∇v +∮∂Ω+v⋅σn(u)−∫Ω+σ(u):∇v =∮∂Ωv⋅σn(u)−∫Ω∖Γσ(u):∇v (7)

The first identity results from integration by parts. The second identity follows from a rearrangement of the boundary terms and

 ∫Γv+⋅σ+n(u)+v−⋅σ−n(u) =∫Γv+⋅σ+ν(u)−v−⋅σ−ν(u) =∫Γ(v+−v−)12(σ+ν(u)+σ−ν(u)) +∫Γ12(v++v−)(σ+ν(u)−σ−ν(u)) (8)

In the weak formulation, the admissible displacement fields will be insufficiently regular to ensure the existence of the tractions . Hence, the terms involving these tractions in (2.2) must be eliminated by means of the boundary conditions and auxiliary conditions on . The traction term on the Neumann boundary can be eliminated by means of (5e). To remove the traction term on , we stipulate that  vanishes on . The traction average in the final term of (2.2) is eliminated by requiring that  be continuous. The traction jump in the final term is deleted by means of the traction-continuity condition (5c).

Summarizing, we find that a solution  of (5) satisfies

 aΓ(u,v)=l(v) (9)

for all sufficiently smooth functions that vanish on , where

 aΓ(u,v) =∫Ω∖Γσ(u):∇v (10a) l(v) =∫Ωv⋅f (10b)

Note that is assumed to be smooth on and, in particular, that it is continuous across the fault .

To furnish a functional setting for the weak formulation of (5) based on (9), we denote by the usual Sobolev space of square-integrable functions from  into with square-integrable distributional derivatives of order , equipped with the inner product

 (u,v)k,Ω=∑|α|≤k∫ΩDαu⋅Dαv

and the corresponding norm and semi-norm . For the square-integrable functions and the corresponding norm and inner-product, we introduce the condensed notation , and . We denote by the subspace of of functions that vanish on . To accommodate the discontinuity corresponding to the dislocation, we introduce the lift operator , which assigns to any suitable slip a function in such that . A precise specification of conditions on the slip distribution is given in Section 4. The weak formulation of (5) based on (9) writes

Weak formulation: given the lift , find such that

 aΓ(u,v)=l(v)∀v∈H10,D(Ω). (11)

The bilinear form and linear form in (11), are the extensions (by continuity) of the corresponding forms in (10). The treatment of the dislocation by means of a lift operator in (11) is analogous to the treatment of inhomogeneous Dirichlet boundary conditions in weak formulations; see, e.g., [7, p.113 ]. Let us note that in the weak formulation (11), we have identified with . This identification is unambiguous, on account of a one-to-one correspondence between the functions in these spaces.

To analyze the weak formulation (11), and to prepare the presentation of the Weakly-enforced-Slip method in Section 3, we note that (11) is to be interpreted in the following manner: find with such that

 a(¯u,v)=l(v)−aΓ(ℓb,v)∀v∈H10,D(Ω). (12)

where the bilinear form ,

 a(u,v)=∫Ωσ(u):∇v, (13)

corresponds to the restriction of to . Indeed, for all pairs , the function is Lebesgue integrable on , and because the manifold corresponds to a set of -Lebesgue measure zero, the integrals of on and on  coincide. It is important to note that the restriction of the bilinear form  to in (13) is independent of the fault . The function in (11) is referred to as the continuous complement of the solution  with respect to the jump lift  and, indeed, it resides in .

For the assumed linear-elastic behavior according to (1) and (2), it follows straightforwardly that

 |aΓ(u,v)| ≤¯¯cA|u|1,Ω∖Γ|v|1,Ω∖Γ (14a) |aΓ(u,u)| ≥c–A|u|21,Ω∖Γ (14b)

for all , where and denote the continuity and strong-positivity constants of the elasticity tensor, respectively, and represents the usual -seminorm. By virtue of Poincaré’s inequality (see, e.g., [2, Theorem 5.3.5]) there exists a bounded positive constant such that

 ∥u∥1,Ω∖Γ≤CP|u|1,Ω∖Γ∀u∈H10,D(Ω∖Γ) (15)

Hence, the -norm and -semi norm are equivalent on . Equations (14) and (15) imply that the bilinear forms  and, accordingly, , are continuous and coercive. Moreover, it is easily verified that the linear form in the right member of (12) is continuous. Problem (12) therefore complies with the conditions of the Lax-Milgram lemma (see, for instance, [2, Theorem 2.7.7]) and, hence, it is well posed.

It is interesting to note that (12) allows for a physical interpretation that is very close to Volterra’s classical construction for dislocations, popularly known as the ‘Volterra knife’: to make a cut in the material, displace the two sides and hold them while welding the seam, and finally release the sides so the material assumes its state of self-stressed equilibrium. The initial cut and displacement is represented by the lift , which is not in equilibrium and is hence maintained by an external load. The addition of represents the transition to a state of equilibrium, by removing the external load but leaving the displacement intact.

### 2.3 Finite element approximation

Galerkin finite-element approximation methods for Volterra’s dislocation problem (11) are generally based on a restriction of the weak formulation to a suitable finite dimensional subspace. The general structure of finite-element methods can be found in many textbooks, for instance, [11, 22, 7, 4]. We present here the main concepts and definitions for the ensuing exposition.

The approximation spaces in finite-element methods are generally subordinate to a mesh , viz., a cover of the domain by non-overlapping element domains . The subscript indicates the dependence of the mesh on a resolution parameter, for instance, the diameter of the largest element in the mesh. In general, we impose some auxiliary conditions on the mesh, such as shape-regularity of the elements and conditions on the connectivity between elements; see, for instance, [7, 5] for further details. A finite-element approximation space subordinate to  can then be defined, for instance, as the subspace of vector-valued continuous element-wise polynomials of degree  which vanish on :

 Vph={vh∈C0(¯¯¯¯Ω,RN):(vh)i|κ∈Pp for all κ∈Th,i=1,2,…,N,vh|D=0} (16)

with the -variate polynomials of degree . Below, our interest is generally restricted to the -dependence of the approximation space and, accordingly, we will suppress . The finite-element approximation of (11) based on an approximation space  writes: find with subject to

 a(¯uh,vh)=l(vh)−aΓ(ℓb,vh)∀vh∈Vh. (17)

We refer the right-most term in the right member of (17) as the lift term.

Approximation properties of the Finite Element Method are generally investigated on the basis of a sequence of meshes , parametrized by a decreasing sequence of mesh parameters with as only accumulation point. A sequence of meshes is called quasi uniform if there exist positive constants and , independent of , such that for all and all . Standard interpolation theory in Sobolev spaces (see, for instance, [7, 2]) conveys that a sequence of approximation spaces of the form (16) based on quasi-uniform meshes possesses the following approximation property: there exists a positive constant §§§We use to denote a generic positive constant, of which the value and connotation may change from one instance to the next, even within a single chain of expressions. independent of  such that for all , all and both , it holds that

 infvh∈Vh∥v−vh∥m,Ω≤Chl+1−m|v|l+1,Ω∀v∈Hk+2(Ω)∩H10,D(Ω) (18)

with . The estimate (18) imparts that for all sufficiently smooth , the -norm of the best approximation in  in that norm decays as as .

The lift in (17) is in principle arbitrary. However, the use of an arbitrary lift carries severe algorithmic complexity, as one has to explicitly construct the lift and evaluate integrals involving products of (gradients of) the lift and finite-element shape functions. Moreover, the evaluation of these integrals by a suitable numerical integration scheme generally leads to a high computational complexity, because the fault is allowed to intersect elements, and there are no efficient quadrature schemes to integrate the discontinuous function that arises. Therefore, in practice, it is convenient to integrate the lift in the finite-element setting. Provided that that the fault coincides with element edges, a lift is then constructed in the broken approximation space:

 ^Vh={u∈C0(¯¯¯¯Ω∖Γ):u|Ω±∈(Vh)|Ω±} (19)

It is to be noted that the slip does not generally reside in and, accordingly, is to be replaced by a suitable interpolant . Moreover, if the fault does not coincide with element edges, then it is to be replaced by an approximation subject to this condition. The aforementioned approach corresponds to the split-node method by Melosh [14], where the adjective ‘split’ refers to the discontinuities between the elements in and contiguous to . The split-node approach is analogous to the standard treatment of Dirichlet boundary conditions; see, for instance, [7, §3.2.2]. The split-node approach bypasses the aforementioned complications of an arbitrary-lift approach and the evaluation of the lift term comes essentially free of charge as part of the regular stiffness matrix. A fundamental disadvantage of the split-node approach, however, is that it requires that the fault coincides with element edges, which connects the fault geometry to the geometry and, generally, the topology of the mesh. As a result, computational primitives such as the stiffness matrix and preconditioners for the stiffness matrix, which are contingent on the mesh, cannot be reused for analyses of alternative fault geometries. This is a prohibitive restriction if many dislocation geometries have to be considered, for instance, in inverse problems.

## 3 The Weakly-enforced Slip Method

In section 2.3 we substantiated that the treatment of the lift term in the split-node approach, which is the natural counterpart of the standard treatment of Dirichlet boundary conditions in finite-element approximations, is unsuitable if many fault geometries have to be analysed, on account of the inherent dependence of the finite-element mesh on the fault geometry.

In this section we propose a new and fundamentally different treatment of the lift term that retains mesh independence: the Weakly-enforced Slip Method (WSM). Below, we present two different formal derivations of the WSM formulation. The derivation in Section 3.1 relies on a limit procedure. In Section 3.2, we derive the WSM formulation on the basis of Nitsche’s variational principle for enforcing Dirichlet-type boundary conditions [15].

### 3.1 Collapsing the lift

Our aim is to derive a tractable finite-element approximation of Volterra’s dislocation problem (11), in which the finite-element space and the bilinear form and, accordingly, the stiffness matrix are independent of the fault geometry.

In principle, the lift-based Galerkin formulation (17) already exhibits the appropriate form. However, as elaborated in Section 2.3, the corresponding finite-element formulation is intractable for general lift operators. One can infer, however, that the complications engendered by a general lift operator can be avoided by collapsing the support of the lift on the fault. The integration of a discontinuous function in  then reduces to the integration of a smooth function on the dislocation. Numerical evaluation of integrals on the dislocation is feasible given a parametrization of the fault. Moreover, the intricate explicit construction of a lift is obviated, and only the slip distribution itself is required, which is presented as part of the problem specification. A further advantage of collapsing the lift is that of localization: instead of having to evaluate the lift term for all shape functions of which the support intersects with the support of the lift, only the shape functions of which the support intersects with the dislocation have to be considered.

To derive the lift term corresponding to a collapsed lift, we consider a symmetric lift as illustrated in Fig. 2, with compact support in an -neighborhood of the fault:

 Γε:={x∈Ω:dist(x,Γ)<ε}. (20)

By virtue of the compact support of in , the lift term of (17) evaluates to

 aΓ(ℓεb,v)=∫Γε∖Γσ(v):∇ℓεb=∫Γb⋅{σν(v)}−∫Γε∖Γℓεb⋅divσ(v). (21)

The first identity follows from the symmetry of the bilinear form in (10a). The second identity results from integration-by-parts and and . Let us note that the second identity is formal in the sense that it requires more regularity of  than is actually provided by . We shall momentarily ignore this aspect, but it manifests itself in the analysis of the approximation properties of the WSM formulation in Section 4. Without loss of generality, we can assume to be bounded independent of . The second term in the final expression in (21) then vanishes if collapses on . Therefore, formally passing to the limit in (21), we obtain

 aΓ(ℓεb,v)ε→+0−−−→∫Γb⋅{σν(v)} (22)

According to (22), the lift term reduces to an integral on in the limit of collapsing the support of the lift onto the fault. The WSM formulation corresponds to replacing the lift term in the right member of (17) by the limit functional according to (22):

Weakly-enforced Slip Method: given a slip distribution , find such that

 a(uh,vh)=l(vh)−∫Γb⋅{σν(vh)}∀vh∈Vh. (23)

The nomenclature Weakly-enforced Slip Method serves to indicate that in (23) the slip discontinuity is weakly enforced in the right-hand side only, and does not appear in the approximation space.

It is to be noted that although the WSM formulation (23) is derived from the lift-based formulation (17) by collapsing the lift, in contrast to (17) we do not add a lift to the continuous complement . Because , WSM thus yields a continuous approximation to the discontinuous solution of the Volterra dislocation problem (11). This implies that the approximation near the dislocation will inevitably be inaccurate. We will however show in Section 4 that away from the dislocation, the error in the WSM approximation converges optimally under mesh refinement.

### 3.2 Alternative derivation via Nitsche’s variational principle

To further elucidate the WSM formulation, we present in this section an alternative derivation of (23) based on Nitsche’s Variationsprinzip [15]. Nitsche presented in [15] a variational principle for weakly imposing Dirichlet boundary conditions in finite-element approximations of elliptic problems, i.e., without incorporating such essential boundary conditions in the approximation space. Nitsche’s variational principle can be extended to the Volterra dislocation problem (11) to weakly impose the slip discontinuity. To specify this extension, we consider a suitable broken space which encapsulates the broken approximation space and contains the solution  to the Volterra dislocation problem (11). We define the quadratic functional :

 (24)

for some suitable constant , generally dependent on . Let denote either the broken approximation space  or the continuous approximation space  and consider the approximation according to:

 ˇuh:=arginfvh∈ˇVhJ(u−vh) (25)

Equation (25) implies that satisfies the Kuhn-Tucker optimality condition for all , where denotes the Fréchet derivative of at . For according to (24), the optimality condition implies that satisfies:

 aΓ(ˇuh,vh)−∫Γ[[ˇuh]]⋅{σν(vh)}−∫Γ[[vh]]⋅{σν(ˇuh)}+ψ∫Γ[[ˇuh]]⋅[[vh]] (26) =l(vh)−∫Γb⋅{σν(vh)}+ψ∫Γb⋅[[v]]h∀vh∈ˇVh.

The final identity follows by invoking integration-by-parts on , a rearrangement of terms and the strong formulation of Volterra’s dislocation problem in (5).

If the broken approximation space  is inserted for , the optimality condition (26) can be reinterpreted as a symmetric-interior-penalty (SIP) discontinuous-Galerkin-type formulation; see, for instance, [5, Sec. 4.2]. In contrast to standard discontinuous Galerkin formulations, the slip terms in the right-hand side, i.e., the terms containing in the ultimate expression in (26), enforce the jump discontinuity at the fault. Convergence results for this formulation can be established in a similar manner as in [15]. For suitable stabilization parameters , the functional in (24) is equivalent to and (25) implies quasi-optimal convergence of .

If the continuous approximation space  is inserted for , the terms containing  and vanish, and we obtain the WSM formulation (23). Hence, WSM can indeed be interpreted as an extension of Nitsche’s variational principle to the Volterra dislocation problem with continuous approximation spaces. Furthermore, in view of , WSM can also be regarded as a SIP discontinuous Galerkin formulation, based on a continuous subspace. One can infer that the WSM approximation retains the quasi-optimal approximation property in . However, since the continuous approximation spaces applied in WSM are not dense in , the immediate significance of this quasi-optimality for the approximation properties of WSM is limited.

## 4 Approximation properties of WSM

An analysis of the approximation properties of the Weakly-enforced Slip Method is non-trivial, owing to the fact that in WSM one considers approximations in -conforming subspaces, while the solution itself resides in , and the embedding of in is non-dense. Essentially, we attempt to approximate a discontinuous function by a continuous one and, in doing so, we incur an error that does not vanish under mesh refinement. Standard techniques to assess global approximation properties, on all of , viz., Céa’s lemma or the Strang lemmas [7, Lems. 2.25-27], therefore provide only partial information; see Section 4.2.

To provide a foundation for analyzing the approximation provided by WSM, we we first recall some aspects of traces and tractions in Section 4.1. Section 4.2 investigates the global approximation properties of WSM, i.e., on the entire domain. Section 4.3 establishes the local approximation properties of WSM, i.e, on the domain excluding a neighborhood of the fault.

### 4.1 Traces and tractions

To enable an analysis of the approximation behavior of WSM, some elementary aspects of trace theory are required. For a comprehensive overview, we refer to [18]. To make the theory applicable to the Volterra dislocation problem, we must impose auxiliary smoothness conditions on the elasticity tensor. In particular, we assume:

 Aijkl∈C1(¯¯¯¯Ω) (27)

It is to be noted that (27) implies that the elasticity tensor is continuous on the domain, including the boundary, and across the fault, including the dislocation. The continuity on the subdomains and  ensures that tractions are well defined. The continuity across the fault is required to establish global convergence of the WSM approximation in the -norm and local convergence in the -norm; see Sections 4.2 and 4.3 below. Let us note that in tectonophysics the elasticity tensor is generally assumed to be uniformly constant in the domain.

Let denote an arbitrary connected domain with Lipschitz boundary. In particular, recalling the partition of into the complementary subsets , we envisage . We denote by the restriction of a function in to the boundary . By virtue of the density of in , the operator can be extended to a linear continuous trace operator, denoted by , from  into . The image of is denoted by . Considering a subset , we denote by the composition of the trace operator and the restriction to . The image of restricted to the class of functions that vanish on is indicated by :

 H1/20(ϰ)={γϰu:u∈H10,∂ω∖ϰ(ω)}. (28)

The space can be endowed with the norm:

 ∥λ∥1/2,ϰ:=inf{∥u∥1,ω:u∈H10,∂ω∖ϰ(ω),γϰu=λ}. (29)

with the obvious extension to . There exist continuous right inverses

 γ−1:H1/2(∂ω)→H1(ω),γ−1ϰ:H1/20(ϰ)→H10,∂ω∖ϰ(ω), (30)

of and . Such a right inverse is called a lifting (or lift) of the trace. It is to be noted that lift operators are generally non-unique.

We denote by the traction of a function in on , where denotes the exterior unit normal vector on . We define

 H1divσ(ω):={u∈H1(ω):divσ(u)∈L2(ω)}. (31)

Applying index notation for transparency, the chain rule yields

 ∂jσij(u)=(∂jAijkl)ϵkl(u)+Aijkl(∂jϵkl(u)) (32)

Therefore, the condition provides a meaningful condition on  if . For , this auxiliary condition on the elasticity tensor is satisfied under the standing assumption (27). The vector space  is a Hilbert space under the inner-product associated with the norm . The traction  can be extended to a bounded linear operator, denoted by , from into , the dual space of . For each , the functional acts on functions in by means of the following duality pairing:

 ⟨σn(u),λ⟩=∫ωdivσ(u)⋅γ−1(λ)+∫ωσ(u):∇γ−1(λ) (33)

One may note that for functions in , Equation (33) corresponds to a standard integration-by-parts identity. Continuity of the operator thus defined follows from the sequence of bounds:

 ∣∣⟨σn(u),λ⟩∣∣ (34) ≤(∥∥divσ(u)∥∥2ω+∥∥σ(u)∥∥2ω)1/2(∥∥γ−1(λ)∥∥2ω+∣∣γ−1(λ)∣∣1,ω)1/2 ≤(1+¯¯c2A)1/2(∥u∥21,ω+∥divσ(u)∥2ω)1/2∥∥γ−1(λ)∥∥1,ω

and the continuity of the lifting of the trace from into . The dual space is a Banach space under the norm

 ∥v∥−1/2,∂ω=supλ∈H1/2(∂ω)⟨v,λ⟩∥λ∥1/2,∂ω.

The restriction of the traction to a subset of the boundary can be extended to a bounded linear operator  from into . The functional acts on functions in  via the duality pairing:

 ⟨σn,ϰ(u),λ⟩=∫ωdivσ(u)⋅γ−1ϰ(λ)+∫ωσ(u):∇γ−1ϰ(λ) (35)

Continuity of the operator thus defined follows in a similar manner as in (34).

For non-rupturing faults, it holds that and the above definitions apply without revisions. The slip discontinuity (5b) and the traction discontinuity (5c) are then to be understood in the sense of traces and tractions outlined above. However, for rupturing faults, i.e., if the dislocation intersects with the boundary of the domain, then , and further consideration is required. We can accommodate in the space:

 ˜H1/2(Γ):={γΓu:u∈H1(ω)} (36)

with . The space is a Banach space under the norm

 ˜∥λ∥1/2,Γ=inf{∥u∥1,ω:u∈H1(ω),γΓu=λ} (37)

The principal complication pertaining to rupturing faults, is that the corresponding slip vectors cannot be lifted into , as traces of functions in do not admit the discontinuity that would otherwise arise at the intersection of and . Hence, we cannot use (35) to define an extension of the restriction of the traction, , to a bounded linear operator from  into . However, there exists a continuous right inverse of the operator , for instance,

 ˜γ−1Γ(λ)=arginf{|u|1,ω:u∈H1(ω),γΓu=λ} (38)

Let us note that the image of the lift operator corresponds to a harmonic function subject to inhomogeneous Dirichlet conditions on with data and a homogeneous Neumann condition on . The lift operator can be modified to include homogeneous Dirichlet conditions on in the codomain, if necessary. The lift operator enables us to extend to a continuous linear operator:

 ˜σn,Γ:{u∈H1divσ(ω):∥∥σn,∂ω∖Γ(u)∥∥−1/2,∂ω∖Γ=0}→[˜H1/2(Γ)]′ (39)

The functional acts on functions in via the duality pairing:

 ⟨˜σn,Γ(u),λ⟩=∫ωdivσ(u)⋅˜γ−1Γ(λ)+∫ωσ(u):∇˜γ−1Γ(λ) (40)

Essentially, in (39) and (40), we have defined the extension of the restriction of the traction to the dislocation, , by restricting the domain of the extended operator to functions for which the traction vanishes on . This restriction in the definition is consistent with the standing assumption that an intersection of the dislocation with the boundary of the domain can only occur at Neumann boundaries.

In the analysis below, we restrict ourselves to non-rupturing faults. The analysis in Section 4.2 however extends to rupturing faults by replacing the spaces and trace and traction operators for non-rupturing faults with those for rupturing faults.

### 4.2 Global approximation properties of WSM

To assess the global approximation properties of WSM, we first construct an upper bound on the functional in the right member of the WSM formulation. It is to be noted that, in general, . Hence, the functional does not admit an interpretation as a duality pairing according to (35). Because is piecewise polynomial, however, an upper bound can be constructed on the basis of inverse and trace inequalities. We refer to [5] for a comprehensive treatment of this subject. Inverse and trace inequalities can generally be derived under suitable (sufficient) regularity conditions on the finite-element mesh; see [5, Chap. 1]. A detailed treatment of the conditions underlying inverse and trace inequalities is beyond the scope of this work. Instead, we shall directly assume that a suitable discrete trace inequality holds. To formulate the assumption, for a given sequence of partitions and for each , we denote by a dense cover of the fault by means of open intersections of the fault with element interiors or with element faces, i.e.,

 Sh ={s⊂Γ:s=Γ∩κ≠∅ for some % κ∈Th} (41) =∪{s⊂Γ:s=int(Γ∩∂κ0∩∂κ1)≠∅ for some κ0,κ1∈Th,κ0≠κ1}

See Figure 3 for an illustration. The separate treatment of element boundaries serves to ensure that subsets of that coincide with element faces are separately included in .

To each segment , we associate a pair of contiguous elements such that and . If (resp. ) for some then and will be distinct (resp. identical). We assume that the following discrete trace inequality holds for all , all and all element-wise polynomial functions  of degree at most :

 (diam(κ))1/2∥vh∥s≤CΓ∥vh∥κ, (42)

for both , for some independent of ; cf. [5, Lemma 1.46]. The constant is allowed to increase with the polynomial order .

###### Lemma 1 (Continuity of the WSM linear form).

Consider a manifold , a slip distribution and a sequence of partitions such that for all , the discrete trace inequality (42) holds for all and all element-wise polynomial functions on  and

 ∥b∥Th,Γ:=(∑s∈Shh−1s∥b∥2s)1/2<∞ (43)

with the harmonic average of the diameters of the elements adjacent to ,

 1hs=1diam(κ+s)+1diam(κ−s) (44)

Then for all , the linear form is continuous and

 ∣∣(b,{σν(v)})Γ∣∣≤2−1/2¯¯cACΓM1/2∥b∥Th,Γ∥v∥1,Ω (45)

with the continuity constant of the elasticity tensor in (3), the constant in the discrete trace inequality (42) and  the maximum multiplicity of the multiset .

Remark. The maximum multiplicity  in Lemma 1 indicates the maximum number of occurrences of any one element in connection to any segment  as a member of the set . For instance, in Figure 3, the element has multiplicity . One can infer that is bounded by the maximum number of faces of any element in the mesh, increased by for interior segments.

###### Proof.

We first separate the integral on  into a sum of contributions from the segments and apply (3) and the Cauchy-Schwarz inequality to obtain

 ∣∣(b,{σν(vh)})Γ∣∣ =∣