Note
Go to the end to download the full example code
Airy Isostatic Moho#
According to the Airy hypothesis of isostasy, topography above sea level is
supported by a thickening of the crust (a root) while oceanic basins are
supported by a thinning of the crust (an anti-root). Function
harmonica.isostatic_moho_airy
computes the depth to crust-mantle
interface (the Moho) according to Airy isostasy. The function takes the depth
to the crystalline basement and optionally any layers on top of it. Each layer
is defined by its thickness and its density. In addition, one must assume
a value for the reference thickness of the continental crust in order to
convert the root/anti-root thickness into Moho depth. The function contains
common default values for the reference thickness and crust, mantle
[TurcotteSchubert2014].
We’ll use our sample topography data
(harmonica.datasets.fetch_topography_earth
) to calculate the Airy
isostatic Moho depth of Africa.
![isostatic moho airy](../_images/sphx_glr_isostatic_moho_airy_001.png)
Topography/bathymetry grid:
<xarray.Dataset> Size: 2MB
Dimensions: (longitude: 481, latitude: 511)
Coordinates:
* longitude (longitude) float64 4kB -20.0 -19.83 -19.67 ... 59.67 59.83 60.0
* latitude (latitude) float64 4kB -40.0 -39.83 -39.67 ... 44.67 44.83 45.0
Data variables:
topography (latitude, longitude) float64 2MB -3.528e+03 -3.51e+03 ... 33.0
Moho depth grid:
<xarray.DataArray 'moho_depth' (latitude: 511, longitude: 481)> Size: 2MB
array([[17510.88, 17574.6 , 17712.66, ..., 11715.9 , 11602.62, 11443.32],
[17833.02, 17942.76, 18165.78, ..., 11609.7 , 11581.38, 11461.02],
[17684.34, 17797.62, 18034.8 , ..., 11503.5 , 11570.76, 11521.2 ],
...,
[15167.4 , 15202.8 , 14937.3 , ..., 30151.2 , 30184.8 , 30184.8 ],
[15347.94, 15503.7 , 15266.52, ..., 30145.6 , 30184.8 , 30184.8 ],
[15362.1 , 15517.86, 15252.36, ..., 30151.2 , 30173.6 , 30184.8 ]])
Coordinates:
* longitude (longitude) float64 4kB -20.0 -19.83 -19.67 ... 59.67 59.83 60.0
* latitude (latitude) float64 4kB -40.0 -39.83 -39.67 ... 44.67 44.83 45.0
Attributes:
isostasy: Airy
density_crust: 2800.0
density_mantle: 3300.0
density_water: 1030
import ensaio
import numpy as np
import pygmt
import xarray as xr
import harmonica as hm
# Load the elevation model and cut out the portion of the data corresponding to
# Africa
fname = ensaio.fetch_earth_topography(version=1)
data = xr.load_dataset(fname)
region = (-20, 60, -40, 45)
data_africa = data.sel(latitude=slice(*region[2:]), longitude=slice(*region[:2]))
print("Topography/bathymetry grid:")
print(data_africa)
# Calculate the water thickness
oceans = np.array(data_africa.topography < 0)
water_thickness = data_africa.topography * oceans * -1
water_density = 1030
# Calculate the isostatic Moho depth using the default values for densities and
# reference Moho with water load. We neglect the effect of sediment here, so
# basement elevation refers to topography.
moho = hm.isostatic_moho_airy(
basement=data_africa.topography,
layers={"water": (water_thickness, water_density)},
)
print("\nMoho depth grid:")
print(moho)
# Draw the maps
fig = pygmt.Figure()
pygmt.grd2cpt(grid=moho, cmap="viridis", reverse=True, continuous=True)
title = "Airy isostatic Moho depth of Africa"
fig.grdimage(
region=region,
projection="Y20/0/10c",
frame=["ag", f"+t{title}"],
grid=moho,
cmap=True,
)
fig.coast(shorelines="0.5p,black", resolution="crude")
fig.colorbar(cmap=True, frame=["a10000f2500", "x+lmeters"])
fig.show()
Total running time of the script: (0 minutes 0.401 seconds)