.. _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 pygmt

   fig = pygmt.Figure()
   fig.grdimage(
       gravity,
       projection="W20c",
       cmap="viridis",
       shading="+a45+nt0.2",
   )
   fig.basemap(frame=["af", "WEsn"])
   fig.colorbar(
       position="JCB+w10c",
       frame=["af", 'y+l"mGal"', 'x+l"observed gravity"'],
   )
   fig.coast(shorelines=True, resolution="c", area_thresh=1e4)
   fig.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::

   fig = pygmt.Figure()
   fig.grdimage(
       normal_gravity,
       projection="W20c",
       cmap="viridis",
       shading="+a45+nt0.2",
   )
   fig.basemap(frame=["af", "WEsn"])
   fig.colorbar(
       position="JCB+w10c",
       frame=["af", 'y+l"mGal"', 'x+l"normal gravity"'],
   )
   fig.coast(shorelines=True, resolution="c", area_thresh=1e4)
   fig.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)

   fig = pygmt.Figure()
   pygmt.makecpt(series=[-maxabs, maxabs], cmap="polar+h")
   fig.grdimage(
       gravity_disturbance,
       projection="W20c",
       cmap=True,
       shading="+a45+nt0.2",
   )
   fig.basemap(frame=["af", "WEsn"])
   fig.colorbar(
       position="JCB+w10c",
       frame=["af", 'y+l"mGal"', 'x+l"gravity disturbance"'],
   )
   fig.coast(shorelines=True, resolution="c", area_thresh=1e4)
   fig.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*.