Non-planar fault models: complexities induced by crustal layering in transcurrent faulting processes

Summary

1 Introduction

Surface maps of faults usually show complex geometrical patterns: fault bending and branching along the strike direction is the rule rather than the exception, suggesting a fractal distribution of faults (e.g. Shaw 1981; Aviles 1987; Okubo & Aki 1987). The irregular fault geometry is considered to be responsible for the generation and preservation of stress heterogeneities (e.g. Andrews 1989). These features are generally interpreted as the result of processes of near-surface segmentation. Many field observations support this interpretation: recent geological and geophysical studies of several stepover zones along the San Andreas Fault system in California reveal that the complex nature of the fault zone at the surface masks a much simpler geometry at depths (>5 km) associated with large earthquakes, as determined by accurately relocated hypocentres (Graymer 2007). In particular, gravity, magnetic and seismicity studies suggest that the Hayward and Calaveras faults in Northern California merge at 3-5 km depth (Ponce 2004); in the Parkfield region (Central California), the complexity shown by surface traces and the discordance of the location and geometry of the mainfault trace relative to the seismicity appears limited to the upper few kilometres (Simpson 2006). It is not a coincidence, probably, that in the fault region a sharp transition of seismic wave velocity is observed at comparable depth (Lin 2010).

Simple at depth extensional stepovers have been also proposed for the North Anatolian Fault section activated during the 1999 Ismit (Turkey) earthquake (Aochi & Madariaga 2003), for the Nojima-Rokko Fault system of the 1995 Kobe (Japan) earthquake (Wald 1996) and for the complex fault system associated with the 1992 Landers, California, earthquake (Felzer & Beroza 1999). In Fig. 1 the system of fractures is shown, which was left in the South Iceland Seismic Zone (SISZ) after the M6.6 earthquake of 2000 June 17. Lines of fractures parallel to the strike direction, but distant hundreds of metres from the fault trace, have been identified and clearly cannot be justified only in terms of the stress released by the vertically dipping main fault, well costrained by the distribution of aftershocks and by the fault plane solution (Stefánsson 2011). The post-seismic microearthquake activity is sharply confined below 3 km depth, so that the bending/branching point is inferred to be shallower and probably connected with major discontinuities in the velocity profile (Vogfjörd 2002; Antonioli 2006).

Figure 1

Aftershocks (dots) and pattern of surface fractures (double lines) in the SISZ earthquake of 2000 June 17 (Hjaltadóttir & Vogfjörd 2005).

Open in new tabDownload slide

When surface fractures appear laterally on one side of the deep fault strike, they can be ascribed to a bending process (Fig. 2a) while a branching process can be invoked when two or more parallel fractures are found on both sides (Fig. 2b).

Figure 2

Schematic representation of a bending process (a) and a branching process (b). The fault sections are shown in grey, the slip is shown by the large arrows. Double lines denote surface fractures, the dashed line is the surface strike of the deep section.

Open in new tabDownload slide

Non-planar fault geometries have been investigated by many authors since geometrical complexities, such as bending (change of strike and/or dip), branching and stepovers, may play an important role in nucleation, propagation and termination processes of an earthquake rupture. The change in strike (bending process) of a strike-slip fault was studied by Duan & Oglesby (2005) to examine the long-term effects of a non-planar fault geometry. Kame (2003) considered the propagation of a mode-II rupture on a pre-existing branched fault to study the effects of pre-stress state and rupture velocity on dynamic branching. Also many experimental studies were performed to study the influence of bending and branching on shear ruptures (e.g. Rousseau & Rosakis 2003), with the purpose of simulating earthquake rupture processes. An important question is why and how such complexities originate since quasi-static crack solutions in a homogeneous medium show that an isolated model fault never bends if the fault extends in the direction of maximum shear traction beyond the tip. Kame & Yamashita (1999) showed that fault bending may take place even in a homogeneous medium if the rupture velocity approaches the shear wave speed, due to dynamic stresses deflecting the direction of maximum shear traction. Employing the same technique, Ando (2004) describe the interaction between fault segments, resulting in cohalescence or repulsion depending on their offset distance. Even if dynamic effects may produce fault bending in a homogeneous medium, non-planar fault geometries must be necessarily taken into account when the faulted medium is heterogeneous, as pointed out by Rybicki & Yamashita (2008), because the fault may bend abruptly at the interface between different media, where the orientation of principal stresses changes discontinuously. In particular, an explanation in terms of elastic heterogeneities allows preserving the fault bending geometry during several earthquake cycles.

The effects of crustal heterogeneities on transcurrent faulting processes have been investigated by many authors. A review of dislocation models in layered media is provided in a beautiful book by Segall (2010). A numerical code for solving point-dislocation problems in elastic and viscoelastic media is provided by Wang (2006). Rybicki (1971) studied the deformation field due to a strike-slip fault in a layered medium, showing that aftershock activity is affected by the residual stress field left after a major earthquake. Rybicki & Yamashita (1998) have shown that low-rigidity media may represent a barrier for fracture propagation in the case of faulting in vertically inhomogeneous media. Analogous results were obtained by Kame (2008) studying the quasi-static growth of strike-slip faults in presence of layered media and Hirano & Yamashita (2011) providing a convenient mathematical framework to analyse faults located along or across a material interface. The presence of a shallow soft layer was invoked by Belardinelli (2000) to interpret the complex pattern of surface fractures observed in the SISZ.

Dislocation models assume that the displacement discontinuity is assigned (and typically is piecewise constant) over the fault surface. On the contrary, crack models assume that the stress drop is assigned as a bounded function prescribed over the fault surface. The stress drop must be bounded because unbounded tractions cannot be sustained over a faulted surface. Due to the singular nature of the integral equation describing crack equilibrium, bounded tractions over the failure surface are generally accompanied by stress singularities beyond crack tips, even in a homogeneous medium. If a crack meets the welded interface between two different elastic media, a further singularity appears at the intersection (e.g. Erdogan 1973).

Bonafede (2002) have employed crack models to study theoretically the deformation and the stress fields induced by a planar strike-slip fault cutting the interface between different layers. A stress drop discontinuity condition was derived as a necessary condition that has considerable influence on the style of faulting in layered media. In fact many relevant geophysical cases are incompatible with the development of a planar fault and therefore complexities including fault bending, fault branching and the development of en echelon system of fractures must be involved in transcurrent fault processes within layered media. The aim of this work is to study fault bending and fault branching using the mathematical framework of crack models in 2-D antiplane configuration. This approach allows finding admissible solutions taking into account that singular stress values are unacceptable over the crack surface even if the dislocation density is found to be generally singular at the interface.

2 Preliminaries

Dislocation models are generally employed to model the deformation and stress fields after faulting events. The simplest dislocation model (hereafter termed ‘elementary dislocation’ accounts for a constant discontinuity Δu of the displacement field (Burgers vector), prescribed over a half-plane terminating at the dislocation line (e.g. Landau & Lifshitz 1967). It is well known that elementary dislocations face severe physical problems: the stress field in proximity of the dislocation line has a non-integrable singularity (of order -1) both within the faulted half-plane and beyond, the energy required to generate an elementary dislocation is infinite, continuity and impenetrability of matter are violated. If Δu is not constant and vanishes approaching the dislocation line, the order of singularity of the stress field changes.

Physical considerations lead us to require that when a fracture takes place in an unfractured medium the released tractions should be bounded over the dislocation surface, while square-integrable singularities may be acceptable outside. These features are provided by crack models, in which the traction Δσ released over the dislocation surface is assigned and the displacement discontinuity Δu is found as a function of the position over the faulted surface, after solving a singular integral equation. In a homogeneous elastic unbounded medium with rigidity μ in an antiplane configuration [in which the only non-vanishing component of displacement is u_y(x, z)], the integral equation is simply

(1)

where the dislocation surface is the strip -ℓ < z < ℓ in the plane zy, and the integral must be evaluated in the ‘principal value’ sense. It is convenient to normalize the crack surface to the interval (-1, 1) and so we introduce the dimensionless variable A6) by . In order that finite values of Δσ are obtained over the crack surface -1 < ξ < 1, the dislocation density distribution ρ(ξ′) must possess a singularity of order close to crack tips (e.g. Muskhelishvili 1953):

where R(ξ′) is a regular function with R(-1) ≠ 0 and R(1) ≠ 0. An efficient method of solution of (1) consists in expanding ρ(ξ′) in Chebyshev polynomials of the first kind T_n(ξ′) = cosθ′nξ′ [with ξ= cosθ′, see e.g. Abramowitz & Stegun (1964)], exploiting the integral property

(2)

where (with ξ= cosθ) are Chebyshev polynomials of the second kind.

Once the regular factor R(ξ′) of the dislocation density distribution is expanded in T_n polynomials

eq. (1) reduces to

and ±_n, n= 1, 2, … may be obtained, if Δσ is assigned, exploiting the orthogonality relation of U_n polynomials

This efficient method of solution was presented in Bilby & Eshelby (1968) (see also Bonafede 1985).

From eq. (2) it appears that the coefficient α₀ is left undetermined by the integral equation, since the term n= 0 does not provide any stress release over the crack surface (the term n= 0 is solution of the homogeneous version of the integral eq. 1). In simple problems, α₀ must vanish if we assume that the crack slip vanishes at both ends. Some dislocation problems however require that a crack may be open at one end: an example is provided by a crack opening at the free surface of a half-space (see case B1 in Bonafede & Rivalta 1999b) and in such a case a supplementary condition must be provided to fix α₀, which may be obtained from the requirement that tractions on the free surface vanish at the intersection with the crack tip. More generally, one might conceive the problem in which one crack is split into two interacting sections: then both sections must be open at the merging point and each of the two dislocation distributions ρ₁ and ρ₂ will contain an undetermined coefficient for n= 0. Two supplementary conditions are then needed: one is simply provided by the request that slip is continuous passing from one section to the other, while the second can be obtained imposing that the singular traction produced by section 1 onto its own crack surface cancels with the contribution produced by section 2 and vice versa. Even if this case may seem only of academic interest, its solution is required to address the more complex problems in which two crack sections are embedded in different media with rigidity µ₁ and µ₂ (see e.g. Bonafede 2002) and/or are not coplanar, which is the target of this paper.

Of course, when a layered medium is considered (µ₁≠ µ₂), the integral eq. (1) must be generalized to account for the welded interface between the two different media. In the case of a vertically dipping crack with the upper tip (ξ=-1) touching the interface between two half-spaces, the integral equation takes the following form (see case B in Bonafede 2002, with ):

(3)

With respect to the crack problem in a homogeneous medium, a second integral kernel appears that is unbounded only if ξ, ξ′ go to -1 simultaneously. This is a generalized Cauchy kernel which changes the singular nature of the dislocation density. In this case the dislocation density distribution may be separated in a singular term and a regular term R(ξ′)

where the order of singularity b close to the tip ξ′=-1 may be determined employing asymptotic methods (see e.g. Erdogan 1973). Examples of the method may be found in Morelli (1987) and Bonafede & Rivalta (1999b). For the sake of completeness, the integral kernels for tensile and for edge dislocations across two welded half-spaces are provided by Bonafede & Rivalta (1999a) and by Rivalta (2002), respectively.

3 Fault Bending Model

We consider an elastic half-space in z < 0 with rigidity µ₁ welded to an elastic half-space in z > 0 with rigidity µ₂. The origin of the reference system is placed at the intersection between the crack and the interface, with the z-axis pointing downwards. We consider an antiplane strain configuration in which the only non-vanishing component of the displacement field is u_y(x, z), which is independent of the y-coordinate.

It is convenient to split the crack into two interacting sections, each embedded in one medium and both open at the interface. The section in the upper layer has length 2l₁, it is inclined by an angle with respect to the negative z-axis and its unit normal is The section embedded in the lower half-space is normal to the interface and its length is 2l₂. A sketch of the model is presented in Fig. 3.

Figure 3

Bending model: the fault sections are shown as thick segments and the slip is normal to the plane of the figure.

Open in new tabDownload slide

Consider the equilibrium equations for the two crack sections

where Δ₁(x, z) is the stress drop in a point of the upper crack section, (r, α) being the polar coordinates of the point (x, z), and Δ₂(0, z) is the stress drop in a point of the lower (vertical) section. The stress drop on each crack section is equal to the sum of two contributions:

• (i)

  are the tractions induced on section 1 due to crack slip on section 1 (I) and on section 2 (II), respectively;

• (ii)

  are the stresses induced on section 2 due to crack slip on section 2 (II) and on section 1 (I), respectively.

On section 2 the stress component released by crack slip is simply σ_xy, while on section 1 the released component is σ_ny, which represents the y component of the traction relative to a surface having as normal versor; we have

(4)

To set up a crack model we employ a continuous distribution ρ₁ of elementary dislocations with dislocation lines parallel to the y-axis, contained in the interval 0 < r₀ < 2l₁ and a distribution ρ₂ with dislocation lines in the interval 0 < z₀ < 2l₂. In this case ρ is split into two subdomains in this way

where ρ₁, ρ₂ are defined through the displacement discontinuity over the crack plane:

(6)

(7)

From (6) and (7) we obtain

where C₁, C₂ are constants that can be fixed imposing the conditions of crack closure

(8)

(9)

The equilibrium equations can be rewritten in the form

(10)

where the elementary solutions and are singular when r=r₀ and z=z₀, is singular only when r=z₀= 0 and is singular only in z=r₀= 0, so that integrals must be evaluated in the principle-value sense.

The integral kernels , for an elementary screw dislocation of arbitrary dip embedded in the upper layer are given by Singh & Rani (1994), while the kernels can be obtained using the analytical solutions proper to a vertical screw dislocation (Bonafede 2002). Once the stress drops Δσ₁, Δσ₂ are assigned, the equilibrium equations become a system of coupled integral equations for the unknown functions ρ₁, ρ₂.

Once the dislocation density ρ is known, the solution of a crack problem is simply obtained by a superposition of elementary solutions; in our case, if f denotes any elementary component of displacement or stress, the crack solution f^c is given by

(11)

if we introduce the following dimensionless variables

(12)

the system of integral equations can be rewritten as follows:

(13)

(14)

where

(15)

The singular integral kernels

(16)

(17)

(18)

(19)

are generalized Cauchy kernels. The presence of the free surface introduces further terms in the integral equations which are regular, so that the singular nature of the dislocation density is not affected by their presence.

3.1 Supplementary conditions

Since both crack sections are open at the interface, two supplementary conditions must be supplied in order that the problem may yield a unique solution. The first supplementary condition is easily derived from the continuity condition of crack slip at the welded interface between the media.

(20)

A second condition is needed to provide finite stress values over the crack surface near the interface (z= 0) and so an asymptotic study was performed (Appendix A) in order that cancellation of singular terms can be guaranteed on each crack section.

The dislocation density distribution has the following form (Muskhelishvili 1953; Erdogan 1973):

(21)

(22)

where and are regular functions [with R₁(± 1), R₂(± 1) ≠ 0]. With respect to the case of a planar fault, a singularity of order ω now appears at the interface and its value is dependent from model parameters α, m, as shown in Fig. 4, while the relative lengths of the two crack sections have no influence on its value [the ratio disappears in the derivation of eq. (A16)]. Once the inclination angle α is fixed, Fig. 4(a) shows that an increase of the rigidity contrast determines a significant increase of the singular degree ω only for m≲ 5. Instead, however we fix the rigidity contrast, the increase of α always determines a significant increase of ω, as shown in Fig. 4(b). We can also observe that the singularity degree ω is always greater than (see Fig. 4b), and this consideration will be useful in the next section to deploy a method of solution.

Figure 4

Dependence of the singularity degree ω from the rigidity contrast m=µ₁/µ₂ and the dip angle δ = 90° -α. Note that the dislocation density is bounded at the interface (ω= 0) only if the two crack sections are coplanar (vertical). The dislocation density must be singular at the bending point even if medium is homogeneous (m= 1).

Open in new tabDownload slide

If the two media have the same rigidity (m= 1), it is possible to explicit the dependence of the singularity on α

and if , it is easy to verify that .

In the homogeneous case, ω takes lower absolute values with respect to the heterogeneous case (m > 1) since only two singular terms must be balanced on both crack sections. The third singular terms are only present if m≠ 1 (Γ≠ 0) and they describe how the rigidity contrast affects the stress release which each crack section produces on itself (second term on the right-hand side of eqs 13 and 14).

Once the value of ω is obtained from the asymptotic study, we find from system (A15) the following condition relating R₁(-1), R₂(-1)

(23)

if α= 0, then ω= 0 and the previous condition reproduces correctly the condition found for a vertical planar crack by Bonafede (2002) (apart from the different sign due to the different convention adopted for the orientation of the upper crack section). We can observe that, when bending takes place, the relation (23) tells us that the ratio between R₁(-1) and R₂(-1) is not only dependent on the rigidity contrast m, but also on the other model parameters α, l₁, l₂.

3.2 Method of solution

The singular nature of the dislocation density suggests expanding the regular factors in Jacobi polynomials , which are the orthogonal polynomials for the weight function

(24)

(25)

However no convenient analytical expression like (2) is available for the singular integrals in the system of integral eqs (13)-(14) and so we are forced to search an approximate solution.

An approximate solution for the systems (13) and (14) can be obtained if we rewrite the dislocation density in the following way:

(26)

where the exponent , since we know that the singularity degree is always in the interval (as shown in Fig. 4b, Section 1.1); in this case can be expanded in a uniformly convergent series of Chebyshev polynomials even if the correct order of singularity of ρ(ξ′) near ξ=-1 cannot be reproduced when the series is truncated to a finite index. After these considerations, the dislocation distributions ρ₁, ρ₂ can be written as

(27)

(28)

where T_n are Chebyshev polynomials of the first kind and

(29)

(30)

Inserting the expressions (27) and (28) for and in the integral equations, the singular Cauchy kernels can be disposed of analytically.

(31)

(32)

In both the previous equations, the first term on the right-hand side was computed employing the integral property (2). Multiplying eq. (31) by and integrating on the interval -1, 1 in the variable ξ₁ we obtain

(33)

In the first term on the right-hand side we have used the orthogonality relation for Chebyshev polynomials of the second kind. S₁₁(k, n) and S₁₂(k, n) are double integrals with the following expression:

(34)

(35)

In the same manner, if we multiply eq. (32) by and then integrate on the interval -1, +1 in the variable ξ₂ we obtain

(36)

where the terms S₂₁(k, n) and S₂₂(k, n) have the following expressions:

(37)

and are, respectively, the coefficients of the expansion of the stress drops Δσ₁ and Δσ₂ in U_k polynomials. If Δσ₁ and Δσ₂ are constants, then only and differ from zero. The system of eqs (33) and (36) may yield a unique solution if we consider the supplementary conditions introduced in Section 3.1. If we introduce in (20) the expansions of ρ₁ and ρ₂ in Chebychev polynomials and we use the orthogonality relation between these polynomials

(38)

the following relation between the coefficients γ₀ and β₀ can be deduced:

(39)

The second supplementary condition between and coefficients can be obtained considering the ratio between ρ₁ and ρ₂ and letting t→-1

(40)

being T_n(-1) = (-1)ⁿ and , as expressed in eq. (23).

To provide an approximate solution to the crack problem, the infinite sums in (33) and (36) are truncated to a finite index N. More specifically, this means that in (33) and (36) we take

• (i)

  S_ij(k, n) = 0 for k > N- 1 ∧n > N;

• (ii)

  for k > N- 1.

Let i=k+ 1, j=n+ 1, we define the following matrix (2N) x (2N+ 2):

(41)

whose elements are the (N) x (N+ 1) matrices

In other terms, the matrix ϒ_ij, with i= 1, 2, …, 2N and j= 1, 2, …, 2N+ 2, is

Furthermore, let η be a column vector with 2N+ 2 components and

(42)

and let _χ be the 2N column vector containing the N components of Δσ₁ and the N components of Δσ₂

(43)

If Δσ₁ and Δσ₂ are constants, only the first and the (N+ 1)th components of _χ do not vanish.

The solution of the truncated problem is given by the following system of 2N equations for the 2N+ 2 unknowns η_j:

(44)

with the supplementary conditions (accounting for continuity of slip at the interface and for regularity of tractions over the crack surface)

(45)

(46)

3.3 Generalized stress drop condition

Before we can use the method of solution proposed, it is necessary to investigate the relation between the stress drops that must be assigned on the two interacting sections of the crack. As shown in Appendix B, Δσ₁ and Δσ₂ must be connected by the following relation:

(47)

where

(48)

This relation generalizes the condition corresponding to a planar crack surface normal to the interface (Bonafede 2002)

(49)

since the second term on right-hand side of (47) vanishes when α= 0.

To understand (47), let us assume that the layered medium before faulting is deformed starting from a reference configuration, where stress and strain both vanish, according to an antiplane displacement field U⁽¹⁾_y(x, z) (in the upper half-space) and U⁽²⁾_y(x, z) (in the lower half-space); the welded boundary conditions require that the traction σ_zy and the displacement U_y must be continuous across the interface z= 0; then, employing the elastic constitutive relation , we obtain

(50)

Before faulting occurs, the traction T⁽¹⁾_y over the shallow inclined fault section and the traction T⁽²⁾_y over the deep vertical fault section are, respectively,

Quantities evaluated after faulting will be denoted by primed variables. Boundary conditions similar to (50) hold in the primed configuration as well, so that the quantities and must be continuous across the interface and the residual tractions acting on the fault sections are

where τ₁ and τ₂ are residual tractions on sections 1 and 2, respectively (e.g. due to friction between fault faces). The traction released over the fault surface (i.e. the stress drop) is

Employing the continuity rules (50) we obtain

(51)

This equation, which is obtained taking into account only the welded boundary conditions, is now compared with eq. (47) obtained from crack equilibrium so that we must obtain

where is the zy-component of released stress at the bottom of the inclined fault section. The first term on the right-hand side vanishes due to the first condition (50) in the primed configuration, so that µ₁F(α, m) can be always interpreted as at the intersection between the crack and the interface. Once µ₁, µ₂, α, Δσ₁ and Δσ₂ are fixed, the stress change is determined; alternatively, if is assigned, α is determined. For instance, in the limit case , we have ; vice versa, if , then . Since has no influence on the deep fault section, a configuration can be always conceived in which ; then F(α, m) = 0 and

(52)

We shall assume in the following that no zy-component of stress is released. The resulting condition (52) shows that, once Δσ₁ and Δσ₂ are fixed by fault rheology, the bending angle α must be chosen appropriately in order that the welded boundary conditions may be fulfilled. A real solution for α can be found only if : if the equality sign holds, then α= 0 (the fault crosses the interface without bending); if Δσ₁= 0 (the initial stress in the upper layer equals the residual frictional level, as expected e.g. in an elasto-plastic medium with Mohr-Coulomb rheology), then and the upper layer unwelds from the deeper layer.

3.4 Numerical results

Following the previous analytic scheme, numerical computations have been performed assuming:

• (i)

  ν₁=ν₂= 0.25 and µ₂= 3µ₁= 30 GPa,

• (ii)

  2l₁= 2l₂= 8 km for the length of crack sections,

• (iii)

  a stepwise constant stress drop profile with Δσ₂= 3 MPa and .

The non-vanishing stress components σ_xy and σ_zy induced by shear cracks corresponding cases

1. 1

   α= 0,

2. 2

   ,

3. 3

   ,

are shown in Fig. 5. If we compare cases (2) and (3) with the case of a planar shear crack (Fig. 5a), we observe an asymmetric stress release with respect to the negative z-axis; the stress release is concentrated in the footwall block, while a minor release of σ_xy is observed in the hanging block, which is accordingly candidate to host new fractures.

Figure 5

Non-vanishing stress components σ_xy, σ_zy induced by shear crack with: (a) α= 0°, (b) α= 22.5°α and (c) α= 45°. A mild singularity is present (out of the crack surface), which can be barely appreciated only when α is large (insets in the lower panels).

Open in new tabDownload slide

From the asymptotic study we know that ρ₁ and ρ₂ are singular at ξ=-1 if α≠ 0 and this implies that also the stress components σ_xy and σ_zy are singular, but from Figs 5(a) and (b) we see that these singularities affect only a limited region around the bending point.

To check the accuracy of the approximate solution found, we show in Fig. 6 the stress components σ_xy, σ_zy computed according to the displacement discontinuity method (Crouch & Starfied 1983). This numerical technique does not involve a study of the singular nature (assuming that the order of singularity is equal to -1 at the ends of each boundary element) and provides accurate results only at distances from the fault much greater than the length of the dislocation elements. Therefore we can check the correctness of the results found employing the crack model. The comparison of stress maps shown in Figs 5 and 6 is satisfactory (away from the fault surface, of course), confirming the crack model approach. It must be stressed that the boundary element method cannot provide indications on fault bending, because singular stress values are always obtained over the fault surface, which must be assumed known a priori, while the crack model provides such indications since singular tractions are not admissible over the crack surface and the bending angle α is strictly linked to the stress drop discontinuity condition.

Figure 6

Non-vanishing stress components σ_xy, σ_zy computed according to the displacement discontinuity method: (a) α= 0°, (b) α= 22.5° and (c) α= 45°.

Open in new tabDownload slide

4 Fault Branching Model

The bending crack model is clearly symmetric if we change the sign of the bending angle α. We have already noted that more stress is released in the footwall block than in the hanging block and this leads to envisage that the hanging block will host an inclined fault section accompanying future deep events. Therefore we consider in the following a branching crack model in which the crack is split into three interacting sections, all open at the interface, but with sections 1a and 1b embedded in the upper layer and section 2 embedded in the lower half-space (Fig. 7). The sections 1a and 1b have length 2l_1a and 2l_1b, respectively; both sections are inclined with respect to the negative z-axis ( and represents the normal versor of section 1a, while is the normal versor of section 1b. As in the previous model, section 2 is normal to the interface, with length 2l₂.

Figure 7

Branching model: fault sections are shown as thick lines. Slip is normal to the plane of the figure.

Open in new tabDownload slide

Now we have to consider three equilibrium equations, one for each crack section

where Δσ_1a(x_1a, z_1a) is the stress drop in a point of section 1a, being (r_1a, α_1a) the polar coordinates of the point (x_1a, z_1a), Δσ_1b(x_1b, z_1b) is the stress drop in a point of section 1b, being (r_1b, α_1b) the polar coordinates of the point (x_1b, z_1b), and Δσ₂(0, z) is the stress drop in a point of the vertical section. The stress drop on each crack section is equal to the sum of three contributions:

• (i)

  are the tractions induced on section 1a due to crack slip on section 1a (Ia) and 1b (Ib), respectively, while is the traction induced on section 1a due to crack slip on section 2 (II);

• (ii)

  are the tractions induced on section 1b due to crack slip on section 1a (Ia) and 1b (Ib), respectively, while is the traction induced on section 1b due to crack slip on section 2 (II);

• (iii)

  is the stress induced on section 2 due to crack slip on section 2 (II), while are the stresses induced on section 2 due to crack slip on section 1a (Ia) and on section 1b (Ib), respectively.

To define a crack model we have to employ suitable superpositions of elementary solutions. We use a continuous distribution ρ_1a of elementary dislocations with dislocation lines contained in the interval , a distribution ρ_1b with dislocation lines in the interval and a distribution ρ₂ with dislocation lines in the interval 0 < z₀ < 2l₂; the dislocation density distribution ρ is so split into three subdomains

(53)

The equilibrium equations can be rewritten in the following form:

(54a)

(54b)

(54c)

where the integrals must be evaluated in the principle-value sense. Once the stress drops Δσ_1a, Δσ_1b and Δσ₂ are assigned, the equilibrium equations become a system of coupled integral equations for the unknown functions ρ_1a, ρ_1b and ρ₂. Once the dislocation density ρ is known, the solution of a crack problem is simply obtained by a superposition of elementary solutions. In our case, if f denotes any elementary component of displacement or stress, the crack solution is given by

(55)

if we introduce the dimensionless variables

(56)

and we limit to consider the case of two half-spaces welded at the interface (H→+∞), then the system of integral equation can be written in the following way:

(57)

(58)

(59)

The eqs (57) and (58) describe the stress drop on section 1a and on section 1b, respectively, and therefore these equations must have the same structure. The singular kernels and are equal to the integral kernels and in eqs. (16) and (17) if we put α=α_1a in case of section 1a and α=α_1b in case of section 1b. Instead are generalized Cauchy kernels that describe the stress drop induced on section 1a due to slip on section 1b.

(60)

(61)

and describe the stress drop induced on section 1b due to slip on section 1a and have the same expressions of the integral kernels and if we exchange the suffix 1a and 1b in the eqs (60) and (61).

In eq. (59), the singular kernel has been already introduced in eq. (18), while and are equal to the integral kernel presented in eq. (19) if we put α=α_1a and α=α_1b, respectively.

4.1 Method of solution

To solve the system of integral equations it is necessary to identify the singular nature of the dislocation density at crack tips. If we perform an asymptotic study similar to that described for the fault bending model (see  Appendix A), we find that and b_1a=b_1b=b₂=ω are necessary conditions to have finite stress release over each crack section (all the three dislocation densities must have the same singularity degree at the interface). To evaluate ω, the condition is imposed that the traction is bounded over the crack surface at the intersection: accordingly, compensation of singular terms provided by each crack section must be imposed and this condition provides three homogeneous linear equations for R_1a(-1), R_1b(-1) and R₂(-1), with coefficients a_ij depending on ω and on α_1a, α_1b, m, l_1a, l_1b and l₂.

(62)

Surprisingly enough, the determinant of this linear system vanishes identically (), so that no condition is extracted for ω. Therefore, any value of ω within the allowable interval) - 1, 0] can be considered as acceptable and only the ratios R_1a(-1)/R₂(-1), R_1b(-1)/R₂(-1) can be evaluated

(63)

Since we must impose that R_1a(-1), R_1b, (-1), R₂(-1) ≠ 0, we can conclude that

• (i)

  the angles α_1a, α_1b must be different (α_1a≠ α_1b) in order that the denominators in (63) do not vanish,

• (ii)

  the values ω=ω_1a and ω=ω_1b for which a single fracture develops in the upper medium with angle α=α_1a and α=α_1b, respectively, must be excluded since R_1a(-1), R_1b(-1) must not vanish. In fact the numerators in (63) reproduce the second factor in eq. (A17), the condition that identifies the singular degree of the dislocation density at ξ=-1 in presence of a single fracture in the upper medium (fault bending model).

To obtain information about the singular nature of the dislocation density at the interface, we are forced to adopt a different approach. We choose to use the displacement discontinuity method, since this numerical technique can be applied without any previous knowledge about the singular nature at crack tips. To apply this method, we need to know which are the stress drops to assign over the crack sections and therefore we must compare, as done for the fault bending model, the finite contributes to stress drop at the branching point.

4.1.1 Stress drop conditions

If we expand ρ_1a, ρ_1b and ρ₂ in Jacobi polynomials

(64)

(65)

(66)

and we compare the finite contributes to stress drop following the procedure (with minor tricks) shown in Appendix B, it is possible to derive the following two stress drop conditions

(67)

(68)

where the function has this form

(69)

The function G(α_1a, α_1b) has the following representation:

(70)

with

(71)

if we assume -1 ≠ωω < 0, while an alternative representation must be adopted if ω= 0,

(72)

with

(73)

since the relation (70) is no more acceptable being undefined.

The relations (67) and (68) are analogues to the condition (47) derived in the case of the fault bending model: considerations similar to those made in Section 3.3 allow interpreting as the zy-component of released stress at the bottom of the two inclined fault sections. If we assume that only σ_xy is released then the stress conditions derived reduce to

(74)

It must be observed that if α_1b=-α_1a and l_1a=l_1b the symmetry of the configuration implies that F(α_1a, α_1b, m) = 0, that instead we prescribe by assumption in case of non-symmetric configurations. These stress drop conditions can be fulfilled if the stress drop in the upper medium is lower than required for a planar through-going surface. Once the stress drop values are assigned, a vertical dipping strike-slip fault at depth may cross the interface and fault branching can take place within the sedimentary layer, provided that the shallower sections are suitably inclined as required by (74).

Once the stress drop conditions have been derived, then we can apply the displacement discontinuity method. If we consider Fig. 8(a), no singularity seems to be present at the crack tip ξ=-1 (the origin of the reference system) and if we perform a more detailed calculation with higher resolution in proximity of the interface (Fig. 8b) this conclusion is confirmed. Thus, there is no need to assume a singularity for ρ at the branching point and in the following we will consider ω= 0.

Figure 8

Stress component σ_xy computed according to the displacement discontinuity method (with α_1a= 30°, α_1b=-45°): (a) in a region of 20 km x 20 km (with N= 20 boundary elements for each crack section), (b) in a region of 200 m x 200 m (in this case a non-uniform distribution of boundary elements was used, with higher density close the origin).

Open in new tabDownload slide

In this case the form of the singular factors suggests expanding R(ξ′) in Jacobi polynomials , which are orthogonal if the weight function w(ξ′) = (1 -ξ′)^-1/2 is employed and the relations (64)-(66) reduce to

(75a)

(75b)

(75c)

Each crack section is open at the interface, therefore we must supply three supplementary conditions in order that the problem may yield a unique solution.

4.2 Supplementary conditions

The first supplementary condition is easily derived from the continuity condition of crack slip at the interface between the media.

(76)

Let us consider the dislocation density ρ₂; Jacobi polynomials in (75c) can be expressed in terms of Legendre polynomials using interrelations among orthogonal polynomials

(77)

where are ultraspherical polynomials that, if , reproduce Legendre polynomials P_k(t). The normalized dislocation density can then be rewritten in the following form:

(78)

where coefficients differ from β_k in eq. (75c) by a factor of . The presence of the singular factor t₂^-1 in eq. (78) should be of no concern, since ρ₂ always appears through integrals when observable quantities are considered and . For instance, if we compute the displacement discontinuity Δu₂, we obtain

(79)

where, in the last equality, the closure condition Δu₂(ξ₂= 1) = 0 has been employed to fix the integration constant C. The slip amplitude at the interface ξ₂=-1 can be easily computed from eq. (79), using the property that P_2k(t) polynomials are even functions and are orthogonal over (-1, 1) to P₀= 1 if k≠ 0.

(80)

Employing the same method, the dislocation densities can be written through expansions similar to (78).

(81)

(82)

and the crack slip at the interface for section 1a and 1b are, respectively,

(83)

(84)

So the first supplementary condition assumes this form:

(85)

The other two supplementary conditions are needed to provide finite stress values over the crack plane near the interface (z= 0) and if we put ω= 0 in eq. (63) we have

(86)

From eq. (86) we see that both sections in the upper layer must be inclined in order that R_1a(-1), R_1b(-1) do not vanish. This result is consistent with the stress drop conditions derived in Section 4.1.1: if the condition (eq. 49) for the development of a vertical planar crack cannot be satisfied, then the crack can cross the interface if the couple of sections in the upper medium are suitably inclined.

Since P_2k(1) = 1 and considering eqs (78), (81) and (82) we obtain a relation between the coefficients and between

(87)

(88)

4.3 Approximate solution: truncated problem

The following infinite systems

(89)

(90)

(91)

can be obtained following these steps

• (i)

  substitution of (78), (81) and (82) in the integral equations and

• (ii)
  use of the integral relation

  

  (92)

  valid for Jacobi polynomials to perform analytically the first integration, where Q_k(t) are Legendre function of the second type (e.g. Gradshteyn & Ryzhik 1965);
• (iii)

  each equation must be multiplied by and then integrated on the interval [-1, 1] in the variable ξ to obtain S_ij(k, n) coefficients.

and are the coefficients of the expansion in U_k polynomials of the stress drops Δσ_1a, Δσ_1b and Δσ₂, respectively. If Δσ_1a, Δσ_1b and Δσ₂ are constants then and are the only coefficients different from zero.

To provide an approximate solution to the crack problem, the infinite sums in (87)-(91) are truncated to a finite index N. More specifically, this means that in (89)- (91) we take

• (i)

  S_ij(k, n) = 0 for k > N- 1 ∧n > N and

• (ii)

  for k > N- 1.

Let i=k+ 1, j=n+ 1, we define the following matrix (3N) x (3N+ 3):

(93)

whose elements are the (N) x (N+ 1) matrices

In other terms, the matrix ϒ_ij, with i= 1, 2, …, 3N and j= 1, 2, …, 3N+ 3, is

Furthermore, let η be a column vector with 3N+ 3 components

(94)

and let χ be the 3N column vector containing the N components of Δσ_1a, the N components of Δσ_1b and the N components of Δσ₂

(95)

If Δσ_1a, Δσ_1b and Δσ₂ are constants, only the first, the (N+ 1)th and the (2N+ 1)th components of χ do not vanish.

The solution of the truncated problem is given by the following system of 3N equations for the 3N+ 3 unknowns η_j

(96)

with the supplementary conditions derived from (85),(87) and (88)

(97)

(98)

(99)

4.4 Numerical results

The computations have been performed assuming:

• (i)

  ν₁=ν₂= 0.25 and µ₂= 3µ₁= 30 GPa;

• (ii)

  2l₁= 2l₂= 8 km for the length of crack sections;

• (iii)

  a stepwise constant stress drop profile with Δσ₂= 3 MPa and

  • 1

    ,

  • 2

    

In Fig. 9, we show the stress drop and the crack slip corresponding to the following three cases:

Figure 9

Stress drop and crack slip on crack sections in three different geometrical configurations: (a,b) α_1a= 30°, α_1b=-30°, (c,d) α_1a= 30°, α_1b=-15°, (e,f) α_1a= 30°, α_1b= 15°.

Open in new tabDownload slide

(I) ;

(II) ;

(III) .

The approximate solutions reproduce with great accuracy the stress drop profile if the infinite sums in (78), (81) and (82) are truncated to N= 10. In case I (Fig. 9a) the stress drop over crack section 1a and over crack section 1b have the same value and furthermore (Figs 10b and c) are identical since the stress drop on section 1b have the same value in both cases.

Figure 10

Stress component σ_xy induced by a fault branching process (with α_1a= 30°, α_1b=-30°) calculated by the method of solution proposed in Section 4 (a) and, for comparison, with a boundary element approach (b) in which 20 elements for each sections have been considered.

Open in new tabDownload slide

As expected, the crack slip on the inclined sections has the same profile when the sections are symmetric with respect to the negative z-axis (Fig. 9b), while in cases II and III the crack slip is significantly affected by the relative position of sections 1a and 1b (Figs 9d and f). First, we can note that the crack slip is split between section 1a and section 1b so that greater slip is observed on the less inclined section, since the other section feels more the higher rigidity of the lower medium (Fig. 9d). If both sections are in the same quadrant (Fig. 9e), the effect is amplified so that the crack slip has very different trends on section 1a and on section 1b.

To evaluate the accuracy of the approximate solution found, we show in Fig. 10 the stress component σ_xy computed according to the displacement discontinuity method. The comparison is satisfactory since the stress map obtained with the boundary element approach is reproduced with great accuracy apart from the very proximity to the crack surface.

Observing the stress maps corresponding to cases I-III (Fig. 11), we note that when the inclined sections are in opposite quadrants, case I (Figs 11a and b) and case II (Figs 11c and d), the release of the stress component σ_xy affects a wider region, while in case III (Figs 11e and f) the pattern of stress release is similar to that obtained in the case of a single inclined section, where no significant stress release is observed in the opposite quadrant. Further we can observe that the presence of two inclined sections provides a very variable pattern for the stress component σ_zy, since this component changes sign across each crack section (Figs 11b, d and f).

Figure 11

Non-vanishing stress components σ_xy, σ_zy induced by branching shear cracks with: (a,b) α_1a= 30°, α_1b=-30°, (c,d) α_1a= 30°, α_1b=-15°, (e,f) α_1a= 30°, α_1b= 15°.

Open in new tabDownload slide

The method of solution proposed involves the expansion of the dislocation density in Legendre polynomials since we have assumed ω= 0. The asymptotic study shows that the singularity degree can assume all the values in the interval] -1, 0]. Therefore a solution can be derived considering another value for ω and in Fig. 12 we show the stress drop computed according to the expansion of the dislocation density in Chebyshev polynomials, the method used to obtain numerical solutions for the fault bending model. As expected, this method of solution has a slower convergence, since now for N= 10 (Fig. 12a) the stress drop is very variable and even increasing the order of truncation to N= 100 (Fig. 12b) the accuracy decreases approaching the interface.

Figure 12

Stress drop on crack sections computed according to the expansion of the dislocation density in Chebyshev polynomials with truncation order N = 10 (a) and N = 100 (b).

Open in new tabDownload slide

5 Conclusions

We have presented two models of non-planar cracks cutting across a rigidity discontinuity, which may be suitable to describe strike-slip faulting in the presence of a shallow soft layer. If we consider the fault bending model, a singularity appears in the dislocation density at the bending point. An asymptotic study shows that the dislocation density has an infinite singularity ˜r^ω at the interface with , depending on the dip angle δ = 90°-α of the upper crack section and on the rigidity contrast m. Once this singularity is properly accounted for, crack equilibrium requires that the stress drop values on either sides of the interface must satisfy eq. (51) to match the welded boundary conditions. Then, if the stress drop is assigned over both faults sections (according to fault rheology), a vertically dipping fault can cross the interface if the upper section is suitably inclined.

Considering the fault branching model, the study of the asymptotic behaviour leaves the singular degree of the dislocation density undetermined, with ω∈] -1, 0]. This indeterminacy has only a mathematical character and it can be interpreted in terms of a further degree of freedom introduced by the presence of the third crack surface. From a physical point of view, the simplest explanation is that no singularity is required at the interface and this guess is supported by the numerical results obtained employing the boundary element method. The stress drop on the three fault sections must satisfy two generalized conditions (67) and (68).

Large rigidity contrasts are typically found in the upper crust and the driving stress of strike-slip faults is strongly discontinuous across such layer interfaces. According to the previous models, we may assess that the presence of crustal layering calls for the formation of geometrical complexities in transcurrent fault regions.

Several field observations regarding the complex pattern of surface fractures connected with major strike-slip earthquakes may be explained in terms of such a mechanism. For instance, if we consider the case of the M_s= 6.6 earthquake of 2000 June 17 in Iceland (Fig. 1), the aftershock activity (dots) shows very well the NS orientation of the vertical main fault, confirming the fault plane solution, while lines of surface fractures (double lines) are observed laterally shifted up to ˜1 km from the deep fault strike.

Considering the fault geometry, since major rigidity discontinuities in the SISZ are found at 3 km and 1.5 km (Antonioli 2006), α can be inferred to be α˜18° or α˜ 30°, respectively, and the stress drop on the shallow section is negligible with respect to the deep section, suggesting that faulting propagated passively to the surface, as suggested also by the lack of aftershocks in the shallower 3 km (Stefànsson 2011). Similar considerations can be proposed also for the San Andreas Fault sections mentioned in the Introduction.

In the bending model, the inclination of the upper section determines an asymmetric stress release with respect to the negative z-axis, therefore the region with minor release is candidate to host new fractures. The stress singularity which is present close to the bending point (see insets in Fig. 5c) may be responsible for the localized microearthquake activity along the slip direction detected on strike-slip faults (e.g. Rubin 1999).

If the generalized stress drop condition (51) imposes low values of the inclination angle, it is reasonable to expect that only a bending process takes place since the singularity degree at the bending point is low (Fig. 3b). Instead, when very low stress drop values are applied on the upper section, the stress drop condition requires higher values for α and the higher values of the singularity should favour a branching process. In Fig. 1 we can observe that multiple surface fractures appear when their distance from the deep fault strike increases.

A limit case of fault branching takes place when a low level of the driving stress () in the upper layer does not allow the propagation of the main fracture across the interface and the σ_zy stress component in excess of the frictional threshold is released; then the unwelding of the interface takes place . In Ferrari & Bonafede (2011) a parametric study of this process was performed employing a boundary element method in a layered medium bounded by a free surface to study the stress redistribution following the unwelding of near surface layers.

It must be mentioned that the models presented in this paper apply to a strictly antiplane strain configuration, assuming laterally uniform properties within each layer. Even in this framework, complexities of fault regions might be generated according to other mechanisms, such as dynamic effects mentioned in the Introduction, the presence of pre-existent weakness planes or strength heterogeneities of shallow layers. Finally, the study of stepover zones, where a planar fault section bifurcates into two separate sections near the surface, requires intrinsically 3-D fault models and cannot be addressed employing the present semi-analytic asymptotic approach.

Acknowledgments

This paper originates from several discussions with our partners in the framework of the UE project PREPARED, focusing on the SISZ (FP5RTD). We thank Eleonora Rivalta and an anonymous reviewer for constructive criticism on a previous version of this paper.

Appendices

Appendix A: Asymptothic Study

The stress drops The stress dropsσ₁ and The stress dropsσ₂ must be bounded, so it is necessary to study the singular behaviour of the dislocation density at crack tips. To this end, we generalize to the system of coupled eqs (13) and (14) the method proposed by Erdogan (1973). To determine the asymptotic behaviour of the integrals, we introduce the following expressions for ρ₁ and ρ₂:

(A1)

(A2)

where R₁(± 1), R₂(± 1) ≠ 0, -1 < a₁, b₁, a₂, b₂≠ ω0 and t denotes either or .

Following Muskhelishvili (1953), we consider the function Φ(ς) of the complex variable ς

(A3)

The asymptotic behaviour of Φ(ς) in the neighbourhood of ς=±1 is assigned by the following expression:

(A4)

where the function Φ₀(ς) is bounded everywhere, with the exception of points ς=±1, where however its order of infinity is less than a_k or b_k (a_k≠≤ 0, b_k≠≤ 0), depending on the point considered.

A1 Singularities of ρ at crack tips

The asymptotic evaluation for the eqs (13) and (14) is obtained summing the asymptotic evaluation of each term

(A5)

(A6)

where

(A7)

and and are functions whose behaviour is similar to that of Φ₀(ς).

A2 Singularity at ξ = 1

To study the singularity of ρ₁ in the crack tip embedded in the upper media, we multiply both sides of eq. (A5) by and letting → 1, we find

(A8)

In the same manner we can study the singularity of ρ₂; multiplying both sides of eq. (A6) by and considering the limit ξ→ 1, we obtain

(A9)

In these points the singular behaviour of the dislocation density corresponds to that observed for cracks embedded in homogeneous media.

A3 Singularity at the interface

Multiplying both sides of eq. (A5) by A6) by and letting ξ-1

(A10)

we deduce that the following inequality must be fulfilled:

(A11)

Multiplying both sides of eq. (A6) by A6) by and considering the limit ξ→1,

(A12)

we find another inequality that must be satisfied

(A13)

The two inequalities found imply that

(A14)

To have a bounded stress drop on the crack surface near the interface, the following two conditions must be fulfilled:

(A15)

This linear system [R₁(-1), R₂(-1) unknowns] admit non-vanishing solutions if its determinant is equal to zero and so we find the following equation for ω:

(A16)

The previous equation can be expressed as the product of two factors

(A17)

that exchange their identities if we substitute α with -α and simply reflect the symmetry of the system with respect to the z-axes.

Appendix B: Stress Drop Condition

If we expand the dislocations densities using Jacobi polynomials

(B1)

(B2)

the integrals present in the system of eqs (13) and (14) can be evaluated using the following integral relation:

(B3)

where A6) by are Jacobi functions of second type that admit this representation

(B4)

that can be derived using the connections between various solutions of the hypergeometric equations (see section 2.9 in Erdelyi 1953) from the representation given in eq. (7.125) by Erdogan (1973). In eq. (B4), A6) by represent the gamma function while F(a, b; c; x) is the hypergeometric function

The Gauss hypergeometric function F(a, b; c; x) is an analytical solution in a neighbourhood of zero of the hypergeometric differential equation and it is uniquely determined for non-integer c by the normalization condition F(a, b; c; 0) = 1.

The stress drop on sections 1 and 2 can be evaluated summing only the second terms of the expression (B4):

(B5)

(B6)

since if ω is solution of (A16) then the compensation of singular terms is obtained. The previous equations can be rewritten in a more compact way

(B7)

(B8)

where

(B9)

(B10)

From eq. (B8) we find that

(B11)

if we insert this expression in eq. (B7) we can rewrite the stress drop on section 1 in this way

(B12)

It is to verify that

(B13)

(B14)

(B15)

(B16)

and so we can finally write the relation between Δσ₁(-1) and Δσ₂(-1) in this form:

(B17)

where

(B18)

References