Effect of Earth Models on Coulomb Stress Change Caused by Surface Load

Abstract

Based on the elastic spherical Earth models, this study presents the radial Love numbers and Green’s functions of the surface load problem, then evaluates the effect of different Earth models on the Coulomb stress changes due to elastic deformation by surface load. First of all, the radial load Love numbers are introduced and verified by comparing them with the results of an analytical method. Second, the formulae of the radial Green’s functions of the internal displacements and strain tensors are obtained by summing the Love numbers. The results are compared with Green’s function of half-space to verify the correctness. Then, the stress tensors are obtained by the constitutive relation. Finally, the Coulomb stress changes due to elastic deformation around the Wenchuan Mw7.9 earthquake by the impoundment load of the Zipingpu reservoir are calculated based on five Earth models. The results show that the elastic deformations caused by the impoundment load do not increase the seismicity rate in the study area. At the epicenter of the great event, the maximum effect due to the Earth models is about 38%. The stratification of the Earth is the main effect factor on the Coulomb stress change, and it should not be ignored in actual studies. The biggest difference of Coulomb stress changes based on the different stratified Earth models is close to 15%, which also cannot be ignored.

Appendices

Appendix A: Load Love Numbers for a Homogeneous and Incompressible Earth

The formulae with the unknown coefficients ($C_{{1,{\text{homo}}}}$, $C_{{2,{\text{homo}}}}$, $C_{{3,{\text{homo}}}}$) of the surface load problem for a homogeneous and incompressible Earth have already been given by Wu and Peltier (1982)

$$\begin{array}{ll} {C_{{1,{\text{homo}}}} = } & {\frac{{2\rho _{0} g_{0} \left( {n + 1} \right)\left( {n - 1} \right)\left( {2n + 3} \right)a^{{ - \left( {n + 1} \right)}} }}{{3m_{e} \beta }}} \\ {C_{{2,{\text{homo}}}} = } & { - \frac{{\rho _{0} g_{0} n^{2} \left( {n + 2} \right)a^{{ - \left( {n - 1} \right)}} }}{{3m_{e} \beta }}} \\ {C_{{3,{\text{homo}}}} = } & { - \frac{{g_{0} \mu _{0} \left( {2n^{2} + 4n + 3} \right)a^{{ - \left( {n + 1} \right)}} }}{{m_{e} \beta }},} \\ \end{array}$$

(A1)

where $m_{e}$ is the mass of the Earth, $g_{0}$ is the gravitational acceleration of the surface of the homogeneous Earth while $\rho_{0}$ and $\mu_{0}$ are the homogeneous density and shear modulus. The term $\beta$ is defined as

$$\beta = a^{ - 1} \left[ {n\rho_{0} g_{0} + \mu_{0} a^{ - 1} \left( {2n^{2} + 4n + 3} \right)} \right].$$

Herein, $a$ is the radius of the Earth, and $n$ is the degree of Legendre function.

The coefficients of the three linear independent solutions for the homogeneous and incompressible Earth model are list below:

$$\begin{gathered} y_{{1,{\text{homo}}}} = \left( {\frac{{nr^{n + 1} }}{{2\left( {2n + 3} \right)}},\frac{{\left( {n + 3} \right)r^{n + 1} }}{{2\left( {2n + 3} \right)\left( {n + 1} \right)}},\frac{{\rho_{0} \xi nr^{n + 2} + 2\mu_{0} r^{n} \left( {n^{2} - n - 3} \right)}}{{2\left( {2n + 3} \right)}},\frac{{2\mu_{0} n\left( {n + 2} \right)r^{n} }}{{2\left( {2n + 3} \right)\left( {n + 1} \right)}},0,\frac{{3\xi nr^{n + 1} }}{{2\left( {2n + 3} \right)}}} \right)^{t} \hfill \\ y_{{2,{\text{homo}}}} = \left( {r^{n - 1} ,\frac{{r^{n - 1} }}{n},\rho_{0} \xi r^{n} + 2\mu_{0} \left( {n - 1} \right)r^{n - 2} ,\frac{{2\mu_{0} \left( {n - 1} \right)}}{n}r^{n - 2} ,0,3\xi r^{n - 1} } \right)^{t} \hfill \\ y_{{3,{\text{homo}}}} = (0,0,\rho_{0} r^{n} ,0,r^{n} ,\left( {2n + 1} \right)r^{n - 1} )^{t} . \hfill \\ \end{gathered}$$

(A2)

where $\xi$ can be expressed as

$$\xi = \frac{{4\pi G\rho_{0} }}{3}.$$

$G$ is the gravitational constant.

The solutions for the interior of Earth can be obtained by the equation

$$\begin{array}{*{20}c} {Y\left( r \right) = \mathop \sum \limits_{i = 1}^{3} C_{{i,{\text{homo}}}} y_{{i,{\text{homo}}}} .} \\ \end{array}$$

(A3)

Then, $U_{n} \left( r \right)$, $V_{n} \left( r \right)$ and $\phi_{n} \left( r \right)$ can be solved as:

$$\begin{array}{ll} {U_{n} \left( r \right) = } & {\frac{{\rho _{0} g_{0} \left[ {n\left( {n + 1} \right)\left( {n - 1} \right)r^{{n + 1}} - n^{2} \left( {n + 2} \right)r^{{n - 1}} a^{2} } \right]}}{{3m_{e} \beta a^{{n + 1}} }}} \\ {V_{n} \left( r \right) = } & {\frac{{\rho _{0} g_{0} \left[ {\left( {n - 1} \right)\left( {n + 3} \right)r^{{n + 1}} - n\left( {n + 2} \right)r^{{n - 1}} a^{2} } \right]}}{{3m_{e} \beta a^{{n + 1}} }}} \\ {\Phi _{n} \left( r \right) = } & { - \frac{{g_{0} \mu _{0} \left( {2n^{2} + 4n + 3} \right)r^{n} }}{{m_{e} \beta a^{{n + 1}} }}.} \\ \end{array}$$

(A4)

The relational equation between $\left( {U_{n} \left( r \right),V_{n} \left( r \right),\Phi_{n} \left( r \right)} \right)$ and the load Love numbers can be combined as:

$$\begin{array}{*{20}c} {\left[ {\begin{array}{*{20}c} {h_{n} \left( r \right)} \\ {l_{n} \left( r \right)} \\ {k_{n} + (r/a)^{n} } \\ \end{array} } \right] = \frac{{m_{e} }}{{ag_{0} }}\left[ {\begin{array}{*{20}c} {U_{n} \left( r \right)g_{0} } \\ {V_{n} \left( r \right)g_{0} } \\ { - \Phi_{n} \left( r \right)} \\ \end{array} } \right],} \\ \end{array}$$

(A5)

The radial load Love numbers of the homogeneous and incompressible Earth are as below:

$$\begin{array}{ll} {\begin{array}{ll} {h_{n} = } & {\frac{{\rho _{0} g_{0} \left[ {n\left( {n + 1} \right)\left( {n - 1} \right)r - n^{2} \left( {n + 2} \right)r^{{ - 1}} a^{2} } \right]}}{{\left( {3\beta a^{2} } \right)\left( {\frac{r}{a}} \right)^{n} }}} \\ {l_{n} = } & {\frac{{\rho _{0} g_{0} \left[ {n\left( {n - 1} \right)\left( {n + 3} \right)r - n\left( {n + 2} \right)r^{{ - 1}} a^{2} } \right]}}{{\left( {3\beta a^{2} } \right)\left( {\frac{r}{a}} \right)^{n} }}} \\ {k_{n} = } & {\left[ { - 1 + \frac{{\mu \left( {2n^{2} + 4n + 3} \right)}}{{\beta a^{2} }}} \right]\left( {\frac{r}{a}} \right)^{n} } \\ \end{array} } \\ \end{array}$$

(A6)

If $r$ equals $a$, then the Eq. (A6) becomes the equations for surface Love numbers.

Appendix B: Comparison with the Previous Studies

See Figs. 8, 9.

Fig. 8

Comparison of the surface load Love numbers between Wang et al. (1996) and this study based on the PREM

Fig. 9

Comparison of the surface Green’s functions between Wang et al. (2012) and this study for radial (a), tangential (b), and strain element (c)

Appendix C: Formulae of Legendre Sums

$$\begin{gathered} \mathop \sum \limits_{n = 0}^{\infty } P_{n} \left( {\cos \theta } \right) = \frac{1}{{2\sin \left( {\theta /2} \right)}} \hfill \\ \mathop \sum \limits_{n = 0}^{\infty } nP_{n} \left( {\cos \theta } \right) = - \frac{1}{{4\sin \left( {\theta /2} \right)}} \hfill \\ \mathop \sum \limits_{n = 0}^{\infty } n^{2} P_{n} \left( {\cos \theta } \right) = - \frac{{\cos^{2} \left( {\theta /2} \right)}}{{8\sin^{3} \left( {\theta /2} \right)}} \hfill \\ \mathop \sum \limits_{n = 0}^{\infty } \frac{{\partial P_{n} \left( {\cos \theta } \right)}}{\partial \theta } = - \frac{{\cos \left( {\theta /2} \right)}}{{4\sin^{2} \left( {\theta /2} \right)}} \hfill \\ \mathop \sum \limits_{n = 0}^{\infty } \frac{1}{n}\frac{{\partial P_{n} \left( {\cos \theta } \right)}}{\partial \theta } = - \frac{{\cos \left( {\theta /2} \right)\left[ {1 + 2\sin \left( {\theta /2} \right)} \right]}}{{2\sin \left( {\theta /2} \right)\left[ {1 + \sin \left( {\theta /2} \right)} \right]}} \hfill \\ \mathop \sum \limits_{n = 0}^{\infty } \frac{1}{n}\frac{{\partial^{2} P_{n} \left( {\cos \theta } \right)}}{{\partial \theta^{2} }} = \frac{{1 + \sin \left( {\theta /2} \right) + \sin^{2} \left( {\theta /2} \right)}}{{4\sin^{2} \left( {\theta /2} \right)\left[ {1 + \sin \left( {\theta /2} \right)} \right]}} \hfill \\ \end{gathered}$$

(C1)

Appendix D: Half-Space Green’s Function

The surface load study of half-space is to calculate the response of a non-gravitating, elastic substance by surface pressure. Since the Green’s function of this study is obtained from Boussinesq (1885), the deformation of the half-space due to the surface load is called the Boussinesq problem. The formulae for the displacement of the Green’s functions introduced by Farrell (1972) was used while the relationship between the displacements and the strain tensor is considered to derive the Green’s functions of the strain tensor.

The displacement the Green’s function for a point load can be expressed as

$$\begin{array}{ll} {u\left( {z,x} \right) = } & { - \frac{{g_{0} }}{{4\pi \mu R}}\left( {\frac{\sigma }{\eta } + \frac{{z^{2} }}{{R^{2} }}} \right)} \\ {v\left( {z,x} \right) = } & { - \frac{{g_{0} }}{{4\pi \eta x}}\left( {1 + \frac{z}{R} + \frac{{\eta x^{2} z}}{{\mu R^{3} }}} \right)} \\ {\sigma = } & {\lambda + 2\mu } \\ {\eta = } & {\lambda + \mu } \\ {R = } & {x^{2} + z^{2} ,} \\ \end{array}$$

(D1)

where $z$ is the depth of the calculation point and $x$ is the horizontal distance from the load to the calculation point while $u$ and $v$ denote the vertical and tangential displacements. Since the formulae in Farrell (1972) consider surface pressure, the surface gravitational acceleration ($g_{0} )$ of the Earth model is included in the formulae to unify the definition of Green’s functions between the Boussinesq problem and the spherical one described in the body text.

Because the Boussinesq problem is axially symmetrical about the load source, the spatial symmetry strain formulae for cylinder coordinate are as below:

$$\begin{array}{ll} { \in_{xx} = } & {\frac{\partial v}{{\partial x}}} \\ { \in_{zz} = } & {\frac{\partial u}{{\partial z}}} \\ { \in_{yy} = } & \frac{v}{x} \\ { \in_{xz} = } & {\frac{\partial v}{{\partial z}} + \frac{\partial u}{{\partial x}},} \\ \end{array}$$

(D2)

where $y$ denotes the tangential direction perpendicular to $x$.

According to the displacement formulae in Eq. D1, the formulae of the partial derivative can be obtained as:

$$\begin{array}{ll} {\frac{{\partial u}}{{\partial z}} = } & {\frac{{gz}}{{4\pi \mu R^{3} }}\left( {\frac{\sigma }{\eta } + \frac{{z^{2} - 2x^{2} }}{{R^{2} }}} \right)} \\ {\frac{{\partial v}}{{\partial x}} = } & {\frac{g}{{4\pi \eta x^{2} }}\left( {1 + \frac{{z\left( {R^{2} + x^{2} } \right)}}{{R^{3} }} + \frac{{\eta x^{2} z\left( {3x^{2} - R^{2} } \right)}}{{\mu R^{5} }}} \right)} \\ {\frac{{\partial u}}{{\partial x}} = } & {\frac{{gr}}{{4\pi \mu R^{3} }}\left( {\frac{\sigma }{\eta } + \frac{{3z^{2} }}{{R^{2} }}} \right)} \\ {\frac{{\partial v}}{{\partial z}} = } & {\frac{{gr}}{{4\pi \mu R^{3} }}\left( { - \frac{\sigma }{\eta } + \frac{{3z^{2} }}{{R^{2} }}} \right).} \\ \end{array}$$

(D3)

Substituting Eq. D3 into Eq. D2, the formulae of the Green’s functions of the strain tensor is acquired as:

$$\begin{array}{ll} {\begin{array}{ll} { \in _{{xx}} \left( {z,x} \right) = } & {\frac{g}{{4\pi \eta x^{2} }}\left( {1 + \frac{{z\left( {R^{2} + x^{2} } \right)}}{{R^{3} }} + \frac{{\eta x^{2} z\left( {3x^{2} - R^{2} } \right)}}{{\mu R^{5} }}} \right)} \\ { \in _{{zz}} \left( {z,x} \right) = } & {\frac{{gz}}{{4\pi \mu R^{3} }}\left( {\frac{\sigma }{\eta } + \frac{{z^{2} - 2x^{2} }}{{R^{2} }}} \right)} \\ { \in _{{yy}} \left( {z,x} \right) = } & { - \frac{g}{{4\pi \eta x^{2} }}\left( {1 + \frac{z}{R} + \frac{{\eta x^{2} z}}{{\mu R^{3} }}} \right)} \\ { \in _{{xz}} \left( {z,x} \right) = } & {\frac{{3gxz^{2} }}{{4\pi \mu R^{5} }}.} \\ \end{array} } \\ \end{array}$$

(D4)

So far, the formulae of Green’s functions of the displacements and strain tensor for the Boussinesq problem are obtained.

Appendix E: Constitutive equation

The relationship equations between the stress tensor ($\tau_{\theta \theta }$, $\tau_{\lambda \lambda }$, $\tau_{rr}$, $\tau_{r\theta }$) and strain tensor ($\in_{\theta \theta }$, $\in_{\lambda \lambda }$, $\in_{rr}$, $\in_{r\theta }$) are as below:

$$\begin{array}{ll} {\tau _{{\theta \theta }} = } & {\lambda \left( r \right)K + 2\mu \in _{{\theta \theta }} } \\ {\tau _{{\lambda \lambda }} = } & {\lambda \left( r \right)K + 2\mu \in _{{\lambda \lambda }} \lambda \left( r \right)K + 2\mu \in _{{\lambda \lambda }} } \\ {\tau _{{rr}} = } & {\lambda \left( r \right)K + 2\mu \in _{{rr}} } \\ {\tau _{{r\theta }} = } & {2\mu \left( r \right) \in _{{r\theta }} } \\ {K = } & { \in _{{\theta \theta }} + \in _{{\lambda \lambda }} + \in _{{rr}} .} \\ \end{array}$$

(E1)

where r is the radial distance from the center of the Earth, which is equal to the difference between the radius of the Earth $a$ and the point of interest depth $d$. $\lambda$ and $\mu$ are the Lamé parameters.