*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]:

$$\sigma \begin{pmatrix} \mathfrak{z}\_{\text{xx}} \\ \mathfrak{z}\_{\text{xy}} \end{pmatrix} = \frac{\mathcal{G}}{4\pi \mathcal{G} (R - D\_{\text{Lth}})} \sum\_{n=2}^{\infty} \frac{2n + 1}{n + 1} \left( \frac{R}{R - D\_{\text{Lth}}} \right)^{n+1} \sum\_{m=-n}^{n} t\_{nm} \begin{pmatrix} \frac{\partial Y\_{nm}(\theta, \lambda)}{\partial \theta} \\ \frac{\partial Y\_{nm}(\theta, \lambda)}{\sin \theta \partial \lambda} \end{pmatrix} \tag{21}$$

where *τzx* and *τzy* are the shear-stresses at the base of the lithosphere toward north and east, respectively. *D*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:

a. the viscosity of mantle is constant.


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–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 *τzx* and *τ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.

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

*<sup>τ</sup>xy* <sup>¼</sup> *<sup>μ</sup>*<sup>~</sup> *r*sin *θ*

constants see [7].

*3.2.1 Earthquakes*

**Figure 4.**

**13**

seismic data. For the coefficients *K<sup>i</sup>*

[24] gravity model in October 2018.

**3.2 Time-variable gravity field**

X∞ *n*¼2

*DOI: http://dx.doi.org/10.5772/intechopen.97459*

X*n m*¼�*n*

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

*tnmK*<sup>5</sup> *n ∂*2

*Ynm*ð Þ *θ*, *λ*

where *μ*~ and ~*λ* are the elasticity coefficients, which can be determined from

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)

**Figure 4** shows the elements of the stress tensor for an earthquake at the depth

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 varia-

Earthquakes are the result of energy extractions in the solid Earth based on different reasons. Whether an Earthquake is detectable by time-variable gravity

*(a) τxx (b) τyy (c) τzz (d) τxz (e) τyz and (f) τxy [MPa], the star is the earthquake epicentre and the small dots*

*are the distribution of seismic points, [18] with permission from the publisher.*

tions like Earthquakes, post-glacial rebound, ground water variations.

*<sup>∂</sup>θ∂<sup>λ</sup>* � cot *<sup>θ</sup>*

*<sup>∂</sup>Ynm*ð Þ *<sup>θ</sup>*, *<sup>λ</sup> ∂λ* � � (27)

*<sup>n</sup>*, *i* = 1,2, … ,5, which are derived from the

**Figure 3.** *(a) τzx and (b) τzy,[MPa], [18] with permission from the publisher.*

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\_{\text{xx}} = \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\_{nm} Y\_{nm}(\theta, \lambda) \tag{22}$$

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

$$\tau\_{\eta\eta} = \frac{1}{r} \sum\_{n=2}^{\text{so}} \sum\_{m=-n}^{n} t\_{nm} \left\{ \left( \bar{\lambda} K\_n^1 + 2 \bar{\mu} K\_n^2 \right) Y\_{nm}(\theta, \lambda) + 2 \bar{\mu} K\_n^5 \left[ \frac{1}{\sin^2 \theta} \frac{\partial^2 Y\_{nm}(\theta, \lambda)}{\partial \lambda^2} + \cot \theta \frac{\partial Y\_{nm}(\theta, \lambda)}{\partial \theta} \right] \right\} \tag{24}$$

$$\tau\_{\text{xx}} = \frac{\tilde{\mu}}{r} \sum\_{n=2}^{\infty} \sum\_{m=-n}^{n} t\_{nm} K\_n^4 \frac{\partial Y\_{nm}(\theta, \lambda)}{\partial \theta} \tag{25}$$

$$\tau\_{yx} = \frac{\tilde{\mu}}{r \sin \theta} \sum\_{n=2}^{\infty} \sum\_{m=-n}^{n} t\_{nm} K\_n^4 \frac{\partial Y\_{nm}(\theta, \lambda)}{\partial \lambda},\tag{26}$$

*The Earth's Gravity Field Role in Geodesy and Large-Scale Geophysics DOI: http://dx.doi.org/10.5772/intechopen.97459*

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

where *μ*~ and ~*λ* are the elasticity coefficients, which can be determined from seismic data. For the coefficients *K<sup>i</sup> <sup>n</sup>*, *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.

## **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

**Figure 4.**

*(a) τxx (b) τyy (c) τzz (d) τxz (e) τyz and (f) τxy [MPa], the star is the earthquake epicentre and the small dots are the distribution of seismic points, [18] with permission from the publisher.*

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 Zare-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.

### *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{b} \mathbf{b} \mathbf{B}^{-1} - \mathbf{I} \right) \tag{28}$$

where.

2 6 4

*Vxx*ð Þ *t*<sup>1</sup> *Vxy*ð Þ *t*<sup>1</sup> *Vxz*ð Þ *t*<sup>1</sup> *Vxy*ð Þ *t*<sup>1</sup> *Vyy*ð Þ *t*<sup>1</sup> *Vyz*ð Þ *t*<sup>1</sup> *Vxz*ð Þ *t*<sup>1</sup> *Vyz*ð Þ *t*<sup>1</sup> *Vzz*ð Þ *t*<sup>1</sup>

*DOI: http://dx.doi.org/10.5772/intechopen.97459*

two epochs of *t*<sup>1</sup> before and *t*<sup>2</sup> after deformation.

3 7 5

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

and **b** ¼

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

<sup>Δ</sup> <sup>¼</sup> *<sup>λ</sup>*eig

*<sup>γ</sup>* <sup>¼</sup> *<sup>λ</sup>*eig

Earthquake will show a high value at the Earthquake epicentre; see [25].

by the USGS, and is not related to the colourbar present for the map.

*The position of the eastern Turkey earthquake epicentre detected by the gravity strain approach and*

In fact, **B** and **b** are the gravitational tensor in the local north-oriented frame at

max þ *λ*

max � *λ*

The map of the maximum shear strain over the area experiencing a shallow

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

eig

eig

min the largest and smaller eigenvalue of the gravity strain

2 6 4

*Vxx*ð Þ *t*<sup>2</sup> *Vxy*ð Þ *t*<sup>2</sup> *Vxz*ð Þ *t*<sup>2</sup> *Vxy*ð Þ *t*<sup>2</sup> *Vyy*ð Þ *t*<sup>2</sup> *Vyz*ð Þ *t*<sup>2</sup> *Vxz*ð Þ *t*<sup>2</sup> *Vyz*ð Þ *t*<sup>2</sup> *Vzz*ð Þ *t*<sup>2</sup>

min (30)

min (31)

3 7 <sup>5</sup>*:* (29)

**B** ¼

respectively

where *λ*eig

tensor.

**Figure 6.**

**15**

*USGS, [27].*

max and *λ*

eig

#### **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.*

*The Earth's Gravity Field Role in Geodesy and Large-Scale Geophysics DOI: http://dx.doi.org/10.5772/intechopen.97459*

where.

$$\mathbf{B} = \begin{bmatrix} V\_{xx}(t\_1) & V\_{xy}(t\_1) & V\_{xx}(t\_1) \\ V\_{xy}(t\_1) & V\_{yy}(t\_1) & V\_{yx}(t\_1) \\ V\_{xx}(t\_1) & V\_{yx}(t\_1) & V\_{xx}(t\_1) \end{bmatrix} \text{and } \mathbf{b} = \begin{bmatrix} V\_{xx}(t\_2) & V\_{xy}(t\_2) & V\_{xx}(t\_2) \\ V\_{xy}(t\_2) & V\_{yy}(t\_2) & V\_{yx}(t\_2) \\ V\_{xx}(t\_2) & V\_{yx}(t\_2) & V\_{xx}(t\_2) \end{bmatrix}. \tag{29}$$

In fact, **B** and **b** are the gravitational tensor in the local north-oriented frame at two epochs of *t*<sup>1</sup> before and *t*<sup>2</sup> after deformation.

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

$$
\Delta = \lambda\_{\text{max}}^{\text{eig}} + \lambda\_{\text{min}}^{\text{eig}} \tag{30}
$$

$$
\lambda = \lambda\_{\text{max}}^{\text{eig}} - \lambda\_{\text{min}}^{\text{eig}} \tag{31}
$$

where *λ*eig max and *λ* eig min 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].*

#### *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}(\theta,\lambda) = \frac{\gamma}{4\pi GR} \sum\_{n=0}^{\infty} \frac{2n+1}{\kappa\_n''} \sum\_{m=-n}^{n} \Delta \dot{N}\_{nm} Y\_{nm}(\theta,\lambda) \tag{32}$$

where Δ*N*\_ *nm* the spherical harmonic coefficients of the geoid rate and

$$
\kappa\_n'' = \rho^\mathbb{C} + \Gamma\_n^{-1} \Delta \rho. \tag{33}
$$

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

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]:

2*n* þ 1

X*n m*¼�*n*

*NnmYnm*ð Þ *θ*, *λ* (34)

*<sup>n</sup>*ð Þ 2*n* þ 4 þ 3*=n*

case of using Eq. (34), it will be (6.0 � 0.3) � <sup>10</sup><sup>21</sup> Pa over Fennoscandia [7].

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

*The global ground water storage (WGS) rate determined from 15 years of GRACE gravity models and*

The mean viscosity of the upper mantle will be (5.0 � 0.2) � <sup>10</sup><sup>21</sup> Pa, and in the

maximum around the centre of Gulf of Bohemia.

*DOI: http://dx.doi.org/10.5772/intechopen.97459*

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

<sup>~</sup>*<sup>η</sup>* ¼ � *<sup>γ</sup>*<sup>2</sup>*ρ<sup>m</sup>* 4*πG*\_

*3.2.4 Monitoring hydrological signals*

**Figure 8.**

**17**

*GLDAS [27].*

*h*ð Þ *θ*, *λ*

where *ρ<sup>m</sup>* is density of the upper mantle.

X 23

*n*¼10

*κ*00

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

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

*The Earth's Gravity Field Role in Geodesy and Large-Scale Geophysics DOI: http://dx.doi.org/10.5772/intechopen.97459*

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]:

$$\ddot{\eta} = -\frac{\eta^2 \rho\_m}{4\pi G \dot{h}(\theta, \lambda)} \sum\_{n=10}^{23} \frac{2n+1}{\kappa\_n''(2n+4+3/n)} \sum\_{m=-n}^{n} N\_{nm} Y\_{nm}(\theta, \lambda) \tag{34}$$

where *ρ<sup>m</sup>* is density of the upper mantle.

The mean viscosity of the upper mantle will be (5.0 � 0.2) � <sup>10</sup><sup>21</sup> Pa, and in the case of using Eq. (34), it will be (6.0 � 0.3) � <sup>10</sup><sup>21</sup> 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

#### **Figure 8.**

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

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\_{nm}^{\text{GWS}} = \sum\_{n=2}^{\infty} \sum\_{m=-n}^{n} \left( \frac{1}{4\pi GR\gamma} \frac{2n+1}{1+k\_n} \delta v\_{nm} - \frac{1}{\rho\_w} \left( \delta \rho\_{nm}^{\text{SM}} + \delta \rho\_{nm}^{\text{SWE}} + \delta \rho\_{nm}^{\text{CAN}} \right) \right) Y\_{nm}(\theta, \lambda) \tag{35}$$

where *ρ<sup>w</sup>* is the density of water, *δvnm* is the changes of the gravitational potential, *kn* is the Love numbers, *δρ*SM *nm*, *δρ*SWE *nm* and *δρ*CAN *nm* 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.

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.
