.. _gravity_disturbance:

Gravity Disturbance
===================

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

.. math::

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

Where :math:`\delta g` is the **gravity disturbance**, :math:`g` the **measured
gravity**, :math:`\gamma` is the **normal gravity** and :math:`\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** :math:`g` is defined as the module of the gravity
    acceleration, i.e. the module of the gradient of Earth's gravity potential
    :math:`W`, which includes both the gravitational potential :math:`V` (due
    to the gravitational attraction of the masses that composes the Earth) and
    the centrifugal potential :math:`\Phi` due to Earth's rotation

    .. math::

       W = V + \Phi

.. tip::

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

    .. math::

       U = V_\text{ell} + \Phi

.. seealso::

   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
:meth:`boule.Ellipsoid.normal_gravity` method from :mod:`boule` and then use it
to compute gravity disturbances.

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

.. jupyter-execute::

   import ensaio
   import xarray as xr

   fname = ensaio.fetch_earth_gravity(version=1)
   gravity = xr.load_dataarray(fname)
   print(gravity)

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:

.. jupyter-execute::

   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 :mod:`boule` and use the
:meth:`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.

.. jupyter-execute::

   import boule as bl

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

And plot it:

.. jupyter-execute::

   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:

.. jupyter-execute::

   gravity_disturbance = gravity - normal_gravity
   print(gravity_disturbance)

And plot it:

.. jupyter-execute::

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