Gravity Disturbance
Gravity Disturbance¶
Gravity disturbances are the differences between the measured gravity and a reference (normal) gravity produced by an ellipsoid:
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
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\):
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()
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(),
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()
The gravity disturbances can be interpreted as the gravitational effect of every anomalous mass, i.e. that is not contained in the normal Earth.