# 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 [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)
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()


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()


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(),