Fault shape effect on SH waves using finite element method

In this paper, we analyse how the combination of fault zone shape and material properties affects the propagation of seismic waves in a two-dimensional domain. We focus on SH wave propagation through several faults with different thicknesses and bending radii, but the theory is easily generalized to the three-dimensional case. We show how the density of energy released is mostly a function of the radius and does not depend on the velocity inside a fault zone.

meters to kilometres in width, which is bounded by major faults within which subordinate faults may be arranged variably or systematically. Single fault zones are marked by fault gouge, breccias, or mylonites" (Allaby 2013). It is well known that a FZ is able to trap and propagate low-velocity waves and/or guide head waves, modifying substantially the damage pattern and the seismic hazard of a strong earthquake. These effects have been observed in many different real situations, see, for example, Rovelli et al. (2002), Li et al. (2004), Li and Leary (1990), Li et al. (1994), and Lewis et al. (2005). More recently some evidence of this phenomena has been observed, for example during the earthquakes (the largest being a 6.1 Mw event) that hit the city of L'Aquila, in 2009 (Avallone et al. 2014;Calderoni et al. 2010). In particular, in Avallone et al. (2014) the authors hypothesize the existence of a low-velocity FZ able to trap the waves and amplify the signal. A comparison between recorded data and an analytical model is also presented.
Due to the well established impact of the fault zones on the damage distribution, the last decades have seen more effort spent modelling wave propagation through a FZ and, although analytical solutions are available in just a few simple cases such as in Ben-Zion and Aki (1990), numerical results are numerous and quite complete. The common methods that are used for the numerical solution of earthquake models in general, and the ones with faults in particular, are described in Igel (2017) and Hori (2011). They include spectral element method (SEM), finite difference method (FDM) (Erickson et al. 2017;Moczo et al. 2014), boundary element method (BEM) and finite volume method or a coupling of several methods (O'Reilly et al. 2015). Finite difference method provides simple ways of discretizing partial differential equations which correspond to the seismological problems, but it still has a number of drawbacks. First of all, FDM is incapable of treating complex geometries adequately. Besides, when using FDM, one should choose small grid sizes to avoid numerical dispersion (Moczo et al. 2000;Liu and Sen 2009b). Although different techniques help to partially eliminate many drawbacks of FDM, this method is not suitable for solving certain seismological problems.
The important aspect of the spectral element method (Komatitsch and Tromp 2002) is the fact that it can combine the approximations of different order in different subdomains, in particular outside and inside the fault. This is done by introducing the discontinuous Galerkin methods or mortar element methods. As a result, on the basis of the spectral element method various software has been created which help model earthquakes and fault zones for 2D and 3D models.
Unlike the SEM and FDM, BEM builds the approximate solution to seismological problems based on integral representation of the solution. The method consists in building the solution to certain boundary integral equations. The main advantage of this method is that it helps to decrease the dimension of the problem. Unfortunately, BEM is not efficient for solving nonlinear or non homogeneous problems. Moreover, the final matrix of the problem is in general full and non symmetric. There are techniques which reduce the number of drawbacks for non stationary problems, though the construction of efficient methods for simulating time-dependent realistic earthquake problems remains an open problem (Álvarez Rubio et al. 2004).
In this paper, we study how seismic waves, generated by a localized source, propagate through a FZ and what happens when the fault has a bend. We only consider the case of horizontally polarized shear waves (SH waves), the general case including pressure waves (P-waves) and vertically polarized shear waves (SV waves) will be discussed in a further paper, where a more general 3D case will also be analysed.
The propagation of SH waves in a homogeneous isotropic linear elastic medium is described by the following equation in a (x 1 , x 2 , x 3 ) Cartesian space.
where u is the component of the displacement in the x 3 direction, V SH = √ μ/ρ is the shear wave velocity, and S is the body force, modelled as a delta-source where x = (x 1 , x 2 ), x s is the location of the source and is a moment-rate time variation and τ is a smoothness parameter related to the frequency which controls the amplitude of the oscillations (Igel 2017).
We consider a 2D model where the fault bends along the strike. Our computational domain is shown in Fig. 1. It can be read as a section of a possible 3D domain parallel to the surface. In this sense, the orientation of the fault is different from much of the literature (see Ben-Zion et al. 2003, Fig. 5). However, the main objective of this paper is to evaluate how the trapped energy is released by a FZ due to its bending. It is not important, at this stage, if it happens on the surface or at depth.
Moreover, this simple 2D model also allows a quantitative analysis (for example in the calculation of the parameter α introduced at the end of this section), that becomes much more complex in a 3D setting. The same problem will be addressed in a more realistic geological setting with 3D physically based simulations (preliminary results are available in Di Michele et al. 2021).
The energy lost as a consequence of the bend has been the subject of intensive study, especially in the field of optics, in order to avoid, or at least reduce the signal lost due to a bend in a waveguide.
The mechanism surrounding this phenomena is a simple application of Snell's law, which says that at the interface between two different materials, having With the label S we identify the source location, with the labels LS C , LS R , LS L we identify respectively the monitored points on the center, right and left side of the fault zone. For the curved fault the radius is labeled as r. The labels l in and l out refer to the cross-section of the fault reflectionindex n 1 and n 2 respectively, the incident and refraction angle are such that where V is the wave propagation velocity. There exists a critical angle for which the incoming waves are not transmitted but all the signal remains trapped inside the guide. For all θ i < arcsin(V 2 /V 1 ) the signal is partially transmitted in the medium labeled as 2. If the guide bends such that the incident angle exceeds the critical angle then all of the signal is reflected within the guide. Obviously in geology it is not so easy to predict and/or evaluate the quantity of energy which will be released, due to the uncertainties in the geometric and mechanical parameters and the source structure. However it could be important to understand how the trapped energy is released as a consequence of the fault geometry and of the mechanical properties of the surrounding media.
During and after an earthquake, as a seismic wave propagates in a real media a certain amount of its carried energy is lost through being converted into heat. This effect is usually modelled introducing an attenuation coefficient named α: where f is the frequency, V is the wave velocity inside the medium and Q is the so-called Q factor defined as where E is the peak strain energy and ΔE is the energy stored per cycle. The energy (density) E is the sum of the kinetic energy K and potential energy W , namely : where and In (9) τ ij and e ij represent the stress and the strain tensors, respectively (Shearer 2009). These quantities can be written as functions of the displacement and its derivative as follows According to our knowledge a specific theory of the bend loss for trapped seismic energy has not been yet developed. However, this phenomenon can also be described using Snell's law and it looks reasonable to assume that a particular geometric configuration, such as the bending of a FZ, can be associated with an extra energy release as for an optical cable. For this reason we introduce a bending attenuation constant α r , that is a complex function of the geometric and mechanical parameters of the problem. α r can be modelled in many different ways, depending on the approximation adopted to solve the propagation of waves in guiding structure (see, for example, Gloge 1976;Marcuse 1971Marcuse , 1976 and reference therein).
Here, due to the complexity of the FZ structure, we adopt the simplest approach by assuming where r is the radius of curvature of the waveguide and C 1 and C 2 are constants. In the optical framework these constants depend on many parameters such as magnetic permeability of the free space, propagation constant, refractive indexes and many others. In this context due to the lack of a suitable theoretical model we calculate C 1 and C 2 experimentally at the end of Section 3. In Section 3 the values of these constants will be estimated using the numerical results in some simple configurations. A third factor producing the seismic wave attenuation is the scattering, due to the presence of obstacles such as fractures and fluid filled pores that mostly constitute the fault zones. These irregularities and their distribution strongly contribute to the energy dissipation and to the seismic wave propagation. Many papers in the literature account for these effects, such as Jahnke et al. (2002) and Gulley et al. (2017ba, b) and many others. However, in our model, we decide to neglect scattering attenuation, and focus on the geometric contribution to the FZ shape. The contribution of the anisotropies, both in the fault zone and in the bedrock, will be included in a forthcoming paper still in preparation.

Numerical method
Let us consider (1) supplemented with the absorbing boundary conditions where n is a unit outer normal vector to the domain Ω.
As for the initial conditions, we assume that both the displacement and velocity at the starting time t = 0 are equal to 0, namely In order to obtain the weak formulation of problem (1) with boundary conditions (12) we first multiply equation (1) by the regular test function v, integrate over Ω, use Green's formula and then boundary conditions (12), which results in where Let us divide the domain Ω into triangles Ω K , K = 1, ..., n e and choose the approximate solution . Substituting the approximate solution into (14) and choosing v = φ i (x) one arrives to the semi-discrete system where Let For the time discretization we use the leapfrog method of the form with the first timestep It should be noted that the choice of the leapfrog timestepping method is motivated by the fact that this method is explicit, and therefore, finding the solution on each timestep does not require solving a linear system (in the case of the mass lumping described below that leads to diagonal mass matrix) unlike in case of implicit methods. This is especially useful since the large number of finite elements in our problem would make it too time consuming to solve the linear system at each timestep (Igel 2017).
Besides, as we are dealing with a hyperbolic equation, the use of explicit over implicit methods is typical because of the CFL condition (Maggio and Quarteroni 1994;Wendroff 1968;Strikwerda 2004).
Let us now introduce the spectral element nodes and basis functions that lead to a diagonal mass matrix M as shown in Cohen et al. (2001). In the case of quadratic elements for the reference triangle the nodes are defined as S i -the vertices of the triangle, M ithe midpoints of the triangle sides, G -its centroid. The basis functions w G , w S i 2i and w M i 2i in this case are described by where λ i are the barycentric coordinates and p S i 2i , p M i 2i are the standard P 2 basis functions on triangular elements.
Unfortunately, in the case of triangular finite elements, the matrix M+ Δt 2 B is non-diagonal. In order to avoid inverting the matrix M + Δt 2 B we proceed in the following way. Let us divide the finite element nodes in two sets: the ones lying inside the domain Ω and the ones lying on the boundary ∂Ω. On each timestep we find the approximate solution separately for these two sets. It can be seen that, for the points inside the domain, the corresponding matrix of the linear system is diagonal, so this part can be solved directly. As for the system that corresponds to the points lying on the boundary ∂Ω, the matrix of the related linear system is non-diagonal, but the dimension of this part is much smaller than of the original system. In other words, the space dimension of this part is reduced by one, since in this case we consider only the boundary nodes. As a result, if the original problem is two-dimensional, this part of the matrix M + Δt 2 B can be inverted without major difficulties.

Results and discussion
We consider a two-dimensional (2D) domain Ω having sides of 10 × 20 km, containing a curved or a straight fault as in Fig. 1. The bottom left corner of the domain corresponds to the origin of the Cartesian coordinates, that is ((0, 0) km). On all the boundaries absorbing conditions are prescribed as in (12).
The geometric parameters of the problem are the thickness of the fault h, the radius of the curvature r and the velocity ratio β, that is the velocity reduction inside the fault zone with respect to surrounding media. From the geological literature we know that h ranges between and hundreds of meters, whereas β is in the range 0.2-0.6 (Huang and Ampuero 2011;Li and Vidale 1996;Kuwahara and Ito 2002). Although β and h are, in general, variable within a real FZ structure, here we will assume that both are constant. From here we label all the physical and geometric parameters related to the fault zone with a subscript F Z and the parameters related with surrounding media with a subscript O.
Within this paper we consider a constant value for velocity in the surrounding media (V O = 1000 m/s), and different values of the velocity inside the fault zone, namely V F Z = 400, 600, 800 m/s. In other words we consider a velocity reduction of 60%, 40% and 20% with respect to the surrounding rock, that means β = 0.6, 0.4, 0.2 respectively. We also point out that in our framework the parameter choice refers to a superficial layer on the crustal model.
Three different values of the fault thickness are also considered, h = 0.25, 0.5, 1 km. We remark that a smaller thickness is theoretically possible but adds numerical difficulties as more finite elements and smaller timestep may be required to get favourable results.
Concerning the radius of curvature r, this parameter is more difficult to estimate because the fault zones are composed of a sequence of cracks and their structure is often very complicated and not well known. Here we consider r = 1 km and r = 4 km. The source, located at (5.0, 3.5) km, is modelled as a delta-source using (2) and (3). Obviously it does not represent a real earthquake, but it is able to generate well trapped waves inside our structure, as we will show later. We remark that there are two parameters which characterize the source, namely M 0 and τ . The first controls the amplitude of the incoming waves, the second their frequencies, which we choose as τ = 0.05 s. The behaviour of the source together with its fast Fourier transform is given in Fig. 2a. Observing Fig. 2b it is clear that the source contains all frequencies, with most of the spectrum in the range [0-5] Hz. Low frequency f means a large wavelength λ ( f = V /λ).
The wavelengths corresponding to f = 5 Hz for the V F Z considered within this paper are λ = 80 m for V F Z = 400 m/s, λ = 120 m for V F Z = 600 m/s, and λ = 160 m for V F Z = 800 m/s. Therefore, all the wavelengths taken into account are small enough to see a fault zone having thickness in the range 0.25-1 km and we expect trapped waves to be observed.
Finally we remark that due to the linearity of the problem, M 0 plays the role of a scale factor. Here we set M 0 = 10 16 Nm.
The signal, the displacement u in this case, is detected at three seismograph locations (LS). We name the LS point located inside the fault LS C , whereas LS R and LS L are on the right and left side of the fault respectively, at distance (in the normal direction) of 500 m from the fault zone boundary. Some related test with distances equal to 50 m and 150 m will also be discussed at the end of this Section.
As we want to focus on the effects of the fault shape on the trapped waves, we locate the point source at the beginning of the fault far enough from the bend and boundary of the domain, in order to allow the wave to be well trapped and to avoid reflection from the domain boundaries.
The source point is located exactly in the middle of the fault although it is known that the maximum capability of the fault to trap waves can be observed for sources located at the fault boundaries. This is because in our model it is very important to have a symmetric source to observe the effect of the curvature. In other words, for a straight fault no differences can be observed for the displacement recorded by LS L and LS R due to the symmetry of the problem and all the differences recorded for the curved shape can be attributed mostly to the curvature itself.
However, for the sake of completeness we have also investigated the effect of the source position, the results of which are displayed in the Appendix.
First of all we verify that, as expected, the waves are well trapped before they reach the curve. We compare the displacement observed at LS C for a straight fault with the signal detected at the same location in a plain domain. We analyse in detail the case corresponding to V F Z = 600 m/s. The displacement has been recorded for all the three values of the thickness at LS C and LS R (see Fig. 3). For the straight fault, after the first peak arrives at LS C , large oscillations in the displacement can be observed for many seconds, for all values of the FZ thickness. On the contrary, in the homogeneous case after the first waves arrive, the displacement decays quite quickly. A similar behaviour of the displacement can be observed at LS R , the oscillations have a similar frequency but smaller amplitude compared to those inside the fault, this is due to the  Fig. 14a and b in the Appendix. Finally, we observe that when decreasing the velocity inside the FZ the trapped waves become less compact as observed in Li and Vidale (1996).
As mentioned before, in our framework, the simulation parameters refer to the crustal layer.  Fig. 4, where we can observe a displacement behaviour qualitatively similar to the one already described above for lower wave speeds.
Having verified that the fault zone, as we designed it, is able to trap waves, we start the parametric analysis of our model focusing on the effects of the A decrease in thickness leads to an increase in the frequency of the trapped waves, as expected. The displacement, recorded at the LS R and LS L , for each of the velocities are plotted in Fig. 5 for a radius of curvature r = 1 km and r = 4 km, and for the straight fault, namely r ∞ . The first peak has a similar amplitude of approximately 0.6 cm. However for r ∞ the oscillation decays in amplitude quite fast, whereas for a bending fault the oscillation remains almost constant in amplitude for 10 s. On the opposite location, namely on the right side, for all three radii we observe similar behaviour. Similar effects can be observed for the other two velocities taken into account (see Fig. 15a and b in the Appendix), except the case V F Z = 400 m/s and r = 4 km where the bending effect is smaller. Indeed according to the definition of the critical angle (4), θ c decreases as V F Z decreases, and a larger amount of energy remains trapped, reducing the effect of the bending.
As noted in the previous section, the displacement in the direction of the LS L is larger than that on the right side of the domain. Therefore, a strong release of energy from the curved faults in the left direction can be expected.
To quantify this effect we plot the integral kinetic energy normalized with respect to the density ρ, namely I K (t) = t 0u 2 . We studied the amount of normalized kinetic energy I K (t) transmitted by the fault to LS L and LS R for three different values of fault thickness h = 0.25 km, h = 0.5 km and h = 1 km. All the combinations of the fault zone velocity and radius of curvature are considered. We discuss, in detail, the case corresponding to the fault thickness of 0.5 km, plotted in Fig. 6 (plots related to the other two cases are reported in Appendix Figs. 16 and 17).
In the case r = 1 km we get I K | LS L > I K | LS R for all the noted velocities. At the beginning of the curve, a transient time for which I K | LS R > I K | LS L is observed. This effect is due to the first peak arrival time that is smaller for LS R because this recorder is closest to the source. Similar observations can be done also for large curvature, that is r = 4 km, except for V F Z = 400 m/s, where the behaviour is completely different and I K | LS R > I K | LS L almost everywhere. As we have observed before in this case the effect of the bending is reduced and the effect of the Similar behaviour can be observed in the Appendix for the other fault thicknesses (h), as displayed in Fig. 16, for h = 0.25 km, and Fig. 17, for h = 1 km.
Let I K (T p ) be the plateau value for the function I K (t). We remark that the plateau is reached at a different time T p for each FZ velocity we select. I K (t) represents the kinetic total energy (normalized by the density) released in a certain point. In Fig. 7 the plateau value is plotted as a function of the radius of curvature, for all the velocities V F Z and the FZ thickness.
The behaviour of the energy and of its cumulative function depends on many parameters, such as geometry and mechanical properties of the domain as well as source frequency spectrum. However some common threads can be observed in Fig. 7  anddecays as the radius of curvature increases, c) I K (T p ) LS L > I K (T p ) LS R . Concerning the last point, we remark that all the data presented within this paper refer to 0.5 km distance between the LS points and the fault boundaries, although we expect that the curvature effect will increase when that distance is reduced. Therefore, we compare the cumulative energy released by the fault at LS R and LS L for different LS-fault boundary distances (500, 150 and 50 m). In particular we consider V F Z = 600 m/s and h = 0.5 km. Results are displayed in Fig. 9. We see that as the distance increases the energy released is increased and the fault bending effect becomes more pronounced. According to Fig. 7 (and Figs. 18 and 19 in the Appendix) the trigger parameter for the energy lost is the curvature radius of the fault, the ratio β between velocity inside and outside the FZ does not play any role. To verify this effect in Fig. 8 ( and Figs. 20 and 21 in the Appendix) we plot I K (T p ) as a function of the velocity inside the FZ. Also in this case some common threads can be observed: a) I K (T p ) LS R and I K (T p ) LS C decay quite fast, although some oscillation is noticeable for small velocities at LS R , b) I K (T p ) LS L looks velocity independent and remains almost constant in a wide range of velocities, c) I K (T p ) LS L > I K (T p ) LS R for all velocities if r = 1 km (Fig. 9).
Finally, to verify the formula (11) we compute numerically the parameter α r , as where l out E out dl out dt are the integrals of incoming energy E in and outgoing energy E out (given by (7)) that pass through the cross-sections l in and l out respectively (Fig. 1). In other words, I E in and I E out are calculated as the line integrals over time along the corresponding cross-section using numerical quadrature. The bending attenuation constant α r is then obtained as the ratio between the integrals of total energy lost per unit fault length through the curved part of the fault and the integral of incoming total energy I E in . Since the length of the curved part of the fault is equal to πr 2 , in order to obtain the energy lost per unit fault length of the curved part, we have to normalize it by the factor of πr 2 . The exponential behaviour is well fitted by the numerical data at least for the higher velocities 600 m/s and 800 m/s (see Fig. 10). The constants C 1 and C 2 are estimated, for τ = 0.05 s, using a SCILAB function for a data fitting. Test

Conclusions and further studies
In the paper we have investigated the effect of the fault geometry on the SH waves propagation. Contrary to the available literature (for example Ben-Zion 1989, 2003 where the 2-D vertical section is analysed, in this study the fault is located at a certain depth in the horizontal plane to better appreciate the effect of the bending. In other words we analyse a horizontal section of the fault zone instead of the usual vertical one. We have investigated bending effects while varying fault zone thickness and found (see Figs. 5 and 15 in the Appendix) that the fault zones are able to guide trapped waves through curved geometries according to the stated aim of this paper. In particular we have shown that a greater amount of energy is released in the direction of the negative curvature of the fault. According to our study, the velocity inside the fault does not affect the energy release in a strong way (see Figs. 7c and 18c and 19c in the Appendix).
However it strongly depends on the radius of curvature (see Figs. 8c and 20c and 21c in the Appendix). We treated only simple 2D geometric cases, a deeper understanding of the phenomena could be obtained from an extension into 3D and the addition of realistic geological features inside the studied domain. This work was performed as a preliminary study into potential fault guide effects, which have been proven, before undertaking more complicated tests on a more realistic domain. A new set of tests will be performed by using SPEED (SPectral Elements in Elastodynamics with Discontinuous Galerkin-http://speed.mox. polimi.it/). SPEED is an open-source code able to simulate seismic events in a three-dimensional reconstructed domain (Mazzieri et al. 2013;Antonietti et al. 2012). Using this tool we will be able to include complex fault zones and we will evaluate in more detail the effects of the fault zone shape on wave propagation.

Appendix
A.1 Derivation of the energy conservation law for seismic waves The energy density carried by seismic waves can be expressed as a sum between the kinetic energy density (K) and the potential energy density (W ), as E = K + W while K, as usual, can be written as K = 1 2 ρu 2 , the derivation of W (strain energy functional), is classically based on mechanical and thermodynamic considerations (Aki and Richards 2002). Therefore, for convenience of the reader, in the following, we sketch the derivation of W (see Aki and Richards 2002 for more details). In our case we consider the system to be adiabatic, namely that there is no heat exchange with the exterior. Let B an elastic body having volume V and boundary surface S. The variation of the total energy (kinetic + potential) per "infinitesimal time"Ė, coincides with the variation of the related mechanical energy.
In view of Gauss's divergence theorem it is not difficult to show thaṫ where u denotes the displacement vector and τ ij and e ij represent the stress and the strain tensors, respectively. As mentioned previously the kinetic contribution to the energy increase follows the standard law of classical mechanics; therefore, In view of the expressions (28)-(29) the internal energy U for adiabatic process, can be written in differential form as hence U = W. Then, for instance in the Hookean case, after some simple calculations it follows that A.2 Source position and bending angles effects Within this work all of the tests use a source located at the center of the fault zone at the same distance with respect to the left and right fault boundary. In this way the system is symmetric except for the fault bending. To evaluate the effect of the source we consider four different sources at different positions as in Fig. 11. The results are displayed in Fig. 12, for two radii of curvature r = 1 km ( Fig. 12a and b) and r = 4 km ( Fig. 12c and d). The fault thickness is 1 km. For both LS L and LS R the behaviour of the recorded signal is quite close to each other, for S 3 this presents remarkable differences due to the fact that the source is located outside the fault in a material having different mechanical properties. Finally we discuss the effect of angle of bending that in Section 3 is assumed to be 90 • . Here we consider and compare 3 different values of the bending angle, namely 90 • , 60 • and 45 • , to better approximate the real fault zone total bending angle (Fig. 13). According to this preliminary results for the selected frequency range, the effects of the bending values remain remarkable for angles greater then 60 • , but different set of mechanical and physical parameters must be considered to better understand the phenomenon not only with respect to the main curvature of the FZ, but also with respect to the often more pronounced local curvatures.

A.3 Other tests
In this section we collect the figures mentioned in the main text and moved here for greater readability. All the pictures are cited in Section 3 and carefully described in the captions Funding Open access funding provided by Gran Sasso Science Institute -GSSI within the CRUI-CARE Agreement. This work was partially supported by the GSSI "Centre for Urban Informatics and Modelling" (CUIM). Italian Government (Presidenza del Consiglio dei Ministri) CUIM project (delibera CIPE n.70/2017). Andriy Styahar was supported by funding from MathMods and InterMaths projects of the University of L'Aquila.

Conflict of interest
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommons.org/licenses/by/4.0/.