Merge branch 'develop' into feature/bypass_flow

tests/models/hydrology/test_below_ground.py

@@ -83,3 +83,30 @@ def test_update_soil_moisture():
     )
 
     np.testing.assert_allclose(result, exp_result, rtol=0.001)
+
+
+def test_soil_moisture_to_matric_potential():
+    """Test."""
+
+    from virtual_rainforest.models.hydrology.below_ground import (
+        soil_moisture_to_matric_potential,
+    )
+    from virtual_rainforest.models.hydrology.constants import HydroConsts
+
+    soil_moisture = np.array([[0.2, 0.2, 0.2], [0.5, 0.5, 0.5]])
+
+    result = soil_moisture_to_matric_potential(
+        soil_moisture=soil_moisture,
+        soil_moisture_capacity=HydroConsts.soil_moisture_capacity,
+        soil_moisture_residual=HydroConsts.soil_moisture_residual,
+        nonlinearily_parameter=HydroConsts.nonlinearily_parameter,
+        alpha=HydroConsts.alpha,
+    )
+    exp_result = np.array(
+        [
+            [26.457513, 26.457513, 26.457513],
+            [5.773503, 5.773503, 5.773503],
+        ]
+    )
+
+    np.testing.assert_allclose(result, exp_result)

tests/models/hydrology/test_hydrology_model.py

@@ -408,6 +408,20 @@ def test_setup(
             dim="layers",
         ).assign_coords(model.data["layer_heights"].coords)
 
+        exp_matric_pot = xr.concat(
+            [
+                DataArray(
+                    np.full((13, 3), np.nan),
+                    dims=["layers", "cell_id"],
+                ),
+                DataArray(
+                    [[5.633922, 5.633922, 5.633922], [6.846069, 6.846069, 6.846069]],
+                    dims=["layers", "cell_id"],
+                ),
+            ],
+            dim="layers",
+        ).assign_coords(model.data["layer_heights"].coords)
+
         exp_surf_prec = DataArray(
             [177.113493, 177.113493, 177.113493],
             dims=["cell_id"],
@@ -463,6 +477,7 @@ def test_setup(
         np.testing.assert_allclose(
             model.data["surface_runoff_accumulated"], exp_runoff_acc
         )
+        np.testing.assert_allclose(model.data["matric_potential"], exp_matric_pot)
 
 
 def test_calculate_layer_thickness():

virtual_rainforest/models/hydrology/below_ground.py

@@ -165,3 +165,52 @@ def update_soil_moisture(
         )
 
     return soil_moisture_updated
+
+
+def soil_moisture_to_matric_potential(
+    soil_moisture: NDArray[np.float32],
+    soil_moisture_capacity: Union[float, NDArray[np.float32]],
+    soil_moisture_residual: Union[float, NDArray[np.float32]],
+    nonlinearily_parameter: Union[float, NDArray[np.float32]],
+    alpha: Union[float, NDArray[np.float32]],
+) -> NDArray[np.float32]:
+    r"""Convert soil moisture to matric potential using van Genuchten/Mualem model.
+
+    The soil water content is converted to matric potential as follows:
+
+    :math:`S = \frac{\Theta-\Theta_{r}}{\Theta_{s}-\Theta_{r}} = (1+(\alpha*\Phi)^n)^-m`
+
+    where :math:`S` is the effective saturation (dimensionless), :math:`\Theta_{r}` and
+    :math:`\Theta_{s}` are the residual and saturated moisture content or soil moisture
+    capacity, respectively. `math`:\Phi` is the soil water matric potential and
+    :math:`m`, :math:`n`, and :math:`\alpha` are shape parameters of the water retention
+    curve.
+
+    Args:
+        soil_moisture: Volumetric relative water content in top soil, [unitless]
+        soil_moisture_capacity: soil moisture capacity, [unitless]
+        soil_moisture_residual: residual soil moisture, [unitless]
+        nonlinearily_parameter: dimensionless parameter n in van Genuchten model that
+            describes the degree of nonlinearity of the relationship between the
+            volumetric water content and the soil matric potential.
+        alpha: dimensionless parameter alpha in van Genuchten model that corresponds
+            approximately to the inverse of the air-entry value, [kPa-1]
+
+    Returns:
+        soil water matric potential, [kPa]
+    """
+
+    shape_parameter = 1 - 1 / nonlinearily_parameter
+
+    # Calculate soil effective saturation in rel. vol. water content for each layer:
+    effective_saturation = (soil_moisture - soil_moisture_residual) / (
+        soil_moisture_capacity - soil_moisture_residual
+    )
+
+    # Solve for phi
+    effective_saturation_m = effective_saturation ** (-1 / shape_parameter)
+    effective_saturation_1 = effective_saturation_m - 1
+    effective_saturation_n = effective_saturation_1 ** (1 / nonlinearily_parameter)
+    soil_matric_potential = effective_saturation_n / alpha
+
+    return soil_matric_potential

virtual_rainforest/models/hydrology/constants.py

@@ -94,6 +94,10 @@ class HydroConsts:
     """Ground water storage capacity in relative volumetric water content. This might be
     replaced with the implementation of below ground horizontal flow."""
 
+    alpha: float = 0.3
+    """Dimensionless parameter alpha in van Genuchten model that corresponds
+    approximately to the inverse of the air-entry value, [kPa-1]"""
+
     infiltration_shape_parameter: float = 1.0
     """Empirical shape parameter that affects how much of the water available for
     infiltration goes directly to groundwater via preferential bypass flow. A value of

virtual_rainforest/models/hydrology/hydrology_model.py

@@ -90,6 +90,7 @@ class HydrologyModel(BaseModel):
         "soil_evaporation",
         "stream_flow",  # P-ET; later surface_runoff_acc + below_ground_acc
         "surface_runoff_accumulated",
+        "matric_potential",
     )
     """Variables updated by the hydrology model."""
 
@@ -354,6 +355,7 @@ def update(self, time_index: int) -> None:
         * soil_evaporation, [mm]
         * stream_flow, [mm/timestep], currently simply P-ET
         * surface_runoff_accumulated, [mm]
+        * matric_potential, [kPa]
         """
         # Determine number of days, currently only 30 days (=1 month)
         if self.update_interval != Quantity("1 month"):
@@ -537,6 +539,16 @@ def update(self, time_index: int) -> None:
                 soil_moisture_updated / soil_layer_thickness
             )
 
+            # Convert soil moisture to matric potential
+            matric_potential = below_ground.soil_moisture_to_matric_potential(
+                soil_moisture=soil_moisture_updated / soil_layer_thickness,
+                soil_moisture_capacity=self.constants.soil_moisture_capacity,
+                soil_moisture_residual=self.constants.soil_moisture_residual,
+                nonlinearily_parameter=self.constants.nonlinearily_parameter,
+                alpha=self.constants.alpha,
+            )
+            daily_lists["matric_potential"].append(matric_potential)
+
             # update soil_moisture_mm for next day
             soil_moisture_mm = soil_moisture_updated
 
@@ -557,28 +569,28 @@ def update(self, time_index: int) -> None:
             coords={"cell_id": self.data.grid.cell_id},
         )
 
-        # Return mean soil moisture, [-], and add atmospheric layers (nan)
-        soil_hydrology["soil_moisture"] = DataArray(
-            np.concatenate(
-                (
-                    np.full(
-                        (
-                            len(self.layer_roles) - self.layer_roles.count("soil"),
-                            self.data.grid.n_cells,
+        # Return mean soil moisture, [-], and matric potential, [kPa], and add
+        # atmospheric layers (nan)
+        for var in ["soil_moisture", "matric_potential"]:
+            soil_hydrology[var] = DataArray(
+                np.concatenate(
+                    (
+                        np.full(
+                            (
+                                len(self.layer_roles) - self.layer_roles.count("soil"),
+                                self.data.grid.n_cells,
+                            ),
+                            np.nan,
+                        ),
+                        np.mean(
+                            np.stack(daily_lists[var], axis=0),
+                            axis=0,
                         ),
-                        np.nan,
-                    ),
-                    np.mean(
-                        np.stack(daily_lists["soil_moisture"], axis=0),
-                        axis=0,
                     ),
                 ),
-            ),
-            dims=self.data["soil_moisture"].dims,
-            coords=self.data["soil_moisture"].coords,
-        )
-
-        # TODO Convert to matric potential
+                dims=self.data["layer_heights"].dims,
+                coords=self.data["layer_heights"].coords,
+            )
 
         # Calculate accumulated surface runoff for model time step
         # Get the runoff created by SPLASH or initial data set