Merge pull request #318 from ImperialCollegeLondon/bugfix/soil_moistu…

…re_positive

soil evaporation constraint

tests/models/hydrology/test_above_ground.py

@@ -33,19 +33,21 @@ def test_calculate_soil_evaporation(wind, dens_air, latvap):
     )
 
     result = calculate_soil_evaporation(
-        temperature=np.array([10.0, 20.0, 30.0]),
+        temperature=np.array([20.0, 20.0, 30.0]),
         wind_speed=wind,
-        relative_humidity=np.array([70, 80, 90]),
+        relative_humidity=np.array([80, 80, 90]),
         atmospheric_pressure=np.array([90, 90, 90]),
-        soil_moisture=np.array([0.1, 0.5, 0.9]),
+        soil_moisture=np.array([0.01, 0.1, 0.5]),
+        soil_moisture_residual=0.1,
+        soil_moisture_capacity=0.9,
         celsius_to_kelvin=HydroConsts.celsius_to_kelvin,
         density_air=dens_air,
         latent_heat_vapourisation=latvap,
         gas_constant_water_vapour=HydroConsts.gas_constant_water_vapour,
         heat_transfer_coefficient=HydroConsts.heat_transfer_coefficient,
     )
 
-    exp_result = np.array([1.523354, 3.86474, 4.473155])
+    exp_result = np.array([0.060855, 0.060855, 4.473155])
     np.testing.assert_allclose(result, exp_result, rtol=0.01)
 
 

virtual_rainforest/models/hydrology/above_ground.py

@@ -19,6 +19,8 @@ def calculate_soil_evaporation(
     relative_humidity: NDArray[np.float32],
     atmospheric_pressure: NDArray[np.float32],
     soil_moisture: NDArray[np.float32],
+    soil_moisture_residual: Union[float, NDArray[np.float32]],
+    soil_moisture_capacity: Union[float, NDArray[np.float32]],
     wind_speed: Union[float, NDArray[np.float32]],
     celsius_to_kelvin: float,
     density_air: Union[float, NDArray[np.float32]],
@@ -35,8 +37,8 @@ def calculate_soil_evaporation(
 
     :math:`E_{g} = \frac{\rho_{air}}{R_{a}} * (\alpha * q_{sat}(T_{s}) - q_{g})`
 
-    where :math:`\Theta` is the top soil moisture (relative volumetric water content),
-    :math:`E_{g}` is the evaporation flux (W m-2), :math:`\rho_{air}` is the
+    where :math:`\Theta` is the available top soil moisture (relative volumetric water
+    content), :math:`E_{g}` is the evaporation flux (W m-2), :math:`\rho_{air}` is the
     density of air (kg m-3), :math:`R_{a}` is the aerodynamic resistance (unitless),
     :math:`q_{sat}(T_{s})` (unitless) is the saturated specific humidity, and
     :math:`q_{g}` is the surface specific humidity (unitless); see Mahfouf (1991).
@@ -49,6 +51,8 @@ def calculate_soil_evaporation(
         relative_humidity: relative humidity at reference height, []
         atmospheric_pressure: atmospheric pressure at reference height, [kPa]
         soil_moisture: Volumetric relative water content, [unitless]
+        soil_moisture_residual: residual soil moisture, [unitless]
+        soil_moisture_capacity: soil moisture capacity, [unitless]
         wind_speed: wind speed at reference height, [m s-1]
         celsius_to_kelvin: factor to convert teperature from Celsius to Kelvin
         density_air: density if air, [kg m-3]
@@ -63,8 +67,15 @@ def calculate_soil_evaporation(
     # Convert temperature to Kelvin
     temperature_k = temperature + celsius_to_kelvin
 
+    # Available soil moisture
+    soil_moisture_free = np.clip(
+        (soil_moisture - soil_moisture_residual),
+        0.0,
+        (soil_moisture_capacity - soil_moisture_residual),
+    )
+
     # Estimate alpha using the Barton (1979) equation
-    barton_ratio = (1.8 * soil_moisture) / (soil_moisture + 0.3)
+    barton_ratio = (1.8 * soil_moisture_free) / (soil_moisture_free + 0.3)
     alpha = np.where(barton_ratio > 1, 1, barton_ratio)
 
     saturation_vapour_pressure = 0.6112 * np.exp(

virtual_rainforest/models/hydrology/hydrology_model.py

@@ -406,6 +406,13 @@ def update(self, time_index: int) -> None:
             * soil_layer_thickness
         ).to_numpy()
 
+        top_soil_moisture_capacity_mm = (
+            self.constants.soil_moisture_capacity * soil_layer_thickness[0]
+        )
+        top_soil_moisture_residual_mm = (
+            self.constants.soil_moisture_residual * soil_layer_thickness[0]
+        )
+
         # Create lists for output variables to store daily data
         daily_lists: dict = {name: [] for name in self.vars_updated}
 
@@ -443,15 +450,19 @@ def update(self, time_index: int) -> None:
             soil_moisture_infiltrated = np.clip(
                 soil_moisture_mm[0] + precipitation_surface,
                 0,
-                (self.constants.soil_moisture_capacity * soil_layer_thickness[0]),
+                top_soil_moisture_capacity_mm,
             )
 
             # Calculate daily soil evaporation, [mm]
+            top_soil_moisture_vol = soil_moisture_infiltrated / soil_layer_thickness[0]
+
             soil_evaporation = above_ground.calculate_soil_evaporation(
                 temperature=subcanopy_temperature,
                 relative_humidity=subcanopy_humidity,
                 atmospheric_pressure=subcanopy_pressure,
-                soil_moisture=soil_moisture_infiltrated / soil_layer_thickness[0],
+                soil_moisture=top_soil_moisture_vol,
+                soil_moisture_residual=self.constants.soil_moisture_residual,
+                soil_moisture_capacity=self.constants.soil_moisture_capacity,
                 wind_speed=0.1,  # m/s TODO wind_speed in data object
                 celsius_to_kelvin=self.constants.celsius_to_kelvin,
                 density_air=self.constants.density_air,
@@ -465,7 +476,11 @@ def update(self, time_index: int) -> None:
             soil_moisture_evap: NDArray[np.float32] = np.concatenate(
                 (
                     np.expand_dims(
-                        (soil_moisture_infiltrated - soil_evaporation),
+                        np.clip(
+                            (soil_moisture_infiltrated - soil_evaporation),
+                            top_soil_moisture_residual_mm,
+                            top_soil_moisture_capacity_mm,
+                        ),
                         axis=0,
                     ),
                     soil_moisture_mm[1:],