INFLUENCE OF HILLS AND VALLEYS ON FREE SURFACE SLIP AMPLIFICATION IN DYNAMIC RUPTURE ON A DIPPING FAULT

1. INTRODUCTION

Slip at the intersection of a fault with a free surface determines the amount of distance by which a structure across the fault gets torn apart and is separated from the portions on each side. Higher slip leads to a larger offset and in case of dipping faults can introduce both vertical and horizontal offsets. The past studies exploring the effect of free surface on rupture dynamics were confined to flat terrain for simplifying the numerical simulations [Chen, Zhang, 2006; Zhang, Chen, 2006a, 2006b; Day et al., 2008; Kaneko, Lapusta, 2010; Xu et al., 2015]. However, the undulations of the earth’s surface differ from one place to another and exhibit an irregular pattern of topography. [Ely et al., 2010] demonstrated that the surface topography has a significant effect on rupture dynamics by modelling the southern San Andrea Fault rupture. [Zhang et al., 2016] exclusively studied the effect of topography on rupture propagation and sub-shear to super-shear transition by comparing rupture dynamics on faults with different topographic features in their models. [Huang et al., 2018] investigated the effect of canyon-shaped and hill-shaped topography on near-fault ground motion using curved-grid-finite-difference models (CGFDM) and concluded that the canyon-shaped topography has a stronger effect on near-fault ground motion compared to the hill-shaped topography. Understanding whether the presence of such topographic features as hills or valleys in case of faults reaching the free surface is a boon or curse would be useful since the construction of any structure across the fault boundary is inevitable due to either lack of land elsewhere or habitat resource availability. Our preliminary results depict that topography in the vicinity of the fault can change the magnitude of slip along the fault. Fences along some bordering countries could also be found among these structures but the damage to a fence may not necessarily have a profound impact on the economy of any country.

Earthquake modelling can range from a simple double couple point source to spontaneous rupture propagation on finite faults governed by frictional contact laws. Many researchers in the past years adopted numerous numerical and analytical methods to solve both point source and spontaneous rupture problems [Zhang, Chen, 2006a, 2006b; Andrews, 1976; Madariaga, 1976; Das, Aki, 1977; Fukuyama, Madariaga, 1998; Aochi et al., 2000; Oglesby et al., 2000a, 2000b; Day et al., 2005; Benjemaa et al., 2007, 2009; Kaneko et al., 2008; Ely et al., 2009; Hok, Fukuyama, 2011; Tago et al., 2012; Zhang et al., 2014; Parla, Somala, 2022a, 2022b; Parla et al., 2022; Somala et al., 2022]. While an earthquake modelled as a spontaneous rupture is dependent on the initial conditions which are often insufficiently known, the spontaneous rupture propagation popularly termed as a dynamic rupture is by far the best possible way of modelling an earthquake. [Hok, Fukuyama, 2011] have modelled the dynamic rupture propagation on shallow dipping fault using the boundary integral equation method (BIEM) to study the effect of free surface on rupture characteristics by introducing virtual free-surface elements. Kinematic ruptures, based on a statistical relation for fault slip obtained from kinematic inversions of earthquake data, will just give as much slip as we prescribe at the junction of a topographic feature and a fault, and are thus considered as if they were not used to understand differences in a slip near the free surface due to topographic feature, while keeping all other conditions the same.

The Southern California Earthquake Center (SCEC) has put together several benchmarks for validation of numerical implementations of dynamic rupture by the researchers [Harris et al., 2009]. Dynamic rupture simulations in this paper are performed using PyLith [Aagaard et al., 2013], a code that is benchmarked with the SCEC dynamic rupture benchmarks. This study addresses the free surface slip amplification owing to the influence of topographic features of a certain class by choosing a TPV10 benchmark from the SCEC dynamic rupture benchmarks. In this study, we constrain ourselves to 2D cases.

2. METHODOLOGY

PyLith, the software chosen to run dynamic rupture simulations in this study, is based on the finite element method of converting strong form of elasto-dynamic equation to weak. For brevity, we outline the solution method for quasi-static problems using only index notation. In quasi-static problems, with neglect of the inertial terms, the reduced weak form of equation will be:

We solve this equation through formulation of a linear algebraic system of equations (Au=b), involving the residual (r=b–Au) and Jacobian (A). The residual is simply:

Combining the expressions for the increment in stresses and making use of the symmetry of the weighting functions, the Jacobian matrix will be:

where σ_ij – stress tensor field, V – volume bounded by surface S, T_i – traction vector field, u – displacement, t – time.

The detailed solution method for quasi-static problem can be referred to PyLith manual [Aagaard et al., 2010] (section 2.3). The domain is discretized into triangular elements and the field variable (displacement in this case) is considered to be varying linearly within the element. Integrals are used as approximation to summation over all the elements and summation over function values evaluated at Gauss quadrature points multiplied by corresponding weights and the Jacobian. The weak form so obtained is solved by minimizing the residual in an iterative fashion after discretizing in time using Newmark’s explicit time-stepping scheme [Aagaard et al., 2010]. While this gives the equations for the bulk domain in terms of displacement, the contact friction conditions are imposed on the fault such that friction is governed by slip-weakening law and interpenetration is avoided. The fault is replaced by zero volume cohesive elements which have additional degrees of freedom in terms of Lagrange multipliers [Aagaard et al., 2010]. These Lagrange multipliers are the energy conjugate of displacements, the tractions. Excess shear tractions over and above the frictional strength get converted into the slip in the direction of resultant shear traction. This constraint is solved together with the weak form discretized in time to obtain displacement field and the spatiotemporal evolution of slip on the fault.

3. REFERENCE CONFIGURATION

The material properties of the model and the initial traction parameters on fault are taken as per 2D version of TPV10 SCEC (Southern California Earthquake Centre) benchmark model. The 2D version of TPV10 SCEC dynamic rupture benchmark [Harris et al., 2009], taken as the reference point, is shown in Fig. 1. The domain is a uniform half-space, and the medium is characterized by the P-wave velocity of 5.72 km/s and S-wave velocity of 3.3 km/s. Mass per unit volume of the material in the medium is 2.7 g/cc. The 2D fault reveals a mechanism of thrust fault rupture with 60° dip and 90° rake angles, reaching the free surface with a length of 15 km. Plane strain assumption is used with a vertical cross-section of the TPV10 benchmark to solve its 2D version. The rupture process on fault is governed by slip-weakening friction law [Ida, 1972]. The linear slip-weakening friction model produces shear tractions equal to the cohesive stress plus a contribution proportional to the fault normal traction that decreases from static to dynamic value as slip progresses. The friction stress on fault will follow the equation:

where µ_s is static coefficient of friction, µ_d is dynamic coefficient of friction, d₀ is slip weakening parameter, d is cumulative slip, T_c is cohesive stress, T_n is normal stress, and T_f is shear traction on fault.

Nucleation is carried out on a 3 km patch located at 12 km down-dip having static and dynamic friction coefficients of 0.76 and 0.448, respectively, and a critical slip-weakening distance of 0.5 m.

Cohesion on the fault is chosen to be 0.2 MPa. Meshing is done with 200m triangular elements. Failure follows the linear slip-weakening law [Ida, 1972]. Initial tractions are assumed to follow a depth-dependent behavior where normal traction increases as 7378 Pa/m along with down-dip distance and the ratio of shear traction to normal traction is 0.55 except on a 3 km nucleation where the shear traction is 67.99 MPa. Shear traction exceeds frictional strength only on the nucleation path where rupture starts and propagates spontaneously.

4. TOPOGRAPHIC MODELS AND RESULTS

4.1. Semi-circular topographic feature

To understand the free surface slip amplification due to a pure upward undulation (hill) and a pure downward undulation (valley) near the fault, we first consider four cases as shown in Fig. 2. The shape of a valley as well as that of a hill is taken as a 1 km radius semi-circle at first and then as a 2 km radius semi-circle. The topographic features considered in this study (including hill and valley) are having dimension as shown in Table 1. The final slip as a function of horizontal distance for these four cases is also plotted in Fig. 3, along with the final slip for the reference case. It can be seen from Fig. 3 that the slip is amplified by a valley like a feature near the intersection of a free surface with the fault. More slip is observed with an increase in the semi-circular valley feature. For a 1 km radius semi-circular valley, peak slip is still at the nucleation patch similar to the reference configuration. However, in case of a 2 km radius semi-circular valley, peak slip occurs near the free surface and is almost twice that of the reference configuration. This could be due to the narrow corner formed near the intersection of a free surface with the fault and a potentially greater surface area of the semi-circular feature on the path taken by the waves to travel from the rupture propagation and to be reflected. Bump feature, especially the one that is 2 km in radius, caused a clear reduction of slip at the free surface. Peak slip, in this case, occurs near the nucleation patch, similar to that of reference configuration, and is lower for overall up-dip regime beyond the nucleation patch when compared with the reference configuration. The bump feature with a 1 km radius shows a mixed pattern. While there is a reduction in the slip, compared to that of reference configuration, just beyond the nucleation patch there is an up-dip direction – a reverse trend at the intersection of a free surface with the fault. Slip near the free surface for this case is found to be more than that observed for the levelled ground reference case.

4.2. Gaussian topographic feature

Though the semi-circular features in topography are well described mathematically, naturally observed hills and valleys are more Gaussian-shaped (Fig. 4). All Gaussian profiles referred to herein are normal distributions. The equivalent of the radius of a semicircle, considered here, is half-width of Gaussian function covering 99.7 % confidence interval. Such a confidence interval is the same as six times the standard deviation (σ). So, the radius referred to in section 4.1 is 3σ. However, the largest amplification is observed for σ=1/3 valley in case of Gaussian topography, unlike that of semi-circular topography where the largest amplification was observed for a 2 km radius valley. Using the normalized Gaussian profile, this difference can be explained by the fact that the σ=1/3 Gaussian valley, exposed to the waves travelling from the up-dip rupture propagation, has a larger surface area than σ=2/3 Gaussian valley (Fig. 5).

5. EFFECT OF TOPOGRAPHY-TO-FAULT DISTANCE

5.1. Semi-circular hill

To understand the slip amplification (or de-amplification) on the fault surface due to upward undulation (hill) ocurring at different distances, we use four models as shown in Fig. 6. This section addresses the effect of distance from hill to fault by considering a 1 km radius semi-circular hill located at 0 km, 2 km and 4 km distance from the tip of the fault. The variation in final slip with horizontal distance for all three cases is shown in Fig. 7, along with the final slip for the reference case. It can be observed in Fig. 7 that the free surface slip is decreasing if a hill-like feature is located far from the fault. It is also observed that peak slip occurs at nucleation patch for all four configurations. In comparison to a hill located at 0 km from the fault, the slip at the intersection of fault with a free surface is considerably less in the reference model, as well as in comparison to 2 km and 4 km distant hill-like feature models.

Unlike purely downward undulation (valley) cases, in topographic environments as those in Fig. 6, there are no potential edge-like features to cause wave reflections and slip amplification on fault surface due to the interference of edge-reflected waves. And this could be the main reason for little-observed amplification of slip for configurations shown in Fig. 6. Further, a trend in variation of slip value after nucleation patch for the case of a 4 km distant hill is the same as that for the reference case, thus depicting that the effect of topography will be insignificant if hill-like features form far from the intersection of fault with surface. As discussed, the hill-like feature with a 1 km radius at 0 km from the intersection of fault with surface shows a mixed pattern consisting of slip variation and horizontal distance patterns.

5.2. Semi-circular valley

As shown in Fig. 8, the study of three configurations with a 1 km radius semi-circular downward undulation (valley) located at 0 km, 2 km, and 4 km from the intersection of fault with a free surface addresses the effect of distance on slip amplification. The final slip values from the numerical simulations for all three cases are compared with the reference case. The final slip on fault as function of a horizontal distance for the three configurations is shown in Fig. 9 along with the reference case. Unlike the case of semi-circular upward undulation (hill), the presence of semi-circular downward undulation (valley) amplifies the slip on fault for all three different valley locations. In Fig. 9 it can be observed that the slip amplification decreases linearly with increasing distance from valley to fault. The case with a valley-like feature at 0 km from the fault exhibits the maximum slip amplification as compared to the cases with 2- and 4-km distant valleys.

Furthermore, the effect of fault-to-feature distance on slip amplification in downward undulation (valley) is more prominent than in upward undulation (hill). Except for valley located at 0 km from the fault, all the cases showed peak slip near nucleation patch on fault. The effect of interaction of edge-reflected waves is insignificant if the valley-like features are more distant and result in minimum amplification of slip on fault. The gap between the slip values of reference case and those of valley cases at the free surface intersection is larger than that at nucleation patch: in other words, the amplification of slip at the free surface intersection is more enhanced than at nucleation patch. Reason for this could be the reflected waves interacting with slip propagation at the free surface intersection rather than with slip propagation at nucleation patch which is located at a greater distance from the reflection-causing feature.

5.3. Gaussian hill

Since natural profiles of hills and valleys are more Gaussian-shaped and located at different distances in the proximity of fault, the effect of distance on slip amplification was studied with three models of Gaussian upward undulation (hill) with a half-width of 1 km and a 99.7 % confidence interval, located at 0 km, 2 km and 4 km from the fault as shown in Fig. 10. Fig. 11 shows the the plotted changes in fault slip values with horizontal distance for all the cases. As can be seen from Fig. 11, for the case of the hill 0 km distant from the fault, the amount of final slip on fault at the intersection with a free surface is larger than that for the reference case and two other cases. The slip values at nucleation patch are nearly the same for all the configurations and the peak slip occurred at nucleation patch in all four cases. Fig. 11 shows that the effect of hill-to-fault distance on slip amplification is insignificant if the distance is beyond a certain limit. As can be seen, slip as a function of horizontal distance is similar for reference case and for the case of the hill 4 km distant from the intersection of fault with a free surface.

5.4. Gaussian valley

Gaussian-shaped downward undulation (valley) 0 km, 2 km and 4km from the intersection of fault with a free surface, as shown in Fig. 12, was studied in order to understand the effect of distance on slip amplification. The final slip amounts from the numerical simulations are compared with reference case and shown in Fig. 13 as a function of horizontal distance. Unlike upward undulation, the Gaussian valley-like topography amplifies the slip on fault in all three cases, and the maximum amplification can be observed for the valley located at 0 km distance. As shown in Fig. 13, the amplification is less for the case of maximum fault-to-valley distance. Except for the valley located 0 km from the fault, peak slip near nucleation patch on fault is observed in all cases. It is clear, that the effect of distance on slip amplification is larger for the cases of downward Gaussian undulation than for the cases of upward undulation. For the case of a 4 km distant valley, the slip amplification started at a horizontal distance of 1 km from nucleation patch, whereas for other cases the amplification started immediately from the nucleation patch in up-dip direction. Except in the case of a 0 km distant valley, the peak slip values are observed near the free surface intersection.

6. CONCLUSIONS

A valley or a depression in the footwall of a dipping fault reaching free surface can be a curse, as there could be an amplification of slip near the intersection of a free surface with fault. Higher slip here is used as a proxy for larger damage in evaluating the curse or boon due to topographic feature. The curse could possibly be worsened by excavating deeper in the footwall, particularly with a profile that is near vertical on the side of a trough closer to the fault. Presence of a hill on footwall near the intersection of fault with a free surface can be a boon or a curse depending on the dimensions and location of the hill, while boon intensity is insignificant compared to the curse a valley can be. In the case of a valley, the effect of fault-to-hill or fault-to-valley distance on slip amplification is larger than that in the case of a hill. The major limitations of the study include the simplicity of the relief models which only have a hill or a valley at a time, whereas in reality both positive and negative relief forms do exist together in the fault zone. This paper can be extended to evaluate the effect of multiple topographic features as boon or curse for other dip angles and different initial conditions.

7. ACKNOWLEDGEMENTS

The authors thank Indian Institute of Tropical Meteorology (IITM) for giving access to Aaditya HPC Supercomputing resources to run PyLith.

8. CONTRIBUTION OF THE AUTHORS

Both authors made an equivalent contribution to this article, read and approved the final manuscript.

9. DISCLOSURE

Both authors declare that they have no conflicts of interest relevant to this manuscript.