From 007b89a790b16f597503a82d76a19767cb901b4a Mon Sep 17 00:00:00 2001 From: vgro Date: Wed, 27 Sep 2023 13:10:14 +0100 Subject: [PATCH 1/2] soil evaporation constraint --- tests/models/hydrology/test_above_ground.py | 3 ++- .../models/hydrology/above_ground.py | 11 +++++++--- .../models/hydrology/hydrology_model.py | 20 ++++++++++++++++--- 3 files changed, 27 insertions(+), 7 deletions(-) diff --git a/tests/models/hydrology/test_above_ground.py b/tests/models/hydrology/test_above_ground.py index 7f8ba8966..2a3ed1592 100644 --- a/tests/models/hydrology/test_above_ground.py +++ b/tests/models/hydrology/test_above_ground.py @@ -38,6 +38,7 @@ def test_calculate_soil_evaporation(wind, dens_air, latvap): relative_humidity=np.array([70, 80, 90]), atmospheric_pressure=np.array([90, 90, 90]), soil_moisture=np.array([0.1, 0.5, 0.9]), + soil_moisture_residual=0.1, celsius_to_kelvin=HydroConsts.celsius_to_kelvin, density_air=dens_air, latent_heat_vapourisation=latvap, @@ -45,7 +46,7 @@ def test_calculate_soil_evaporation(wind, dens_air, latvap): heat_transfer_coefficient=HydroConsts.heat_transfer_coefficient, ) - exp_result = np.array([1.523354, 3.86474, 4.473155]) + exp_result = np.array([0.053332, 3.86474, 4.473155]) np.testing.assert_allclose(result, exp_result, rtol=0.01) diff --git a/virtual_rainforest/models/hydrology/above_ground.py b/virtual_rainforest/models/hydrology/above_ground.py index 8fb6e56c2..fb654cbb0 100644 --- a/virtual_rainforest/models/hydrology/above_ground.py +++ b/virtual_rainforest/models/hydrology/above_ground.py @@ -19,6 +19,7 @@ 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]], wind_speed: Union[float, NDArray[np.float32]], celsius_to_kelvin: float, density_air: Union[float, NDArray[np.float32]], @@ -35,8 +36,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 +50,7 @@ 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] 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 +65,11 @@ def calculate_soil_evaporation( # Convert temperature to Kelvin temperature_k = temperature + celsius_to_kelvin + # Available soil moisture + soil_moisture_free = soil_moisture - 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( diff --git a/virtual_rainforest/models/hydrology/hydrology_model.py b/virtual_rainforest/models/hydrology/hydrology_model.py index 5725061bf..25a0e64ca 100644 --- a/virtual_rainforest/models/hydrology/hydrology_model.py +++ b/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,18 @@ 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, 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 +475,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:], From 9ff37270dbd2270fe0f2330d3452a3a5ca64dde8 Mon Sep 17 00:00:00 2001 From: vgro Date: Thu, 28 Sep 2023 09:11:29 +0100 Subject: [PATCH 2/2] test includes case where soil moisture < residual --- tests/models/hydrology/test_above_ground.py | 9 +++++---- virtual_rainforest/models/hydrology/above_ground.py | 8 +++++++- virtual_rainforest/models/hydrology/hydrology_model.py | 1 + 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/tests/models/hydrology/test_above_ground.py b/tests/models/hydrology/test_above_ground.py index 2a3ed1592..894a63927 100644 --- a/tests/models/hydrology/test_above_ground.py +++ b/tests/models/hydrology/test_above_ground.py @@ -33,12 +33,13 @@ 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, @@ -46,7 +47,7 @@ def test_calculate_soil_evaporation(wind, dens_air, latvap): heat_transfer_coefficient=HydroConsts.heat_transfer_coefficient, ) - exp_result = np.array([0.053332, 3.86474, 4.473155]) + exp_result = np.array([0.060855, 0.060855, 4.473155]) np.testing.assert_allclose(result, exp_result, rtol=0.01) diff --git a/virtual_rainforest/models/hydrology/above_ground.py b/virtual_rainforest/models/hydrology/above_ground.py index fb654cbb0..b484b4eaa 100644 --- a/virtual_rainforest/models/hydrology/above_ground.py +++ b/virtual_rainforest/models/hydrology/above_ground.py @@ -20,6 +20,7 @@ def calculate_soil_evaporation( 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]], @@ -51,6 +52,7 @@ def calculate_soil_evaporation( 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] @@ -66,7 +68,11 @@ def calculate_soil_evaporation( temperature_k = temperature + celsius_to_kelvin # Available soil moisture - soil_moisture_free = soil_moisture - soil_moisture_residual + 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_free) / (soil_moisture_free + 0.3) diff --git a/virtual_rainforest/models/hydrology/hydrology_model.py b/virtual_rainforest/models/hydrology/hydrology_model.py index 25a0e64ca..0cfc472c1 100644 --- a/virtual_rainforest/models/hydrology/hydrology_model.py +++ b/virtual_rainforest/models/hydrology/hydrology_model.py @@ -462,6 +462,7 @@ def update(self, time_index: int) -> None: atmospheric_pressure=subcanopy_pressure, 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,