Gravity disturbances#

One of the main uses of Boule in geophysics is the calculation of gravity disturbances. Our closed-form implementation of normal gravity allows accurate gravity disturbance calculation at any height without the errors associated with approximate free-air corrections.

As an example, lets calculate a gravity disturbance grid for the Earth using a global gravity from the EIGEN-6C4 spherical harmonic model at a height of 10 km.

First off, we import the required packages to do this:

import boule as bl
import ensaio
import xarray as xr  # For loading the grid data
import pygmt         # For plotting nice maps

Next, we can download and cache the gravity grid using ensaio and load it with xarray:

fname = ensaio.fetch_earth_gravity(version=1)
observed_gravity = xr.load_dataarray(fname)
observed_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...

Let’s plot this data on a map using pygmt to see what it looks like:

fig = pygmt.Figure()
fig.grdimage(
    observed_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()
../_images/gravity_disturbance_2_0.png

Now we can calculate the WGS84 normal gravity at the same points as the observed gravity using boule.Ellipsoid.normal_gravity:

normal_gravity = bl.WGS84.normal_gravity(
    observed_gravity.latitude, observed_gravity.height,
)
normal_gravity

<xarray.DataArray (latitude: 1081, longitude: 2161)>
array([[980142.33509235, 980142.33509235, 980142.33509235, ...,
        980142.33509235, 980142.33509235, 980142.33509235],
       [980142.29097894, 980142.29097894, 980142.29097894, ...,
        980142.29097894, 980142.29097894, 980142.29097894],
       [980142.15864029, 980142.15864029, 980142.15864029, ...,
        980142.15864029, 980142.15864029, 980142.15864029],
       ...,
       [980142.15864035, 980142.15864035, 980142.15864035, ...,
        980142.15864035, 980142.15864035, 980142.15864035],
       [980142.29097894, 980142.29097894, 980142.29097894, ...,
        980142.29097894, 980142.29097894, 980142.29097894],
       [980142.33509235, 980142.33509235, 980142.33509235, ...,
        980142.33509235, 980142.33509235, 980142.33509235]])
Coordinates:
  * latitude   (latitude) float64 -90.0 -89.83 -89.67 -89.5 ... 89.67 89.83 90.0
  * longitude  (longitude) float64 -180.0 -179.8 -179.7 ... 179.7 179.8 180.0
    height     (latitude, longitude) float32 1e+04 1e+04 1e+04 ... 1e+04 1e+04
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()
../_images/gravity_disturbance_4_2.png

We can arrive at a grid of gravity disturbances by subtracting normal gravity from the data grid:

disturbance = observed_gravity - normal_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

Finally, we can display the disturbance in a nice global map:

fig = pygmt.Figure()
fig.grdimage(
    disturbance, projection="W20c", cmap="polar+h", 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()
../_images/gravity_disturbance_6_2.png

Gravity disturbances can be interpreted geophysically as the gravitational effect of anomalous masses, i.e. those that are not present in the normal (ellipsoidal) Earth. The disturbance clearly highlights all of the major subduction zones and large oceanic island chains, all of which are not in local isostatic equilibrium.

Tip

We used the WGS84 ellipsoid here but the workflow is the same for any other ellipsoid. Checkout Available ellipsoids for options.