The Earth’s Gravity Field Role in Geodesy and Large-Scale Geophysics

Abstract

The Earth gravity field is a signature of the Earth’s mass heterogeneities and structures and applied in Geodesy and Geophysics for different purposes. One of the main goals of Geodesy is to determine the physical shape of the Earth, geoid, as a reference for heights, but Geophysics aims to understand the Earth’s interior. In this chapter, the general principles of geoid determination using the well-known methods of Remove-Compute-Restore, Stokes-Helmert and least-squares modification of Stokes’ formula with additive corrections are shortly discussed. Later, some Geophysical applications like modelling the Mohorovičić discontinuity and density contrast between crust and uppermantle, elastic thickness, ocean depth, sediment and ice thicknesses, sub-lithospheric and lithospheric stress, Earthquakes and epicentres, post-glacial rebound, groundwater storage are discussed. The goal of this chapter is to briefly present the roll of gravity in these subjects.

1. Introduction

The Earth’s gravity field reflects of the Earth’s interior and is an interesting subject in Geodesy and Geophysics with various applications. Geodesy aims to determine three types of the shape and size of the Earth, the Earth’s surface, geoid as the physical shape, reference ellipsoid as the mathematical one. Physical Geodesy deals with determination of the physical shape of the Earth or the geoid, which is a reference for heights, from gravimetric data. In this chapter, short descriptions of three known methods of geoid determination such as Remove-Compute-Restore (RCR) [1], Stokes-Helmert (SH) [2] and least-squares modification of the Stokes formula with addition corrections (LSMSF) [3] are presented.

In Geophysics, understanding the Earth’s physics, dynamics and interior geometry is of interest using such data. Gravity measurements can be analysed over small or large area depending on the geophysical purpose. For instance, in exploration Geophysics they are used to detect or discover near surface resources and for such a goal precision and accuracy of these data should be high. Here, such applications are named small-scale Geophysics. However, understanding or studying the deep Earth’s interior physics, dynamics or geometry does not require high spatial resolutions and long wavelength portions of the gravity data are more suitable. In addition, large areas are considered for such purpose and therefore, here, such subjects are called large-scale Geophysics. Some of these large-scale phenomena are modelling the Mohorovičić discontinuity, elastic thickness of the lithosphere, sub-lithospheric/lithospheric stress, and thickness of ocean water, sediments, and ice; land uplift, mantle viscosity and groundwater storage; and post-seismic studies of Earthquakes, detecting the epicentre points of shallow Earthquakes, which are briefly presented in this chapter.

2. Geoid determination as a purpose in Geodesy

The geoid is a reference surface for heights and if this reference is not enough precise and accurate, all determined heights from it will be unreliable. Having a precise geoid model simplifies the lengthy and costly work of levelling by simply using a global navigation satellite systems (GNSS) receiver, the height above the geoid can be determined. However, to reach to this goal, deep knowledge about the Earth’s gravity field and the Physical Geodesy theories, skills in numerical modelling and precise data in all frequency bands are required.

The main task of Physical Geodesy is to develop theories and methods to model a precise geoid. Different approaches have been developed toward this goal. As known, the surface terrestrial gravity data are sensitive to high frequencies and near surface mass variations, but their low frequencies of the signal are weak unlike satellite-only Earth gravity models (EGMs) having better qualities only in the low frequency band. In geoid modelling approaches the terrestrial gravity data are used for recovering the high frequencies of the geoid and the satellite EGMs for the lower. Generally, in geoid modelling the following issues should also be considered:

1. The effects of topographic and atmospheric masses.

2. Downward continuation of gravity data.

3. Conversion of gravity anomalies to geoid/co-geoid.

The differences between the geoid modelling methods are related to how mathematically these issues are handled. In the following the three methods of RCR, SH and LSMSF are shortly presented.

2.1 Remove-compute-restore approach

In the RCR scheme, the low and high frequencies of terrestrial gravity data are removed by an EGM and topographic heights. Because the low frequencies are global and for converting of the gravity data with a regional coverage to the geoid, the low frequencies of geoid cannot be recovered well. In addition, removing the effect of topographic terrain makes the gravity data smoother, and simplifies the computations. After computing the geoid excluding the low and high frequencies, the removed frequencies are restored to it to complete all frequencies of the geoid. Mathematically, the idea is presented by:

$$N=\frac{R}{4\mathit{\pi \gamma }}\underset{\sigma }{\iint }S\left(\psi \right)\left({\left(\mathrm{\Delta }g-\mathrm{\Delta }{g}^{\mathrm{RTM}}\right)}^{\ast }-\sum _{n=2}^L\mathrm{\Delta }{g}_n^{\mathrm{EGM}}\right)\mathrm{d}\sigma +\sum _{n=2}^L{N}_n^{\mathrm{EGM}}+N^{\mathrm{RTM}}$$E1

where

$$\left.\begin{array}{l}\mathrm{\Delta }{g}_n^{\mathrm{EGM}}\\ N_n^{\mathrm{EGM}}\end{array}\right\}=\sum _{m=-n}^n\left.\begin{array}{l}\mathrm{\Delta }{g}_{\mathit{nm}}^{\mathrm{EGM}}\\ N_{\mathit{nm}}^{\mathrm{EGM}}\end{array}\right\}{Y}_{\mathit{nm}}\left(\theta ,\lambda \right)$$E2

N is the geoid height, $S\left(\psi \right)$ the Stokes function [4] converting gravity to geoid, $\mathrm{\Delta }g$ the gravity anomaly, $\mathrm{\Delta }{g}^{\mathrm{RTM}}$ the residual terrain effect (RTM) on the anomaly, ${\left(\mathrm{\Delta }g-\mathrm{\Delta }{g}^{\mathrm{RTM}}\right)}^{\ast }$ means the downward continued $\mathrm{\Delta }g$- $\mathrm{\Delta }{g}^{\mathrm{RTM}}$,${\sigma }_0$ the integration domain, $\mathrm{d}\sigma$ the surface integration element, $\mathrm{\Delta }{g}_{\mathit{nm}}^{\mathrm{EGM}}$ and $N_{\mathit{nm}}^{\mathrm{EGM}}$ are, respectively, the spherical harmonic coefficients of $\mathrm{\Delta }g$ and N of degree n and order m, derived from an EGM, limited to the maximum degree L, $Y_{\mathit{nm}}\left(\theta ,\lambda \right)$ the spherical harmonic with arguments of co-latitude $\theta$ and longitude $\lambda$,$N^{\mathrm{RTM}}$ the restored RTM effect on the computed geoid.

The first term on the right-hand side (rhs) of Eq. (1), is the Stokes integral, which converts the gravity anomalies to geoid height, and since the long and short frequencies of the anomalies are removed, the solution of this integral is the geoid height excluding these frequencies. In fact, the addition of the second and third terms of Eq. (1) is to restore these frequencies back to the computed geoid height. In the following, some issues regarding the RCR method is presented and discussed:

1. Three types of data are used in Eq. (1), terrestrial gravity data, EGM and RTM with own error properties. According to the error propagation law, the error of the reduced gravity anomalies is the square root of summation of variances of the terrestrial data, EGM and RTM, which is surely larger than the error of gravity data. If the discretisation error of the Stokes integral is assumed small, in the restoring step, the errors of the EGM and RTM will be propagated to the final solution again.

2. Most portion of the geoid signal comes from its low frequencies, and by removing it by an EGM, this part of signal is assumed as known. By increasing the maximum degree of it, more portion of the geoid signal removed and restored. This means that the sensitivity to the terrestrial data will be reduced. The high frequencies of the geoid comes from topographic masses, and by considering it known as well, the main task will be to recover the medium wavelengths from gravity data.

3. After the removal step, gravity data are converted to a medium frequency geoid height, using Stokes formula, least-squares collocation (LSC) [5] or Fast Fourier Transform (FFT); see e.g. [1]. The low frequencies are restored by the same EGM applied limited to the same maximum degree, and the high frequencies from the RTM effect.

4. Downward continuation (DWC) of the gravity data should be performed before applying the Stokes formula or FFT. However, by using LSC, the conversion of the gravity data to the geoid heights and DWC can be done in one-step simultaneously.

2.2 Stokes-Helmert approach

The Stokes-Helmert (SH) method was proposed by Vanicek and Martinec [2] and developed further by Martinec [6]. Theoretically, the gravity data on a spherical surface are needed to numerically solve the Stokes integral for computing a geoid height. In addition, this integral is the solution of the gravimetric boundary-problem, the Laplace equation, with gravity anomalies at the boundary, the geoid. This means that this solution is theoretically valid where there is no mass outside the boundary surface. However, in practice, the gravity data are collected at the Earth’s surface, on topographic and under atmospheric masses. The presence of such masses violates the theory. Therefore, the gravitational effects on the gravity data should be removed to fulfil the Laplace equation. The result will be a no-topography and no-atmosphere computational space, or the Helmert space. After removing these effects the gravity data still remain above the boundary and need to be continued downward. By solving the Stokes integral numerically, these continued data are converted to a surface similar to geoid, known as co-geoid. The next step will be to convert this co-geoid by restoring the effects of topographic and atmospheric masses. The principle of SH method is:

$$N=\frac{R}{4\mathit{\pi \gamma }}\underset{{\sigma }_0}{\iint }{S}^L\left(\psi \right)\left(\mathrm{\Delta }{g}^L-\mathrm{\Delta }{g}_{\mathrm{dir}}^{\mathrm{TAe}}\right)\mathrm{d}\sigma +\frac{R}{2\gamma }\sum _{n=2}^L\frac{2}{n-1}\mathrm{\Delta }{g}_n^{\mathrm{EGM}}+N_{\mathrm{Ind}}^{\mathrm{TAe}}$$E3

where $S^L\left(\psi \right)$ is the modified Stokes function, $\mathrm{\Delta }{g}_{\mathrm{dir}}^{\mathrm{TAe}}$ the joint direct effect of topographic and atmospheric masses as well as the ellipticity of the Earth, which should be removed from the gravity data. $\mathrm{\Delta }{g}^L$ means the gravity anomalies excluding the frequencies to degree L. $N_{\mathrm{Ind}}^{\mathrm{TAe}}$ is the joint indirect effect of the removed masses and ellipticty.

The SH method has the following properties:

1. The topographic and atmospheric effects are removed from the gravity anomalies directly. The topographic effect (TE) over mountainous areas is considerably larger than the gravity anomalies, then a compensation/condensation mechanism is required to reduce their values to the order of the terrestrial ones. The anomalies will be smoother and the Stokes integral can be solved numerically with a better precision. The same mechanism is used for restoring the TE as Indirect TE. This means that it has no effect on the resulted geoid height, because of being added and subtracted during the process.

2. The TE and atmospheric effect (AE) are computed by taking the radial derivative of their respective gravitational potentials. This means the effect on the gravity, and not the gravity anomaly which is used in the Stokes formula. There are two terms in the fundamental equation of Physical Geodesy as the definition of gravity anomaly. The first term is the radial derivative of the disturbing potential and the second is 2/R times of the potential. Considering the TE and AE on the gravity solely by taking the radial derivatives of their potentials, means ignorance of the second term. Therefore, the restoration is done in two step, primary indirect effects, in which the removed effected are restored, and the secondary indirect effect when the effect of the missing second term is restored. If the fundamental equation of Physical Geodesy is applied for computing the direct and indirect TE and AE, this secondary effect will not be needed.

3. The Stokes-integral is modified meaning that its kernel function is changed in such a way that the contribution of the anomalies outside the integration cap is minimised. The effect of the truncation of the integration domain will be restored after integrating the reduced anomalies.

4. The TE and AE have their own error properties, in addition by removing the long wavelengths of the anomalies by an EGM, Therefore, the reduced anomalies contaminate larger stochastic error than the measured ones. The errors of the gravity anomalies and EGMs are not considered in the solution.

5. The reduced gravity anomalies in the Helmert space need to be continued downward to see level prior to integrating them. To do so, inverse solution of the Poisson integral is applied; see [6], which is an ill-posed problem and complicated when the resolution of the anomalies is high.

2.3 Least-squares modification of stokes formula with additive corrections

Unlike the RCR and SH methods, neither the TE and AE nor the long wavelength portion of the anomalies are removed from the gravity anomalies. Instead, the terrestrial anomalies and EGM are spectrally weighted which means that the Stokes integral is modified in such a way their errors and the truncation error of the Stokes integral outside the integration cap are minimised in a least-squares sense. This method is called least-squares modification of Stokes formula [3]. In this method, the terrestrial anomalies are integrated directly by the modified Stokes formula to estimate a geoid model. Later the total TE and AE, DWC and ellipsoidal corrections will be added to the modelled geoid to make it precise. This method can be mathematically presented by:

$$N=\frac{R}{4\mathit{\pi \gamma }}\underset{{\sigma }_0}{\iint }{S}^L\left(\psi \right)\mathrm{\Delta }g\mathrm{d}\sigma +\frac{R}{2\gamma }\sum _{n=2}^L{b}_n\mathrm{\Delta }{g}_n^{\mathrm{EGM}}+N^{\mathrm{T}}+N^{\mathrm{A}}+N^e+N^{\mathrm{DWC}}$$E4

where $b_n$ is a parameter depending on type of modification, $N^{\mathrm{T}}$ and $N^{\mathrm{A}}$ are, respectively, the total TE and AE. $N^e$ is the ellipsoidal correction, $N^{\mathrm{DWC}}$ the DWC effect on the geoid.

The properties of this method are:

1. The measured gravity anomalies are used in the modified Stokes integral. However, gridded anomalies are not at the boundary surface, which is not theoretically corrected, also no mass should exist outside the geoid when applying the stokes formula.

2. This method considers the errors of the terrestrial data, EGM and truncation of the integral formula and modify the integral in an optimal way, meaning that the quality of the data play important role in geoid modelling. The data contributes to solution according to their precision.

3. Because of neglecting the TE and AE on the gravity data, results of the modified Stokes integral will contain biases. However, the total TE and AE will be removed. In fact, the gravitational potential of these masses are computed for points at the surface of the Earth. Later they are continued downward to the boundary, and subtracted from the indirect gravitational potential of the points under the masses at the boundary. Such a potential will be converted simply to correction to geoid using the Bruns formula. Note that no compensation or condensation mechanism is required in this method

4. The DWC process is done directly on the potential, the gravity data converted to the potential and continued downward analytically. Therefore, no integral equation needs to be solved numerically.

5. The effect of ignoring the ellipticity of the Earth will be considered as an extra correction to the geoid directly.

3. Gravity field and large-scale Geophysics

In Geophysics, the gravimetric data are used for different purposes; e.g. in exploration and prospecting for detecting near surface sources, or studying the Earth’s deep interior, which are named here the large-scale Geophysics. The Earth gravity field is determined in two ways. If the temporal variations of the gravity is considered the time-variable gravity field can be determined, otherwise, the static field. In this section, some of the well-known applications of static and time-variable gravity data in large-scale Geophysics are presented and discussed.

3.1 Static gravity field and large-scale Geophysics

A static gravity field reflects the physics of the Earth’s interior, which is not fully known. Therefore, different assumptions are used to extract the desired information from the gravity field. Here, the use of the static gravity data and modelling of crustal structure, elastic thickness and rigidity, ice thickness, bathymetry, sediment basement, lithospheric and sub-lithospheric stresses due to mantle convection are presented briefly.

3.1.1 Determination of Moho depth

One the assumptions about the Earth’s interior is Isostasy, which is a state of equilibrium between the crust and upper mantle. Aity-Heiskanen, in which the mountains have roots beneath to keep them in isostatic balance, and Partt-Hayford theory, which states that the mountains loads are compensated by density variations inside the crust are two known models of Isostasy. The gravimetric isostasy mean that the isostatic gravity anomaly ($\mathrm{\Delta }{g}^{\mathrm{I}}$) should be zero to have the crust in isostatic equilibrium. The mathematical description of the gravimetric isostasy is [4]:

$$\mathrm{\Delta }{g}^{\mathrm{I}}=\mathrm{\Delta }g-\mathrm{\Delta }{g}^{\text{TBSCI}}+\mathrm{\Delta }{g}^{\mathrm{CMP}}=0$$E5

where $\mathrm{\Delta }g$ is gravity anomaly, $\mathrm{\Delta }{g}^{\text{TBSCI}}$ total effect of the topographic and bathymetric masses, sediments, crustal crystalline and ice on $\mathrm{\Delta }g$ and finally, $\mathrm{\Delta }{g}^{\mathrm{CMP}}$ the compensation effect on $\mathrm{\Delta }g$. Eq. (5) means that there are some compensation attraction, which is equal to the gravitational difference between the effect of loads on the crust and gravity.

In Eq. (5), when $\mathrm{\Delta }g$ =0, then $\mathrm{\Delta }{g}^{\text{TBSCI}}=\mathrm{\Delta }{g}^{\mathrm{CMP}}$, meaning that the gravimetric isostasy becomes the Airy-Heiskanen model having a local compensation property. The presence of $\mathrm{\Delta }g$ in Eq. (5), makes the compensation mechanism regional and $\mathrm{\Delta }g$ acts as a smoother or regularisation factor of the compensation [7].

Two factors are important for modelling the compensation depth, so-called the Mohorovic discontinuity (Moho), a) the mean compensation depth (${\tilde{D}}_0$) and b) the density contrast ($\mathrm{\Delta }\rho$) between the crust and upper-mantle. If either of $\mathrm{\Delta }\rho$ or ${\tilde{D}}_0$ is known the other one can be estimated from the model. The variation of Moho depth around ${\tilde{D}}_0$ can be determined by; see [7]:

$$\mathrm{\Delta }\tilde{D}=\frac{1}{4\mathit{\pi GR}\mathrm{\Delta }\rho }\sum _{\begin{array}{l}n=0\\ n\ne 1\end{array}}^{\infty }\frac{2n+1}{n-1}{\beta }_n^{\ast }{\mathrm{\Gamma }}_n\sum _{m=-n}^n\left(\mathrm{\Delta }{g}_{\mathit{nm}}^{\text{TBSCI}}-\mathrm{\Delta }{g}_{\mathit{nm}}\right)Y_{\mathit{nm}}\left(\theta ,\lambda \right)$$E6

where G is the Newtonian gravitational constant and

$${\beta }_n^{\ast }=\left\{\begin{array}{ll}1& \text{over oceans}\\ {\left(1+\left(n+2\right)\frac{{\tilde{D}}_0}{2R}\right)}^{-1}& \text{over continents}\end{array}\right.$$E7

$${\mathrm{\Gamma }}_n={\left(\frac{R}{R-{\tilde{D}}_0}\right)}^{n+2}.$$E8

The factor ${\mathrm{\Gamma }}_n$ is a signal amplifier, and when n increases this factor grows. For large values of ${\tilde{D}}_0$, this amplification starts from lower frequencies, therefore, the series in Eq. (6) should be truncated at lower degrees [7].

$\mathrm{\Delta }\rho$ can also be determined from $\mathrm{\Delta }\tilde{D}$, if its available or even the product $\mathrm{\Delta }\rho$ $\mathrm{\Delta }\tilde{D}$; e.g. see [7] in which the GOCE data are constrained to seismic data for determination of $\mathrm{\Delta }\rho$ $\mathrm{\Delta }\tilde{D}$.

CRUST1.0 [8] is a global model having information about the thicknesses and densities of sediments, crustal crystalline, topographic heights and bathymetric depths, and the Moho depths with a spatial resolution of 1^° × 1^°. This means that $\mathrm{\Delta }{g}_{\mathit{nm}}^{\text{TBSCI}}$ can be generated from CRUST1.0. In addition, numerous EGMs have been provided, which applicable for computing $\mathrm{\Delta }{g}_{\mathit{nm}}$. Figure 1a shows the Moho flexure/variation computed based on Eq. (6) the CRUST1.0 model, and EGM08 [9] limited to degree and order 180, corresponding to the resolution 1^° × 1^°. Figure 1b showed the contribution of $\mathrm{\Delta }g$ ranging from −15 to 15 km to the estimated Moho depth.

Figure 1.

(a) Global Moho flexure, and (b) the contribution of the gravity data to the Moho flexure.

3.1.2 Elastic thickness and rigidity

In flexural isostasy [10] the lithospheric is considered as an elastic shell, being flexes under loads. This shell bends based on its own mechanical properties and pressure of the loads. Elastic thickness (Te) is one of the properties of this shell. Admittance and coherence analyses between the topography and gravity anomalies; see [11] are known methods for estimating this elastic thickness. By combining the gravimetric and flexure isostasy models the elastic thickness or rigidity of the lithosphere can be estimated as well [7]. The main assumption of this approach is that the Moho variations derived from the gravimetric and flexural isostasy theories are equal. Therefore, the elastic thickness is estimated such a way that the Moho variation estimated from the gravimetric isostasy becomes closer to that from the flexural isostasy, by this assumption the resulted $\mathrm{\Delta }g$ from the lithospheric properties will be [7]:

$$\mathrm{\Delta }g=\mathrm{\Delta }{g}^{\text{TBSCI}}-4\mathit{\pi GR}\mathrm{\Delta }\rho \sum _{n=0}^{\infty }{\mathrm{\Gamma }}_n^{-1}\frac{n-1}{2n+1}{\beta }_n^{\ast }{C}_n^{-1}\sum _{m=-n}^n{\overline{K}}_{\mathit{nm}}{Y}_{\mathit{nm}}\left(\theta ,\lambda \right)$$E9

with

$${\overline{K}}_{\mathit{nm}}={\left(\overline{\rho }d\right)}_{\mathit{nm}}+{\left({\rho }^{\mathrm{S}}{d}_{\mathrm{S}}\right)}_{\mathit{nm}}+{\left({\rho }^{\mathrm{C}}{d}_{\mathbf{C}}\right)}_{\mathit{nm}}+{\left({\rho }^{\mathrm{I}}{d}_{\mathrm{I}}\right)}_{\mathit{nm}}$$E10

where $\overline{\rho }$ is the density of the topographic masses when the computation point is in continents and the density contrast between the water and topographic masses when it is over ocean, d stands for the topography height or bathymetric depth based on the position of the computation point. ${\rho }^{\mathrm{S}}$ and $d_{\mathrm{S}}$ are, respectively the density and the thickness of sediment layers, ${\rho }^{\mathrm{C}}$ and $d_{\mathbf{C}}$ the corresponding one for crustal crystalline, and ${\rho }^{\mathrm{I}}$ and $d_{\mathrm{I}}$ those of the ice. ${\left(•\right)}_{\mathit{nm}}$ means the spherical harmonic coefficients. $C_n$ is the compensation degree, which is derived from the flexure isostasy model:

$$C_n=\frac{n^2{\left(n+1\right)}^2}{R^{4}g}\tilde{\mathrm{\Theta }}+\mathrm{\Delta }\rho \text{and} \tilde{\mathrm{\Theta }}=\left\{\begin{array}{l}{D}^{\mathrm{Rig}}\text{if flexural rigidity is desired}\\ \frac{{\mathit{ET}}_e^3}{12\left(1-{\nu }^2\right)}\text{if elastic thickness is desired}\end{array}\right.$$E11

g is the gravity attraction, E stands for the Young modulus and v the Poisson ratio. In fact, $C_n$ carries the mechanical information of lithosphere including the elastic thickness.

The gravity anomaly on the left-hand side of Eq. (9), is generated from the lithosphere’s mass and density structures excluding the signal from sub-lithosphere. By comparison of this gravity anomaly and the observed ones excluding the lower degrees, coming from sub-lithosphere say to degree 15 [12] elastic thickness is determined in a trial and error process.

Figure 2 is the map of elastic thickness determined from GOCE gradiometric data over Africa in [13] the same procedure as explained for Eq. (9). The large elastic thickness over the tectonic border in the ocean is not realistic.

Figure 2.

Elastic thickness determined from GOCE data over Africa [13].

3.1.3 Bathymetry

Determining the ocean depths using gravity data is an old subject. Over offshore areas, hydrographic surveying methods are applicable by boats and Echo-sounders, known as traditional methods, modernised today by being equipped by GNSS technologies. However, they are costly and not practicable over oceans. Satellite altimetry data cover oceans sufficiently well and bathymetry can be done with acceptable precision, but the shortcoming is the low quality of them over shallow water. In this section, the focus will be on application of gravimetry over oceans for bathymetry purpose, based on isostasy. The theory and mathematical developments are available in [7], but they are not applied so far. Then the strengths and weaknesses of the method is still unknown.

Satellite altimetry data are the distance between the satellite and the sea surface, which is not fully-coincidence to the geoid. The departure of the sea surface from geoid is called sea surface topography. For determining the geoid from satellite altimetry, the sea surface topography should be known; and for determining the sea surface topography, the geoid is needed. The satellite gravimetry data or gravity models can be used without any involvement with the sea surface topography, but they have low resolutions. If the average depth of ocean $d_0$ is available, the variations of the seafloor topography around it will be [7]:

$$\mathrm{\Delta }d=-\frac{1}{4\mathit{\pi G}\overline{\rho }}\sum _{n=0}^{\infty }\frac{2n+1}{B_n}\sum _{m=-n}^n\left(t_{\mathit{nm}}-v_{\mathit{nm}}^{\mathrm{B}/\mathrm{Iso}}\right)Y_{\mathit{nm}}\left(\theta ,\lambda \right)-\left(\frac{R^{2}A{\delta }_{n0}}{B_0}\right)$$E12

where ${\delta }_{n0}$ stands for the Kronecker delta and

$$v_{\mathit{nm}}^{\mathrm{B}/\mathrm{Iso}}=-\left(v_{\mathit{nm}}^{\mathrm{S}}+v_{\mathit{nm}}^{\mathrm{C}}\right)+4\mathit{\pi GR}\mathrm{\Delta }\rho {\mathrm{\Gamma }}_n^{-1}{\beta }_n^{\ast }{C}_n^{-1}\left({\left({\rho }^{\mathrm{S}}{d}_{\mathrm{S}}\right)}_{\mathit{nm}}+{\left({\rho }^{\mathrm{C}}{d}_{\mathbf{C}}\right)}_{\mathit{nm}}\right)$$E13

$$A=\frac{1}{3R}\left(1-{\left(1-\frac{d_0}{R}\right)}^3\right)-\mathrm{\Delta }\rho \frac{1}{R}{\mathrm{\Gamma }}_0^{-1}{\beta }_0^{\ast }{C}_0^{-1}{d}_0$$E14

$$B_n=R{\left(\frac{R-d_0}{R}\right)}^{n+2}+R\mathrm{\Delta }\rho {\mathrm{\Gamma }}_n^{-1}{\beta }_n^{\ast }{C}_n^{-1}$$E15

$v_{\mathit{nm}}^{\mathrm{S}}$ and $v_{\mathit{nm}}^{\mathrm{C}}$ are gravitational potential of the sediment and crustal crystalline masses. Eq. (13) means the compensated gravitational potentials of sediment and crustal crystalline by the flexure isostasy. A and B_n are the contribution of the mean depth and its flexural compensation.

The important factor in bathymetry using this method is the elastic thickness of the lithosphere over oceans, which can be independently determined with a proper approximation from the age of the oceanic lithosphere by [14]:

$$T_{\mathrm{e}}=2.7\sqrt{t}$$E16

where t is the age of oceanic lithosphere in Ma.

3.1.4 Ice thickness

Determination the thickness of continental ice and its changes is important these days because of global warming. The continental ice is melted and water flow enters oceans and causes the sea level to rise. This is an issue which affect the Earth climate and may have risk of entering water to low land areas. Some satellite missions have been designed and developed for measuring the ice thickness using radar signals directly. This thickness can also be determined indirectly using gravimetry data. By assuming that the ice mass covers the surface of the Earth and it is part of the Earth’s solid topography, its thickness can be determined using the following spherical harmonic expansion [7]:

$$d_{\mathbf{I}}=\frac{1}{4\mathit{\pi GR}}\sum _{n=0}^{\infty }{\left(-\frac{\mathrm{\Delta }{\rho }^{\mathrm{I}}}{2n+1}+\mathrm{\Delta }{\mathit{\rho \rho }}^{\mathrm{I}}{\mathrm{\Gamma }}_n^{-1}{\beta }_n^{\ast }{C}_n^{-1}\right)}^{-1}\sum _{m=-n}^n\left(t_{\mathit{nm}}-v_{\mathit{nm}}^{\mathrm{TSC}/\mathrm{Iso}}\right)Y_{\mathit{nm}}\left(\theta ,\lambda \right)$$E17

where $\mathrm{\Delta }{\rho }^{\mathrm{I}}$ is the density contrast between the upper crust and ice, ${\rho }^{\mathrm{I}}$ stands for the density of ice, and

$$v_{\mathit{nm}}^{\mathrm{TSC}/\mathrm{Iso}}=-\left(v_{\mathit{nm}}^{\mathrm{T}}+v_{\mathit{nm}}^{\mathrm{S}}+v_{\mathit{nm}}^{\mathrm{C}}\right)+4\mathit{\pi Gr}\mathrm{\Delta }\rho {\mathrm{\Gamma }}_n^{-1}{\beta }_n^{\ast }{C}_n^{-1}\left({\left({\rho }^{\mathrm{T}}H\right)}_{\mathit{nm}}+{\left({\rho }^{\mathrm{S}}{d}_{\mathrm{S}}\right)}_{\mathit{nm}}+{\left({\rho }^{\mathrm{C}}{d}_{\mathbf{C}}\right)}_{\mathit{nm}}\right).$$E18

Note that Eq. (17) is based on the linear approximation of the involved binomial terms related to the topographic heights. Such an approximation is good as long as the heights are not large and the maximum degree of the expansion is not high. For example, for a height of 10 km and maximum degree 360, the relative error of this approximation will be 11%, for degree 180 is 4% and when the eight is 5 km for the maximum degree of 360 it will be 4% and less than 1% for 180. Since we have applied isostasy principle to obtain this equation, higher resolution than 180 is not needed, then approximation should be rather fine. One issue is the elastic thickness of lithosphere which is needed to determine the compensation degrees, which should be known from independent sources.

3.1.5 Sediment basement determination

Sediments are located at the surface of the upper-crust resulted from erosion during a long period of time. They are compacted by time meaning that their density will be high at their bottom and low at the surface. Therefore, the process of determining their thickness is not simple because the sediment density is an exponential function increasing by depth. In [7] some of the density contrast models have been presented and the gravitational potential of sediments have been modelled in spherical harmonics series. If we assume an average density for sediments the following approximate formula can be used to determine its thickness

$$d_{\mathrm{S}}=\frac{{\rho }^{\mathrm{S}}}{-4\mathit{\pi GR}}\sum _{n=0}^{\infty }{\left(-\frac{1}{2n+1}+\mathrm{\Delta }\rho {\mathrm{\Gamma }}_n^{-1}{\beta }_n^{\ast }{C}_n^{-1}\right)}^{-1}\sum _{m=-n}^n\left(t_{\mathit{nm}}-v_{\mathit{nm}}^{\mathrm{BC}/\mathrm{Iso}}\right)Y_{\mathit{nm}}\left(\theta ,\lambda \right)$$E19

where $\mathrm{\Delta }{\rho }^{\mathrm{I}}$ is the density contrast between the upper crust and ice, ${\rho }^{\mathrm{I}}$ stands for the density of ice, and

$$v_{\mathit{nm}}^{\mathrm{BC}/\mathrm{Iso}}=-\left(v_{\mathit{nm}}^{\mathrm{B}}+v_{\mathit{nm}}^{\mathrm{C}}\right)+4\mathit{\pi GR}\mathrm{\Delta }\rho {\mathrm{\Gamma }}_n^{-1}{\beta }_n^{\ast }{C}_n^{-1}\left({\left({\rho }^{\mathrm{B}}H\right)}_{\mathit{nm}}+{\left({\rho }^{\mathrm{C}}{d}_{\mathbf{C}}\right)}_{\mathit{nm}}\right).$$E20

One important parameter which should be known for sediment thickness determination using Eq. (19) is the elastic thickness of lithosphere, needed for computing the compensation degree. Over ocean there is a known relation between the lithospheric plate age and its elastic thickness see Eq. (16), but not over continents.

3.1.6 Runcorn’s theory and sub-lithospheric shear stresses

Mantle convection can also be studied by the long wavelength structure of the Earth gravity field. The Navier–Stokes equations of convection can be solved and simplified it in such a way that simple formula for shear stress at the base of lithosphere is obtained see [15]:

$$\left(\begin{array}{l}{\tau }_{\mathit{zx}}\\ {\tau }_{\mathit{zy}}\end{array}\right)=\frac{g}{4\mathit{\pi G}\left(R-D_{\text{Lith}}\right)}\sum _{n=2}^{\infty }\frac{2n+1}{n+1}{\left(\frac{R}{R-D_{\text{Lith}}}\right)}^{n+1}\sum _{m=-n}^n{t}_{\mathit{nm}}\left(\begin{array}{l}\frac{\partial Y_{\mathit{nm}}\left(\theta ,\lambda \right)}{\partial \theta }\\ \frac{\partial Y_{\mathit{nm}}\left(\theta ,\lambda \right)}{sin\theta \partial \lambda }\end{array}\right)$$E21

where ${\tau }_{\mathit{zx}}$ and ${\tau }_{\mathit{zy}}$ are the shear-stresses at the base of the lithosphere toward north and east, respectively. $D_{\text{Lith}}$ is the depth of boundary between lithosphere and mantle.

Eq. (21) is known as Runcorn’s formulae. He assumed that the mantle convection creates only the shear stresses at the base of the lithosphere. Most importantly, he assumed that:

1. the viscosity of mantle is constant.

2. the toroidal flow in the mantle is negligible.

3. the mantle is Newtonian.

Only by these assumptions the simple formula having a direct relation with the gravity data can be obtained. Many believe that the Runcorn simple solution is not realistic and successful, in spite of different efforts for justifying the applicability of this theory [16, 17, 18].

In Eq. (21), the maximum degree of expansion should not be infinity as the mantle convection contributes mainly in low frequencies of the gravity field. In [16] degrees between 13 to 25 are suggested to reduce the contributions from the core and lithosphere. However, in [12] the degrees below 15 are considered as contributions from sub-lithosphere.

Figure 3a and b show the map of the sub-lithospheric shear stresses ${\tau }_{\mathit{zx}}$ and ${\tau }_{\mathit{zy}}$, respectively, using Eq. (21) at the lithospheric depths of Conrad and Lithgow-Bertelloni model [19] over Iran. One issue in applying Eq. (21) is the choice of the maximum degree of expansion based on the lithospheric depth. When the base of the lithosphere is deeper, this degree should be lower and vice versa.

Figure 3.

(a) ${\tau }_{\mathit{zx}}$ and (b) ${\tau }_{\mathit{zy}}$,[MPa], [18] with permission from the publisher.

In [20] a better theory was developed for modelling the mantle convection using the displacement vectors of and tectonic movement. They also use the long wavelength portion of a geoid model in their solution, but the contribution of geoid is not very significant. This could be the reason that Runcorn has simplified the same mathematical models by ignoring the significant parameters and emphasising on the weakest one.

3.1.7 Stress propagation through the lithosphere from its base

By assuming that the lithosphere is an elastic shell, solution of the spherical boundary-value problem of elasticity can be applied for presenting the stress status inside the lithosphere. The stresses at base and top of the lithosphere is considered as the boundary-values. This subject was investigated by in [21] based on the solution of this problem by [22], and developed further by Fu and Huang [17]. The general solution of the boundary-value problem of elasticity is a displacement vector with four constants, which should be determined from the boundary-values. To do so this vector should be converted to general solutions for stress by the known relation between displacement and strain; and stain and stress [7]. The general solutions for stress include also those four constants. The Runcorn formula (21) can be considered as the low boundary-values of stresses, and it is assumed that the stress will disappear at the upper boundary, meaning that there is no stress. By selecting these boundary-values, a system of four equations is constructed and its solution will be those constants. By inserting these constants into the general solutions the stresses at a geocentric distance of r inside the lithospheric shell can be estimated. Also, they can be used in the general solutions of the strain and displacement to determine the strain tensor and displacement vector; for details see [7]. The elements of the stress tensor are:

$${\tau }_{\mathit{zz}}=\frac{1}{r}\sum _{n=2}^{\infty }\left(\tilde{\lambda }{K}_n^1+2\tilde{\mu }{K}_n^2\right)\sum _{m=-n}^n{t}_{\mathit{nm}}{Y}_{\mathit{nm}}\left(\theta ,\lambda \right)$$E22

$${\tau }_{\mathit{xx}}=\frac{1}{r}\sum _{n=2}^{\infty }\sum _{m=-n}^n{t}_{\mathit{nm}}\left\{\left(\tilde{\lambda }{K}_n^1+2\tilde{\mu }{K}_n^3\right)Y_{\mathit{nm}}\left(\theta ,\lambda \right)+2\tilde{\mu }{K}_n^5\frac{{\partial }^2{Y}_{\mathit{nm}}\left(\theta ,\lambda \right)}{\partial {\theta }^2}\right\},$$E23

$${\tau }_{\mathit{yy}}=\frac{1}{r}\sum _{n=2}^{\infty }\sum _{m=-n}^n{t}_{\mathit{nm}}\left\{\left(\tilde{\lambda }{K}_n^1+2\tilde{\mu }{K}_n^3\right)Y_{\mathit{nm}}\left(\theta ,\lambda \right)+2\tilde{\mu }{K}_n^5\left[\frac{1}{{sin}^2\theta }\frac{{\partial }^2{Y}_{\mathit{nm}}\left(\theta ,\lambda \right)}{\partial {\lambda }^2}+cot\theta \frac{\partial Y_{\mathit{nm}}\left(\theta ,\lambda \right)}{\partial \theta }\right]\right\}$$E24

$${\tau }_{\mathit{xz}}=\frac{\tilde{\mu }}{r}\sum _{n=2}^{\infty }\sum _{m=-n}^n{t}_{\mathit{nm}}{K}_n^4\frac{\partial Y_{\mathit{nm}}\left(\theta ,\lambda \right)}{\partial \theta }$$E25

$${\tau }_{\mathit{yz}}=\frac{\tilde{\mu }}{rsin\theta }\sum _{n=2}^{\infty }\sum _{m=-n}^n{t}_{\mathit{nm}}{K}_n^4\frac{\partial Y_{\mathit{nm}}\left(\theta ,\lambda \right)}{\partial \lambda },$$E26

$${\tau }_{\mathit{xy}}=\frac{\tilde{\mu }}{rsin\theta }\sum _{n=2}^{\infty }\sum _{m=-n}^n{t}_{\mathit{nm}}{K}_n^5\left(\frac{{\partial }^2{Y}_{\mathit{nm}}\left(\theta ,\lambda \right)}{\partial \theta \partial \lambda }-cot\theta \frac{\partial Y_{\mathit{nm}}\left(\theta ,\lambda \right)}{\partial \lambda }\right)$$E27

where $\tilde{\mu }$ and $\tilde{\lambda }$ are the elasticity coefficients, which can be determined from seismic data. For the coefficients $K_n^i$, i = 1,2,…,5, which are derived from the constants see [7].

Figure 4 shows the elements of the stress tensor for an earthquake at the depth of 10 km occurred in 25th of November 2018 with the magnitude of 6.3. The earthquake epicentre (34.361° N, 45.744° E) was located near the town Sar-e-Pol Zahab in West Iran close to the border with Iraq. The stress tensor has been determined by the Gravity field and Climate Experiment follow-on (GRACE-FO) [24] gravity model in October 2018.

Figure 4.

(a) ${\tau }_{\mathit{xx}}$ (b) ${\tau }_{\mathit{yy}}$ ₍c) ${\tau }_{\mathit{zz}}$ (d) ${\tau }_{\mathit{xz}}$ (e) ${\tau }_{\mathit{yz}}$ and (f) ${\tau }_{\mathit{xy}}$ [MPa], the star is the earthquake epicentre and the small dots are the distribution of seismic points, [18] with permission from the publisher.

3.2 Time-variable gravity field

The gravity field of the Earth varies in time due to different Geodynamical phenomena. This means that time-variable gravity data can be used for studying them. For example, the satellite missions Gravity Recovery and Climate Experiment (GRACE) [23] and GRACE-FO [24] have been designed and developed for determining temporal variations of the gravity field. Here, some of the phenomena causing such variations like Earthquakes, post-glacial rebound, ground water variations.

3.2.1 Earthquakes

Earthquakes are the result of energy extractions in the solid Earth based on different reasons. Whether an Earthquake is detectable by time-variable gravity data depends on the magnitude of the Earthquake, resolution and sensitivity of gravimetry. The GRACE and GRACE-FO monthly gravity models are applicable for studying the large Earthquakes. Geoid, gravity anomalies/disturbances, gravity gradients, stress, strain and even displacements can be computed from such gravity models before and after the Earthquakes. Changes of each quantities before and after the Earthquake provides information about the effect of the Earthquake on the gravity regime of area. However, one important issue is that the changes due to the non-Earthquake variations, like hydrological signals, should be removed or reduced from the gravity data/models prior to analysing any Earthquake.

Figure 5 shows the map of changes of gravity anomaly before and after the Zar-e-Pol Zahab Earthquake. Positive values are seen over the area and around the Earthquake’s epicentre illustrates by a circle, meaning increase of gravity, whilst the negative values are seen in eastern part of the area or gravity reduction. The black dot are earthquake points.

Figure 5.

Changes of the gravity anomalies before and after the Sar-e-pol Zahab earthquake on 25th November 2018, determined by the GRACE-FO gravity models in December 2018 and January 2019, [μGal]. Black dots are active seismic points and the start the earthquake epicentre. [18] With permission from publisher.

3.2.2 Determination of epicentre of shallow earthquakes

In [25] a connection between the maximum shear strain of the gravity strain tensor and epicentre of shallow Earthquakes were presented. A theory was presented in [26] for determining gravity strain tensor. In order to explain this type of strain, consider the geoid as a deforming surface. The changes of the geoid surface are regarded as a displacement field, and accordingly, this field is converted to a strain tensor, named gravity strain tensor with the following mathematical definition [26]:

$$\mathbf{S}=\frac{1}{2}\left({\mathbf{B}}^{-1}{\mathbf{bbB}}^{-1}-\mathbf{I}\right)$$E28

where.

$$\mathbf{B}=\left[\begin{array}{ccc}{V}_{\mathit{xx}}\left(t_1\right)& V_{\mathit{xy}}\left(t_1\right)& V_{\mathit{xz}}\left(t_1\right)\\ V_{\mathit{xy}}\left(t_1\right)& V_{\mathit{yy}}\left(t_1\right)& V_{\mathit{yz}}\left(t_1\right)\\ V_{\mathit{xz}}\left(t_1\right)& V_{\mathit{yz}}\left(t_1\right)& V_{\mathit{zz}}\left(t_1\right)\end{array}\right]\text{and}\mathbf{b}=\left[\begin{array}{ccc}{V}_{\mathit{xx}}\left(t_2\right)& V_{\mathit{xy}}\left(t_2\right)& V_{\mathit{xz}}\left(t_2\right)\\ V_{\mathit{xy}}\left(t_2\right)& V_{\mathit{yy}}\left(t_2\right)& V_{\mathit{yz}}\left(t_2\right)\\ V_{\mathit{xz}}\left(t_2\right)& V_{\mathit{yz}}\left(t_2\right)& V_{\mathit{zz}}\left(t_2\right)\end{array}\right].$$E29

In fact, B and b are the gravitational tensor in the local north-oriented frame at two epochs of t₁ before and t₂ after deformation.

Dilatation and maximum shear strain of the gravity strain tensor are, respectively

$$\mathrm{\Delta }={\lambda }_{\max }^{\mathrm{eig}}+{\lambda }_{\min }^{\mathrm{eig}}$$E30

$$\gamma ={\lambda }_{\max }^{\mathrm{eig}}-{\lambda }_{\min }^{\mathrm{eig}}$$E31

where ${\lambda }_{\max }^{\mathrm{eig}}$ and ${\lambda }_{\min }^{\mathrm{eig}}$ the largest and smaller eigenvalue of the gravity strain tensor.

The map of the maximum shear strain over the area experiencing a shallow Earthquake will show a high value at the Earthquake epicentre; see [25].

In order to represent an example about the application of this theory, the eastern Turkey Earthquake occurred on 2010-2103-08 at 7:41:41 UTC and depth of 10 km is considered. The position of the earthquake epicentre is 38.709°N and 40.051°E according to the United States Geological Survey (USGS); see Figure 6. In this figure, the map of the maximum shear strain determined from the GRACE monthly gravity models are over the area. The maximum shear strain have been computed from two years of gravity models before and after the Earthquake, and the yellow rectangle shows the approximate position of the Earthquake epicentre. Note that the colour of the circle was chosen for better visualisation of the epicentre reported by the USGS, and is not related to the colourbar present for the map.

Figure 6.

The position of the eastern Turkey earthquake epicentre detected by the gravity strain approach and USGS, [27].

Figure 7.

(a) Geoid trend during 15 years of GRACE mission, (b) land uplift model determined from geoid rate of change.

3.2.3 Post-glacial rebound and mantle viscosity

Mantle is a viscous medium and its viscosity creates an upward force at the base of the lithosphere against bending due to loads. The lithosphere would bend more if it was buoyant over a less viscous medium. In addition, the age of load is an important factor in the lithosphere thrusting into the mantle; older lithosphere is more inside the mantle than the younger one. Both of the lithospheric strength and the mantle’s viscosity keep the lithosphere in an isostatic equilibrium against loads pushing the lithosphere downwards. If these loads are removed, this balance is destructed and the mantle pushes the lithosphere upwards causing the land rise or uplift.

In the ice age period, huge ice masses existed at the surface of the lithosphere, and by the increase of the Earth’s temperature, they were melted and the lithospheric rebound began toward the isostatic equilibrium. This phenomena is called post-glacial rebound, or glacial isostatic adjustment, causing land uplift, which can be monitored by the temporal changes of gravity data. For example, if the geoid rate is determined using time-variable gravity models, the land uplift rate due to this rebound can be computed by [7, 28]:

$$\dot{h}\left(\theta ,\lambda \right)=\frac{\gamma }{4\mathit{\pi GR}}\sum _{n=0}^{\infty }\frac{2n+1}{{\kappa }_n^{″}}\sum _{m=-n}^n\mathrm{\Delta }{\dot{N}}_{\mathit{nm}}{Y}_{\mathit{nm}}\left(\theta ,\lambda \right)$$E32

where $\mathrm{\Delta }{\dot{N}}_{\mathit{nm}}$ the spherical harmonic coefficients of the geoid rate and

$${\kappa }_n^{″}={\rho }^{\mathrm{C}}+{\mathrm{\Gamma }}_n^{-1}\mathrm{\Delta }\rho .$$E33

The effect of hydrological signals should be removed from the time-variable gravity models prior to applying them for determining the geoid rate by a linear regression. Figure 7a is the map of this rate showing variation from −0.6 to 0.4 mm/yr., determined from the GRACE time-variable gravity models during the of the GRACE mission and after removing the hydrological signals using Global Land data Assimilation System (GLDAS) [29] model over Fennoscandia, which is experiencing the post-glacial rebound after the ice age. The geoid rate of change has been computed globally and after performing a spherical harmonic analysis its harmonics have been computed and inserted into Eq. (32) for estimating the land uplift rate; see Figure 7b. The uplift rate varies from −4 to 9 mm/yr. with the maximum around the centre of Gulf of Bohemia.

In [30] methods for determining the mantle viscosity were presented, but the geodetic approach was proposed in [31] and shown that the highest correlation between the land uplift data and geoid is achieved when the geoid is computed from degree 10 to 70. In [29, 32] used the spherical harmonic degrees to 23 instead of 70. In fact, degree 23 is obtained by performing a correlation analysis between the geoid derived to different maximum degrees and the land uplift model determined by the GNSS measurements. In [33] there is a discussion about some frequencies window of the geoid signal affected by the post-glacial rebound and later investigation in [33] it is shown that this frequency window is limited between degrees 10 to 23. If we accept this theory the viscosity of the upper mantle can be determined by [7]:

$$\tilde{\eta }=-\frac{{\gamma }^2{\rho }_m}{4\mathit{\pi G}\dot{h}\left(\theta ,\lambda \right)}\sum _{n=10}^{23}\frac{2n+1}{{\kappa }_n^{″}\left(2n+4+3/n\right)}\sum _{m=-n}^n{N}_{\mathit{nm}}{Y}_{\mathit{nm}}\left(\theta ,\lambda \right)$$E34

where ${\rho }_m$ is density of the upper mantle.

The mean viscosity of the upper mantle will be (5.0 ± 0.2) × 10²¹ Pa, and in the case of using Eq. (34), it will be (6.0 ± 0.3) × 10²¹ Pa over Fennoscandia [7].

3.2.4 Monitoring hydrological signals

Hydrological signals are the main surfaces of fast temporal changes of the gravity field. They come from ground water storage (GWS), snow water equivalent (SWE), solid moisture (SM) and Canopy (CAN). Different models have been presented these signals except for the GWS and the most known one is the GLDAS model [29] which has had good agreement with the temporal variations of the gravity field determined by GRACE. However, the GRACE models provide information about the total water, or equivalent water height, or a summations of SM, SWE, CAN and GWS. Therefore, if one of these hydrological signals is required, it can be determined from a combination of the GRACE and hydrological models; see [7, 34]. In the case where the GWS is needed, it can be computed by:

$$\delta h_{\mathit{nm}}^{\mathrm{GWS}}=\sum _{n=2}^{\infty }\sum _{m=-n}^n\left(\frac{1}{4\mathit{\text{πGRγ}}}\frac{2n+1}{1+k_n}\delta v_{\mathit{nm}}-\frac{1}{{\rho }_w}\left({\mathit{\delta \rho }}_{\mathit{nm}}^{\mathrm{SM}}+{\mathit{\delta \rho }}_{\mathit{nm}}^{\mathrm{SWE}}+{\mathit{\delta \rho }}_{\mathit{nm}}^{\mathrm{CAN}}\right)\right)Y_{\mathit{nm}}\left(\theta ,\lambda \right)$$E35

where ${\rho }_w$ is the density of water, $\delta v_{\mathit{nm}}$ is the changes of the gravitational potential, $k_n$ is the Love numbers, ${\mathit{\delta \rho }}_{\mathit{nm}}^{\mathrm{SM}}$, ${\mathit{\delta \rho }}_{\mathit{nm}}^{\mathrm{SWE}}$ and ${\mathit{\delta \rho }}_{\mathit{nm}}^{\mathrm{CAN}}$ are the spherical harmonic coefficients of the densities of SM, SWE and CAN, respectively.

Figure 8 is the global map of the GWS all over the globe computed by the GRACE gravity models during 15 years, 2002 to 2017. Note that the post-glacial rebound and earthquake signals have not be excluded in the computations.

Figure 8.

The global ground water storage (WGS) rate determined from 15 years of GRACE gravity models and GLDAS [27].

The largest GWS is seen over Hudson Bay in Canada, and the green land. Both of these places are known as active areas for post-glacial rebound. Reduction of GWS is seen in the Middle East and eastern Africa, and Western Australia and increase in Russia, western Africa, eastern Australia.

4. Concluding remarks

The goal of this chapter was to demonstrate applications of the Earth’s gravity field in Geodesy and Global Geophysics. In Geodesy, the main goal is to determine the physical shape of the Earth, or the geoid, and its importance in levelling and height systems were discussed. Philosophies behind three well-known methods of geoid determination, such as Remove-Compute-Restore, Stokes-Helmert and Least-squares modification of the Stokes formula with additive corrections, were discussed.

When the temporal variation of the gravity field is disregarded and the field is considered static, some geophysical subjects can be studies by them to understand the Earth’s interior such as the crustal structure, density contrast between the crust and mantle, sediment basement, ice thickness, and depth of ocean water determination. In addition, the sub-lithospheric stress induced by mantle convection and its propagation through the lithosphere can also be determined using gravity data.

By studying the temporal variations of the gravity field, Geodynamic phenomena can be studies. Post-glacial rebound, determining the land uplift rate and mantle viscosity, studying the earthquakes and their epicentre and also ground water storage mapping are the subject which can be studied by these variations.

One distinction between the application of the gravity field in Geodesy and global geophysics is the resolution of gravity field. The main purpose of Geodesy is to determine the shape of the Earth as precise and accurate as possible, and focus is on recovering the high frequencies of the gravity field, by combining satellite and terrestrial data with mathematical tools with the least approximations. However, in Geophysics due to lack our knowledge about the Earth’s interior structure and dynamic, different assumptions have be made and also the mathematical models are developed based on them. In addition high resolution gravity data do not play a significant role in global Geophysics.

Acknowledgments

The author is thankful to Farzam Fatolazadeh for generating maps of the maximum shear gravity strain and global ground water storage.