Avoid kc becomes zero when NDVI == 1.1 * NDVI_min

doc/source/changelog.rst

@@ -223,6 +223,15 @@ Fixed
   precipitation was being added to every value
   (`#319 <https://github.com/LabSid-USP/RUBEM/issues/319>`__). The golden
   fixtures were regenerated (see ``tests/fixtures/AUDIT.md``).
+- Evapotranspiration: a cell whose NDVI equals ``1.1 * NDVI_min`` takes the
+  ``kc_min`` branch of the crop coefficient, as documented, instead of a
+  crop coefficient of zero
+  (`#320 <https://github.com/LabSid-USP/RUBEM/issues/320>`__). The reference
+  dataset has no cell on the threshold, so the golden fixtures are unchanged.
+- Baseflow: the recession baseflow is limited to the water available in the
+  saturated zone (``TU_S`` of the previous step plus the recharge), so the
+  saturated-zone storage can no longer become negative
+  (`#322 <https://github.com/LabSid-USP/RUBEM/issues/322>`__).
 
 Removed
 ```````

doc/source/overview.rst

@@ -368,28 +368,35 @@ where:
 Baseflow
 `````````
 
-The basic flow is the one that occurs in the saturated soil layer. The base flow is determined by [TERINK2015]_ by groundwater recharge and the recession coefficient. It is calculated only if the moisture content in the saturated zone exceeds a specified threshold.
+The basic flow is the one that occurs in the saturated soil layer. The base flow is determined by [TERINK2015]_ by groundwater recharge and the recession coefficient. It is calculated only if the moisture content in the saturated zone exceeds a specified threshold, Equation :eq:`baseflow-recession`, and it is limited to the water actually available in the saturated zone during the time step, Equation :eq:`baseflow`, so that the saturated zone storage (Equation :eq:`tus`) never becomes negative.
 
 .. math::
-   :label: baseflow
+   :label: baseflow-recession
    :nowrap:
 
     \[
-    BF = \left\{
+    BF_{calc} = \left\{
             \begin{array}{ll}
-            0,  & \mbox{if } TU_S \leq BF_{thresh} \\
-            BF_{T-1} \cdot e^{-\alpha_{GW}} + (1-e^{-\alpha_{GW}}) \cdot REC,  & \mbox{if } TU_S > BF_{thresh} \\
+            0,  & \mbox{if } TU_{S,T-1} \leq BF_{thresh} \\
+            BF_{T-1} \cdot e^{-\alpha_{GW}} + (1-e^{-\alpha_{GW}}) \cdot REC,  & \mbox{if } TU_{S,T-1} > BF_{thresh} \\
             \end{array}
     \right.
     \]
 
+.. math::
+   :label: baseflow
+   :nowrap:
+
+    \[BF = \min\left(BF_{calc},\; TU_{S,T-1} + REC\right)\]
+
 where:
 
 - :math:`BF` – Baseflow (mm);
+- :math:`BF_{calc}` – Baseflow given by the recession equation, before the availability limit (mm);
 - :math:`BF_{T-1}` – Baseflow at the previous time step (mm);
 - :math:`\alpha_{GW}` – Baseflow decay coefficient (-) [parameter to be calibrated];
 - :math:`REC` – Recharge (mm);
-- :math:`TU_S` – Saturated zone moisture content (mm);
+- :math:`TU_{S,T-1}` – Saturated zone moisture content at the previous time step (mm);
 - :math:`BF_{thresh}` – Threshold baseflow, attributed for each watershed (mm).
 
 Water Balance
@@ -420,7 +427,8 @@ where:
 
 - :math:`TU_S` – Saturated zone moisture content (mm);
 - :math:`TU_{S,T-1}` – Saturated zone moisture content at the previous time step (mm);
-- :math:`BF` – Baseflow (mm).
+- :math:`BF` – Baseflow (mm);
+- :math:`REC` – Recharge (mm).
 
 .. math::
    :label: tur

rubem/_dynamic_model.py

@@ -355,11 +355,12 @@ def dynamic(self):
         partial_crop_coef = Interception.get_crop_coef(
             current_ndvi, self.ndvi_min, self.ndvi_max, min_crop_coef, max_crop_coef
         )
-        # If NDVI < 1.1 * NDVI_min, kc = kc_min
-        crop_coef_lt_min_ndvi = pcrfw.scalar(current_ndvi < 1.1 * self.ndvi_min)
-        crop_coef_gt_min_ndvi = pcrfw.scalar(current_ndvi > 1.1 * self.ndvi_min)
-        current_crop_coef = pcr.scalar(
-            (crop_coef_gt_min_ndvi * partial_crop_coef) + (crop_coef_lt_min_ndvi * min_crop_coef)
+        # If NDVI <= 1.1 * NDVI_min, kc = kc_min
+
+        current_crop_coef = pcr.ifthenelse(
+            current_ndvi <= 1.1 * self.ndvi_min,
+            min_crop_coef,
+            partial_crop_coef,
         )
 
         water_stress_coef = pcr.scalar(

rubem/hydrological_processes/_soil.py

@@ -80,31 +80,48 @@ def get_baseflow(
     ) -> Field:
         """Return Baseflow in the pixel [mm].
 
+        The baseflow given by the recession equation is limited to the water
+        available in the saturated zone during the timestep, i.e. the storage
+        of the previous timestep plus the recharge, so that the saturated-zone
+        balance never becomes negative.
+
         :param previous_baseflow: Baseflow at timestep t-1 [mm]
         :type previous_baseflow: Field ``PCRASTER_VALUESCALE=VS_SCALAR``
 
         :param baseflow_recession_coef: Baseflow recession coefficient (Calibrated) [-]
         :type baseflow_recession_coef: float
 
-        :param recharge: Monthly Recharge at timestep t
+        :param recharge: Monthly Recharge at timestep t [mm]
         :type recharge: Field ``PCRASTER_VALUESCALE=VS_SCALAR``
 
-        :param water_cont_sat_zone: Water content at saturated zone [mm]
+        :param water_cont_sat_zone: Water content at saturated zone at timestep t-1 [mm]
         :type water_cont_sat_zone: Field ``PCRASTER_VALUESCALE=VS_SCALAR``
 
         :param threshold_for_baseflow_ocurrence: Threshold for baseflow occurrence [mm]
         :type threshold_for_baseflow_ocurrence: Field ``PCRASTER_VALUESCALE=VS_SCALAR``
 
-        :returns: Monthly Baseflow [mm]
+        :returns: Monthly Baseflow [mm], limited to ``water_cont_sat_zone + recharge``
         :rtype: Field ``PCRASTER_VALUESCALE=VS_SCALAR``
         """
-        # limit condition for base flow
+
+        # Baseflow occurs only above the defined saturated-zone threshold
         cond_lim_for_baseflow = pcr.scalar(water_cont_sat_zone > threshold_for_baseflow_ocurrence)
-        return (
+
+        # Baseflow calculated according to the RUBEM recession equation
+        calculated_baseflow = (
             (previous_baseflow * ((pcr.exp(1)) ** -baseflow_recession_coef))
             + (1 - ((pcr.exp(1)) ** -baseflow_recession_coef)) * recharge
         ) * cond_lim_for_baseflow
 
+        # Water effectively available in the saturated zone
+        available_water = water_cont_sat_zone + recharge
+
+        # Baseflow cannot exceed available water
+        return pcr.min(
+            calculated_baseflow,
+            available_water,
+        )
+
     # First soil layer
     @staticmethod
     def get_actual_soil_moist_cont(

tests/unit/core/test_dynamic_model.py

@@ -6,7 +6,7 @@
 from rubem.configuration.model_configuration import ModelConfiguration
 from rubem.core import DynamicFrameworkWrapper
 from tests.helpers.compare import compare_rasters
-from tests.helpers.synthetic import series_name, write_synthetic_dataset
+from tests.helpers.synthetic import COLS, MISSING, ROWS, series_name, write_synthetic_dataset
 from tests.unit.core.test_core import expected_outputs, run_model
 
 STEP_FLUXES = (
@@ -91,6 +91,56 @@ def test_ndvi_and_landuse_gaps_fall_back_to_the_previous_step(self, tmp_path):
             gap_dir / "out" / interception, distinct_dir / "out" / interception
         ), "the fallback reproduced the second-step NDVI instead of the first"
 
+    @pytest.mark.unit
+    def test_ndvi_at_the_crop_coefficient_threshold_takes_the_minimum_branch(self, tmp_path):
+        """A cell whose NDVI equals 1.1 * NDVI_min gets kc = kc_min (issue #320).
+
+        The documented rule is kc = kc_min whenever NDVI <= 1.1 * NDVI_min, so
+        the first-step evapotranspiration of a run whose NDVI sits exactly on
+        the threshold must reproduce the run whose NDVI equals NDVI_min. Two
+        strict comparisons around the threshold left that cell with kc = 0 and
+        no vegetated-area evapotranspiration at all. A third run above the
+        threshold, whose kc is interpolated, must differ, so that a constant
+        output cannot pass.
+        """
+        import numpy as np
+        import pcraster as pcr
+
+        threshold_dir = tmp_path / "threshold"
+        below_dir = tmp_path / "below"
+        above_dir = tmp_path / "above"
+
+        runs = [
+            (directory, write_synthetic_dataset(str(directory)))
+            for directory in (threshold_dir, below_dir, above_dir)
+        ]
+        # The threshold is the Float32 product PCRaster evaluates from the
+        # ndvi_min raster the model reads: the Python literal 1.1 * 0.2 rounds
+        # to a Float32 below it and would take the kc_min branch with any
+        # comparison operator.
+        ndvi_min = pcr.readmap(str(threshold_dir / "maps" / "ndvi" / "ndvi_min.map"))
+        threshold = float(pcr.pcr2numpy(1.1 * ndvi_min, np.nan)[0, 0])
+        below = float(pcr.pcr2numpy(ndvi_min, np.nan)[0, 0])
+        # Write every step-1 NDVI raster while the synthetic clone is still set;
+        # each model run replaces the PCRaster clone.
+        for (directory, _), value in zip(runs, (threshold, below, 0.6), strict=True):
+            ndvi = pcr.numpy2pcr(
+                pcr.Scalar, np.full((ROWS, COLS), value, dtype=np.float32), MISSING
+            )
+            pcr.report(ndvi, str(directory / "maps" / "ndvi" / series_name("ndvi", 1)))
+        for directory, config in runs:
+            run_model(str(directory), config=config)
+
+        name = series_name("eta", 1)
+        same = compare_rasters(threshold_dir / "out" / name, below_dir / "out" / name)
+        assert same.equal, (
+            f"{name} at NDVI = 1.1 * NDVI_min differs from the NDVI = NDVI_min control:\n"
+            f"{same.report()}"
+        )
+        assert not compare_rasters(threshold_dir / "out" / name, above_dir / "out" / name).equal, (
+            "the threshold run reproduced the interpolated kc of the above-threshold control"
+        )
+
     @pytest.mark.unit
     def test_disabling_tss_produces_no_time_series(self, tmp_path):
         config = write_synthetic_dataset(str(tmp_path))

tests/unit/hydrological_processes/test_module_soil.py

@@ -1,3 +1,5 @@
+import math
+
 import pcraster as pcr
 import pytest
 from pcraster.framework import generalfunctions
@@ -108,6 +110,51 @@ def test_baseflowCalc_cond_false_TUs_lt_EBlim(self):
         expected = 0.0
         assert result == pytest.approx(expected)
 
+    @pytest.mark.unit
+    def test_baseflowCalc_capped_by_available_water(self):
+        """Baseflow is limited to the water available in the saturated zone (issue #322)."""
+        eb_prev = pcr.scalar(10.0)
+        alpha = pcr.scalar(1.0)
+        rec = pcr.scalar(0.5)
+        tus = pcr.scalar(2.0)
+        eb_lim = pcr.scalar(1.0)
+        field = Soil.get_baseflow(eb_prev, alpha, rec, tus, eb_lim)
+        result = generalfunctions.getCellValue(field, 0, 0)
+        # Recession value (about 3.995) exceeds the available water (2.0 + 0.5)
+        expected = 2.5
+        assert result == pytest.approx(expected)
+        # The saturated-zone balance ends at zero instead of going negative
+        balance = Soil.get_actual_water_cont_sat_zone(tus, rec, field)
+        balance_result = generalfunctions.getCellValue(balance, 0, 0)
+        assert balance_result == pytest.approx(0.0, abs=1e-6)
+
+    @pytest.mark.unit
+    def test_baseflowCalc_not_capped_below_available_water(self):
+        eb_prev = pcr.scalar(2.0)
+        alpha = pcr.scalar(1.0)
+        rec = pcr.scalar(1.0)
+        tus = pcr.scalar(5.0)
+        eb_lim = pcr.scalar(1.0)
+        field = Soil.get_baseflow(eb_prev, alpha, rec, tus, eb_lim)
+        result = generalfunctions.getCellValue(field, 0, 0)
+        # Recession value (about 1.368) is below the available water (5.0 + 1.0),
+        # so it is returned unchanged
+        expected = 2.0 * math.exp(-1.0) + (1.0 - math.exp(-1.0)) * 1.0
+        assert result == pytest.approx(expected, rel=1e-5)
+
+    @pytest.mark.unit
+    def test_baseflowCalc_cond_false_water_available(self):
+        eb_prev = pcr.scalar(10.0)
+        alpha = pcr.scalar(1.0)
+        rec = pcr.scalar(0.5)
+        tus = pcr.scalar(0.5)
+        eb_lim = pcr.scalar(1.0)
+        field = Soil.get_baseflow(eb_prev, alpha, rec, tus, eb_lim)
+        result = generalfunctions.getCellValue(field, 0, 0)
+        # Below the threshold no baseflow occurs even though water is available
+        expected = 0.0
+        assert result == pytest.approx(expected)
+
     @pytest.mark.unit
     def test_baseflowCalc_None_values(self):
         with pytest.raises(TypeError):