Elastic deformation resulting from internal mass loading: semi-analytical solution

SUMMARY

1 INTRODUCTION

Various factors, such as ocean tides, fluctuations in atmospheric pressure, ice and snow accumulation in high mountains or polar regions, variations in lake water levels and heavy rainfall, can perturb the Earth's internal stress field, which in turn causes crustal deformation (Farrell 1979; Tregoning & van Dam 2005; Cavalié et al. 2007; Dietrich et al. 2010; Bos et al. 2015; Penna et al. 2015; Heki & Arief 2022).

Alongside these phenomena, which can be approximated by a surface loading problem, there are also natural processes and human activities that can be partially treated as internal mass loading problems. For instance, the volcano and seamount, including the surface mass above the average topographic surface and the internal mass anomaly originating from during penetration of the Earth's crust from deep, generate loading effects in the lithosphere (Harmon et al. 2006; Klügel et al. 2015). In addition, during the eruption process of some active volcanoes, a large amount of mass of magma migrates from the deep interior to the shallow region, which may cause observable surface displacement and gravity field changes (Furuya et al. 2003; Schmincke 2004; Lu & Dzurisin 2010; Carbone & Poland 2012; Bagnardi et al. 2014). Also, the process of obtaining industrial raw materials such as ores, coal seams and natural gas from the crust and pumping groundwater from the basin also involves the movement of the internal mass in the crust (Trešl 1988; Bell et al. 2000; Phien-wej et al. 2006; Litvinenko 2018). Backfilling in mining will significantly reduce the risk of collapse (Guo et al. 2011; Gupta & Paul 2015) and then the ground deformation may be partly attributable to the change of internal density. Although the mechanical mechanisms underlying geophysical processes are complex and involve numerous factors, such as material migration, pore pressure and gravitational collapse, investigating the deformation caused by mass movement is meaningful in certain geophysical contexts.

To quantitatively estimate the influence on gravity field and surface deformation, some theoretical modelling, to a certain extent, is necessary and beneficial. For instance, the classic loading theories assume that the mass distribution is located outside the Earth's crust. This theory model, summarized by Bos & Scherneck (2013) and mainly originally developed by Longman (1962), Longman (1963) and Farrell (1972), provides a widely used tool for modelling the theoretical deformation of the layered Earth under the action of mass transports in global geophysical fluids (Chao et al. 2000; Lu et al. 2018). However, as explained above, mass redistribution does not occur entirely on the surface of the Earth, but at a certain depth range. Investigation of internal mass loading deformation is necessary for this issue.

Then, the following problems naturally arise: how to get the deformation of Green's function generated by internal mass load? What is the difference between surface and internal mass loading problems from the perspective of the displacement field and gravity change? Is it necessary to consider the inner mass loading effect in the general loading problem? In this study, we will put the loading mass in the interior instead of on the surface and partially extend the loading theory in Section 2 to investigate the difference in theoretical deformation between them and try to answer these questions in other sections. An example is also given for the disc-like mass loading to show the practicability of the theoretical results. Finally, there are some discussions and conclusions on the applicable scope of internal load deformation theory.

2 METHODOLOGY

2.1 Boundary value problem for internal mass loading

Early theoretical geophysicists gave the governing equations for the quasi-static deformation of stratified Earth deformation (Love 1911; Backus 1967) in the spherical coordinate system |$( {r,\theta ,\phi } )$|⁠. In this study, we consider an internal point mass |${m}_0$| located at |$( {{r}_s,0,0} )$| below the North Pole. This additional mass introduces a body force per volume |${\boldsymbol{f}} = - {g}_0( r ){\rho }_L{{\boldsymbol{e}}}_r$| assuming bulk density |${\rho }_L = {m}_0\delta ( {{\boldsymbol{r}} - {{\boldsymbol{r}}}_s} )$| due to the self-gravitation of the Earth. Here, the Dirac's delta function is |$\delta ( {{\boldsymbol{r}} - {{\boldsymbol{r}}}_s} ) = \frac{{\delta ( {r - {r}_s} )}}{{{r}^2}}\frac{{\delta ( \theta )}}{{\sin \theta }}\delta ( \phi )$|⁠, where |${{\boldsymbol{e}}}_r$| is the unit outward vector, and |$- {g}_0{{\boldsymbol{e}}}_r$| is the gravitational acceleration in the absence of motion (note here |${g}_0 > 0$|⁠). The loading bulk density can be expansion in a Legendre series as |${\rho }_L = \mathop \sum \limits_n \frac{{2n + 1}}{{4\pi }}\frac{{\delta ( {r - {r}_s} )}}{{{r}^2}}{m}_0{P}_n( {\cos \theta } )$|⁠. Hereafter, the Legendre function may be abbreviated as |${P}_n$|⁠. Then, the body force is expanded as |${\boldsymbol{f}} = \delta ( {r - {r}_s} ){{\boldsymbol{e}}}_r\mathop \sum \limits_n {f}_2{P}_n$| with |${f}_2 = - \frac{{2n + 1}}{{4\pi {r}^2}}{m}_0{g}_0( r )$|⁠.

The elastic deformation of a regular layered earth model under this internal loading force |${\boldsymbol{f}}$| is described by three equations, namely the stress equilibrium equation, Poisson's equation and the constitutive relation (Backus 1967; Takeuchi & Saito 1972; Farrell 1972; Wang 1997; Bos & Scherneck 2013):

Here we set the geopotential as a positive value, and the force is its gradient; |${\rho }_0$| is the initial density before any motion and |${\boldsymbol{I}}$| is the unit tensor; |${\boldsymbol{u}}$| is the displacement; |${\boldsymbol{\tau }}$| is the stress tensor; |$\lambda$| and |$\mu$| are Lamé constants; |$\psi$| is the total potential perturbation including the gravitational effect of the mass source and the perturbation increment of the initial potential field. Note the second term in Poisson's equation is added here, responsible for the loading mass source. If we denote the initial potential of positive value as |${\psi }_0$| then the corresponding gravitational acceleration vector is |${{\boldsymbol{g}}}_0 = \nabla {\psi }_0 = \frac{{{\rm d}{\psi }_0}}{{{\rm d}r}}{{\boldsymbol{e}}}_r = - {g}_0{{\boldsymbol{e}}}_r$| and the Poisson's equation before motion yields |${\nabla }^2{\psi }_0 = \frac{{{{\rm d}}^2}}{{{\rm d}{r}^2}}{\psi }_0 + \frac{2}{r}\frac{{\rm d}}{{{\rm d}r}}{\psi }_0 = - 4\pi G{\rho }_0( r )$| when |$r \ne {r}_s$|⁠.

Using the y-notation used by Takeuchi & Saito (1972) and developed by Alterman et al. (1959), the three equations in SNREI (spherical, non-rotating, elastic and isotropic) earth model can be rewritten as an ordinary differential equation system (Pekeris & Jarosch 1958) only considering the spheroidal competent [the torsional part is not needed here, see Farrell (1972)] in the spherical coordinate system |$( {r,\theta ,\phi } )$| after implementing vector spherical harmonic expansion. The displacement is expressed by the coefficients |${y}_1$| and |${y}_3$| as |${\boldsymbol{u}} = \mathop \sum \limits_n {y}_1( r ){P}_n{{\boldsymbol{e}}}_r + {y}_3( r )\frac{{{\rm d}{P}_n}}{{{\rm d}\theta }}{{\boldsymbol{e}}}_\theta$|⁠; the traction stress is expressed by |${y}_2$| and |${y}_4$| as |${\boldsymbol{\ }}{{\boldsymbol{e}}}_r \cdot {\boldsymbol{\tau }} = \mathop \sum \limits_n {y}_2( r ){P}_n{{\boldsymbol{e}}}_r + {y}_4( r )\frac{{{\rm d}{P}_n}}{{{\rm d}\theta }}{{\boldsymbol{e}}}_\theta$|⁠. The total potential perturbation is expressed as |$\psi = \mathop \sum \limits_n {y}_5( r ){P}_n$|⁠. We use the definition of |${y}_6( r ) = \frac{{{\rm d}{y}_5( r )}}{{{\rm d}r}} - 4\pi G{\rho }_0{y}_1( r ) + \frac{{n + 1}}{r}{y}_5( r )$| following Takeuchi & Saito (1972). The total potential perturbation |${\psi }_i$| inside the Earth for |$r \le a$| is expressed by |${y}_5$|⁠, and |${y}_6$| is called gravitational flux density defined by Takeuchi & Saito (1972). There are differences in the definitions of |${y}_5$| and |${y}_6$| in numerous literatures, which may lead to significant difference in the expression of the differential equations and boundary conditions. The explicit form of matrix |${\boldsymbol{A}}$| used here in Appendix A is identical to that in the previous papers (Okubo & Endo 1986; Sun 1992). To illustrate the total potential perturbation |$\psi$| in the entire space in Section 2, the total potential perturbation outside the Earth expressed by |${\psi }_e( r )$| for |$r \ge a$| need to be introduced. Obviously, here we use piecewise functions to express the total perturbation |$\psi$|⁠.

After expansion, we can get the common ordinary differential equation system (Takeuchi & Saito 1972; Alterman et al. 1959):

where |${\boldsymbol{Y}} = {( {{y}_1,{y}_2,{y}_3,{y}_4,{y}_5,{y}_6} )}^{\rm{T}}$|⁠, coefficient matrix |${\boldsymbol{A}}$| is of order |$6{\rm{\ }} \times {\rm{\ }}6$| (see Appendix A for their values), depending on the spherical harmonic degree, density and Lamé parameters, and the factor |${\boldsymbol{F}}$| of the inhomogeneous part is

Instead of solving the inhomogeneous eq. (1b), one can solve the homogeneous counterpart while considering an added boundary condition at the source, which involves defining some possible y-notation jumps |${\boldsymbol{s}}$| with six components as |${s}_i = {{y}_i( {{r}_s} )} |_ - ^ + = {y}_i( {{r}_s + 0} ) - {y}_i( {{r}_s - 0} )$|⁠. After some derivation which are discussed in detail for seismic deformation in the section 2.2.2 Propagators and sources of the monograph (Kennett 1983), one can gets |${\boldsymbol{s}} = {\boldsymbol{F}}( {{r}_s} ) = {( {0,\frac{{2n + 1}}{{4\pi r_s^2}}{m}_0{g}_0( {{r}_s} ),0,0,0, - \frac{{2n + 1}}{{r_s^2}}G{m}_0} )}^{\rm{T}}$|⁠. As the derivation process is similar for the issues in this paper, we will not repeat it here. It typically incorporates the mass source indirectly by imposing boundary conditions on the y-variables. This is a general method of obtaining source functions in mathematics. Understanding this method requires some knowledge of ordinary differential equation solutions. In this study, we will use an alternative approach to achieve the same discontinuities at the mass source as expressed in eqs. (2) and (7).

Different quasi-static deformation problems of the earth usually have individual boundary conditions, but they share the identical homogeneous equation |${\rm d}{\boldsymbol{Y}}/{\rm d}r = {\rm{\ }}{\boldsymbol{AY}}$| because it expresses the static equilibrium of a self-gravitational layered sphere due to a point force when using the same notations. One can introduce the so-called source function |${\boldsymbol{s}}$| for different static problems, as discussed by Takeuchi & Saito (1972). Then the source function, a boundary condition at the source, is the key to distinguishing a deformation problem and needs special attention. Mathematically, this is an alternative expression of the actual force. In this investigation, our primary focus lies on the internal mass loading problem, which diverges from the well-explored point dislocation problem (e.g. Sabadini et al. 1984; Pollitz 1992; Sun 1992) and surface mass loading problem (e.g. Farrell 1972) in previous studies. The origin of Earth deformation examined in this research is fundamentally distinct from that of the internal point dislocation source. The point dislocation source arises from usual double-couple forces, whereas the internal loading is attributed to a mass source. Compared to surface loading issues, our mass source is internal. This mass source encompasses both a direct pressure effect, resulting from the Earth's self-gravitation and a direct and indirect influence on the gravitational potential distribution (Farrell 1972). Consequently, it is imperative to derive new boundary conditions and address this unique problem. Fortunately, the methodology for resolving the simplified system of ordinary differential equations bears resemblance to previous approaches, and the accumulated experience can be partially applied to solving our problem.

Takeuchi & Saito (1972) presented source functions following from an internal unit force in direction |$( {{\nu }_x,{\nu }_y,{\nu }_z} )$| in rectangular Cartesian co-ordinates. We assume an arbitrary internal mass source |${m}_0$| at |${r}_s$|⁠. The value of |${m}_0$| can be arbitrary. For example, some authors assume that the unit gravitational potential (⁠|$G{m}_0/{r}_s = 1$|⁠) will be generated by |${m}_0$| (Saito 1974, 1978; Sun & Sjoberg 1999), and some authors assume that |${m}_0{\rm{\ }}$| is the Earth's mass (Guo et al. 2004; Zhou & Pan 2021). To facilitate a comprehensive dimensional analysis of the final formula, which describes the deformation caused by a unit point mass (1 kg), we refrain from making specific assumptions regarding |${m}_0$|⁠. Setting the direction to |$( {0,0, - 1} )$| and the magnitude to |${m}_0{g}_0( {{r}_s} )$|⁠, the source functions of |${y}_1, \cdots ,{y}_4$| are:

For source functions of involving changes in the potential |${y}_5$| and |${y}_6$|⁠, their expressions are more complicated, and details on their derivation are presented below. Both the displacement field |${\boldsymbol{u}}( {r,\theta ,\phi } )$| and the additional loading mass |${m}_0$| at r generate potential change (we temporarily set the mass location at r for a general discussion here). On any interface or surface of the earth model, the geopotential change |$\psi ( r )$| can be represented by areal density disturbance |${\rm{\Sigma }}( r )$|⁠:

where |${\rho ( r )} |_ - ^ +$| indicates the density difference between two sides of an interface before deformation, |$\gamma ( r )$| denotes the surface density expansion of point mass |${m}_0$| and |${\Gamma }_n = \frac{{2n + 1}}{{4\pi {r}^2}}{m}_0$| for the areal density. Within a sphere of radius a, the following condition should be satisfied for the potential and its derivative:

After submitting the spherical harmonics, it reduces to:

Considering the definition of |${y}_6( r )$| by Takeuchi & Saito (1972), the second in eq. (5) is rewritten as [Note |${\Gamma }_n$| is defined in eq. (3)]:

It can be simplified as follows considering the continuity of vertical displacement and total potential:

Therefore, the source functions of |${y}_5( r )$| and |${y}_6( r )$| at |${r}_s$| are:

Up to here, eqs. (2) and (7) give the source functions to the internal mass loading problem. In addition to the source function at the source, the ground surface |$r = a$| should comply with the stress-free boundary condition and the potential continuity condition. The traction vector on the surface is zero and it results the normal stress and shear stress at |$r = a$| are zeros. This surface boundary condition leads to |${y}_2( a ) = {y}_4( a ) = 0$|⁠. For the boundary conditions of the total potential perturbation, a more careful derivation is needed. Considering the continuity of the total potential perturbation at the surface and the fact that the total potential perturbation outside the Earth satisfies the Laplace equation, one gets the external total potential perturbation |${\psi }_e( r ) = \mathop \sum \limits_n {y}_5( a ){( {\frac{a}{r}} )}^{n + 1}{P}_n$| for |$r \ge a$|⁠. Like the meaning expressed in second one of eq. (4), but only considering the additional mass layer caused by the displacement, the derivative value of total potential perturbation when crosses the surface boundary has a jump, that is,

Here |$r = {a}^ +$| and |$r = {a}^ -$| denotes the outer and inner boundary, respectively and |${u}_r( a )$| means the vertical displacement. Then, the special definition of |${y}_6( r )$| by Takeuchi & Saito (1972) leads to |${y}_6( a ) = 0$| from eq. (8). Finally, the surface boundary conditions are:

The boundary conditions for a degree one mode (corresponding to the term |$n = 1$| in the spherical harmonic expansion) should be given additional attention (Okubo & Endo 1986; Blewitt 2003). In this case, the consistency relation (Farrell 1972) holds as |${y}_2( r ) + 2{y}_4( r ) + \frac{{g( r )}}{{4\pi G}}{y}_6( r ) = 0$|⁠. For this special case, the centre of mass of the Earth plus load is fixed in space, but there is no constraint on the Earth alone (Farrell 1972). As Okubo & Endo (1986) stated, the origin at the centre of mass of the undeformed Earth could be chosen. This implies that the potential change owing to solid Earth deformation is zero (Okubo & Endo 1986; Martens et al. 2019), that is |${k}_1( a ) = 0$|⁠. Subsequently, this results in a constraint of |${y}_5( a ) = \frac{{G{m}_0}}{{{r}_s}}{( {\frac{{{r}_s}}{a}} )}^2$| [see eq. (13) for a detailed definition of |${k}_1$|]. Afterward, boundary conditions

are used to determine a unique special solution of degree one.

2.2 Solutions for a layered model and loading Love numbers

The equation system of y-numbers can be solved to obtain a particular solution of |${y}_i( r )$| using the above source functions and surface boundary conditions. We assumed a point mass of unit magnitude. The resulting deformation, such as displacement and geopotential changes are called Green's functions. For a homogeneous sphere without a solid inner core and an outer liquid core, the exact analytical solution of |${y}_i$| could be obtained. For a layered spherical earth model, such as PREM, the preliminary reference earth model (Dziewonski & Anderson 1981), a standard numerical scheme such as the Runge‒Kutta integration method on ordinary differential equations could be applied for numerical solutions. We present both the analytical and the numerical solutions in this study. Approaches to solving these differential equations have been discussed in detail in numerous studies (Longman 1962, 1963; Farrell 1972; Saito 1974; Sun 1992; Bos & Scherneck 2013; Pan et al. 2015).

Here, we only briefly describe the particular solution method following Sun (1992) and Farrell (1972). A numerical integration scheme must be applied for a layered spherical earth model. We integrated upward from a homogeneous sphere with a small radius |${r}_0$||$( {0 < {r}_0 \ll a} )$|⁠, usually 100 km. Three independent solutions are in the inner core |$( {r < c} )$| that are regular at the centre. The general case is the solution for degree |$n \ge 2$| mode. These are integrated upward into the inner and outer core boundaries |$( {r = c} )$|⁠. Subsequently, these three solutions are obtained into one through the conditions at the inner and outer core boundaries (Saito 1974). From this combined solution, the initial values are obtained, that is of |${y}_5( c )$| and |${y}_7( c )$| in the outer liquid core and integrated into the core mantle boundary |$( {r = b} )$| through auxiliary variables |${y}_5( r )$| and |${y}_7( r )$| (Saito 1974; Mochizuki 1988). Afterward, three independent solutions |$( {{\boldsymbol{Y}}_1^ \uparrow ,{\boldsymbol{Y}}_2^ \uparrow ,{\boldsymbol{Y}}_3^ \uparrow } )$| could be built in the mantle based on |${y}_5( b )$| and |${y}_7( b )$| and, subsequently, integrated upward to the source |$( {r = {r}_s} )$| awaiting follow-up treatment.

In addition, an integration operation from the surface downwards is required. Taking the matrix |${\boldsymbol{D}}$|⁠, consisting of 3 column vectors of unit length and non-zero elements at (1, 1) and (1, 3) and (1, 5), as the initial value at the surface and integrating it downward to the source location results in three independent solutions |$( {{\boldsymbol{Y}}_1^ \downarrow ,{\boldsymbol{Y}}_2^ \downarrow ,{\boldsymbol{Y}}_3^ \downarrow } )$| in the mantle. We assumed that the combination coefficient of the three solutions after upward integration is |${\boldsymbol{p}}$|⁠, and the combination coefficient of the three solutions after downward integration is |${\boldsymbol{q}}$|⁠. The solutions of the upper and lower parts of the source satisfy the jump [see eqs. (2) and (7)] at the source:

Therefore, the column vector |$[ {{\boldsymbol{q}};{\boldsymbol{p}}} ]$| could be obtained as

and the surface y-numbers satisfying all conditions could be obtained:

Here, ‘,’ in expression means horizontal concatenation of two matrices, and ‘;’ means vertical concatenation.

For the solution of degree one and degree zero, special attention is required. One can refer to Takeuchi & Saito (1972), Saito (1974) and others (Farrell 1972; Dahlen 1974) for a detailed clarification. We will not repeat these contents in detail and only give some key prompts as follows. For degree one mode, only two of the three solutions |$( {{\boldsymbol{Y}}_1^ \uparrow ,{\boldsymbol{Y}}_2^ \uparrow ,{\boldsymbol{Y}}_3^ \uparrow } )$| propagating upward from the mantle are linearly independent, and the first two solutions were chosen. The third solution is |${\boldsymbol{Y}}_3^ \uparrow = {[ {a,0,a,0,ag( r ),0} ]}^{\rm{T}}$| (Farrell 1972; Okubo & Endo 1986). For degree zero mode, the 6 by 6 equation system in the solid core degenerates into two decoupled equations [see p. 240 in Takeuchi & Saito (1972)]. The initial values of the integrations from a small radius |${r}_0$| is given by eq. (11) in Saito (1974) and the general solution in the solid core is obtained from eq. (12) in Saito (1974). Then, the subsequent processing is like the previous discussion of the general case of |$n \ge 2$| mode.

After the above process, one can easily calculate Green's functions, including displacement, potential perturbation, gravity change and others, using y-numbers. This is more straightforward than using the Love number expression (defined later in this study or see Farrell's work). In addition, the difference between Euler's and Lagrange's descriptions should be paying attention when one expresses the change of potential or gravity. Here, we use a Lagrangian description of the gravity change on the deformed surface |$g( a ) = - ( {{{ {\frac{{{\rm d}[ {{\psi }_0( r ) + {\psi }_i( r )} ]}}{{{\rm d}r}}} |}}_{r = a + {u}_r( a )} - {{ {\frac{{{\rm d}{\psi }_0( r )}}{{{\rm d}r}}} |}}_{r = a}} ) = {g}^N + {g}^E$| (positive downward). Using Taylor expansion and ignoring higher-order terms, one can obtain|$\ g( a ) \cong - {u}_r( a ){\ddot{\phi }}_0( a ) - {\dot{\psi }}_i( a ) = - {u}_r( a )[ { - 4\pi G{\rho }_0( a ) + 2{g}_0( a )/a} ] - {\dot{\psi }}_i( a )$|⁠. It contains two contributions: the direct download gravitation of the internal mass layer |${g}^N$| which can be calculated using Newton's formula in eq. (16b), an additional gravitational effect |${g}^E$| in eq. (16a) caused by the deformation field, that is, the volume strain and the displacement of the density interface.

We listed the relationship of internal loading Love numbers and those of surface loading for comparison with Farrell's calculation (Farrell 1972) in the following. Assuming a point mass |${m}_0$| located at |${r}_s$| below the North Pole, it could be expanded harmonically as an areal density distribution over this sphere (Okubo & Saito 1983): |$\sigma = \mathop \sum \limits_{n = 0}^\infty {\sigma }_n{P}_n( {\cos \theta } )$| with |${\sigma }_n = \frac{{2n + 1}}{{4\pi {r}_s^2}}{m}_0$|⁠. The gravitational potential |${\psi }_2( r )$| due to this mass's gravitational effect (Okubo & Saito 1983) is:

Considering total potential perturbation is |$\psi$|⁠, then the geopotential change owing to the deformation field is |${\psi }_1 = \psi - {\psi }_2$| and the corresponding spherical harmonics coefficient is |${\psi }_{1,n} = {\psi }_n - {\psi }_{2,n} = {y}_5( r ) - {\psi }_{2,n}$|⁠.

The Love numbers for internal mass loading |$\{ {{h}_n( r ),{l}_n( r ),{k}_n( r )} \}$| are defined according to Farrell's classical work (Farrell 1972), with a minor adjustment incorporated, for example choosing |${\psi }_{2,n}( {{r}_s} )$| as reference and |${m}_0$| for unit potential by setting |$G{m}_0/{r}_s \equiv 1$| m² s⁻² (in the following formulas, we retain this constant factor for the purpose of distinguishing the dimensions of all physical quantities):

Here |${g}_0( a )$| is the surface mean gravity on the Earth. The internal loading Love number triplets could be obtained. And Green's functions of displacement and potential change due to unit internal mass loading for different r values are listed below for different radius intervals:

Here the factor |$n + 1$| appears in eq. (16b) because of the observation point is above the applied mass layer. And the infinite series summation formula involved in eq. (16b) can be verified through table 2 provided by Tang & Sun (2021). It is clear that when |${r}_s \to a$|⁠, the internal Love numbers and Green's function are identical to the expressions of Farrell (1972). This shows the consistency in the final Love number and Green's function expression form of the internal loading with the surface loading as its location tends to the surface.

2.3 Analytical expression of the displacement in an elastic half-space model

For deformation owing to internal unit force without self-gravitational effect, Okada (1992) listed the closed formulas, that is eq. (1) in his paper for internal displacement in a homogeneous elastic half-space based on classic works (Mindlin 1936; Press 1965). Assuming the elastic half-space with shear modulus |$\mu$| and Poisson's ratio |$\nu$|⁠, and taking the projection point of the source at the surface as the reference point, the vertical displacement |${u}_r$| (positive up) and horizontal displacement |${u}_\theta$| (positive when away from the source) at a distance of R to loading force F at depth d are:

When |$d \to 0$|⁠, the above expressions due to an internal force are identical to the classic solution of Boussinesq's problem: |${u}_r = - \frac{{1 - \nu }}{{2\pi R\mu }}$| and |${u}_\theta = - \frac{{1 - 2\nu }}{{4\pi R\mu }}$|⁠. These analytical formulae in elastic-half space can provide some reference for later calculation and verification, especially when the source tends to the surface from the inside. Here, we have not given the formula for surface gravity changes because it seems that it cannot be easily derived after many attempts after failing to find it in the literature.

3 VERIFICATION AND SIMULATIONS

After describing our calculation method and implementation algorithm, we verify the calculated Love number and Green's function in Appendix. Note that all calculations in this paper are about the deformation of the Earth's surface, not the Earth's interior, but for the source, there are two cases, one at the surface and the other in the interior.

Then, the difference between the deformation caused by internal and surface loadings must be shown. For very shallow sources, simulating the resulting displacements using a half-space model is easier than a spherical model. But for deeper sources, the necessity of the spherical earth model will be gradually demonstrated. For example, many studies have revealed the importance of the curvature of the Earth to seismic deformation (Sun & Okubo 2002). For internal load deformation, the importance of curvature should depend on the source depth and the horizontal distance away from the source. Therefore, we first discuss the difference in loading deformation caused by different shallow source models for the half-space model. After that, we discuss the differences between the half-space and spherical models for the deformation of sources located at different depths.

3.1 Surface deformation Green's functions in half-space model

For the half-space model, the deformation resulting from internal mass loading can be described by the analytical formulation in Section 2.3. We set up several very shallow sources, calculated their loading displacement and gave a visualization in Fig. 1 for source depth within 0–0.1 km and within 0–5 km.

Figure 1.

Deformation generated by several different depths of loading sources in the half-space. The first and third rows are for vertical displacement (⁠|${u}_r$|⁠) and second and fourth rows for horizontal displacement (⁠|${u}_\theta$|⁠), respectively. The second and fourth row regularizes the results of the first and third row, respectively. (a, b, c, d) share the same legend and (e, f, g, h) share the same legend. Lines of different thicknesses and colour in legend for different source depths. The result of a source at 0 km is indicated by the red dashed curve.

Open in new tabDownload slide

For the vertical displacement of different depth sources, its envelope is close to the deformation of the surface source, but the horizontal displacement does not have this characteristic. The shallower the source is, the larger the deformation on the surface will be, and the importance of source depth will be more obvious. The deformation caused by sources at different depths on the surface shows different decay characteristics with the increase in distance. For sources up to 100 m, vertical deformation is not sensitive to it, and horizontal displacement shows a slight sensitivity (as in Fig. 1); the deformation slightly decreases as the depth decreases. Results corresponding to sources deeper than 100 m show that the magnitude of the deformation becomes smaller as the source depth deepens.

Here we do not draw the result of the change in gravity because the simple formula seems unavailable in the half-space model. The direct gravitation effect can be simply calculated from Newton's Law of gravitation, but the gravity effect caused by the volume strain and the surface displacement is rather complex. Fortunately, the impact of different components on gravity can be easily considered within the spherical model, and some result will be shown later.

3.2 Green's functions in half-space model versus spherical earth model

We also calculated Green's function by several loading sources at different depths using the half-space and spherical models PREM with or without considering the self-gravitation. In Fig. 2, we list the normalized vertical and horizontal displacements for unit mass loading at 0.01 km, and 0.1 km, 1.0 km and 10 km. Here, the normalized factor is |${10}^{12}( {a\theta } ),$| as a function of horizontal distance |$a\theta$| from the mass source (here a denotes the Earth radius), as done by Farrell (1972). Note that for the half-space model, when the loading is at the surface, the displacement decays strictly in inverse proportion to the distance (see eq. 17), resulting in the grey dashed horizontal line after normalization in Fig. 2.

Figure 2.

Surface and internal loading displacement normalized Green's functions of the PREM model, a homogeneous sphere model, and a half-space model. (a, e), a source at 0.01 km, (b, f) at 0.1 km, (c, g) at 1.0 km, (d, h) at 10.0 km. Red lines with markers of cycles for the internal loading for the PREM model, blue lines for a homogeneous sphere without self-gravitation, solid black lines for surface loading in PREM, cyan lines for internal loading using the half-space model, magenta lines for surface loading using the half-space model. The unit for unnormalized displacement is m² kg⁻¹. The normalization factor of displacement is |${10}^{12}R = {10}^{12}a\theta$| (a of unit m and |$\theta$| is in radius).

Open in new tabDownload slide

The deformations exhibit differences in the normalized Green's functions between surface and internal loads over a horizontal distance of approximately 0.1 km when the source depth is 0.01 km and approximately a few kilometres when the source depth is 0.1 km, respectively. For source of depth 1.0 and 10 km, several models show differences at about 1–2 km away of source. In the nearfield range, all internal loading and surface loading behave consistently, both in the PREM model and in the homogeneous sphere model without self-gravitation, as well as in the half-space model. For the discussed internal sources, the uniform sphere model and the half-space model give good approximations within the 1 km range. As the source depth increases, the magnitudes of horizontal and vertical displacements decrease in the near field. The displacement Green's function varies inversely with distance for surface loading and exhibits more complex characteristics with distance increases for internal loading.

In the far field, the effects of internal and surface loading in the PREM model tend to be consistent, as shown by the black and red cycled lines. This phenomenon is predictable and understandable as the location of the source becomes less sensitive when viewed from a greater distance. In addition, in the far field, both the uniform sphere model and the half-space model differ very much from the PREM model, which demonstrates the superiority of the PREM model in studying global deformation.

As the source depth decreases, the results of both the PREM and half-space models approach the same value. Physically, the effect of the different depths of the source is more significant on higher-degree deformation than on low-degree deformation. This phenomenon leads to the difference in deformation in the near field and the consistency in the far field. Moreover, the gravity effect in spherical earth models has little impact on nearfield deformation but is essential for far-field deformation. This phenomenon is consistent with the results obtained by Pollitz (1992) in a study on seismic deformations.

3.3 Displacement and gravity in spherical earth model

We investigate the Green's functions of displacements |${u}_{r,\theta }$| and gravity changes |${g}^E$| due to different depth sources in PREM model. In Fig. 3, various regularization factors are used for displacement and gravity, black curves indicate surface loading sources, and other red lines of different thicknesses indicate loading sources at different depths.

Figure 3.

Normalized Green's functions of surface and internal loading displacements and gravity changes in the PREM model. Solid red lines of different thicknesses for the internal loading of different depths and solid black lines for surface loading. The unnormalized displacement and gravity change unit is m kg⁻¹ and (m s⁻²) kg⁻¹, respectively. The normalization factor near the y-axis is |${10}^{12}a\theta$| or |${10}^{18}a\theta$| (a of unit m and |$\theta$| is in radius).

Open in new tabDownload slide

The three physical quantities share a common feature: in the nearfield region, the deformation produced by the internal load decays faster than that of the surface load; the depth of the load in the extreme far-field region does not affect the deformation. Upon closer observation of vertical displacement and gravity changes, it is impossible to distinguish whether the source is located on the surface or in the interior beyond a range of two to three times the depth of the load source, as the corresponding Green's function curves nearly coincide. However, the horizontal displacement shows a unique characteristic; outside the depth range of 30–40 times, the Green's function curves of the internal source and the surface source only tend to coincide.

In addition, one also needs to pay attention to the contribution of the gravitation of the loading mass itself to the surface gravity change. It will significantly depend on the relationship between the observation position and the relative position of the source. Fortunately, most studies can directly calculate this part by integrating the gravitational force.

3.4 An application case of conceived disk-like internal mass loading

Various physical scenarios involve the movement of the mass within the Earth, such as the movement of groundwater and magma, and underground ore mining. For example, the mass change in Lava tube (Sawłowicz 2020) or underground rivers in karst areas (Pivetta et al. 2021) may be of relatively large magnitude, especially during extreme precipitation events. In addition, groundwater extraction also involves a large amount of mass, resulting in a certain loading displacement and poroelastic stress change (Wang & Kümpel 2003). To briefly demonstrate the possible application scenarios of the above theory. Here, we assumed a disk-shaped mass layer source and subsequently calculated the surface displacement field generated by the source when located at the surface or underground. In this example, we will not compare the theoretical calculated values with the observed data, which will be carried out in future work.

We consider two disc-shaped mass layers with different radii. The first mass layer has a radius of 0.5 km, while the second mass layer has a radius of 20 km. Both layers have a height of 1 meter and a density of 1000 kg m⁻³ at their centres. From the centre to the edges of each layer, the density linearly decreases to zero (see Fig. 4). We then visualized its loading displacement, as shown in Fig. 5.

Figure 4.

Schematic diagram of underground mass migration. (a) Extraction of groundwater from underground aquifers, (b) magma formed from deep flows into magma chamber along the original channel and (c) approximating the mass change underground as a conical thin mass layer with height of 1 m.

Open in new tabDownload slide

Figure 5.

Disc loading displacement from internal and surface mass loading. Assuming a disc-shaped mass layer with a height of 1 meter, where the density is 1000 kg m⁻³ at the centre and linearly decreases to zero towards the edges. (a, b) Disc loading of radius 0.5 km, (c, d) Disc loading of radius 20 km. (a, c) for vertical displacement and (b, d) for horizontal displacement. Different styles of thin lines for different depths, and thick red lines for surface loading.

Open in new tabDownload slide

As the distance from the centre of the disc gradually increases, the vertical displacement gradually decreases, the horizontal displacement presents a single peak phenomenon and the largest displacement occurs near the edge of the disc. As regards mass distribution with a radius of 0.5 km, in the mass coverage area, the loading displacements of the 1 km source and 0 km source could reach a difference of two to three times. However, these differences decrease to approximately 10 per cent. We found that when the range of mass distribution expands gradually, the loading displacement caused by internal mass and external mass gradually tends to become consistent, as expected in physics (see Figs. 5c and 5d).

4 DISCUSSION

To show the influence of different source depths on deformation, we calculated the deformation Green's function generated by sources at varying depths from 0.01 to 50 km. When the source of immovable depth is within 10 km, a significant difference occurs in the deformation within 10 km away of the source. Interestingly, the envelope curve of Green's functions at these depths is consistent with the deformation of the surface source. Under the current normalization factor, for the several different depths considered, the deformation of the internal source tends to 0 with decreasing distance. This means that the deformation of the internal source on the surface is smaller than that of the surface source. In addition, horizontal displacement and vertical displacement show slight differences. For horizontal deformation, sources at different depths need larger distances of approximately 100 km to approach the same values.

We also considered a source shallower than 0.01 km for the elastic half-space model. Our results showed that as the source depth tends to 0, the corresponding Green's function tends to the Green's function of the surface source. In this case, the semi-infinite elastic space is a good approximation for studying displacement and deformation problems. However, more efforts are needed to get the final analytical expression when it comes to geoid or gravity-related physical quantities.

The mass source's location and the earth model's structure should be investigated for loading problems. Considering the depth variation of groundwater distribution, using the internal loading Green's function can obtain more accurate calculation results in the near area than Green's functions due to surface loading. Strictly speaking, surface and internal loads approximate this kind of mass load problem, especially when considering the complexity of other geometric or physical properties such as actual terrain, heterogeneity, etc.

When studying crust deformation or gravity changes resulting from long-term pumping ground water and gas extraction, a two-phase theory considering pore-elastic deformation and consolidation deformation (Wang & Kümpel 2003; Schmidt & Bürgmann 2003; Teatini et al. 2011; Candela & Koster 2022) can be a suitable choice. However, the surface and internal loading effects due to the underground mass itself should also be considered (Amos et al. 2014; Hammond et al. 2016). These processes may all simultaneously play important roles, although the specific magnitudes of these effects may depend on the problem case by case. It is necessary to develop a comprehensive theory that accounts for multiple effects in the same study and conduct comparative analyses using observed data for specific issues. Furthermore, the complexity of the underground medium's structure and its mechanical and hydraulic behaviour adds to the complexity of the problem.

However, it is difficult to semi-analytically solve the elastic equations and obtain the theoretical deformation value directly aimed at a 3-D solid earth model comprising the topography, irregular internal interface and even lateral inhomogeneity (Kuriyama & Mizuta 1993; Pan & Amadei 1995; Williams & Wadge 2000; Zhou & Chen 2008; Lan & Zhang 2011). Owing to the computational complexity of 3-D models, this situation cannot be solved concisely, even using numerical techniques such as finite element, finite difference or approximate spectral methods (Tromp & Mitrovica 1999a,b, 2000; Wong & Wu 2019). The most realistic and accessible option is retaining the layered sphere model but considering the loading deformations caused by mass redistribution inside the solid Earth. Regarding the development of the loading theory itself, dealing with the deformation of 3-D models by mixed symbolic-numeric methods appears promising. In this regard, further studying is the following research work.

5 CONCLUSIONS

In this study, we partly extended the classic loading theory (Longman 1962, 1963; Farrell 1972) by solving the boundary value problem due to internal mass loading instead of the surface mass redistribution. We presented the source function and boundary conditions corresponding to the internal mass source, solved the corresponding elastic equation, and finally obtained the surface deformation Green's function generated by sources at arbitrary depth. We also presented the analytical solutions of two uniform sphere models. We overcame the computational difficulty of the Green's functions of shallow sources in the nearfield by using analytical solutions and Kummer transforms. Then, the results of several spherical models were compared and verified against the half-space model results. The consistency of the results proves the correctness of our calculations.

A comparison of sources at different depths indicated that the source depth significantly affects the Green's function calculation in the nearfield. The depth of the source has a greater impact on horizontal displacement compared to vertical displacement and gravity change. In general, the deformation of the internal source in the nearfield is smaller than that of the surface source, and the deformation of the internal source in the far field is consistent with that of the surface source. Away from the source at horizontal observation distances beyond a range of three to five times the depth of the source, the surface loading can approximate the deformation produced by the internal loading. The semi-infinite space model can be a useful approximation when studying internal deformations caused by very shallow internal sources. In such cases, simple analytic expressions for displacements may exist, but obtaining expressions for gravity change can be difficult. The deformation caused by an internal loading cannot be expressed analytically using the layered sphere model, but Love numbers and Green's function still be used to easily calculate physical quantities such as displacement, gravity, geoid and others in the range of different observation distances, considering conditions such as the compressibility and self-gravitation of the Earth. Therefore, we think that for the general mass loading problem, if it involves the change of mass in the vertical direction, it is recommended to use the internal loading Green's function, because it will not greatly increase the complexity of modelling.

To facilitate the study of internal load deformation using our model, we have published the corresponding computer code. An example of an application that illustrates the difference of surface displacements by disc-like lake mass change when located at the surface or underground. We did not conduct a theoretical and observational comparative study on actual observational data due to space limitations. However, our examples could encourage readers to participate in discussing this topic and stimulate more research.

Acknowledgement

We would like to express our sincere gratitude to the two anonymous reviewers, as well as the Editor, Prof. Kosuke Heki. Their constructive and thoughtful comments greatly enhanced the quality of this work. HT would also like to thank Dr. Yang Wang for helpful discussions on the volcanic process during the preparation of this paper. The following foundations support HT, namely National Natural Science Foundation of China (No. 42104002), China Postdoctoral Science Foundation (No. 2020M680649), Fundamental Research Funds for the Central Universities (No. E1E40415X2 and E2ET0411X2) and State Key Laboratory of Geodesy and Earth's Dynamics, Innovation Academy for Precision Measurement Science and Technology, Chinese Academy of Sciences (No. SKLGED2021-4-2). The following foundations support other authors, namely National Natural Science Foundation of China (Nos. 42174097 and 41974093).

CONFLICT OF INTEREST

No financial interests or connections, direct or indirect or other situations might affect this work.

DATA AVAILABILITY

No observational data were utilized in this study. All the information presented in the figures can be replicated using the provided formulas. The internal loading Green's function and the associated Fortran program are freely available in the public repository (https://github.com/UCAStanghe2014/internal-mass-loading-Green-functions) or upon reasonable request from H. Tang ([email protected]).

References

APPENDIX: ANALYTICAL SOLUTION FOR THE HOMOGENEOUS EARTH MODEL AND USES IT TO VERIFY THE NUMERICAL RESULT

A. Earth model with self-gravitation but not compressibility

According to our definition of y-notation, the explicit form of matrix A in eq. (1b) is

For a homogeneous self-gravitational sphere with incompressibility, the analytical solution of |${y}_i( r )$| is obtained by a matrix method. In this case, the coefficient matrix |${\boldsymbol{A}}$| degenerates, as shown:

The general solution |${{\boldsymbol{Y}}}^{\boldsymbol{*}}( r ) = [ {{\boldsymbol{Y}}_L^{\boldsymbol{*}}( r ),{\boldsymbol{Y}}_U^{\boldsymbol{*}}( r )} ]$| of order |$6{\rm{\ }} \times {\rm{\ }}6$| with the first three solutions |${\boldsymbol{Y}}_L^{\boldsymbol{*}}( r )$| are regular at the centre of the Earth, and the others |${\boldsymbol{Y}}_U^{\boldsymbol{*}}( r )$| are regular at infinity:

The propagation matrix could be set as

Subsequently, a solution |${\boldsymbol{Y}}( {{r}_0} )$| and |${\boldsymbol{Y}}( r )$| at different locations is expressed by |${\boldsymbol{Y}}( r ) = {\boldsymbol{\Omega }}( {r,{r}_0} ){\boldsymbol{Y}}( {{r}_0} )$|⁠. We set the particular solution at the free surface as |${\boldsymbol{Y}}( {\rm{a}} ) = {\boldsymbol{Dq}}$|⁠, with unknown vector |${\boldsymbol{q}}$| of order |$3{\rm{\ }} \times {\rm{\ }}1$| and a constant matrix |${\boldsymbol{D}}$| of order |$6{\rm{\ }} \times {\rm{\ }}3$| with non-zero elements:

Subsequently, the solution above the source is

The solution below the source is

with unknown vector |${\boldsymbol{p}}$| of order |$3{\rm{\ }} \times {\rm{\ }}1$|⁠.

Subsequently, considering the definition of the source function, the following is obtained:

The above equations are solved analytically to obtain unknown vectors:

Finally, the surface solution |${y}_i( a )$| is written analytically as|${\boldsymbol{\ Y}}( a ) = {\boldsymbol{Dq}}$|⁠:

with a defined auxiliary variable:|$\ \frac{{{r}_s}}{a} = \delta ,\ \frac{{{\rho }_0}}{{\bar{\rho }}} = \alpha ,\frac{\mu }{{a{g}_0\bar{\rho }}} = \beta ,G = \frac{{3{g}_0}}{{4\pi a\bar{\rho }}}$| here |$\bar{\rho }$| is the mean density of the Earth.

The asymptotic expansion of y-numbers as |${y}_{i\infty }( a ) = {\boldsymbol{N}}( n ){\boldsymbol{y}}_i^{asy}$|⁠:

B. Earth model with compressibility but not self-gravitation

Similarly, for a homogeneous compressible sphere without considering self-gravitation, the coefficient matrix |${\boldsymbol{A}}$| of order |$4{\rm{\ }} \times {\rm{\ }}4$| and |${\boldsymbol{Y}} = {( {{y}_1,{y}_2,{y}_3,{y}_4} )}^{\rm{T}}$| is:

The general solution |${{\boldsymbol{Y}}}^{\boldsymbol{*}}( r ) = [ {{\boldsymbol{Y}}_L^{\boldsymbol{*}}( r ),{\boldsymbol{Y}}_U^{\boldsymbol{*}}( r )} ]$| of order |$4{\rm{\ }} \times {\rm{\ }}4,$| with the first two solutions |${\boldsymbol{Y}}_L^{\boldsymbol{*}}( r )$| are regular at the centre of Earth, and the others |${\boldsymbol{Y}}_U^{\boldsymbol{*}}( r )$| are regular at infinity:

Using a similar approach, that is setting the special solution at the free surface as |${\boldsymbol{Y}}( a ) = {\boldsymbol{Dq}},$| with unknown vector |${\boldsymbol{q}}$| of order |$2{\rm{\ }} \times {\rm{\ }}1$||${\rm{and}}$| a constant matrix |${\boldsymbol{D}}$| of order |$4{\rm{\ }} \times {\rm{\ }}1$| is obtained:

and source function:

results in the surface analytical solution |${y}_i( a )$| as:

with the first Lamé parameter |$\lambda$| expressed as |$\lambda = \frac{{2\nu \mu }}{{1 - 2\nu }}$| by Poisson's ratio |$\nu$|⁠.

Also, the asymptotic expansion of y-numbers as |${y}_{i\infty }( a ) = {\boldsymbol{N}}( n ){\boldsymbol{y}}_i^{asy}$| (here |${\boldsymbol{N}}( n )$| is defined in eq. A13):

C. Analytical Green's functions for a homogeneous sphere

To obtain the closed form of Green's function, the asymptotic expansion of y-numbers is required, as |${y}_{i\infty }( a ) = {\boldsymbol{N}}( n ){\boldsymbol{y}}_i^{asy}$| in above. Subsequently, the closed-form Green's functions of displacement |${u}_r$| and |${u}_\theta$| are:

where |${p}_i$| and |${\rm d}{p}_i$| are the analytical sums for an infinite series involving Legendre functions |$\mathop \sum \limits_{n = 0}^\infty {\boldsymbol{N}}( n ){P}_n( \theta )$| and |$\mathop \sum \limits_{n = 0}^\infty {\boldsymbol{N}}( n )\frac{{{\rm d}{P}_n( \theta )}}{{{\rm d}\theta }}$| in tables 1 and 2, given by Tang & Sun (2021). We call these approximate expressions asymptotic Green's functions here.

The above asymptotic Love numbers and the asymptotic Green's functions play a critical role in accelerating the convergence of the series of Love numbers. Usually, when directly calculating the partial sum of this series involving numerical Love numbers of a layered earth model, convergence difficulties would inevitably appear, particularly when the angular distance |$\theta$| is relatively small. As Farrell (1972) emphasized, it is the terms with n large that dominate the spherical harmonic expansions of the loading response in the near field. Fortunately, the Kummer's transformation could be used to for improving the rate of convergence of the series. The asymptotic Love numbers are deducted from the Love numbers and then the residual series can be safely truncated after hundreds of items are summed. The contribution of the subtracted series to the final Green's functions can be analytically obtained. Two further devices, that is, putting a disk factor and using Euler's transformation can be applied to speed the convergence of the residue series (Farrell 1972; Martens et al. 2019). Different physical quantities have different convergence rates, for example, tilt and strain converge more slowly than displacement. Overall, applying these mathematical skills, we found that it ensures the accuracy of the final Green's function when the angular distance is in the range |${10}^{ - 6} < \theta \le \pi$| (0.00006°–180°). Alternatively, an analytical solution for a half-space model would work well when the convergence difficulty exists at small distances. However, we prefer to express Green's function using a spherical earth model to ensure that the loading deformation is expressed under a consistent earth model.

D. Verification of the Love numbers and Green's functions

The numerical and analytical solutions were compared to verify the correctness of our method. Here, we compare the computational results of the three models for verification. Model 1: an approximate homogeneous earth model with a core. Model 2: Strict homogeneous elastic sphere model without core. Model 3: Elastic half-space model. Numerical methods are used for the first model, and analytical methods are used for the latter two. Because our Fortran program is developed for a layered model that includes a solid core, a liquid outer core and a solid mantle-crust, such as PREM, we can approximate the homogeneous elastic sphere model (Model 2) by a spherical earth model with inner and outer cores (Model 1) for the sake of simplicity. The following calculations will show the consistency of Model 1 and Model 2 in terms of higher order spherical harmonic deformations and the consistency of the three models in nearfield Green's functions.

As an example, a unit loading force is buried at a depth of 5 km both in Model 1. We assumed no self-gravitation in the layered earth model (Model 1) by reducing the gravity to a minimal value and compared |${y}_1$| and|${\rm{\ }}{y}_3$| with their analytical results in a simple sphere (Model 2). In order to enable interested researchers to verify our results, we do not normalize y-numbers here. The y-numbers of each degree are shown in Fig. A1. After approximately 500°, the two results tend to be consistent. The difference in low-degree results indicates the influence of the inner and outer cores of the earth model (Note the intrinsic difference between Model 1 and Model 2). The consistency of higher degree results shows the stability and accuracy of the numerical calculation.

Fig. A2 shows the results of |${y}_1$|⁠, |${y}_3$| and |${y}_5$| for the case considering self-gravitation but ignoring compressibility. Similarly, the numerical result is obtained by using a uniform mantle layer in layered earth model (as Model 1). The analytical results are obtained using Model 2, a sphere model. Again, a good consistency of higher degree y-numbers is obtained and the difference of low degree results is due to the fact that one of the models containing a core.

Figure A1.

Numerical and analytical surface y-numbers for a sphere without gravitation: (a, b) |${y}_1$| and (c, d) |${y}_3$|⁠. The first column indicates the results of the two models: The blue cycles for numerical result of Model 1, and the red lines for analytical result of Model 2. The second column indicates the ratio of Model 1 to Model 2. The dashed horizontal line indicates a value of 1. Here, a unit loading force is buried at a depth of 5 km. Note Model 1 containing a core which is slightly different to Model 2. Refer to the first paragraph of this section for a detailed description of the three models.

Open in new tabDownload slide

Figure A2.

Numerical and analytical surface y-numbers for a sphere without compressibility: (a, b) |${y}_1$|⁠, (c, d) |${y}_3$| and (e, f) |${y}_5$|⁠. The blue cycles indicate numerical results of Model 1, and the red lines indicate analytical results of Model 2. The second column indicates the ratio of Model 1 to Model 2. The dashed horizontal line indicates a value of 1. Here, a unit loading force is buried at a depth of 5 km.

Open in new tabDownload slide

Because the deformation field decays as the distance from the source increases, to clearly show the nearfield and far-field Green's functions in the same plot, we normalized the displacement Green's function by multiplied a factor |${10}^{15}a\theta$|⁠. Fig. A3 shows a good agreement of the corresponding displacement Green's function for three models in the near field. And the difference between the Green's functions of a spherical earth model with a core (Model 1) and a strict homogeneous sphere (Model 2) in the far field can be also illustrated in the figure. The ratio of Model 1 to Model 3 results also shows the gradually increasing difference between the half-space model and the spherical earth model outside a certain range, approximately 5 times the source depth. Because the main purpose of this section is to verify effectiveness of numerical calculation method, we will not further discuss this point here. Please refer to the main text for a discussion of the necessity of using a spherical earth model.

Figure A3.

Numerical and analytical surface displacement Green's function normalized by a factor |${10}^{15}a\theta$| for a sphere without gravitation from a unit loading force measured at a depth of 5 km. (a, b) Vertical displacement and (c, d) horizontal displacement. (b, d) The ratio of the displacement of a spherical model (Model 1 and Model 2) to that of a half-space model (Model 3).

Open in new tabDownload slide