Equivalent formulae of stress Green’s functions for a constant slip rate on a triangular fault

We present an equivalent form of the expressions first obtained by Tada (Geophys J Int 164:653–669, 2006. doi:10.1111/j.1365-246X.2006.03868.x), which represents the transient stress response of an infinite, homogeneous and isotropic medium to a constant slip rate on a triangular fault that continues perpetually after the slip onset. Our results are simpler than Tada’s, and the corresponding codes have a higher running speed.


Introduction
The boundary integral equation method (BIEM) is one of the most effective ways for studying the rupture propagation on faults. The method starts with the representation theorem (Aki and Richards 1980) to establish the relationship between the stress and discontinuity of displacement with a appropriate discrete scheme. Then we can obtain the processes of the rupture propagation combining a specific friction law.
Recently, Aochi (1999), Aochi et al. (2000a) and Tada et al. (2000) acquired explicit expressions for the stress Green's functions representing the transitory response of a 3-D infinite, homogeneous and isotropic medium on arbitrary curved fault. Tada discretized the faults into rectangular-shaped elements with the constant slip rate (Tada 2005) and triangular-shaped elements (Tada 2006) in order to simulate the complex system of the 3-D non-plane faults.
In this paper, we have obtained rigorous expressions of boundary integral equation on a triangular fault in an infinite, homogeneous and isotropic medium, which are equivalent to Tada's work, but we circumvent transformation of coordinates and miscellaneous nomenclature in the Tada's paper and the terminal form of our result is simpler than Tada's. It is easy for us to generate the codes using these formulae, which have a higher running speed compared with codes using Tada's expressions, due to the alternative method calculating the function E(c).

Definition of the problem
We assume the slip occurs on a triangle fault element with its vertices named A, B and C. We assume the medium as a 3-D infinite, homogeneous and isotropic medium. We transform the triangle ABC inside the plane z ¼ 0 by coordinate transformation. The vertices A, B and C are aligned counterclockwise in the order, with the coordinate value Aða 1 ; a 2 ; 0Þ, Bðb 1 ; b 2 ; 0Þ and Cðc 1 ; c 2 ; 0Þ. We intend to calculate the stress response of the medium at the receiver Rðx 1 ; x 2 ; x 3 Þ to a hypothetical slip of the form in which Du i ðn; tÞ is the discontinuity of displacement in the i direction at location n and time t, D i the corresponding discrete discontinuity of displacement in the i direction in the triangular fault ABC for the semi-infinity duration t [ 0, and HðÁÞ the Heaviside function. Throughout the article, dðÁÞ represents Dirac function, d ij the Kronecker symbol, u i ðx; tÞ the displacement in the i direction, r ij ðx; tÞ the stress component in the ij direction, S the whole surface of the fault, l the rigidity, c L the P wave velocity, and c T the S wave velocity. We denote by HðDABCÞ a function which equals 1 when the slip source n lies within ABC and 0 otherwise. We also define the symbols: Roman subscripts are meant to run 1, 2 and 3; and Greek subscripts are meant to run over 1, 2.
P abc represents the circular sum of a 1 , a 2 , b 1 , b 2 , c 1 and c 2 . Our nomenclature is very similar to Tada's.

Calculating the stress Green's functions
Our starting point is the expressions of Tada et al. (2000, eqs. (51)- (54)). Substituting Eq. (1) into these equations, with a series of algebra calculations, we obtain with a ¼ 3 À a and f ¼ 3 À f. a and f need to satisfy Einstein's summation convention. In Eqs. (4)-(8), and EðcÞ in which c represents the velocity of P wave or S wave. Other components of the stress tensor are not mentioned such as r 13 , r 23 and r 21 . They may be obtained by considering the symmetry between the coordinates x 1 and x 2 . We notice that Eqs. (4)-(20) are equivalent to Tada (2006, eqs. (6)- (25)). We can obtain the explicit BIEM if we succeed to calculate the definite integrals.
3.1 Calculating B if ðcÞ, Q if ðcÞ and R 3ijf ðcÞ We will solve these integrals together by three steps.
In order to solve the elementary integral, we use a series of methods of substitution. We omit the dreary process and just list the eventual results we need.
Earthq Sci (2017) where where A, B, C are constant through the calculation.
The expressions of B if ðcÞ, Q if ðcÞ and R 3ijf ðcÞ can be eventually written as

Calculating E(c)
Via the property of d function, we obtain where Dh is the sum of arc angles corresponding to the parts of the circle, centered at point ðx 1 ; x 2 Þ and with a radius r 1 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi c 2 t 2 À x 2 3 p . We use a different method to obtain the Dh based on analytic geometry. We express the point on the circle by n 1 ¼ x 1 þ r 1 cos h n 2 ¼ x 2 þ r 1 sin h & and present the equation of the straight line AB by the symbol f AB ðx; yÞ. So, f AB ðx; yÞ ¼ ða 1 À b 1 Þðy À b 2 Þ þ ðb 2 À a 2 Þðx À b 1 Þ; f BC ðx; yÞ ¼ ðb 1 À c 1 Þðy À c 2 Þ þ ðc 2 À b 2 Þðx À c 1 Þ; f CA ðx; yÞ ¼ ðc 1 À a 1 Þðy À a 2 Þ þ ða 2 À c 2 Þðx À a 1 Þ: So the statement that the point ðn 1 ; n 2 Þ is inside the triangle DABC is equivalent to a system of inequalities Using the universal formulae t ¼ tan h 2 , the first inequality can be transformed as We can obtain another two inequalities similarly. We solve the system of the inequalities in order to acquire the range of t, which can be used by h ¼ 2 arctan t to get the range of Dh: It is essential to point out that the code generated by this method is much quicker than that generated by Tada.

Correctness validation
In order to corroborate our expressions, we should compare the numerical results with Tada's. We choose the same parameter values in Tada (2006). The fault is divided into two triangular faults ABC and BCD, with their vertices located at A(0, 0, 0), B(1, 0, 0), C(1, 1, 0) and D (0, 1, 0). The receiver point R locates at (0.6, 0.3, 0.2) (Fig. 1). The velocity of the P wave is ffiffi ffi 3 p , and that of S wave is 1 (dimensionless value). The slip lasts perpetually after the onset of slip at t ¼ 0.
We generate the corresponding codes using MATLAB to obtain some Green's functions L ij=k . Compared with the figures copied from Tada's paper (Figs. 2b,3b and 4b), it is

Comparison of run times
The fact which we emphasize several times that our codes have a higher speed results from a different method to solve the integral E(c). We transform the integral into a elementary algebraic problem to solve a system of inequalities, and Tada considered the problem from geometric perspective. To get a more effective program, we express the result with a series of set operations of closed intervals. We manage to circumvent complicated nested control flows that will exist according to Tada's method, which will be served as the main reason to acquire a higher speed. The same parameter values in the previous correctness validation will be chosen to calculate the integral E(c). We choose a series of time durations of rupture propagation starting at t ¼ 0 and record the corresponding run times in order to compare the run times of our method and Tada's (Fig. 5). Observing the figure, the obvious advantage of our codes is attested.
We choose L 31=1 which contain E(c) to demonstrate the percentages of time for the calculation of E(c) in the whole calculation with different time durations (Fig. 6). Using our code, the percentage is about 15 times for the calculation of the whole expressions of L ij=k by our method is obviously shorter than that of E(c) using Tada's formulae. The optimization for the calculation of E(c) is valuable. The improvement of the efficiency of the algorithm is conducive to the simulation of rupture propagation which we will concentrate on in the future.
Finally, it is essential to point out that singularity will be generated when the receiver point locates at the vertices of triangular element. In the actual simulations of the rupture propagation, we choose the inner center of each triangular element as the receiver points to circumvent the problem. Throughout this section, we use MATLAB 2015a on the ThinkPad T450 which contains 4 GB internal storage to complete all the numerical examples.

Conclusions
It is presented in this article that rigorous time-domain expressions for the Green's functions represent the transient stress response of an infinite, homogeneous and isotropic medium to a hypothetical constant slip rate on a triangular fault that continues perpetually after the slip onset. We obtain different forms of the expressions of the stress Green's function which can be proved equivalent to the former results by Tada and the corresponding program has a higher speed.
We believe our expressions can be used to simulate the rupture propagation in some actual earthquakes (Hok and Fukuyama 2011;Hu et al. 2017) due to the convenience in handling the complex geometry structures with a higher speed. We also intend to generalize the formulae into halfspace medium considering the effect of free surface (Zhang and Chen 2006) as the strong influence of a free surface usually occurs in some disastrous earthquakes. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.