Tesseroids#
When our region of interest covers several longitude and latitude degrees, utilizing Cartesian coordinates to model geological structures might introduce significant errors: they don’t take into account the curvature of the Earth. Instead, we would need to work in Spherical coordinates. A common approach to forward model bodies in geocentric spherical coordinates is to make use of tesseroids.
A tesseroid (a.k.a spherical prism) is a three dimensional body defined by the volume contained by two longitudinal boundaries, two latitudinal boundaries and the surfaces of two concentric spheres of different radii (see Figure: Tesseroid).
Figure: Tesseroid#
Tesseroid defined by two longitude coordinates (\(\lambda_1\) and \(\lambda_2\)), two latitude coordinates (\(\phi_1\) and \(\phi_2\)) and the surfaces of two concentric spheres of radii \(r_1\) and \(r_2\). This figure is a modified version of [Uieda2015].
Through the harmonica.tesseroid_gravity function we can calculate the
gravitational field of any tesseroid with a given density on any computation
point. Each tesseroid can be represented through a tuple containing its six
boundaries in the following order: west, east, south, north, bottom,
top, where the former four are its longitudinal and latitudinal boundaries in
decimal degrees and the latter two are the two radii given in meters.
These two radii represent the top and bottom surfaces of the tesseroid, and should be given as distances from the center of the Earth. Note this is different from the vertical boundaries used for prisms in Cartesian coordinates, which are given as heights above or below some reference level (e.g., mean sea level or a reference ellipsoid).
Note
The harmonica.tesseroid_gravity numerically computed the
gravitational fields of tesseroids by applying a method that applies the
Gauss-Legendre Quadrature along with a bidimensional adaptive discretization
algorithm. Refer to [Soler2019] for more details.
Lets define a single tesseroid and compute the gravitational potential it generates on a regular grid of computation points located at 10 km above its top boundary.
Get the WGS84 reference ellipsoid from boule so we can obtain its mean
radius:
import boule as bl
ellipsoid = bl.WGS84
mean_radius = ellipsoid.mean_radius
Define the tesseroid and its density (in kg per cubic meters):
tesseroid = (-70, -50, -40, -20, mean_radius - 10e3, mean_radius)
density = 2670
Define a set of computation points located on a regular grid at 100 km above the top surface of the tesseroid:
import verde as vd
coordinates = vd.grid_coordinates(
region=[-80, -40, -50, -10],
shape=(80, 80),
extra_coords=100e3 + mean_radius,
)
Lets compute the downward component of the gravitational acceleration it generates on the computation point:
import harmonica as hm
gravity = hm.tesseroid_gravity(coordinates, tesseroid, density, field="g_z")
Important
The downward component \(g_z\) of the gravitational acceleration computed in spherical coordinates corresponds to \(-g_r\), where \(g_r\) is the radial component.
And finally plot the computed gravitational field
import pygmt
grid = vd.make_xarray_grid(
coordinates, gravity, data_names="gravity", extra_coords_names="extra")
fig = pygmt.Figure()
title = "Downward component of gravitational acceleration"
with pygmt.config(FONT_TITLE="12p"):
fig.grdimage(
region=[-80, -40, -50, -10],
projection="M-60/-30/10c",
grid=grid.gravity,
frame=["a", f"+t{title}"],
cmap="viridis",
)
fig.colorbar(cmap=True, frame=["a200f50", "x+lmGal"])
fig.coast(shorelines="1p,black")
# Plot edges of tesseroid
fig.plot(
x=[tesseroid[0], tesseroid[1], tesseroid[1], tesseroid[0], tesseroid[0]],
y=[tesseroid[2], tesseroid[2], tesseroid[3], tesseroid[3], tesseroid[2]],
pen="1p,red",
label="Tesseroid boundary",
)
fig.legend()
fig.show()
Multiple tesseroids#
We can compute the gravitational field of a set of tesseroids by passing a list
of them, where each tesseroid is defined as mentioned before, and then making
a single call of the harmonica.tesseroid_gravity function.
Lets define a set of four prisms along with their densities:
tesseroids = [
[-70, -65, -40, -35, mean_radius - 100e3, mean_radius],
[-55, -50, -40, -35, mean_radius - 100e3, mean_radius],
[-70, -65, -25, -20, mean_radius - 100e3, mean_radius],
[-55, -50, -25, -20, mean_radius - 100e3, mean_radius],
]
densities = [2670 , 2670, 2670, 2670]
Compute their gravitational effect on a grid of computation points:
coordinates = vd.grid_coordinates(
region=[-80, -40, -50, -10],
shape=(80, 80),
extra_coords=100e3 + mean_radius,
)
gravity = hm.tesseroid_gravity(coordinates, tesseroids, densities, field="g_z")
And plot the results:
grid = vd.make_xarray_grid(
coordinates, gravity, data_names="gravity", extra_coords_names="extra")
fig = pygmt.Figure()
title = "Downward component of gravitational acceleration"
with pygmt.config(FONT_TITLE="12p"):
fig.grdimage(
region=[-80, -40, -50, -10],
projection="M-60/-30/10c",
grid=grid.gravity,
frame=["a", f"+t{title}"],
cmap="viridis",
)
fig.colorbar(cmap=True, frame=["a1000f500", "x+lmGal"])
fig.coast(shorelines="1p,black")
# Plot edges of tesseroids
for i, tesseroid in enumerate(tesseroids):
if i == 0:
label="Tesseroid boundaries"
else:
label=None
fig.plot(
x=[tesseroid[0], tesseroid[1], tesseroid[1], tesseroid[0], tesseroid[0]],
y=[tesseroid[2], tesseroid[2], tesseroid[3], tesseroid[3], tesseroid[2]],
pen="1p,red",
label=label,
)
fig.legend()
fig.show()
Tesseroids with variable density#
The harmonica.tesseroid_gravity is capable of computing the
gravitational effects of tesseroids whose density is defined through
a continuous function of the radial coordinate. This is achieved by the
application of the method introduced in [Soler2021].
To do so we need to define a regular Python function for the density, which
should have a single argument (the radius coordinate) and return the
density of the tesseroids at that radial coordinate.
In addition, we need to decorate the density function with
numba.jit(nopython=True) or numba.njit for short.
Lets compute the gravitational effect of four tesseroids whose densities are
given by a custom linear density function.
Start by defining the tesseroids
tesseroids = (
[-70, -60, -40, -30, mean_radius - 3e3, mean_radius],
[-70, -60, -30, -20, mean_radius - 5e3, mean_radius],
[-60, -50, -40, -30, mean_radius - 7e3, mean_radius],
[-60, -50, -30, -20, mean_radius - 10e3, mean_radius],
)
Then, define a linear density function. We need to use the jit decorator so
Numba can run the forward model efficiently.
from numba import njit
@njit
def density(radius):
"""Linear density function"""
top = mean_radius
bottom = mean_radius - 10e3
density_top = 2670
density_bottom = 3000
slope = (density_top - density_bottom) / (top - bottom)
return slope * (radius - bottom) + density_bottom
Lets create a set of computation points located on a regular grid at 100km above the mean Earth radius:
coordinates = vd.grid_coordinates(
region=[-80, -40, -50, -10],
shape=(80, 80),
extra_coords=100e3 + ellipsoid.mean_radius,
)
And compute the gravitational fields the tesseroids generate:
gravity = hm.tesseroid_gravity(coordinates, tesseroids, density, field="g_z")
Finally, lets plot it:
grid = vd.make_xarray_grid(
coordinates, gravity, data_names="gravity", extra_coords_names="extra")
fig = pygmt.Figure()
title = "Downward component of gravitational acceleration"
with pygmt.config(FONT_TITLE="12p"):
fig.grdimage(
region=[-80, -40, -50, -10],
projection="M-60/-30/10c",
grid=grid.gravity,
frame=["a", f"+t{title}"],
cmap="viridis",
)
fig.colorbar(cmap=True, frame=["a200f100", "x+lmGal"])
fig.coast(shorelines="1p,black")
# Plot edges of tesseroids
for i, tesseroid in enumerate(tesseroids):
if i == 0:
label="Tesseroid boundaries"
else:
label=None
fig.plot(
x=[tesseroid[0], tesseroid[1], tesseroid[1], tesseroid[0], tesseroid[0]],
y=[tesseroid[2], tesseroid[2], tesseroid[3], tesseroid[3], tesseroid[2]],
pen="1p,red",
label=label,
)
fig.legend()
fig.show()
