An Approximate Method to Simulate Post-Seismic Deformations in a Realistic Earth Model

The geodetic observations of static deformations, including gravity perturbations and displacement ﬁelds due to huge earthquakes, are understood and explained using recent dislocation theories. Due to multiple possible mechanisms for the post-seismic phase of earthquakes, the dominant mechanism may change at different spatiotemporal ranges for different earthquake types. Accurate forward and inverse modeling of post-seismic deformations is valuable and needed information for geoscience communities. The existing methods for calculating gravitational viscoelastic relaxation can be improved or simpliﬁed to make them more suitable for more realistic Earth models and/or to overcome the poor convergence performance and/or overﬂow risks during numerical calculations. In this study, a simple and effective method for calculating the post-seismic relaxation deformations is proposed. This method is different from previous methods, such as the normal mode summation and rectangle integration methods. The proposed method consists of a rational functional approximation of the integral kernel and a transformation of the numerical inverse Laplace transform (NILT) into an alternating series summation using the residual theorem. Then the intrinsic oscillation and overﬂow risks are thoroughly suppressed. The accuracy of the calculated Green’s functions can be easily controlled by choosing a suitable parameter. In addition, the proposed method also has applicability in different Earth models with linear rheological proﬁles.


Introduction
An entire earthquake process may contain three major phases: the inter-seismic accumulation, the co-seismic transient deformation, and the long-term post-seismic adjustment. The viscoelastic relaxation deformation, as a major mechanism of the post-seismic phase, can be captured using modern geodetic technologies for megathrust earthquakes (e.g., Wang et al. 2012;Sun et al. 2014;Li et al. 2018;Qiu et al. 2019;Agata et al. 2019). To explain the multiple observations and investigate the physical mechanisms, accurate modeling for such processes is a necessity and attractive to geophysicists and geodesists.
Some dislocation theories based on a flat viscoelastic Earth model have been well developed in recent decades. For instance, post-seismic deformations due to a point source and a finite fault have been studied by Rundle (1982) using a two-layer viscoelastic-gravitational model. The gravity effect, however, was not considered correctly in this and other previous papers. Wang (2005) presented a consistent method to include the gravity effect in viscoelastic Earth models. His group contributed an open source Fortran code PSGRN/PSCMP for calculating co-and post-seismic deformations on a layered half-space Earth model (Wang et al. 2006). \looseness1fgMany scientists have also proposed dislocation theories in spherical Earth models. For example, the normal mode summation method was applied for calculating the global deformations in a stratified, self-gravitating, and incompressible viscoelastic Earth model by Sabadini et al. (1984), Pollitz (1992), Piersanti et al. (1995), and others. However, the continuous structure and compressibility of the Earth cannot be considered simultaneously due to difficulties in finding the innumerable mode. The first improvement was proposed by Tanaka et al. (2006Tanaka et al. ( , 2007 where a closed path integration method in a complex plane was used to bypass this fault in the normal mode method. Once a suitable integration path was determined, the deformation Love numbers (DLNs) (or Green's functions) could be obtained from the complex ones using numerical integration along a rectangle path. An alternative method was also proposed by Spada and Boschi (2006) that modified the normal mode formulation and used the Post-Widder Laplace inversion formula (Gaver 1966). This method is beneficial because it bypasses the numerical Laplace integral and root-finding procedure at the same time.
These two major methods have been successful in simulating post-seismic deformations; however, their distinct features and applicability scope needs to be discussed. The accuracy and reliability of the closed path integration depend on a rectangle path and a sampling scheme along the path. As pointed out by Tanaka et al. (2006), two trade-off conditions need to be considered: a smooth variation of the integrand and low oscillation of the integral kernel. These considerations result in the utilization of a root-finding scheme to locate the positions of the smallest and largest poles for a given viscoelastic Earth model. This may limit its application in inverse problems, such as estimating the viscosity of the asthenosphere. The modified normal mode method proposed by Spada and Boschi (2006) has suffered due to its slow convergence performance of the series summation formula. In addition, it requires a multiple-precision computer system. Therefore, alternative methods to calculate post-seismic deformations in a realistic Earth model are still required. Tang and Sun (2019) confirmed that the approximate inverse Laplace transform method presented by Valsa and Brančik (1998) can be applied to simulate post-seismic deformations. However, the mathematical process of this method is complicated and not easy to follow. In this paper, a straightforward approach is proposed as a simplification of the previous study. The integral kernel is directly approximated as a rational function, and the residual theorem is then applied to obtain a similar formula. This new approach is concise and clear, focusing on the nature of the problem and avoiding lengthy non-essential processes. Moreover, the modified method might be considered as common one that can be extended to deal with other similar problems. That is, the kernel of an integral transform (such as the Laplace and Mellin transforms) can be replaced by approximate functions that result in a series of new algorithms.
In the following sections, this simple and effective algorithm for calculating the post-seismic deformations is presented in Sect. 2. A numerical verification of the method on the Love numbers is given in Sect. 3. Finally, the discussions and conclusions are presented.

Complex DLNs in the Laplace Domain
The focus of this research is on the time-dependent postseismic deformations in a spherically symmetric, viscoelastic, and self-gravitating Earth. There are three equations that govern the geodynamic process: the equation of the equilibrium, the stress-strain relationship of the viscoelastic martial, and Poisson's equation (Takeuchi and Saito 1972). Because the stress-strain relationship is time-dependent, the correspondence principle (McConnell 1965) is utilized to transform the viscoelastic problem into an equivalent elastic problem via the Laplace transform. Then all of the equations have the same form as they are in the elastic case, and the only difference is that the Lamé parameters are functions of the Laplace variable, s, for viscoelastic material. After a vector spherical harmonic expansion of the transformed equations, the equations system is uncoupled into two systems: the spheroidal portion involving a 6-D displacementstress vector, and the toroidal part involving a 2-D one. The boundary conditions of these two systems at the center of the Earth, the source location, the free surface, and the internal discontinuities need to be satisfied. Then, the Runge-Kutta scheme is used to solve the transformed deformation problem to obtain DLNs, similar to what Sun and Okubo (1993) did in treating the co-seismic deformation in the pure elastic Earth model. To transform the DLNs at a radius, r, in the Laplace domain into the time domain is quite difficult compared with other processes. In addition, the full analytical method, such as used by Tang and Sun (2018a), Tang and Sun (2018b), is restricted in this more realistic Earth model. A numerical method might be the only possible choice. Following the definition of the inverse Laplace transform, the timedependent DLNs can be written as: is the DLNs in the time domain; i is the imaginary unit of the complex number; and c is a real number larger than all of the real parts of the singularities of the integrand. Here, 1/s appears because the focus is on the relaxation deformation of the main shock, i.e., the relaxation of a heaviside slip is considered here. The complex DLNs change smoothly along the vertical line Re(s) D c; however, the integral kernel, e st , is an oscillating function along this path. More explicitly, the integral kernel can be written as e st D e ct [cos(Im(s)t) C i sin (Im(s)t)], and it reveals the intrinsic oscillational feature with respect to Im(s)t, which may trouble the inverse Laplace transform.
This problem is avoided using an approximation of the integral kernel (Valsa and Brančik 1998): e st e st 1 e 2.a st/ ; a > Re.s/t: (2) Then the inverse Laplace transform of Y h (r, s) is approximately written as: To calculate the integral of G(r, s), another curve integral is added along a semicircle centered at (c, 0) with an infinite radius to obtain a closed path integral (Fig. 1b): It has been shown that the complex DLNs approach a real constant as the Laplace variable jsj approaches infinity [see Eq. (23) of Cambiotti et al. (2009)]. Thus, it is clear that Y h (r, s) ! 0 as s ! 1. This fact results in G(r, s) ! 0 as s ! 1. Hence, the integration of the second portion of Eq. (4) is zero, and it contributes nothing to the time-dependent Love numbers.
For this closed path integral, the residual theorem can be directly applied. For simplicity, we define the variable P(s) D Y h (r, s)e st and Q(s) D 1 e 2(a st) , then Eq. (4) is rewritten as: where all of the zero points of Q(s) are q k D (a i k)/t and Q ' (s) denotes the derivative of it. After some simple algebraic operations, Eq. (5) can be written as: a=t > c: Here, the identical equation Y h D Y h .s / is used during simplification, and denotes a conjugate operator. Note that the above formula will result in a real value, although it involves a complex series.
The real DLNs can be found by approximately expressing them as a weighted sum of the samples of the complex counterpart at some special points along a vertical path in the complex plane. To speed up the convergence of the series and maintain high accuracy, the first portion of the series with 0 Ä k < k c should be summed directly, and the Euler's transforms can be applied for the other portions with k > k c . Note that the above formula only holds for post-seismic deformations when t > 0. The co-seismic DLNs at t D 0 can be obtained from the complex DLNs by setting a large enough jsj.

The Oscillational Integral During NILT
The oscillational inverse Laplace transform in Eq. (1) is approximately rewritten as an alternating series of DLNs sampled on a special vertical path without any oscillational terms. Hence, the overflow risk due to the integral kernel is completely avoided, and an accurate calculation of the post-seismic deformations can be achieved using the current method. The new method can be applied to calculate postseismic relaxation deformations due to four point sources (labeled as 12, 32, 22, and 33) as defined by Sun and Okubo (1993) in the realistic spherical model with a continuous stratification.
To show this point in a straightforward manner, a complex DLN of k 32 (r D R, s) was used as an example on the surface of the Earth (set r D R, the mean radius) at time t D 5 years with degree n D 20 due to a point dip-slip source with depth of 20 km calculated using this proposed method (the parameter was set to a D 6) and a common numerical integration. The elastic parameters of the PREM (Dziewonski and Anderson 1981) were used. The inner core within a radius r < 1,221.5 km and the crust with a thickness of 50 km were elastic. In addition, the viscosity of the mantle (r > 5,871 km) was set at 10 18 Pa s, and that of the other portion was set at 10 19 Pa s. In Fig. 2, the k 32 (s) approaches the elastic DLN and k 32 (s)/s approaches zero when the norm of the Im(s) is large. The integrand e st k 32 (s)/s then oscillated quickly, which was the trouble source of the inverse Laplace transform. In the following subsection, the solution of this problem will be demonstrated using our proposed method.

Convergence of the Proposed Method
The time-dependent DLN k 32 (t) was computed using the proposed method and is plotted in Fig. 3. A simple rectangular integration along the path Re(s) D a/t D 6/t D 1.2 year 1 shows a very poor convergence performance (Fig. 3a). The final series is truncated at various values of the maximum number of k. The summation of the first kth terms are denoted as S k (Fig. 3b). It was found that S k exhibited regular oscillations around the convergence value ( 1.73E 5, calculated using Euler's accelerated sum in Fig. 3c). In Fig. 3c, the first 21 terms were directly summed and the following 19 terms were accelerated using Euler's transform. It is clear that a satisfactory convergence performance was achieved using approximately 40 terms of the complex DLNs if a suitable convergence acceleration of the series was utilized. The direct integration of the inverse Laplace transform was too rough to be available for post-seismic simulations. The integral kernel along a vertical line was proportional to sin(Im(s)t), and this means the kernel oscillated with the period of T 0 D 2 /t along the integral path. If we assume that the maximum sample step length is T 0 /20 and truncate 10 -4 10 -2 10 0 10 2 10 4 Im(s)  Here, the normalization factor of the DLNs is the same as in Sun and Okubo (1993). The solid lines denote those calculated using the previous method (Tanaka et al. 2006), and the cycles denote the results calculated using the proposed method the integral at a length of 100 periods, then it requires at least 2,000 samples. In fact, a numerical integration cannot achieve a stable value after approximately 2,000 samplings along the vertical line in Fig. 3a. However, only 40 terms resulted in an acceptable accuracy of the result using this proposed method and Euler's transform to accelerate the convergence speed (Fig. 3c).

Verifying Using a Comparison with a Previous Method
The time-dependent DLNs were then calculated using the proposed method to verify them. Only k 32 was used as an example, and the result is shown here (other variables are also verified but are not shown here). Most of the parameters of the Earth model were used, as discussed in the above Sects. 3.1 and 3.2. In addition, the crust was set to 30 km, and the inner core was kept elastic. The degrees of 5 and 20 DLNs of the dip-slip source with depths of 1 km were calculated using Tanaka's method (Tanaka et al. 2006) and our proposed method [Eq. (6)] in the time range from 0 to 50 years. The results are plotted in Fig. 4. The results of the two schemes agree well, and this confirms the correctness of this method.

Discussion and Conclusions
In this study, an approximate approach was presented to calculate the post-seismic deformations due to point sources in a realistic viscoelastic spherical Earth model with the continuous stratification, compressibility, and self-gravitation was proposed. The good agreement of the time-dependent dislocation Love numbers calculated using a previous method (Tanaka et al. 2006) and by our proposed method demonstrated that this new approximation scheme had high accuracy. The idea of the method was to approximately treat the integral kernel (Valsa and Brančik 1998) and to implement the residue theorem in a semicircle. The location of this semicircle, i.e., the parameter, can be adjusted according to a specific problem. Then, the numerical oscillation and overflow risk of some previous methods are avoided by using a skillful treatment of the integral kernel. Essentially, the approximate approach is different in comparison with other published methods. The proposed method is as simple and clear as the Post-Widder method. They have very similar forms of the final expression, i.e., the inversed function is expressed using a weighted sum of the special sampled complex functions. Both of them can be applied for different viscoelastic Earth models because all of the singular points are negative (Plag and Jüttner 1995;Vermeersen and Mitrovica 2000) for a stable Earth model without a heavier layer overlying a lighter one. However, they are different in essence. The Post-Widder method is a differential approach along the positive real axis; while, the current method is an integral approach along the vertical line in complex plane. There are 40 sampling points of the complex function in our method and 16 in the Post-Widder method (Melini et al. 2008) for calculating postseismic deformations. The computational efficiency of the two methods were nearly the same. However, a commercial Fortran compiler supporting 30 significant digits with the IEEE extended-precision format cannot ensure success in the application of the Post-Widder method (Melini et al. 2008). On the contrary, the current method can be easily implemented using any scientific computing language without an extended precision algebraic library. The less terms in the sampling points of the Post-Widder method primarily benefits from the extended library. In addition, an extended floating-number requirement means large computational memory and performance degradation. If considering computational efficiency and memory consumption together, these two methods have similar performances.
Compared with the rectangle integration method or the normal mode method, the current method focused on how to treat the oscillational integral kernel rather than the poles of a specific problem. The rectangle integration method (Tanaka et al. 2006(Tanaka et al. , 2007 treats the inverse Laplace transform bypassing the innumerable poles. However, there are at least three parameters that should be determined in the numerical computation for a given earth model. In contrast, this proposed method only has one adjustable parameter. In fact, this parameter can be fixed at approximately six because the largest pole of the post-seismic deformation will not be larger than zero in a stable Earth model. In this point of view, the current method is suitable for an inverse problem. In addition, the current method can also be used for other linear rheological models, such as the Kelvin-Voigt material, the Maxwell solid, Burgers body, and their combination.