diff --git a/tests/models/litter/conftest.py b/tests/models/litter/conftest.py index bdad8df06..26f972eab 100644 --- a/tests/models/litter/conftest.py +++ b/tests/models/litter/conftest.py @@ -67,12 +67,14 @@ def dummy_litter_data(fixture_core_components): # Vertically structured variables data["soil_temperature"] = lyr_strct.from_template() - data["soil_temperature"][lyr_strct.index_all_soil] = 20 + data["soil_temperature"][lyr_strct.index_topsoil] = 20 + data["soil_temperature"][lyr_strct.index_subsoil] = [19.5, 18.7, 18.7, 17.6] # At present the soil model only uses the top soil layer, so this is the # only one with real test values in data["matric_potential"] = lyr_strct.from_template() data["matric_potential"][lyr_strct.index_topsoil] = [-10.0, -25.0, -100.0, -100.0] + data["matric_potential"][lyr_strct.index_subsoil] = [-11.0, -29.5, -123.0, -154.1] data["air_temperature"] = lyr_strct.from_template() data["air_temperature"][lyr_strct.index_filled_atmosphere] = np.array( diff --git a/tests/models/litter/test_env_factors.py b/tests/models/litter/test_env_factors.py index f8eb4df30..d375155f3 100644 --- a/tests/models/litter/test_env_factors.py +++ b/tests/models/litter/test_env_factors.py @@ -1,6 +1,7 @@ """Test module for litter.env_factors.py.""" import numpy as np +import pytest from virtual_ecosystem.models.litter.constants import LitterConsts @@ -27,21 +28,151 @@ def test_calculate_temperature_effect_on_litter_decomp( assert np.allclose(actual_factor, expected_factor) -def test_calculate_soil_water_effect_on_litter_decomp(): +def test_calculate_soil_water_effect_on_litter_decomp( + dummy_litter_data, fixture_core_components +): """Test that soil moisture effects on decomposition are calculated correctly.""" from virtual_ecosystem.models.litter.env_factors import ( calculate_soil_water_effect_on_litter_decomp, ) - water_potentials = np.array([-10.0, -25.0, -100.0, -400.0]) - - expected_factor = [1.0, 0.88496823, 0.71093190, 0.53689556] + expected_factor = [1.0, 0.88496823, 0.71093190, 0.71093190] actual_factor = calculate_soil_water_effect_on_litter_decomp( - water_potentials, + water_potential=dummy_litter_data["matric_potential"][ + fixture_core_components.layer_structure.index_topsoil_scalar + ], water_potential_halt=LitterConsts.litter_decay_water_potential_halt, water_potential_opt=LitterConsts.litter_decay_water_potential_optimum, moisture_response_curvature=LitterConsts.moisture_response_curvature, ) assert np.allclose(actual_factor, expected_factor) + + +@pytest.mark.parametrize( + "increased_depth,expected_av_temps", + [ + pytest.param( + True, + [18.6319817, 18.498648, 18.498648, 18.315315], + id="increased depth", + ), + pytest.param( + False, + [18.0729725, 18.0729725, 18.0729725, 18.0729725], + id="normal depth", + ), + ], +) +def test_average_temperature_over_microbially_active_layers( + dummy_litter_data, fixture_core_components, increased_depth, expected_av_temps +): + """Check averaging of temperatures over soil layers works correctly.""" + from virtual_ecosystem.models.litter.env_factors import ( + average_temperature_over_microbially_active_layers, + ) + + if increased_depth: + fixture_core_components.layer_structure.soil_layer_active_thickness = np.array( + [0.5, 0.25] + ) + fixture_core_components.layer_structure.max_depth_of_microbial_activity = 0.75 + + actual_av_temps = average_temperature_over_microbially_active_layers( + soil_temperatures=dummy_litter_data["soil_temperature"], + surface_temperature=dummy_litter_data["air_temperature"][ + fixture_core_components.layer_structure.index_surface + ].to_numpy(), + layer_structure=fixture_core_components.layer_structure, + ) + + assert np.allclose(actual_av_temps, expected_av_temps) + + +@pytest.mark.parametrize( + "increased_depth,expected_water_pots", + [ + pytest.param( + True, + [-10.1667, -25.750, -103.8333, -109.0167], + id="increased depth", + ), + pytest.param( + False, + [-10.0, -25.0, -100.0, -100.0], + id="normal depth", + ), + ], +) +def test_average_water_potential_over_microbially_active_layers( + dummy_litter_data, fixture_core_components, increased_depth, expected_water_pots +): + """Check averaging of water potentials over soil layers works correctly.""" + from virtual_ecosystem.models.litter.env_factors import ( + average_water_potential_over_microbially_active_layers, + ) + + if increased_depth: + fixture_core_components.layer_structure.soil_layer_active_thickness = np.array( + [0.5, 0.25] + ) + fixture_core_components.layer_structure.max_depth_of_microbial_activity = 0.75 + + actual_water_pots = average_water_potential_over_microbially_active_layers( + water_potentials=dummy_litter_data["matric_potential"], + layer_structure=fixture_core_components.layer_structure, + ) + + assert np.allclose(actual_water_pots, expected_water_pots) + + +@pytest.mark.parametrize( + "increased_depth,expected_factors", + [ + pytest.param( + True, + { + "temp_above": [0.1878681, 0.1878681, 0.1878681, 0.1878681], + "temp_below": [0.2407699, 0.2377353, 0.2377353, 0.2335993], + "water": [0.9979245, 0.8812574, 0.7062095, 0.7000939], + }, + id="increased depth", + ), + pytest.param( + False, + { + "temp_above": [0.1878681, 0.1878681, 0.1878681, 0.1878681], + "temp_below": [0.2281971, 0.2281971, 0.2281971, 0.2281971], + "water": [1.0, 0.88496823, 0.71093190, 0.71093190], + }, + id="normal depth", + ), + ], +) +def test_calculate_environmental_factors( + dummy_litter_data, fixture_core_components, increased_depth, expected_factors +): + """Check that the calculation of the relevant environmental factors is correct.""" + from virtual_ecosystem.models.litter.env_factors import ( + calculate_environmental_factors, + ) + + if increased_depth: + fixture_core_components.layer_structure.soil_layer_active_thickness = np.array( + [0.5, 0.25] + ) + fixture_core_components.layer_structure.max_depth_of_microbial_activity = 0.75 + + actual_factors = calculate_environmental_factors( + air_temperatures=dummy_litter_data["air_temperature"], + soil_temperatures=dummy_litter_data["soil_temperature"], + water_potentials=dummy_litter_data["matric_potential"], + layer_structure=fixture_core_components.layer_structure, + constants=LitterConsts, + ) + + assert set(expected_factors.keys()) == set(actual_factors.keys()) + + for key in actual_factors.keys(): + assert np.allclose(actual_factors[key], expected_factors[key]) diff --git a/tests/models/litter/test_litter_model.py b/tests/models/litter/test_litter_model.py index 2ac06753a..f0300fd14 100644 --- a/tests/models/litter/test_litter_model.py +++ b/tests/models/litter/test_litter_model.py @@ -311,12 +311,12 @@ def test_update(fixture_litter_model, dummy_litter_data): end_above_meta = [0.32072786, 0.15473132, 0.08523907, 0.08074153] end_above_struct = [0.50470382, 0.25068224, 0.09843778, 0.11163532] end_woody = [4.7745168, 11.89872931, 7.3614112, 7.3314112] - end_below_meta = [0.4090768, 0.37287148, 0.06883228, 0.08315412] - end_below_struct = [0.6066315, 0.31860251, 0.02010566, 0.03038382] + end_below_meta = [0.41087696, 0.37434507, 0.06905624, 0.08337808] + end_below_struct = [0.6066914, 0.31869812, 0.02010607, 0.03038423] end_lignin_above_struct = [0.49790843, 0.10067782, 0.70495536, 0.71045831] end_lignin_woody = [0.49580586, 0.79787834, 0.35224223, 0.35012603] - end_lignin_below_struct = [0.50313604, 0.2658639, 0.7499951, 0.82142894] - c_mineral = [0.02987233, 0.02316114, 0.00786517, 0.00786517] + end_lignin_below_struct = [0.50313573, 0.26585915, 0.7499951, 0.82142798] + c_mineral = [0.02652423, 0.02033658, 0.00746131, 0.00746131] fixture_litter_model.update(time_index=0) diff --git a/tests/models/litter/test_litter_pools.py b/tests/models/litter/test_litter_pools.py index 40dd0d962..5d81c0a0d 100644 --- a/tests/models/litter/test_litter_pools.py +++ b/tests/models/litter/test_litter_pools.py @@ -96,8 +96,8 @@ def test_calculate_decay_rates(dummy_litter_data, fixture_core_components): "metabolic_above": [0.00450883, 0.00225442, 0.00105206, 0.00105206], "structural_above": [1.6742967e-4, 6.1857359e-4, 1.1086908e-5, 1.1086908e-5], "woody": [0.0004832, 0.00027069, 0.0015888, 0.0015888], - "metabolic_below": [0.01092804, 0.00894564, 0.00135959, 0.00135959], - "structural_below": [3.6365995e-04, 5.803657e-04, 2.469074e-06, 2.469074e-06], + "metabolic_below": [0.00912788, 0.00747205, 0.00113563, 0.00113563], + "structural_below": [3.0375501e-4, 4.8476324e-4, 2.0623487e-6, 2.0623487e-6], } actual_decay = calculate_decay_rates( @@ -109,15 +109,10 @@ def test_calculate_decay_rates(dummy_litter_data, fixture_core_components): lignin_above_structural=dummy_litter_data["lignin_above_structural"].to_numpy(), lignin_woody=dummy_litter_data["lignin_woody"].to_numpy(), lignin_below_structural=dummy_litter_data["lignin_below_structural"].to_numpy(), - surface_temp=dummy_litter_data["air_temperature"][ - fixture_core_components.layer_structure.index_surface_scalar - ].to_numpy(), - topsoil_temp=dummy_litter_data["soil_temperature"][ - fixture_core_components.layer_structure.index_topsoil_scalar - ].to_numpy(), - water_potential=dummy_litter_data["matric_potential"][ - fixture_core_components.layer_structure.index_topsoil_scalar - ].to_numpy(), + air_temperatures=dummy_litter_data["air_temperature"], + soil_temperatures=dummy_litter_data["soil_temperature"], + water_potentials=dummy_litter_data["matric_potential"], + layer_structure=fixture_core_components.layer_structure, constants=LitterConsts, ) diff --git a/virtual_ecosystem/models/litter/env_factors.py b/virtual_ecosystem/models/litter/env_factors.py index fc5cc982a..27ba1eaf8 100644 --- a/virtual_ecosystem/models/litter/env_factors.py +++ b/virtual_ecosystem/models/litter/env_factors.py @@ -5,6 +5,91 @@ import numpy as np from numpy.typing import NDArray +from xarray import DataArray + +from virtual_ecosystem.core.core_components import LayerStructure +from virtual_ecosystem.models.litter.constants import LitterConsts + + +def calculate_environmental_factors( + air_temperatures: DataArray, + soil_temperatures: DataArray, + water_potentials: DataArray, + layer_structure: LayerStructure, + constants: LitterConsts, +): + """Calculate the impact of the environment has on litter decay across litter layers. + + For the above ground layer the impact of temperature is calculated, and for the + below ground layer the effect of temperature and soil water potential are both + calculated. + + The relevant above ground temperature is the surface temperature, which can be + easily extracted from the temperature data. It's more complex for the below ground + temperature and the water potential as the relevant values are averages across the + microbially active depth. These are calculated by averaging across the soil layers + with each layer weighted by the proportion of the total microbially active depth it + represents. + + If a shallow microbially active depth is used then below ground litter decomposition + will be exposed to a high degree of environmental variability. This is + representative of the real world, but needs to be kept in mind when comparing to + other models. + + Args: + air_temperatures: Air temperatures, for all above ground layers [C] + soil_temperatures: Soil temperatures, for all soil layers [C] + water_potentials: Water potentials, for all soil layers [kPa] + layer_structure: The LayerStructure instance for the simulation. + constants: Set of constants for the litter model + + Returns: + A dictionary containing three environmental factors, one for the effect of + temperature on above ground litter decay, one for the effect of temperature on + below ground litter decay, and one for the effect of soil water potential on + below ground litter decay. + """ + + temperatures = { + "surface": air_temperatures[layer_structure.index_surface_scalar].to_numpy(), + # TODO - This currently takes uses the surface temperature for the first layer. + # Once we start change the default to use a thin topsoil layer that should be + # used here instead + "below_ground": average_temperature_over_microbially_active_layers( + soil_temperatures=soil_temperatures, + surface_temperature=air_temperatures[ + layer_structure.index_surface_scalar + ].to_numpy(), + layer_structure=layer_structure, + ), + } + water_potential = average_water_potential_over_microbially_active_layers( + water_potentials=water_potentials, layer_structure=layer_structure + ) + + temperature_factors = { + level: calculate_temperature_effect_on_litter_decomp( + temperature=temp, + reference_temp=constants.litter_decomp_reference_temp, + offset_temp=constants.litter_decomp_offset_temp, + temp_response=constants.litter_decomp_temp_response, + ) + for (level, temp) in temperatures.items() + } + + # Calculate the water factor (relevant for below ground layers) + water_factor = calculate_soil_water_effect_on_litter_decomp( + water_potential=water_potential, + water_potential_halt=constants.litter_decay_water_potential_halt, + water_potential_opt=constants.litter_decay_water_potential_optimum, + moisture_response_curvature=constants.moisture_response_curvature, + ) + + return { + "temp_above": temperature_factors["surface"], + "temp_below": temperature_factors["below_ground"], + "water": water_factor, + } def calculate_temperature_effect_on_litter_decomp( @@ -69,3 +154,91 @@ def calculate_soil_water_effect_on_litter_decomp( ) ** moisture_response_curvature return 1 - supression + + +def average_temperature_over_microbially_active_layers( + soil_temperatures: DataArray, + surface_temperature: NDArray[np.float32], + layer_structure: LayerStructure, +) -> NDArray[np.float32]: + """Average soil temperatures over the microbially active layers. + + First the average temperature is found for each layer. Then an average across the + microbially active depth is taken, weighting by how much of the microbially active + depth lies within each layer. + + Args: + soil_temperatures: Soil temperatures to be averaged [C] + surface_temperature: Air temperature just above the soil surface [C] + layer_structure: The LayerStructure instance for the simulation. + + Returns: + The average temperature across the soil depth considered to be microbially + active [C] + """ + + # Find weighting for each layer in the average by dividing the microbially active + # depth in each layer by the total depth of microbial activity + layer_weights = ( + layer_structure.soil_layer_active_thickness + / layer_structure.max_depth_of_microbial_activity + ) + + # Find the average for each layer + layer_averages = np.empty((layer_weights.shape[0], soil_temperatures.shape[1])) + layer_averages[0, :] = ( + surface_temperature + soil_temperatures[layer_structure.index_topsoil] + ) / 2.0 + + for index in range(1, len(layer_structure.soil_layer_active_thickness)): + layer_averages[index, :] = ( + soil_temperatures[layer_structure.index_topsoil_scalar + index - 1] + + soil_temperatures[layer_structure.index_topsoil_scalar + index] + ) / 2.0 + + return np.dot(layer_weights, layer_averages) + + +def average_water_potential_over_microbially_active_layers( + water_potentials: DataArray, + layer_structure: LayerStructure, +) -> NDArray[np.float32]: + """Average water potentials over the microbially active layers. + + The average water potential is found for each layer apart from the top layer. This + is because for the top layer a sensible average can't be taken as water potential is + not defined for the surface layer. In this case, the water potential at the maximum + layer height is just treated as the average of the layer. This is a reasonable + assumption if the first soil layer is shallow. + + These water potentials are then averaged across the microbially active depth, + weighting by how much of the microbially active depth lies within each layer. + + Args: + water_potentials: Soil water potentials to be averaged [kPa] + layer_structure: The LayerStructure instance for the simulation. + + Returns: + The average water potential across the soil depth considered to be microbially + active [kPa] + """ + + # Find weighting for each layer in the average by dividing the microbially active + # depth in each layer by the total depth of microbial activity + layer_weights = ( + layer_structure.soil_layer_active_thickness + / layer_structure.max_depth_of_microbial_activity + ) + + # Find the average for each layer + layer_averages = np.empty((layer_weights.shape[0], water_potentials.shape[1])) + # Top layer cannot be averaged + layer_averages[0, :] = water_potentials[layer_structure.index_topsoil] + + for index in range(1, len(layer_structure.soil_layer_active_thickness)): + layer_averages[index, :] = ( + water_potentials[layer_structure.index_topsoil_scalar + index - 1] + + water_potentials[layer_structure.index_topsoil_scalar + index] + ) / 2.0 + + return np.dot(layer_weights, layer_averages) diff --git a/virtual_ecosystem/models/litter/litter_model.py b/virtual_ecosystem/models/litter/litter_model.py index aad51932b..3e64414e6 100644 --- a/virtual_ecosystem/models/litter/litter_model.py +++ b/virtual_ecosystem/models/litter/litter_model.py @@ -228,15 +228,10 @@ def update(self, time_index: int, **kwargs: Any) -> None: lignin_above_structural=self.data["lignin_above_structural"].to_numpy(), lignin_woody=self.data["lignin_woody"].to_numpy(), lignin_below_structural=self.data["lignin_below_structural"].to_numpy(), - surface_temp=self.data["air_temperature"][ - self.layer_structure.index_surface_scalar - ].to_numpy(), - topsoil_temp=self.data["soil_temperature"][ - self.layer_structure.index_topsoil_scalar - ].to_numpy(), - water_potential=self.data["matric_potential"][ - self.layer_structure.index_topsoil_scalar - ].to_numpy(), + air_temperatures=self.data["air_temperature"], + soil_temperatures=self.data["soil_temperature"], + water_potentials=self.data["matric_potential"], + layer_structure=self.layer_structure, constants=self.model_constants, ) diff --git a/virtual_ecosystem/models/litter/litter_pools.py b/virtual_ecosystem/models/litter/litter_pools.py index 4b8838c4c..0b43b08e7 100644 --- a/virtual_ecosystem/models/litter/litter_pools.py +++ b/virtual_ecosystem/models/litter/litter_pools.py @@ -18,12 +18,13 @@ import numpy as np from numpy.typing import NDArray +from xarray import DataArray from virtual_ecosystem.core.constants import CoreConsts +from virtual_ecosystem.core.core_components import LayerStructure from virtual_ecosystem.models.litter.constants import LitterConsts from virtual_ecosystem.models.litter.env_factors import ( - calculate_soil_water_effect_on_litter_decomp, - calculate_temperature_effect_on_litter_decomp, + calculate_environmental_factors, ) @@ -36,9 +37,10 @@ def calculate_decay_rates( lignin_above_structural: NDArray[np.float32], lignin_woody: NDArray[np.float32], lignin_below_structural: NDArray[np.float32], - surface_temp: NDArray[np.float32], - topsoil_temp: NDArray[np.float32], - water_potential: NDArray[np.float32], + air_temperatures: DataArray, + soil_temperatures: DataArray, + water_potentials: DataArray, + layer_structure: LayerStructure, constants: LitterConsts, ) -> dict[str, NDArray[np.float32]]: """Calculate the decay rate for all five of the litter pools. @@ -54,70 +56,59 @@ def calculate_decay_rates( lignin_woody: Proportion of dead wood pool which is lignin [unitless] lignin_below_structural: Proportion of below ground structural pool which is lignin [unitless] - surface_temp: Temperature of soil surface, which is assumed to be the same - temperature as the above ground litter [C] - topsoil_temp: Temperature of topsoil layer, which is assumed to be the same - temperature as the below ground litter [C] - water_potential: Water potential of the topsoil layer [kPa] + air_temperatures: Air temperatures, for all above ground layers [C] + soil_temperatures: Soil temperatures, for all soil layers [C] + water_potentials: Water potentials, for all soil layers [kPa] + layer_structure: The LayerStructure instance for the simulation. constants: Set of constants for the litter model + Decay rates depend on lignin proportions as well as a range of environmental + factors. These environmental factors are calculated as part of this function. + Returns: A dictionary containing the decay rate for each of the five litter pools. """ - # Calculate temperature factor for the above ground litter layers - temperature_factor_above = calculate_temperature_effect_on_litter_decomp( - temperature=surface_temp, - reference_temp=constants.litter_decomp_reference_temp, - offset_temp=constants.litter_decomp_offset_temp, - temp_response=constants.litter_decomp_temp_response, - ) - # Calculate temperature factor for the below ground litter layers - temperature_factor_below = calculate_temperature_effect_on_litter_decomp( - temperature=topsoil_temp, - reference_temp=constants.litter_decomp_reference_temp, - offset_temp=constants.litter_decomp_offset_temp, - temp_response=constants.litter_decomp_temp_response, - ) - # Calculate the water factor (relevant for below ground layers) - water_factor = calculate_soil_water_effect_on_litter_decomp( - water_potential=water_potential, - water_potential_halt=constants.litter_decay_water_potential_halt, - water_potential_opt=constants.litter_decay_water_potential_optimum, - moisture_response_curvature=constants.moisture_response_curvature, + # Calculate environmental factors + env_factors = calculate_environmental_factors( + air_temperatures=air_temperatures, + soil_temperatures=soil_temperatures, + water_potentials=water_potentials, + layer_structure=layer_structure, + constants=constants, ) # Calculate decay rate for each pool metabolic_above_decay = calculate_litter_decay_metabolic_above( - temperature_factor_above, - above_metabolic, + temperature_factor=env_factors["temp_above"], + litter_pool_above_metabolic=above_metabolic, litter_decay_coefficient=constants.litter_decay_constant_metabolic_above, ) structural_above_decay = calculate_litter_decay_structural_above( - temperature_factor_above, - above_structural, - lignin_above_structural, + temperature_factor=env_factors["temp_above"], + litter_pool_above_structural=above_structural, + lignin_proportion=lignin_above_structural, litter_decay_coefficient=constants.litter_decay_constant_structural_above, lignin_inhibition_factor=constants.lignin_inhibition_factor, ) woody_decay = calculate_litter_decay_woody( - temperature_factor_above, - woody, - lignin_woody, + temperature_factor=env_factors["temp_above"], + litter_pool_woody=woody, + lignin_proportion=lignin_woody, litter_decay_coefficient=constants.litter_decay_constant_woody, lignin_inhibition_factor=constants.lignin_inhibition_factor, ) metabolic_below_decay = calculate_litter_decay_metabolic_below( - temperature_factor_below, - water_factor, - below_metabolic, + temperature_factor=env_factors["temp_below"], + moisture_factor=env_factors["water"], + litter_pool_below_metabolic=below_metabolic, litter_decay_coefficient=constants.litter_decay_constant_metabolic_below, ) structural_below_decay = calculate_litter_decay_structural_below( - temperature_factor_below, - water_factor, - below_structural, - lignin_below_structural, + temperature_factor=env_factors["temp_below"], + moisture_factor=env_factors["water"], + litter_pool_below_structural=below_structural, + lignin_proportion=lignin_below_structural, litter_decay_coefficient=constants.litter_decay_constant_structural_below, lignin_inhibition_factor=constants.lignin_inhibition_factor, )