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

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


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),
)


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 matplotlib.pyplot as plt
import cartopy.crs as ccrs

fig = plt.figure(figsize=(8, 9))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-60))
img = ax.pcolormesh(
coordinates[0], coordinates[1], gravity, transform=ccrs.PlateCarree()
)
plt.colorbar(img, ax=ax, pad=0, aspect=50, orientation="horizontal", label="mGal")
ax.coastlines()
ax.set_title("Downward component of gravitational acceleration")
plt.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 = [
]
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),
)
gravity = hm.tesseroid_gravity(coordinates, tesseroids, densities, field="g_z")


And plot the results:

fig = plt.figure(figsize=(8, 9))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-60))
img = ax.pcolormesh(
coordinates[0], coordinates[1], gravity, transform=ccrs.PlateCarree()
)
plt.colorbar(img, ax=ax, pad=0, aspect=50, orientation="horizontal", label="mGal")
ax.coastlines()
ax.set_title("Downward component of gravitational acceleration")
plt.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 = (
)


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
"""Linear density function"""
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),
)


And compute the gravitational fields the tesseroids generate:

gravity = hm.tesseroid_gravity(coordinates, tesseroids, density, field="g_z")


Finally, lets plot it:

fig = plt.figure(figsize=(8, 9))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-60))
img = ax.pcolormesh(*coordinates[:2], gravity, transform=ccrs.PlateCarree())
plt.colorbar(img, ax=ax, pad=0, aspect=50, orientation="horizontal", label="mGal")
ax.coastlines()
ax.set_title("Downward component of gravitational acceleration")
plt.show()