Gravity Disturbance

Gravity Disturbance

Gravity disturbances are the differences between the measured gravity and a reference (normal) gravity produced by an ellipsoid:

\[\delta g(\mathbf{p}) = g(\mathbf{p}) - \gamma(\mathbf{p})\]

Where \(\delta g\) is the gravity disturbance, \(g\) the measured gravity, \(\gamma\) is the normal gravity and \(\mathbf{p}\) is the observation point where the gravity measurement has been carried out. The disturbances are what allows geoscientists to infer the internal structure of the Earth.

Tip

The measured gravity \(g\) is defined as the module of the gravity acceleration, i.e. the module of the gradient of Earth’s gravity potential \(W\), which includes both the gravitational potential \(V\) (due to the gravitational attraction of the masses that composes the Earth) and the centrifugal potential \(\Phi\) due to Earth’s rotation

\[W = V + \Phi\]

Tip

The normal gravity \(\gamma\) is defined as the gradient of the potential gravity field generated by the reference ellipsoid \(U\) composed by the sum of the gravitational field \(V_\text{ell}\) and the centrifugal potential \(\Phi\):

\[U = V_\text{ell} + \Phi\]

See also

See [Hofmann-WellenhofMoritz2006] for detailed explainations of Earth gravity and the definition of its gravity and gravitational potentials.

We can compute the normal gravity generated by any ellipsoid through the boule.Ellipsoid.normal_gravity method from boule and then use it to compute gravity disturbances.

Lets start by loading a sample gravity dataset for the whole Earth:

import ensaio
import xarray as xr

fname = ensaio.fetch_earth_gravity(version=1)
gravity = xr.load_dataarray(fname)
print(gravity)
<xarray.DataArray 'gravity' (latitude: 1081, longitude: 2161)>
array([[980106.5 , 980106.5 , 980106.5 , ..., 980106.5 , 980106.5 ,
        980106.5 ],
       [980108.25, 980108.25, 980108.25, ..., 980108.25, 980108.25,
        980108.25],
       [980108.8 , 980108.8 , 980108.8 , ..., 980108.75, 980108.75,
        980108.8 ],
       ...,
       [980153.8 , 980153.75, 980153.6 , ..., 980153.94, 980153.8 ,
        980153.8 ],
       [980160.44, 980160.44, 980160.44, ..., 980160.44, 980160.44,
        980160.44],
       [980157.5 , 980157.5 , 980157.5 , ..., 980157.5 , 980157.5 ,
        980157.5 ]], dtype=float32)
Coordinates:
  * longitude  (longitude) float64 -180.0 -179.8 -179.7 ... 179.7 179.8 180.0
  * latitude   (latitude) float64 -90.0 -89.83 -89.67 -89.5 ... 89.67 89.83 90.0
    height     (latitude, longitude) float32 1e+04 1e+04 1e+04 ... 1e+04 1e+04
Attributes:
    Conventions:     CF-1.8
    title:           Gravity acceleration (EIGEN-6C4) at a constant geometric...
    crs:             WGS84
    source:          Generated from the EIGEN-6C4 model by the ICGEM Calculat...
    license:         Creative Commons Attribution 4.0 International Licence
    references:      https://doi.org/10.5880/icgem.2015.1
    long_name:       gravity acceleration
    description:     magnitude of the gravity acceleration vector (gravitatio...
    units:           mGal
    actual_range:    [974748.6 980201.9]
    icgem_metadata:  generating_institute: gfz-potsdam\ngenerating_date: 2021...

These observations are located on a regular grid on geodetic coordinates at the same height of 10 km above the reference ellipsoid. Lets plot it:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Robinson())
pc = gravity.plot.pcolormesh(
    ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False, cmap="viridis"
)
plt.colorbar(
    pc, label="mGal", orientation="horizontal", aspect=50, pad=0.01, shrink=0.5
)
ax.set_title("Gravity of the Earth")
ax.coastlines()
plt.show()
../_images/gravity_disturbance_1_0.png

We can then get the WGS84 ellipsoid defined in boule and use the boule.Ellipsoid.normal_gravity to compute the normal gravity (the gravity acceleration generated by the ellipsoid) on every observation point. This method implements the closed-form formula of [LiGotze2001], which calculates the normal gravity at any latitude and (geometric) height through an analytic solution.

import boule as bl

ellipsoid = bl.WGS84
normal_gravity = ellipsoid.normal_gravity(gravity.latitude, gravity.height)

And plot it:

plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Robinson())
pc = ax.pcolormesh(
    gravity.longitude,
    gravity.latitude,
    normal_gravity,
    transform=ccrs.PlateCarree(),
    cmap="viridis"
)
plt.colorbar(
    pc, label="mGal", orientation="horizontal", aspect=50, pad=0.01, shrink=0.5
)
ax.set_title("Normal gravity of the Earth")
ax.coastlines()
plt.show()
../_images/gravity_disturbance_3_0.png

Now we can compute the gravity disturbance:

gravity_disturbance = gravity - normal_gravity
print(gravity_disturbance)
<xarray.DataArray (latitude: 1081, longitude: 2161)>
array([[-35.83509235, -35.83509235, -35.83509235, ..., -35.83509235,
        -35.83509235, -35.83509235],
       [-34.04097894, -34.04097894, -34.04097894, ..., -34.04097894,
        -34.04097894, -34.04097894],
       [-33.34614029, -33.34614029, -33.34614029, ..., -33.40864029,
        -33.40864029, -33.34614029],
       ...,
       [ 11.65385965,  11.59135965,  11.46635965, ...,  11.77885965,
         11.65385965,  11.65385965],
       [ 18.14652106,  18.14652106,  18.14652106, ...,  18.14652106,
         18.14652106,  18.14652106],
       [ 15.16490765,  15.16490765,  15.16490765, ...,  15.16490765,
         15.16490765,  15.16490765]])
Coordinates:
  * longitude  (longitude) float64 -180.0 -179.8 -179.7 ... 179.7 179.8 180.0
  * latitude   (latitude) float64 -90.0 -89.83 -89.67 -89.5 ... 89.67 89.83 90.0
    height     (latitude, longitude) float32 1e+04 1e+04 1e+04 ... 1e+04 1e+04

And plot it:

import verde as vd

maxabs = vd.maxabs(gravity_disturbance)

plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Robinson())
pc = gravity_disturbance.plot.pcolormesh(
    ax=ax,
    transform=ccrs.PlateCarree(),
    add_colorbar=False,
    cmap="seismic",
    vmin=-maxabs,
    vmax=maxabs,
)
plt.colorbar(
    pc, label="mGal", orientation="horizontal", aspect=50, pad=0.01, shrink=0.5
)
ax.set_title("Gravity disturbance of the Earth")
ax.coastlines()
plt.show()
../_images/gravity_disturbance_5_0.png

The gravity disturbances can be interpreted as the gravitational effect of every anomalous mass, i.e. that is not contained in the normal Earth.