# 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 xarray as xr
[docs]def isostatic_moho_airy(
basement,
layers=None,
density_crust=2.8e3,
density_mantle=3.3e3,
reference_depth=30e3,
):
r"""
Calculate the isostatic Moho depth using Airy's hypothesis.
Take the height of the crystalline basement and optional additional layers
located on top of it. Each layer must be specified through its vertical
thickness and its corresponding density. Return the depth to the
Mohorovicic discontinuity (crust-mantle interface).
Parameters
----------
basement : float or array
Height of the crystalline basement in meters.
It usually refer to topography and bathymetry height without sediment
cover.
When considering sedimentary basins, it refers to crystalline basement
(topography/bathymetry minus sediment thickness).
layers : dict (optional)
Dictionary that contains information about the thickness and density of
the layers located above the ``basement``.
For each layer, a single item should be created: its key will be the
layer name as a ``str`` and its values must be tuples containing the
layer thickness in meters and the layer density (in :math:`kg/m^3`)
in that order.
Thicknesses and densities can be floats or arrays.
If ``None``, no layers will be considered.
Default as ``None``.
density_crust : float or array (optional)
Density of the crust in :math:`kg/m^3`.
density_mantle : float or array (optional)
Mantle density in :math:`kg/m^3`.
reference_depth : float or array (optional)
The reference Moho depth (:math:`H`) in meters.
Returns
-------
moho_depth : float or array
The isostatic Moho depth in meters.
Notes
-----
According to the Airy hypothesis of isostasy, rock equivalent topography
above sea level is supported by a thickening of the crust (a root) while
rock equivalent topography below sea level is supported by a thinning of
the crust (an anti-root). This assumption is usually
.. figure:: ../../_static/figures/airy-isostasy-moho.svg
:align: center
:width: 400px
*Schematic of isostatic compensation following the Airy hypothesis.*
The relationship between the rock equivalent topography (:math:`r_{et}`)
and the root thickness (:math:`r`) is governed by mass balance relations
and can be found in classic textbooks like [TurcotteSchubert2014]_ and
[Hofmann-WellenhofMoritz2006]_.
Compress all layers' mass above basement (:math:`h`) into *rock equivalent
topography* [Balmino1973]_ :
.. math ::
r_{et} = h + \sum\limits_{i=1}^N \frac{\rho_{i}}{\rho_{c}} t_{i}
Based on rock equivalent topography, the root is calculated as:
.. math ::
r = \frac{\rho_{c}}{\rho_m - \rho_{c}} r_{et}
in which :math:`r_{et}` is the rock equivalent topography , :math:`\rho_m`
is the density of the mantle, 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.
Examples
--------
Simple model of continental topography with a sedimentary basin on top
>>> # Define crystalline basement height (in meters)
>>> basement = 1200
>>> # Define a layer of sediments with a thickness of 200m
>>> sediments_thickness = 200
>>> sediments_density = 2300
>>> # Get depth of the Moho following Airy's isostatic hypothesis
>>> moho_depth = isostatic_moho_airy(
... basement,
... layers={"sediments": (sediments_thickness, sediments_density)}
... )
>>> moho_depth
37640.0
Simple model of oceanic sedimentary basin
>>> # Define bathymetry (in meters)
>>> bathymetry = -3000
>>> # Define a layer of sediments with a thickness of 400m
>>> sediments_thickness = 400
>>> sediments_density = 2200
>>> # Define a layer for the oceanic water
>>> water_thickness = abs(bathymetry)
>>> water_density = 1040
>>> # Get depth of the Moho following Airy's isostatic hypothesis
>>> moho_depth = isostatic_moho_airy(
... bathymetry - sediments_thickness,
... layers={
... "sediments": (sediments_thickness, sediments_density),
... "water": (water_thickness, water_density),
... }
... )
>>> moho_depth
18960.0
"""
# Compute equivalent topography for the layers (if any)
layers_equivalent_topography = 0
if layers is not None:
for thickness, density in layers.values():
layers_equivalent_topography += thickness * density
layers_equivalent_topography /= density_crust
# Calculate rock equivalent topography
rock_equivalent_topography = basement + layers_equivalent_topography
# Calculate Moho depth
scale = density_crust / (density_mantle - density_crust)
moho = rock_equivalent_topography * scale + reference_depth
# Add attributes to the xr.DataArray
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)
if layers is not None:
for name, (_, density) in layers.items():
moho.attrs[f"density_{name}"] = density
return moho