Rename function to isostatic_moho_airy

1. Rename function to isostatic_moho_airy
2. Test density array

examples/isostasy_airy.py

@@ -15,8 +15,7 @@
 (the Moho) according to Airy isostasy. 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, and water densities
-[TurcotteSchubert2014]_.
+values for the reference thickness and crust, mantle [TurcotteSchubert2014]_.
 
 We'll use our sample topography data
 (:func:`harmonica.datasets.fetch_topography_earth`) to calculate the Airy
@@ -35,9 +34,15 @@
 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
-moho = hm.isostasy_airy(data_africa.topography)
+# reference Moho with water load
+moho = hm.isostasy_airy(data_africa.topography, layer_thickness=(water_thickness), 
+       layer_density=(water_density))
 print("\nMoho depth grid:")
 print(moho)
 

harmonica/isostasy.py

@@ -8,50 +8,52 @@
 Function to calculate the thickness of the roots and antiroots assuming the
 Airy isostatic hypothesis.
 """
-import numpy as np
 import xarray as xr
 
 
 def isostasy_airy(
-    topography,
+    basement_elevation,
+    layer_thickness=(),
+    layer_density=(),
     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.
+    Calculate the isostatic Moho depth from rock equvalent 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
+    According to the Airy hypothesis of isostasy, rock equvalent topography
+    above sea level is supported by a thickening of the crust (a root) while
+    rock equvalent topography below sea level is supported by a thinning of
+    the crust (an anti-root). This assumption is usually
 
-    .. figure:: ../../_static/figures/airy-isostasy.svg
+    .. figure:: ../../_static/figures/airy-isostasy.png
         :align: center
         :width: 400px
 
         *Schematic of isostatic compensation following the Airy hypothesis.*
 
-    The relationship between the topographic/bathymetric heights (:math:`h`)
+    The relationship between the rock equvalent topography (:math:`ret`)
     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):
+    Compress all layers's mass above basement (:math:`h`) into rock equvalent
+    topography [Balmino_etal1973]_ :
 
     .. math ::
 
-        r = \frac{\rho_{c}}{\rho_m - \rho_{c}} h
+        ret = h + \frac{\rho_{i}}{\rho_{c}}th_{i} + ..
 
-    while on the oceans (negative topographic heights):
+    Based on rock equvalent topography, the root is calculated as:
 
     .. math ::
-        r = \frac{\rho_{c} - \rho_w}{\rho_m - \rho_{c}} h
+        r = \frac{\rho_{c}}{\rho_m - \rho_{c}} ret
 
-    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.
+    in which :math:`ret` is the rock equvalent 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
@@ -60,10 +62,18 @@ def isostasy_airy(
 
     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.
+    basement_elevation : array or :class:`xarray.DataArray`
+        Basement elevation in meters. It usually refer to topography height
+        and bathymetry depth without sediment cover. When considering
+        sedimentary layer, it refer to crystalline basement. It is usually
+        prudent to use floating point values instead of integers to avoid
+        integer division errors.
+    layer_thickness : tuple contains `xarray.DataArray`
+        Layer thickness in meters. It refer to all layers above
+        topography/basement, including ice, water, and sediment.
+    layer_density : tuple contains `xarray.DataArray` or float
+        Layer density in :math:`kg/m^3`. It refer to all layers above
+        topography/basement, including ice, water, and sediment.
     density_crust : float
         Density of the crust in :math:`kg/m^3`.
     density_mantle : float
@@ -77,21 +87,30 @@ def isostasy_airy(
     -------
     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
+
+    # Define scale factor to calculate Airy root
+    scale = density_crust / (density_mantle - density_crust)
+
+    # For the case: multi-layers above topography/basement
+    if type(layer_thickness) is tuple:
+        # Define total mass of layers abrove topography/basement
+        layer_mass = tuple(
+            ele1 * ele2 for ele1, ele2 in zip(layer_thickness, layer_density)
+        )
+        # Calculate equvalent topography
+        rock_equavalent_topography = (
+            basement_elevation + sum(list(layer_mass)) / density_crust
+        )
+    else:
+        # For the case: a single layer or no layer above topography/basement
+        layer_mass = layer_thickness * layer_density
+        rock_equavalent_topography = basement_elevation + layer_mass / density_crust
+    # Calculate Moho depth
+    moho = rock_equavalent_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

harmonica/tests/test_isostasy.py

@@ -29,12 +29,15 @@ def test_isostasy_airy_zero_topography():
 def test_isostasy_airy():
     "Use a simple integer topography to check the calculations"
     topography = np.array([-2, -1, 0, 1, 2, 3])
+    thickness_water = np.array([2, 1, 0, 0, 0, 0])
+    density_water = 0.5
     true_root = np.array([-0.5, -0.25, 0, 0.5, 1, 1.5])
     root = isostasy_airy(
         topography,
+        layer_thickness=(thickness_water),
+        layer_density=(density_water),
         density_crust=1,
         density_mantle=3,
-        density_water=0.5,
         reference_depth=0,
     )
     npt.assert_equal(root, true_root)
@@ -45,12 +48,17 @@ def test_isostasy_airy_dataarray():
     topography = xr.DataArray(
         np.array([-2, -1, 0, 1, 2, 3]), coords=(np.arange(6),), dims=["something"]
     )
+    thickness_water = xr.DataArray(
+        np.array([2, 1, 0, 0, 0, 0]), coords=(np.arange(6),), dims=["something"]
+    )
+    density_water = 0.5
     true_root = np.array([-0.5, -0.25, 0, 0.5, 1, 1.5])
     root = isostasy_airy(
         topography,
+        layer_thickness=(thickness_water),
+        layer_density=(density_water),
         density_crust=1,
         density_mantle=3,
-        density_water=0.5,
         reference_depth=0,
     )
     assert isinstance(root, xr.DataArray)