# Source code for harmonica.isostasy

```
# Copyright (c) 2018 The Harmonica Developers.
# Distributed under the terms of the BSD 3-Clause License.
# SPDX-License-Identifier: BSD-3-Clause
#
# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
#
"""
Function to calculate the thickness of the roots and antiroots assuming the
Airy isostatic hypothesis.
"""
import numpy as np
import xarray as xr
[docs]def isostasy_airy(
topography,
density_crust=2.8e3,
density_mantle=3.3e3,
density_water=1e3,
reference_depth=30e3,
):
r"""
Calculate the isostatic Moho depth from topography using Airy's hypothesis.
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). This assumption is
usually
.. figure:: ../../_static/figures/airy-isostasy.svg
:align: center
:width: 400px
*Schematic of isostatic compensation following the Airy hypothesis.*
The relationship between the topographic/bathymetric heights (:math:`h`)
and the root thickness (:math:`r`) is governed by mass balance relations
and can be found in classic textbooks like [TurcotteSchubert2014]_ and
[Hofmann-WellenhofMoritz2006]_.
On the continents (positive topographic heights):
.. math ::
r = \frac{\rho_{c}}{\rho_m - \rho_{c}} h
while on the oceans (negative topographic heights):
.. math ::
r = \frac{\rho_{c} - \rho_w}{\rho_m - \rho_{c}} h
in which :math:`h` is the topography/bathymetry, :math:`\rho_m` is the
density of the mantle, :math:`\rho_w` is the density of the water, and
:math:`\rho_{c}` is the density of the crust.
The computed root thicknesses will be added to the given reference Moho
depth (:math:`H`) to arrive at the isostatic Moho depth. Use
``reference_depth=0`` if you want the values of the root thicknesses
instead.
Parameters
----------
topography : array or :class:`xarray.DataArray`
Topography height and bathymetry depth in meters. It is usually prudent
to use floating point values instead of integers to avoid integer
division errors.
density_crust : float
Density of the crust in :math:`kg/m^3`.
density_mantle : float
Mantle density in :math:`kg/m^3`.
density_water : float
Water density in :math:`kg/m^3`.
reference_depth : float
The reference Moho depth (:math:`H`) in meters.
Returns
-------
moho_depth : array or :class:`xarray.DataArray`
The isostatic Moho depth in meters.
"""
# Need to cast to array to make sure numpy indexing works as expected for
# 1D DataArray topography
oceans = np.array(topography < 0)
continent = np.logical_not(oceans)
scale = np.full(topography.shape, np.nan, dtype="float")
scale[continent] = density_crust / (density_mantle - density_crust)
scale[oceans] = (density_crust - density_water) / (density_mantle - density_crust)
moho = topography * scale + reference_depth
if isinstance(moho, xr.DataArray):
moho.name = "moho_depth"
moho.attrs["isostasy"] = "Airy"
moho.attrs["density_crust"] = str(density_crust)
moho.attrs["density_mantle"] = str(density_mantle)
moho.attrs["density_water"] = str(density_water)
moho.attrs["reference_depth"] = str(reference_depth)
return moho
```