diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 1f4ecb1f..d041b9ec 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -223,6 +223,15 @@ Fixed precipitation was being added to every value (`#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 `__). 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 `__). Removed ``````` diff --git a/doc/source/overview.rst b/doc/source/overview.rst index cb9401d2..992138a8 100644 --- a/doc/source/overview.rst +++ b/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 diff --git a/rubem/_dynamic_model.py b/rubem/_dynamic_model.py index a98ed6ba..77f8cf7c 100644 --- a/rubem/_dynamic_model.py +++ b/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( diff --git a/rubem/hydrological_processes/_soil.py b/rubem/hydrological_processes/_soil.py index 91223ab2..294e1c26 100644 --- a/rubem/hydrological_processes/_soil.py +++ b/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( diff --git a/tests/unit/core/test_dynamic_model.py b/tests/unit/core/test_dynamic_model.py index fc0afce1..d28d7887 100644 --- a/tests/unit/core/test_dynamic_model.py +++ b/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)) diff --git a/tests/unit/hydrological_processes/test_module_soil.py b/tests/unit/hydrological_processes/test_module_soil.py index 27eb7fa6..bf30ed1d 100644 --- a/tests/unit/hydrological_processes/test_module_soil.py +++ b/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):