diff --git a/tests/conftest.py b/tests/conftest.py index f9db7c36e..7bd1956ad 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -139,7 +139,10 @@ def dummy_carbon_data(layer_roles_fixture): DataArray(np.full((13, 4), np.nan), dims=["layers", "cell_id"]), # At present the soil model only uses the top soil layer, so this is the # only one with real test values in - DataArray([[0.5, 0.7, 0.6, 0.2]], dims=["layers", "cell_id"]), + DataArray( + [[0.472467929, 0.399900047, 0.256053640, 0.153616897]], + dims=["layers", "cell_id"], + ), DataArray(np.full((1, 4), np.nan), dims=["layers", "cell_id"]), ], dim="layers", @@ -151,6 +154,23 @@ def dummy_carbon_data(layer_roles_fixture): "cell_id": data.grid.cell_id, } ) + # TODO - Eventually this should replace the dummy soil moisture entirely + data["matric_potential"] = xr.concat( + [ + DataArray(np.full((13, 4), np.nan), dims=["layers", "cell_id"]), + # At present the soil model only uses the top soil layer, so this is the + # only one with real test values in + DataArray([[-3.0, -10.0, -250.0, -10000.0]], dims=["layers", "cell_id"]), + DataArray(np.full((1, 4), np.nan), dims=["layers", "cell_id"]), + ], + dim="layers", + ).assign_coords( + { + "layers": np.arange(0, 15), + "layer_roles": ("layers", layer_roles_fixture), + "cell_id": data.grid.cell_id, + } + ) data["soil_temperature"] = xr.concat( [ DataArray(np.full((13, 4), np.nan), dims=["dim_0", "cell_id"]), diff --git a/tests/models/litter/conftest.py b/tests/models/litter/conftest.py index 9d9d1c5cb..fa6e2586c 100644 --- a/tests/models/litter/conftest.py +++ b/tests/models/litter/conftest.py @@ -68,12 +68,15 @@ def dummy_litter_data(layer_roles_fixture): ) # The layer dependant data has to be handled separately - data["soil_moisture"] = concat( + data["matric_potential"] = concat( [ DataArray(np.full((13, 3), np.nan), dims=["layers", "cell_id"]), # At present the soil model only uses the top soil layer, so this is the # only one with real test values in - DataArray([[0.25, 0.45, 0.3]], dims=["layers", "cell_id"]), + DataArray( + [[-10.0, -25.0, -100.0]], + dims=["layers", "cell_id"], + ), DataArray(np.full((1, 3), np.nan), dims=["layers", "cell_id"]), ], dim="layers", diff --git a/tests/models/litter/test_litter_model.py b/tests/models/litter/test_litter_model.py index 6d2797657..790e5964b 100644 --- a/tests/models/litter/test_litter_model.py +++ b/tests/models/litter/test_litter_model.py @@ -441,12 +441,12 @@ def test_update(litter_model_fixture, dummy_litter_data): end_above_meta = [0.29587973, 0.14851276, 0.07041856] end_above_struct = [0.50055126, 0.25010012, 0.0907076] end_woody = [4.702103, 11.802315, 7.300997] - end_below_meta = [0.394145, 0.35923, 0.069006] - end_below_struct = [0.60027118, 0.30975403, 0.02047743] + end_below_meta = [0.38949196, 0.36147436, 0.06906041] + end_below_struct = [0.60011634, 0.30989963, 0.02047753] end_lignin_above_struct = [0.4996410, 0.1004310, 0.6964345] end_lignin_woody = [0.49989001, 0.79989045, 0.34998229] end_lignin_below_struct = [0.499760108, 0.249922519, 0.737107757] - c_mineral = [0.0212182, 0.02746286, 0.00796359] + c_mineral = [0.02987233, 0.02316114, 0.00786517] litter_model_fixture.update(time_index=0) @@ -468,23 +468,3 @@ def test_update(litter_model_fixture, dummy_litter_data): dummy_litter_data["lignin_below_structural"], end_lignin_below_struct ) assert np.allclose(dummy_litter_data["litter_C_mineralisation_rate"], c_mineral) - - -def test_convert_soil_moisture_to_water_potential( - dummy_litter_data, top_soil_layer_index -): - """Test that function to convert soil moisture to a water potential works.""" - from virtual_rainforest.models.litter.litter_model import ( - convert_soil_moisture_to_water_potential, - ) - - expected_potentials = [-297.14104, -4.2647655, -79.666189] - - actual_potentials = convert_soil_moisture_to_water_potential( - dummy_litter_data["soil_moisture"][top_soil_layer_index].to_numpy(), - air_entry_water_potential=LitterConsts.air_entry_water_potential, - water_retention_curvature=LitterConsts.water_retention_curvature, - saturated_water_content=LitterConsts.saturated_water_content, - ) - - assert np.allclose(actual_potentials, expected_potentials) diff --git a/tests/models/litter/test_litter_pools.py b/tests/models/litter/test_litter_pools.py index 3632989b7..447cc86de 100644 --- a/tests/models/litter/test_litter_pools.py +++ b/tests/models/litter/test_litter_pools.py @@ -14,30 +14,25 @@ def temp_and_water_factors( dummy_litter_data, surface_layer_index, top_soil_layer_index ): """Temperature and water factors for the various litter layers.""" - from virtual_rainforest.models.litter.litter_model import ( - convert_soil_moisture_to_water_potential, - ) from virtual_rainforest.models.litter.litter_pools import ( calculate_environmental_factors, ) - water_potentials = convert_soil_moisture_to_water_potential( - dummy_litter_data["soil_moisture"][top_soil_layer_index], - air_entry_water_potential=LitterConsts.air_entry_water_potential, - water_retention_curvature=LitterConsts.water_retention_curvature, - saturated_water_content=LitterConsts.saturated_water_content, - ) - environmental_factors = calculate_environmental_factors( surface_temp=dummy_litter_data["air_temperature"][surface_layer_index], topsoil_temp=dummy_litter_data["soil_temperature"][top_soil_layer_index], - water_potential=water_potentials, + water_potential=dummy_litter_data["matric_potential"][top_soil_layer_index], constants=LitterConsts, ) return environmental_factors +# TODO - Compare the below +# [-297.1410435034187, -4.264765510307134, -79.66618999943468] +# [-10.0, -25.0, -100.0] + + @pytest.fixture def decay_rates(dummy_litter_data, temp_and_water_factors): """Decay rates for the various litter pools.""" @@ -59,9 +54,6 @@ def test_calculate_environmental_factors( calculate_environmental_factors, ) - # TODO - hardcoding this in, but eventually this should be added to the data object. - water_potentials = [-10.0, -25.0, -100.0] - expected_water_factors = [1.0, 0.88496823, 0.71093190] expected_temp_above_factors = [0.1878681, 0.1878681, 0.1878681] expected_temp_below_factors = [0.2732009, 0.2732009, 0.2732009] @@ -69,7 +61,7 @@ def test_calculate_environmental_factors( environmental_factors = calculate_environmental_factors( surface_temp=dummy_litter_data["air_temperature"][surface_layer_index], topsoil_temp=dummy_litter_data["soil_temperature"][top_soil_layer_index], - water_potential=water_potentials, + water_potential=dummy_litter_data["matric_potential"][top_soil_layer_index], constants=LitterConsts, ) @@ -104,7 +96,7 @@ def test_calculate_moisture_effect_on_litter_decomp(top_soil_layer_index): calculate_moisture_effect_on_litter_decomp, ) - water_potentials = [-10.0, -25.0, -100.0, -400.0] + water_potentials = np.array([-10.0, -25.0, -100.0, -400.0]) expected_factor = [1.0, 0.88496823, 0.71093190, 0.53689556] @@ -139,9 +131,6 @@ def test_calculate_change_in_litter_variables( dummy_litter_data, surface_layer_index, top_soil_layer_index ): """Test that litter pool update calculation is correct.""" - from virtual_rainforest.models.litter.litter_model import ( - convert_soil_moisture_to_water_potential, - ) from virtual_rainforest.models.litter.litter_pools import ( calculate_change_in_litter_variables, ) @@ -150,22 +139,14 @@ def test_calculate_change_in_litter_variables( "litter_pool_above_metabolic": [0.29587973, 0.14851276, 0.07041856], "litter_pool_above_structural": [0.50055126, 0.25010012, 0.0907076], "litter_pool_woody": [4.702103, 11.802315, 7.300997], - "litter_pool_below_metabolic": [0.394145, 0.35923, 0.069006], - "litter_pool_below_structural": [0.60027118, 0.30975403, 0.02047743], + "litter_pool_below_metabolic": [0.38949196, 0.36147436, 0.06906041], + "litter_pool_below_structural": [0.60011634, 0.30989963, 0.02047753], "lignin_above_structural": [0.4996410, 0.1004310, 0.6964345], "lignin_woody": [0.49989001, 0.79989045, 0.34998229], "lignin_below_structural": [0.499760108, 0.249922519, 0.737107757], - "litter_C_mineralisation_rate": [0.0212182, 0.02746286, 0.00796359], + "litter_C_mineralisation_rate": [0.02987233, 0.02316114, 0.00786517], } - # Calculate water potential - water_potential = convert_soil_moisture_to_water_potential( - dummy_litter_data["soil_moisture"][top_soil_layer_index].to_numpy(), - air_entry_water_potential=LitterConsts.air_entry_water_potential, - water_retention_curvature=LitterConsts.water_retention_curvature, - saturated_water_content=LitterConsts.saturated_water_content, - ) - result = calculate_change_in_litter_variables( surface_temp=dummy_litter_data["air_temperature"][ surface_layer_index @@ -173,7 +154,9 @@ def test_calculate_change_in_litter_variables( topsoil_temp=dummy_litter_data["soil_temperature"][ top_soil_layer_index ].to_numpy(), - water_potential=water_potential, + water_potential=dummy_litter_data["matric_potential"][ + top_soil_layer_index + ].to_numpy(), above_metabolic=dummy_litter_data["litter_pool_above_metabolic"].to_numpy(), above_structural=dummy_litter_data["litter_pool_above_structural"].to_numpy(), woody=dummy_litter_data["litter_pool_woody"].to_numpy(), @@ -200,8 +183,8 @@ def test_calculate_decay_rates(dummy_litter_data, temp_and_water_factors): "metabolic_above": [0.00450883, 0.00225442, 0.00105206], "structural_above": [1.67429665e-4, 6.18573593e-4, 1.10869077e-5], "woody": [0.0004832, 0.00027069, 0.0015888], - "metabolic_below": [0.00627503, 0.01118989, 0.00141417], - "structural_below": [2.08818451e-4, 7.25965456e-4, 2.56818870e-6], + "metabolic_below": [0.01092804, 0.00894564, 0.00135959], + "structural_below": [3.63659952e-04, 5.80365659e-04, 2.46907410e-06], } actual_decay = calculate_decay_rates( @@ -362,7 +345,7 @@ def test_calculate_litter_decay_metabolic_below( calculate_litter_decay_metabolic_below, ) - expected_decay = [0.00627503, 0.01118989, 0.00141417] + expected_decay = [0.01092804, 0.00894564, 0.00135959] actual_decay = calculate_litter_decay_metabolic_below( temperature_factor=temp_and_water_factors["temp_below"], @@ -382,7 +365,7 @@ def test_calculate_litter_decay_structural_below( calculate_litter_decay_structural_below, ) - expected_decay = [2.08818451e-04, 7.25965456e-04, 2.56818870e-06] + expected_decay = [3.63659952e-04, 5.80365659e-04, 2.46907410e-06] actual_decay = calculate_litter_decay_structural_below( temperature_factor=temp_and_water_factors["temp_below"], diff --git a/tests/models/soil/test_carbon.py b/tests/models/soil/test_carbon.py index a700f2259..8c7da068a 100644 --- a/tests/models/soil/test_carbon.py +++ b/tests/models/soil/test_carbon.py @@ -14,41 +14,48 @@ @pytest.fixture -def moist_temp_scalars(dummy_carbon_data, top_soil_layer_index): - """Combined moisture and temperature scalars based on dummy carbon data.""" - from virtual_rainforest.models.soil.carbon import ( - convert_moisture_to_scalar, - convert_temperature_to_scalar, - ) +def moist_scalars(dummy_carbon_data, top_soil_layer_index): + """Moisture scalars based on dummy carbon data.""" + from virtual_rainforest.models.soil.carbon import convert_moisture_to_scalar moist_scalars = convert_moisture_to_scalar( - dummy_carbon_data["soil_moisture"][top_soil_layer_index], + np.array([0.5, 0.7, 0.6, 0.2]), SoilConsts.moisture_scalar_coefficient, SoilConsts.moisture_scalar_exponent, ) - temp_scalars = convert_temperature_to_scalar( - dummy_carbon_data["soil_temperature"][top_soil_layer_index], - SoilConsts.temp_scalar_coefficient_1, - SoilConsts.temp_scalar_coefficient_2, - SoilConsts.temp_scalar_coefficient_3, - SoilConsts.temp_scalar_coefficient_4, - SoilConsts.temp_scalar_reference_temp, + + return moist_scalars + + +@pytest.fixture +def water_factors(dummy_carbon_data, top_soil_layer_index): + """Water potential factors based on dummy carbon data.""" + from virtual_rainforest.models.soil.carbon import ( + calculate_water_potential_impact_on_microbes, + ) + + moist_scalars = calculate_water_potential_impact_on_microbes( + water_potential=dummy_carbon_data["matric_potential"][top_soil_layer_index], + water_potential_halt=SoilConsts.soil_microbe_water_potential_halt, + water_potential_opt=SoilConsts.soil_microbe_water_potential_optimum, + moisture_response_curvature=SoilConsts.moisture_response_curvature, ) - return moist_scalars * temp_scalars + return moist_scalars def test_top_soil_data_extraction(dummy_carbon_data, top_soil_layer_index): """Test that top soil data can be extracted from the data object correctly.""" top_soil_temps = [35.0, 37.5, 40.0, 25.0] - top_soil_moistures = [0.5, 0.7, 0.6, 0.2] + top_soil_water_potentials = [-3.0, -10.0, -250.0, -10000.0] assert np.allclose( dummy_carbon_data["soil_temperature"][top_soil_layer_index], top_soil_temps ) assert np.allclose( - dummy_carbon_data["soil_moisture"][top_soil_layer_index], top_soil_moistures + dummy_carbon_data["matric_potential"][top_soil_layer_index], + top_soil_water_potentials, ) @@ -58,10 +65,10 @@ def test_calculate_soil_carbon_updates(dummy_carbon_data, top_soil_layer_index): from virtual_rainforest.models.soil.carbon import calculate_soil_carbon_updates change_in_pools = { - "soil_c_pool_lmwc": [0.00509349, 0.02038796, 0.71353723, 0.00125521], - "soil_c_pool_maom": [0.13088391, 0.05654771, -0.39962841, 0.00533357], - "soil_c_pool_microbe": [-0.33131188, -0.16636299, -0.76078599, -0.01275669], - "soil_c_pool_pom": [-0.00168464, -0.00936813, -0.00873008, 0.00403342], + "soil_c_pool_lmwc": [0.25847352, 0.59685734, 0.97122869, 0.02112077], + "soil_c_pool_maom": [0.10296645, 0.04445693, -0.31401747, 0.00422143], + "soil_c_pool_microbe": [-0.16002809, -0.07621711, -0.36452355, -0.01140071], + "soil_c_pool_pom": [-0.25377786, -0.58704913, -0.41245236, -0.01566563], } # Make order of pools object @@ -70,17 +77,20 @@ def test_calculate_soil_carbon_updates(dummy_carbon_data, top_soil_layer_index): pool_order[pool] = np.array([]) delta_pools = calculate_soil_carbon_updates( - dummy_carbon_data["soil_c_pool_lmwc"].to_numpy(), - dummy_carbon_data["soil_c_pool_maom"].to_numpy(), - dummy_carbon_data["soil_c_pool_microbe"].to_numpy(), - dummy_carbon_data["soil_c_pool_pom"].to_numpy(), - dummy_carbon_data["pH"], - dummy_carbon_data["bulk_density"], - dummy_carbon_data["soil_moisture"][top_soil_layer_index], - dummy_carbon_data["soil_temperature"][top_soil_layer_index], - dummy_carbon_data["percent_clay"], - dummy_carbon_data["litter_C_mineralisation_rate"], - pool_order, + soil_c_pool_lmwc=dummy_carbon_data["soil_c_pool_lmwc"].to_numpy(), + soil_c_pool_maom=dummy_carbon_data["soil_c_pool_maom"].to_numpy(), + soil_c_pool_microbe=dummy_carbon_data["soil_c_pool_microbe"].to_numpy(), + soil_c_pool_pom=dummy_carbon_data["soil_c_pool_pom"].to_numpy(), + pH=dummy_carbon_data["pH"], + bulk_density=dummy_carbon_data["bulk_density"], + soil_moisture=np.array([0.5, 0.7, 0.6, 0.2]), + soil_water_potential=dummy_carbon_data["matric_potential"][ + top_soil_layer_index + ].to_numpy(), + soil_temp=dummy_carbon_data["soil_temperature"][top_soil_layer_index], + percent_clay=dummy_carbon_data["percent_clay"], + mineralisation_rate=dummy_carbon_data["litter_C_mineralisation_rate"], + delta_pools_ordered=pool_order, constants=SoilConsts, ) @@ -90,12 +100,12 @@ def test_calculate_soil_carbon_updates(dummy_carbon_data, top_soil_layer_index): assert np.allclose(delta_pools[i * 4 : (i + 1) * 4], change_in_pools[pool]) -def test_calculate_mineral_association(dummy_carbon_data, moist_temp_scalars): +def test_calculate_mineral_association(dummy_carbon_data, moist_scalars): """Test that mineral_association runs and generates the correct values.""" from virtual_rainforest.models.soil.carbon import calculate_mineral_association - output_l_to_m = [-7.35822655e-03, -1.27716013e-02, -7.16245859e-01, 3.29436494e-05] + output_l_to_m = [-5.78872135e-03, -1.00408341e-02, -5.62807109e-01, 2.60743689e-05] # Then calculate mineral association rate lmwc_to_maom = calculate_mineral_association( @@ -103,7 +113,7 @@ def test_calculate_mineral_association(dummy_carbon_data, moist_temp_scalars): dummy_carbon_data["soil_c_pool_maom"], dummy_carbon_data["pH"], dummy_carbon_data["bulk_density"], - moist_temp_scalars, + moist_scalars, dummy_carbon_data["percent_clay"], constants=SoilConsts, ) @@ -199,22 +209,30 @@ def test_calculate_binding_coefficient(dummy_carbon_data): assert np.allclose(binding_coefs, output_coefs) -def test_convert_temperature_to_scalar(dummy_carbon_data, top_soil_layer_index): - """Test that scalar_temperature runs and generates the correct value.""" - from virtual_rainforest.models.soil.carbon import convert_temperature_to_scalar - - output_scalars = [1.27113, 1.27196, 1.27263, 1.26344] +@pytest.mark.parametrize( + "activation_energy,expected_factors", + [ + (30000.0, [2.57153601, 2.82565326, 3.10021393, 1.73629781]), + (45000.0, [4.12371761, 4.74983258, 5.45867825, 2.28789625]), + (57000.0, [6.01680536, 7.19657491, 8.58309980, 2.85289648]), + ], +) +def calculate_temperature_effect_on_microbes( + dummy_carbon_data, top_soil_layer_index, activation_energy, expected_factors +): + """Test function to calculate microbial temperature response.""" + from virtual_rainforest.models.soil.carbon import ( + calculate_temperature_effect_on_microbes, + ) - temp_scalar = convert_temperature_to_scalar( - dummy_carbon_data["soil_temperature"][top_soil_layer_index], - SoilConsts.temp_scalar_coefficient_1, - SoilConsts.temp_scalar_coefficient_2, - SoilConsts.temp_scalar_coefficient_3, - SoilConsts.temp_scalar_coefficient_4, - SoilConsts.temp_scalar_reference_temp, + actual_factors = calculate_temperature_effect_on_microbes( + soil_temperature=dummy_carbon_data["soil_temperature"][top_soil_layer_index], + activation_energy=activation_energy, + reference_temperature=SoilConsts.arrhenius_reference_temp, + gas_constant=SoilConsts.universal_gas_constant, ) - assert np.allclose(temp_scalar, output_scalars) + assert np.allclose(expected_factors, actual_factors) @pytest.mark.parametrize( @@ -257,7 +275,7 @@ def test_convert_moisture_to_scalar( ) else: moist_scalar = convert_moisture_to_scalar( - dummy_carbon_data["soil_moisture"][top_soil_layer_index], + np.array([0.5, 0.7, 0.6, 0.2]), SoilConsts.moisture_scalar_coefficient, SoilConsts.moisture_scalar_exponent, ) @@ -267,30 +285,52 @@ def test_convert_moisture_to_scalar( log_check(caplog, expected_log_entries) -def test_calculate_maintenance_respiration(dummy_carbon_data, moist_temp_scalars): +def test_calculate_water_potential_impact_on_microbes(): + """Test the calculation of the impact of soil water on microbial rates.""" + from virtual_rainforest.models.soil.carbon import ( + calculate_water_potential_impact_on_microbes, + ) + + water_potentials = np.array([-3.0, -10.0, -250.0, -10000.0]) + + expected_factor = [1.0, 0.94414168, 0.62176357, 0.07747536] + + actual_factor = calculate_water_potential_impact_on_microbes( + water_potential=water_potentials, + water_potential_halt=SoilConsts.soil_microbe_water_potential_halt, + water_potential_opt=SoilConsts.soil_microbe_water_potential_optimum, + moisture_response_curvature=SoilConsts.moisture_response_curvature, + ) + + assert np.allclose(actual_factor, expected_factor) + + +def test_calculate_maintenance_respiration( + dummy_carbon_data, moist_scalars, top_soil_layer_index +): """Check maintenance respiration cost calculates correctly.""" from virtual_rainforest.models.soil.carbon import calculate_maintenance_respiration - expected_resps = [0.19906823, 0.0998193, 0.45592854, 0.00763283] + expected_resps = [0.05443078, 0.02298407, 0.12012258, 0.00722288] main_resps = calculate_maintenance_respiration( - dummy_carbon_data["soil_c_pool_microbe"], - moist_temp_scalars, - SoilConsts.microbial_turnover_rate, + soil_c_pool_microbe=dummy_carbon_data["soil_c_pool_microbe"], + soil_temp=dummy_carbon_data["soil_temperature"][top_soil_layer_index], + constants=SoilConsts, ) assert np.allclose(main_resps, expected_resps) -def test_calculate_necromass_adsorption(dummy_carbon_data, moist_temp_scalars): +def test_calculate_necromass_adsorption(dummy_carbon_data, moist_scalars): """Check maintenance respiration cost calculates correctly.""" from virtual_rainforest.models.soil.carbon import calculate_necromass_adsorption - expected_adsorps = [0.13824183, 0.06931897, 0.31661708, 0.00530057] + expected_adsorps = [0.10875517, 0.05449776, 0.24878964, 0.00419536] actual_adsorps = calculate_necromass_adsorption( dummy_carbon_data["soil_c_pool_microbe"], - moist_temp_scalars, + moist_scalars, SoilConsts.necromass_adsorption_rate, ) @@ -359,49 +399,52 @@ def test_calculate_pom_decomposition_saturation(dummy_carbon_data): def test_calculate_microbial_carbon_uptake( - dummy_carbon_data, top_soil_layer_index, moist_temp_scalars + dummy_carbon_data, top_soil_layer_index, water_factors ): """Check microbial carbon uptake calculates correctly.""" from virtual_rainforest.models.soil.carbon import calculate_microbial_carbon_uptake - expected_uptake = [0.00599894, 0.00277614, 0.01176059, 0.00017683] + expected_uptake = [3.15786124e-3, 1.26472838e-3, 4.38868085e-3, 1.75270792e-5] actual_uptake = calculate_microbial_carbon_uptake( - dummy_carbon_data["soil_c_pool_lmwc"], - dummy_carbon_data["soil_c_pool_microbe"], - moist_temp_scalars, - dummy_carbon_data["soil_temperature"][top_soil_layer_index], + soil_c_pool_lmwc=dummy_carbon_data["soil_c_pool_lmwc"], + soil_c_pool_microbe=dummy_carbon_data["soil_c_pool_microbe"], + water_factor=water_factors, + soil_temp=dummy_carbon_data["soil_temperature"][top_soil_layer_index], constants=SoilConsts, ) assert np.allclose(actual_uptake, expected_uptake) -def test_calculate_labile_carbon_leaching(dummy_carbon_data, moist_temp_scalars): +def test_calculate_labile_carbon_leaching(dummy_carbon_data, moist_scalars): """Check leaching of labile carbon is calculated correctly.""" from virtual_rainforest.models.soil.carbon import calculate_labile_carbon_leaching - expected_leaching = [7.15045537e-05, 3.61665981e-05, 1.68115460e-04, 1.59018704e-06] + expected_leaching = [5.62526764e-05, 2.84336164e-05, 1.32100695e-04, 1.25860748e-06] actual_leaching = calculate_labile_carbon_leaching( dummy_carbon_data["soil_c_pool_lmwc"], - moist_temp_scalars, + moist_scalars, SoilConsts.leaching_rate_labile_carbon, ) assert np.allclose(actual_leaching, expected_leaching) -def test_calculate_pom_decomposition(dummy_carbon_data, moist_temp_scalars): +def test_calculate_pom_decomposition( + dummy_carbon_data, top_soil_layer_index, water_factors +): """Check that particulate organic matter decomposition is calculated correctly.""" from virtual_rainforest.models.soil.carbon import calculate_pom_decomposition - expected_decomp = [0.0038056940, 0.0104286084, 0.0092200655, 0.0014665616] + expected_decomp = [0.25589892, 0.58810966, 0.41294236, 0.02116563] actual_decomp = calculate_pom_decomposition( dummy_carbon_data["soil_c_pool_pom"], dummy_carbon_data["soil_c_pool_microbe"], - moist_temp_scalars, + water_factor=water_factors, + soil_temp=dummy_carbon_data["soil_temperature"][top_soil_layer_index], constants=SoilConsts, ) diff --git a/tests/models/soil/test_soil_model.py b/tests/models/soil/test_soil_model.py index ac8688620..ef8a8f4aa 100644 --- a/tests/models/soil/test_soil_model.py +++ b/tests/models/soil/test_soil_model.py @@ -222,7 +222,7 @@ def test_soil_model_initialization( }, }, pint.Quantity("12 hours"), - 0.01, + 0.2, does_not_raise(), ( ( @@ -369,6 +369,7 @@ def test_update(mocker, soil_model_fixture, dummy_carbon_data): end_lmwc = [0.04980117, 0.01999411, 0.09992829, 0.00499986] end_maom = [2.50019883, 1.70000589, 4.50007171, 0.50000014] end_microbe = [5.8, 2.3, 11.3, 1.0] + end_pom = [0.25, 2.34, 0.746, 0.3467] mock_integrate = mocker.patch.object(soil_model_fixture, "integrate") @@ -377,6 +378,7 @@ def test_update(mocker, soil_model_fixture, dummy_carbon_data): soil_c_pool_lmwc=DataArray(end_lmwc, dims="cell_id"), soil_c_pool_maom=DataArray(end_maom, dims="cell_id"), soil_c_pool_microbe=DataArray(end_microbe, dims="cell_id"), + soil_c_pool_pom=DataArray(end_pom, dims="cell_id"), ) ) @@ -388,6 +390,8 @@ def test_update(mocker, soil_model_fixture, dummy_carbon_data): # Check that data fixture has been updated correctly assert np.allclose(dummy_carbon_data["soil_c_pool_lmwc"], end_lmwc) assert np.allclose(dummy_carbon_data["soil_c_pool_maom"], end_maom) + assert np.allclose(dummy_carbon_data["soil_c_pool_microbe"], end_microbe) + assert np.allclose(dummy_carbon_data["soil_c_pool_pom"], end_pom) @pytest.mark.parametrize( @@ -399,17 +403,17 @@ def test_update(mocker, soil_model_fixture, dummy_carbon_data): Dataset( data_vars=dict( lmwc=DataArray( - [0.05283173, 0.03076934, 1.87751275, 0.00561615], dims="cell_id" + [0.13122264, 0.30432201, 0.50601794, 0.01541629], dims="cell_id" ), maom=DataArray( - [2.56409185, 1.72670677, 2.84168443, 0.50266489], dims="cell_id" + [2.54647825, 1.71347325, 4.32317632, 0.50157362], dims="cell_id" ), microbe=DataArray( - [5.63681284, 2.21869505, 10.96048678, 0.99364838], + [5.72558633, 2.27805994, 11.21111949, 0.99495352], dims="cell_id", ), pom=DataArray( - [0.09916254, 0.99531802, 0.69563758, 0.35201611], dims="cell_id" + [0.02068157, 0.71317448, 0.50004912, 0.34220329], dims="cell_id" ), ) ), @@ -475,6 +479,7 @@ def test_order_independance(dummy_carbon_data, soil_model_fixture): "pH", "bulk_density", "soil_moisture", + "matric_potential", "soil_temperature", "percent_clay", "litter_C_mineralisation_rate", @@ -516,24 +521,26 @@ def test_construct_full_soil_model(dummy_carbon_data, top_soil_layer_index): from virtual_rainforest.models.soil.soil_model import construct_full_soil_model - delta_pools = [ - 0.00509349, - 0.02038796, - 0.71353723, - 0.00125521, - 0.13088391, - 0.05654771, - -0.39962841, - 0.00533357, - -0.33131188, - -0.16636299, - -0.76078599, - -0.01275669, - -0.00168464, - -0.00936813, - -0.00873008, - 0.00403342, - ] + delta_pools = ( + [ + 0.25809707, + 0.59264789, + 0.56851013, + 0.02112901, + 0.0962045, + 0.02576618, + -0.08926844, + 0.0029497, + -0.15288599, + -0.05330495, + -0.18645947, + -0.01013683, + -0.25377786, + -0.58704913, + -0.41245236, + -0.01566563, + ], + ) # make pools pools = np.concatenate( @@ -557,7 +564,7 @@ def test_construct_full_soil_model(dummy_carbon_data, top_soil_layer_index): dummy_carbon_data, 4, top_soil_layer_index, - delta_pools_ordered, + delta_pools_ordered=delta_pools_ordered, constants=SoilConsts, ) diff --git a/virtual_rainforest/models/litter/constants.py b/virtual_rainforest/models/litter/constants.py index ca2fc4490..4d402dc34 100644 --- a/virtual_rainforest/models/litter/constants.py +++ b/virtual_rainforest/models/litter/constants.py @@ -162,38 +162,6 @@ class LitterConsts: being. """ - air_entry_water_potential: float = -3.815 - """Water potential at which soil pores begin to aerate [kPa]. - - The constant is used to estimate soil water potential from soil moisture. As this - estimation is a stopgap this constant probably shouldn't become a core constant. The - value is the average across all soil types found in - :cite:t:`cosby_statistical_1984`. In future, this could be calculated based on soil - texture. - - TODO: This estimation of water potential will either be moved to or superseded by - the hydrology model. So these constants should **not** remain here in the long term. - """ - - water_retention_curvature: float = -7.22 - """Curvature of the water retention curve. - - The value is the average across all soil types found in - :cite:t:`cosby_statistical_1984`; see documentation for - :attr:`air_entry_water_potential` for further details. - """ - - saturated_water_content: float = 0.457 - """Relative water content at which the soil is completely saturated [unitless]. - - The value is the average across all soil types found in - :cite:t:`cosby_statistical_1984`; see documentation for - :attr:`air_entry_water_potential` for further details. - - TODO: This constant already exists within the hydrology model, once the hydrology - model handles soil water potential, this constant **must** be deleted. - """ - lignin_inhibition_factor: float = -5.0 """Exponential factor expressing the extent that lignin inhibits litter breakdown. diff --git a/virtual_rainforest/models/litter/litter_model.py b/virtual_rainforest/models/litter/litter_model.py index 94d9f20da..fe0cf870e 100644 --- a/virtual_rainforest/models/litter/litter_model.py +++ b/virtual_rainforest/models/litter/litter_model.py @@ -36,7 +36,6 @@ class instance. If errors crop here when converting the information from the con from typing import Any import numpy as np -from numpy.typing import NDArray from pint import Quantity from xarray import DataArray @@ -97,6 +96,7 @@ class LitterModel(BaseModel): "lignin_above_structural", "lignin_woody", "lignin_below_structural", + "litter_C_mineralisation_rate", ) """Variables updated by the litter model.""" @@ -208,16 +208,6 @@ def update(self, time_index: int, **kwargs: Any) -> None: time_index: The index representing the current time step in the data object. """ - # Estimate water potentials based on soil moistures - water_potential = convert_soil_moisture_to_water_potential( - soil_moisture=self.data["soil_moisture"][ - self.top_soil_layer_index - ].to_numpy(), - air_entry_water_potential=self.constants.air_entry_water_potential, - water_retention_curvature=self.constants.water_retention_curvature, - saturated_water_content=self.constants.saturated_water_content, - ) - # Find change in litter variables using the function updated_variables = calculate_change_in_litter_variables( surface_temp=self.data["air_temperature"][ @@ -226,7 +216,9 @@ def update(self, time_index: int, **kwargs: Any) -> None: topsoil_temp=self.data["soil_temperature"][ self.top_soil_layer_index ].to_numpy(), - water_potential=water_potential, + water_potential=self.data["matric_potential"][ + self.top_soil_layer_index + ].to_numpy(), constants=self.constants, update_interval=self.update_interval.to("day").magnitude, above_metabolic=self.data["litter_pool_above_metabolic"].to_numpy(), @@ -252,35 +244,3 @@ def update(self, time_index: int, **kwargs: Any) -> None: def cleanup(self) -> None: """Placeholder function for litter model cleanup.""" - - -def convert_soil_moisture_to_water_potential( - soil_moisture: NDArray[np.float32], - air_entry_water_potential: float, - water_retention_curvature: float, - saturated_water_content: float, -) -> NDArray[np.float32]: - """Convert soil moisture into an estimate of water potential. - - This function provides a coarse estimate of soil water potential. It is taken from - :cite:t:`campbell_simple_1974`. - - TODO - This is a stopgap solution until we decide on a systematic way of handling - water potentials across the relevant models (`soil`, `litter`, `plants` and - `hydrology`). - - Args: - soil_moisture: Volumetric relative water content [unitless] - air_entry_water_potential: Water potential at which soil pores begin to aerate - [kPa] - water_retention_curvature: Curvature of water retention curve [unitless] - saturated_water_content: The relative water content at which the soil is fully - saturated [unitless]. - - Returns: - An estimate of the water potential of the soil [kPa] - """ - - return air_entry_water_potential * ( - (soil_moisture / saturated_water_content) ** water_retention_curvature - ) diff --git a/virtual_rainforest/models/litter/litter_pools.py b/virtual_rainforest/models/litter/litter_pools.py index 0b1e55001..130389431 100644 --- a/virtual_rainforest/models/litter/litter_pools.py +++ b/virtual_rainforest/models/litter/litter_pools.py @@ -497,10 +497,13 @@ def calculate_moisture_effect_on_litter_decomp( decomposition [unitless] """ + # TODO - Need to make sure that this function is properly defined for a plausible + # range of matric potentials. + # Calculate how much moisture suppresses microbial activity supression = ( - (np.log10(np.abs(water_potential)) - np.log10(abs(water_potential_opt))) - / (np.log10(abs(water_potential_halt)) - np.log10(abs(water_potential_opt))) + (np.log10(-water_potential) - np.log10(-water_potential_opt)) + / (np.log10(-water_potential_halt) - np.log10(-water_potential_opt)) ) ** moisture_response_curvature return 1 - supression diff --git a/virtual_rainforest/models/soil/carbon.py b/virtual_rainforest/models/soil/carbon.py index 49f657772..f1ef4d2a4 100644 --- a/virtual_rainforest/models/soil/carbon.py +++ b/virtual_rainforest/models/soil/carbon.py @@ -6,10 +6,14 @@ import numpy as np from numpy.typing import NDArray +from scipy.constants import convert_temperature, gas_constant from virtual_rainforest.core.logger import LOGGER from virtual_rainforest.models.soil.constants import SoilConsts +# TODO - Once enzymes are added, temperature dependence of saturation constants should +# be added. + def calculate_soil_carbon_updates( soil_c_pool_lmwc: NDArray[np.float32], @@ -19,6 +23,7 @@ def calculate_soil_carbon_updates( pH: NDArray[np.float32], bulk_density: NDArray[np.float32], soil_moisture: NDArray[np.float32], + soil_water_potential: NDArray[np.float32], soil_temp: NDArray[np.float32], percent_clay: NDArray[np.float32], mineralisation_rate: NDArray[np.float32], @@ -39,6 +44,7 @@ def calculate_soil_carbon_updates( pH: pH values for each soil grid cell bulk_density: bulk density values for each soil grid cell [kg m^-3] soil_moisture: relative water content for each soil grid cell [unitless] + soil_water_potential: Soil water potential for each grid cell [kPa] soil_temp: soil temperature for each soil grid cell [degrees C] percent_clay: Percentage clay for each soil grid cell mineralisation_rate: Amount of litter mineralised into POM pool [kg C m^-3 @@ -51,46 +57,63 @@ def calculate_soil_carbon_updates( A vector containing net changes to each pool. Order [lmwc, maom]. """ - # Find scalar factors that multiple rates - temp_scalar = convert_temperature_to_scalar( - soil_temp, - constants.temp_scalar_coefficient_1, - constants.temp_scalar_coefficient_2, - constants.temp_scalar_coefficient_3, - constants.temp_scalar_coefficient_4, - constants.temp_scalar_reference_temp, - ) + # TODO - At present we have two different factors that capture the impact of soil + # water. water_factor is based on soil water potential and applies to the microbial + # processes, whereas moist_scalar is based on soil water content and applies to the + # chemical processes. In future, all processes will be microbially mediated, meaning + # that moist_scalar can be removed. + + # Find the impact of soil water content on chemical soil processes moist_scalar = convert_moisture_to_scalar( soil_moisture, constants.moisture_scalar_coefficient, constants.moisture_scalar_exponent, ) - moist_temp_scalar = moist_scalar * temp_scalar - + # Find the impact of soil water potential on the biochemical soil processes + water_factor = calculate_water_potential_impact_on_microbes( + water_potential=soil_water_potential, + water_potential_halt=constants.soil_microbe_water_potential_halt, + water_potential_opt=constants.soil_microbe_water_potential_optimum, + moisture_response_curvature=constants.moisture_response_curvature, + ) # Calculate transfers between pools lmwc_to_maom = calculate_mineral_association( - soil_c_pool_lmwc, - soil_c_pool_maom, - pH, - bulk_density, - moist_temp_scalar, - percent_clay, - constants, + soil_c_pool_lmwc=soil_c_pool_lmwc, + soil_c_pool_maom=soil_c_pool_maom, + pH=pH, + bulk_density=bulk_density, + moisture_scalar=moist_scalar, + percent_clay=percent_clay, + constants=constants, ) microbial_uptake = calculate_microbial_carbon_uptake( - soil_c_pool_lmwc, soil_c_pool_microbe, moist_temp_scalar, soil_temp, constants + soil_c_pool_lmwc=soil_c_pool_lmwc, + soil_c_pool_microbe=soil_c_pool_microbe, + water_factor=water_factor, + soil_temp=soil_temp, + constants=constants, ) microbial_respiration = calculate_maintenance_respiration( - soil_c_pool_microbe, moist_temp_scalar, constants.microbial_turnover_rate + soil_c_pool_microbe=soil_c_pool_microbe, + soil_temp=soil_temp, + constants=constants, ) necromass_adsorption = calculate_necromass_adsorption( - soil_c_pool_microbe, moist_temp_scalar, constants.necromass_adsorption_rate + soil_c_pool_microbe=soil_c_pool_microbe, + moisture_scalar=moist_scalar, + necromass_adsorption_rate=constants.necromass_adsorption_rate, ) labile_carbon_leaching = calculate_labile_carbon_leaching( - soil_c_pool_lmwc, moist_temp_scalar, constants.leaching_rate_labile_carbon + soil_c_pool_lmwc=soil_c_pool_lmwc, + moisture_scalar=moist_scalar, + leaching_rate=constants.leaching_rate_labile_carbon, ) pom_decomposition_to_lmwc = calculate_pom_decomposition( - soil_c_pool_pom, soil_c_pool_microbe, moist_temp_scalar, constants + soil_c_pool_pom=soil_c_pool_pom, + soil_c_pool_microbe=soil_c_pool_microbe, + water_factor=water_factor, + soil_temp=soil_temp, + constants=constants, ) # Determine net changes to the pools @@ -117,7 +140,7 @@ def calculate_mineral_association( soil_c_pool_maom: NDArray[np.float32], pH: NDArray[np.float32], bulk_density: NDArray[np.float32], - moist_temp_scalar: NDArray[np.float32], + moisture_scalar: NDArray[np.float32], percent_clay: NDArray[np.float32], constants: SoilConsts, ) -> NDArray[np.float32]: @@ -135,8 +158,8 @@ def calculate_mineral_association( soil_c_pool_maom: Mineral associated organic matter pool [kg C m^-3] pH: pH values for each soil grid cell bulk_density: bulk density values for each soil grid cell [kg m^-3] - moist_temp_scalar: A scalar capturing the combined impact of soil moisture and - temperature on process rates + moisture_scalar: A scalar capturing the impact of soil moisture and on process + rates [unitless] percent_clay: Percentage clay for each soil grid cell constants: Set of constants for the soil model. @@ -153,9 +176,7 @@ def calculate_mineral_association( ) equib_maom = calculate_equilibrium_maom(pH, Q_max, soil_c_pool_lmwc, constants) - return ( - moist_temp_scalar * soil_c_pool_lmwc * (equib_maom - soil_c_pool_maom) / Q_max - ) + return moisture_scalar * soil_c_pool_lmwc * (equib_maom - soil_c_pool_maom) / Q_max def calculate_max_sorption_capacity( @@ -251,48 +272,80 @@ def calculate_binding_coefficient( return 10.0 ** (binding_with_ph_slope * pH + binding_with_ph_intercept) -def convert_temperature_to_scalar( - soil_temp: NDArray[np.float32], - temp_scalar_coefficient_1: float, - temp_scalar_coefficient_2: float, - temp_scalar_coefficient_3: float, - temp_scalar_coefficient_4: float, - temp_scalar_reference_temp: float, +def calculate_temperature_effect_on_microbes( + soil_temperature: NDArray[np.float32], + activation_energy: float, + reference_temperature: float, ) -> NDArray[np.float32]: - """Convert soil temperature into a factor to multiply rates by. + """Calculate the effect that temperature has on microbial metabolic rates. - This form is used in :cite:t:`abramoff_millennial_2018` to minimise differences with - the CENTURY model. We very likely want to define our own functional form here. I'm - also a bit unsure how this form was even obtained, so further work here is very - needed. + This uses a standard Arrhenius equation to calculate the impact of temperature. + + This function takes temperatures in Celsius as inputs and converts them into Kelvin + for use in the Arrhenius equation. TODO - review this after we have decided how to + handle these conversions in general. Args: - soil_temp: soil temperature for each soil grid cell [degrees C] - temp_scalar_coefficient_1: Unclear exactly what this parameter is [degrees C] - temp_scalar_coefficient_2: Unclear exactly what this parameter is [unclear] - temp_scalar_coefficient_3: Unclear exactly what this parameter is [unclear] - temp_scalar_coefficient_4: Unclear exactly what this parameter is [unclear] - temp_scalar_reference_temp: Reference temperature for temperature scalar - [degrees C] + soil_temperature: The temperature of the soil [C] + activation_energy: Energy of activation [J mol^-1] + soil_temperature: The reference temperature of the Arrhenius equation [C] Returns: - A scalar that captures the impact of soil temperature on process rates + A multiplicative factor capturing the effect of temperature on microbial rates """ - # This expression is drawn from Abramoff et al. (2018) - numerator = temp_scalar_coefficient_2 + ( - temp_scalar_coefficient_3 / np.pi - ) * np.arctan(np.pi * (soil_temp - temp_scalar_coefficient_1)) - - denominator = temp_scalar_coefficient_2 + ( - temp_scalar_coefficient_3 / np.pi - ) * np.arctan( - np.pi - * temp_scalar_coefficient_4 - * (temp_scalar_reference_temp - temp_scalar_coefficient_1) + # Convert the temperatures to Kelvin + soil_temp_in_kelvin = convert_temperature( + soil_temperature, old_scale="Celsius", new_scale="Kelvin" + ) + ref_temp_in_kelvin = convert_temperature( + reference_temperature, old_scale="Celsius", new_scale="Kelvin" ) - return np.divide(numerator, denominator) + return np.exp( + (-activation_energy / gas_constant) + * ((1 / (soil_temp_in_kelvin)) - (1 / (ref_temp_in_kelvin))) + ) + + +def calculate_water_potential_impact_on_microbes( + water_potential: NDArray[np.float32], + water_potential_halt: float, + water_potential_opt: float, + moisture_response_curvature: float, +) -> NDArray[np.float32]: + """Calculate the effect that soil water potential has on microbial rates. + + This function only returns valid output for soil water potentials that are less than + the optimal water potential. + + Args: + water_potential: Soil water potential [kPa] + water_potential_halt: Water potential at which all microbial activity stops + [kPa] + water_potential_opt: Optimal water potential for microbial activity [kPa] + moisture_response_curvature: Parameter controlling the curvature of the moisture + response function [unitless] + + Returns: + A multiplicative factor capturing the impact of moisture on soil microbe rates + decomposition [unitless] + """ + + # If the water potential is greater than the optimal then the function produces NaNs + # so the simulation should be interrupted + if any(water_potential > water_potential_opt): + err = ValueError("Water potential greater than minimum value") + LOGGER.critical(err) + raise err + + # Calculate how much moisture suppresses microbial activity + supression = ( + (np.log10(-water_potential) - np.log10(-water_potential_opt)) + / (np.log10(-water_potential_halt) - np.log10(-water_potential_opt)) + ) ** moisture_response_curvature + + return 1 - supression def convert_moisture_to_scalar( @@ -333,42 +386,48 @@ def convert_moisture_to_scalar( def calculate_maintenance_respiration( soil_c_pool_microbe: NDArray[np.float32], - moist_temp_scalar: NDArray[np.float32], - microbial_turnover_rate: float, + soil_temp: NDArray[np.float32], + constants: SoilConsts, ) -> NDArray[np.float32]: """Calculate the maintenance respiration of the microbial pool. Args: soil_c_pool_microbe: Microbial biomass (carbon) pool [kg C m^-3] - moist_scalar: A scalar capturing the combined impact of soil moisture and - temperature on process rates - microbial_turnover_rate: Rate of microbial biomass turnover [day^-1] + soil_temp: soil temperature for each soil grid cell [degrees C] + constants: Set of constants for the soil model. Returns: - Total respiration for all microbial biomass + Total maintenance respiration for all microbial biomass """ - return microbial_turnover_rate * moist_temp_scalar * soil_c_pool_microbe + temp_factor = calculate_temperature_effect_on_microbes( + soil_temperature=soil_temp, + activation_energy=constants.activation_energy_microbial_turnover, + reference_temperature=constants.arrhenius_reference_temp, + ) + + return constants.microbial_turnover_rate * temp_factor * soil_c_pool_microbe def calculate_necromass_adsorption( soil_c_pool_microbe: NDArray[np.float32], - moist_temp_scalar: NDArray[np.float32], + moisture_scalar: NDArray[np.float32], necromass_adsorption_rate: float, ) -> NDArray[np.float32]: """Calculate adsorption of microbial necromass to soil minerals. Args: soil_c_pool_microbe: Microbial biomass (carbon) pool [kg C m^-3] - moist_temp_scalar: A scalar capturing the combined impact of soil moisture and - temperature on process rates + moisture_scalar: A scalar capturing the impact of soil moisture on process rates + [unitless] necromass_adsorption_rate: Rate at which necromass is adsorbed by soil minerals + [day^-1] Returns: Adsorption of microbial biomass to mineral associated organic matter (MAOM) """ - return necromass_adsorption_rate * moist_temp_scalar * soil_c_pool_microbe + return necromass_adsorption_rate * moisture_scalar * soil_c_pool_microbe def calculate_carbon_use_efficiency( @@ -379,6 +438,8 @@ def calculate_carbon_use_efficiency( ) -> NDArray[np.float32]: """Calculate the (temperature dependant) carbon use efficiency. + TODO - This should be adapted to use an Arrhenius function at some point. + Args: soil_temp: soil temperature for each soil grid cell [degrees C] reference_cue: Carbon use efficiency at reference temp [unitless] @@ -468,7 +529,7 @@ def calculate_pom_decomposition_saturation( def calculate_microbial_carbon_uptake( soil_c_pool_lmwc: NDArray[np.float32], soil_c_pool_microbe: NDArray[np.float32], - moist_temp_scalar: NDArray[np.float32], + water_factor: NDArray[np.float32], soil_temp: NDArray[np.float32], constants: SoilConsts, ) -> NDArray[np.float32]: @@ -477,8 +538,8 @@ def calculate_microbial_carbon_uptake( Args: soil_c_pool_lmwc: Low molecular weight carbon pool [kg C m^-3] soil_c_pool_microbe: Microbial biomass (carbon) pool [kg C m^-3] - moist_temp_scalar: A scalar capturing the combined impact of soil moisture and - temperature on process rates + water_factor: A factor capturing the impact of soil water potential on microbial + rates [unitless] soil_temp: soil temperature for each soil grid cell [degrees C] constants: Set of constants for the soil model. @@ -496,6 +557,11 @@ def calculate_microbial_carbon_uptake( microbial_saturation = calculate_microbial_saturation( soil_c_pool_microbe, constants.half_sat_microbial_activity ) + temp_factor = calculate_temperature_effect_on_microbes( + soil_temperature=soil_temp, + activation_energy=constants.activation_energy_microbial_uptake, + reference_temperature=constants.arrhenius_reference_temp, + ) # TODO - the quantities calculated above can be used to calculate the carbon # respired instead of being uptaken. This isn't currently of interest, but will be @@ -503,7 +569,8 @@ def calculate_microbial_carbon_uptake( return ( constants.max_uptake_rate_labile_C - * moist_temp_scalar + * water_factor + * temp_factor * soil_c_pool_lmwc * microbial_saturation * carbon_use_efficency @@ -512,7 +579,7 @@ def calculate_microbial_carbon_uptake( def calculate_labile_carbon_leaching( soil_c_pool_lmwc: NDArray[np.float32], - moist_temp_scalar: NDArray[np.float32], + moisture_scalar: NDArray[np.float32], leaching_rate: float, ) -> NDArray[np.float32]: """Calculate rate at which labile carbon is leached. @@ -522,21 +589,22 @@ def calculate_labile_carbon_leaching( Args: soil_c_pool_lmwc: Low molecular weight carbon pool [kg C m^-3] - moist_temp_scalar: A scalar capturing the combined impact of soil moisture and - temperature on process rates + moisture_scalar: A scalar capturing the impact of soil moisture on process rates + [unitless] leaching_rate: The rate at which labile carbon leaches from the soil [day^-1] Returns: The amount of labile carbon leached """ - return leaching_rate * moist_temp_scalar * soil_c_pool_lmwc + return leaching_rate * moisture_scalar * soil_c_pool_lmwc def calculate_pom_decomposition( soil_c_pool_pom: NDArray[np.float32], soil_c_pool_microbe: NDArray[np.float32], - moist_temp_scalar: NDArray[np.float32], + water_factor: NDArray[np.float32], + soil_temp: NDArray[np.float32], constants: SoilConsts, ) -> NDArray[np.float32]: """Calculate decomposition of particulate organic matter into labile carbon (LMWC). @@ -547,8 +615,9 @@ def calculate_pom_decomposition( Args: soil_c_pool_pom: Particulate organic matter pool [kg C m^-3] soil_c_pool_microbe: Microbial biomass (carbon) pool [kg C m^-3] - moist_temp_scalar: A scalar capturing the combined impact of soil moisture and - temperature on process rates + water_factor: A factor capturing the impact of soil water potential on microbial + rates [unitless] + soil_temp: soil temperature for each soil grid cell [degrees C] constants: Set of constants for the soil model. Returns: @@ -564,9 +633,17 @@ def calculate_pom_decomposition( soil_c_pool_pom, constants.half_sat_pom_decomposition ) + # Calculate the impact of temperature on the rate + temp_factor = calculate_temperature_effect_on_microbes( + soil_temperature=soil_temp, + activation_energy=constants.activation_energy_pom_decomp, + reference_temperature=constants.arrhenius_reference_temp, + ) + return ( constants.max_decomp_rate_pom * saturation_with_pom * saturation_with_biomass - * moist_temp_scalar + * water_factor + * temp_factor ) diff --git a/virtual_rainforest/models/soil/constants.py b/virtual_rainforest/models/soil/constants.py index b9640e22c..1ec32dca0 100644 --- a/virtual_rainforest/models/soil/constants.py +++ b/virtual_rainforest/models/soil/constants.py @@ -7,9 +7,6 @@ # TODO - Need to figure out a sensible area to volume conversion -# TODO - Need to either work out what the temp and moisture scalars are, or find an -# alternative way of capturing temperature and moisture effects - @dataclass(frozen=True) class SoilConsts: @@ -56,86 +53,125 @@ class SoilConsts: Units of [(RWC)^-1] """ - temp_scalar_coefficient_1: float = 15.4 - """Used in :cite:t:`abramoff_millennial_2018`, can't find original source. + reference_cue: float = 0.6 + """Carbon use efficiency of community at the reference temperature [no units]. - Unclear exactly what this parameter is [degrees C] + Default value taken from :cite:t:`abramoff_millennial_2018`. """ - temp_scalar_coefficient_2: float = 11.75 - """Used in :cite:t:`abramoff_millennial_2018`, can't find original source. + cue_reference_temp: float = 15.0 + """Reference temperature for carbon use efficiency [degrees C]. - Unclear exactly what this parameter is [units unclear] + Default value taken from :cite:t:`abramoff_millennial_2018`. """ - temp_scalar_coefficient_3: float = 29.7 - """Used in :cite:t:`abramoff_millennial_2018`, can't find original source. + cue_with_temperature: float = 0.012 + """Change in carbon use efficiency with increasing temperature [degree C^-1]. - Unclear exactly what this parameter is [units unclear] + Default value taken from :cite:t:`abramoff_millennial_2018`. """ - temp_scalar_coefficient_4: float = 0.031 - """Used in :cite:t:`abramoff_millennial_2018`, can't find original source. + necromass_adsorption_rate: float = 0.025 + """Rate at which necromass is adsorbed by soil minerals [day^-1]. - Unclear exactly what this parameter is [units unclear] + Taken from :cite:t:`abramoff_millennial_2018`, where it was obtained by calibration. """ - temp_scalar_reference_temp: float = 30.0 - """Used in :cite:t:`abramoff_millennial_2018`, can't find original source. + half_sat_microbial_activity: float = 0.0072 + """Half saturation constant for microbial activity (with increasing biomass). - Reference temperature [degrees C] + Units of [kg C m^-2]. """ - reference_cue: float = 0.6 - """Carbon use efficiency of community at the reference temperature [no units]. + half_sat_microbial_pom_mineralisation: float = 0.012 + """Half saturation constant for microbial POM mineralisation [kg C m^-2].""" - Default value taken from :cite:t:`abramoff_millennial_2018`. + leaching_rate_labile_carbon: float = 1.5e-3 + """Leaching rate for labile carbon (lmwc) [day^-1].""" + + half_sat_pom_decomposition: float = 0.150 + """Half saturation constant for POM decomposition to LMWC [kg C m^-2].""" + + soil_microbe_water_potential_optimum: float = -3.0 + """The water potential at which soil microbial rates are maximised [kPa]. + + Value is taken from :cite:t`moyano_responses_2013`. """ - cue_reference_temp: float = 15.0 - """Reference temperature for carbon use efficiency [degrees C]. + soil_microbe_water_potential_halt: float = -15800.0 + """The water potential at which soil microbial activity stops entirely [kPa]. - Default value taken from :cite:t:`abramoff_millennial_2018`. + Value is taken from :cite:t`moyano_responses_2013`. """ - cue_with_temperature: float = 0.012 - """Change in carbon use efficiency with increasing temperature [degree C^-1]. + moisture_response_curvature: float = 1.47 + """Curvature of the soil microbial moisture response function [unitless]. - Default value taken from :cite:t:`abramoff_millennial_2018`. + Value is taken from :cite:t`moyano_responses_2013`. """ - microbial_turnover_rate: float = 0.036 - """Microbial turnover rate [day^-1], this isn't a constant but often treated as one. + arrhenius_reference_temp: float = 12.0 + """Reference temperature for the Arrhenius equation [C]. + + This is the reference temperature used in :cite:t:`wang_development_2013`, which is + the source of the activation energies and corresponding rates. """ - max_uptake_rate_labile_C: float = 0.35 - """Maximum (theoretical) rate at which microbes can take up labile carbon [day^-1]. + # TODO - Split this and the following into 2 constants once fungi are introduced + max_uptake_rate_labile_C: float = 0.04 + """Maximum rate at the reference temperature of labile carbon uptake [day^-1]. + + The reference temperature is given + by :attr:`arrhenius_reference_temp`, and the corresponding activation energy is + given by :attr:`activation_energy_microbial_uptake`. + + TODO - Source of this constant is not completely clear, investigate this further + once fungi are added. """ - necromass_adsorption_rate: float = 0.025 - """Rate at which necromass is adsorbed by soil minerals [day^-1]. + activation_energy_microbial_uptake = 47000 + """Activation energy for microbial uptake of low molecular weight carbon [J K^-1]. - Taken from :cite:t:`abramoff_millennial_2018`, where it was obtained by calibration. + Value taken from :cite:t:`wang_development_2013`. The maximum labile carbon uptake + rate that this activation energy corresponds to is given by + :attr:`max_uptake_rate_labile_C`. """ - half_sat_microbial_activity: float = 0.0072 - """Half saturation constant for microbial activity (with increasing biomass). + # TODO - Add another set of constants once we start tracking lignin + max_decomp_rate_pom: float = 0.2 + """Maximum rate for particulate organic matter break down (at reference temp). - Units of [kg C m^-2]. + Units of [day^-1]. The reference temperature is given by + :attr:`arrhenius_reference_temp`, and the corresponding activation energy is given + by :attr:`activation_energy_pom_decomp`. + + TODO - Once enzymes are included this should take the value of 200.0 [day^-1]. + + TODO - Source of this constant is not completely clear, investigate this further + once lignin chemistry is added. """ - half_sat_microbial_pom_mineralisation: float = 0.012 - """Half saturation constant for microbial POM mineralisation [kg C m^-2].""" + activation_energy_pom_decomp = 37000 + """Activation energy for decomposition of particulate organic matter [J K^-1]. - max_decomp_rate_pom: float = 0.01 - """Maximum (theoretical) rate for particulate organic matter break down. + Taken from :cite:t:`wang_development_2013`. + """ - Units of [kg C m^-2 day^-1]. Taken from :cite:t:`abramoff_millennial_2018`, where it - was obtained by calibration. + # TODO - Split this and the following into 2 constants once fungi are introduced + microbial_turnover_rate: float = 0.005 + """Microbial turnover rate at reference temperature [day^-1]. + + The reference temperature is given by :attr:`arrhenius_reference_temp`, and the + corresponding activation energy is given by + :attr:`activation_energy_microbial_turnover`. + + TODO - Source of this constant is not completely clear, investigate this further + once fungi are added. """ - leaching_rate_labile_carbon: float = 1.5e-3 - """Leaching rate for labile carbon (lmwc) [day^-1].""" + activation_energy_microbial_turnover = 20000 + """Activation energy for microbial maintenance turnover rate [J K^-1]. - half_sat_pom_decomposition: float = 0.150 - """Half saturation constant for POM decomposition to LMWC [kg C m^-2].""" + Value taken from :cite:t:`wang_development_2013`. The microbial turnover rate that + this activation energy corresponds to is given by :attr:`microbial_turnover_rate`. + """ diff --git a/virtual_rainforest/models/soil/soil_model.py b/virtual_rainforest/models/soil/soil_model.py index 400136e4b..5855635eb 100644 --- a/virtual_rainforest/models/soil/soil_model.py +++ b/virtual_rainforest/models/soil/soil_model.py @@ -288,6 +288,7 @@ def construct_full_soil_model( pH=data["pH"].to_numpy(), bulk_density=data["bulk_density"].to_numpy(), soil_moisture=data["soil_moisture"][top_soil_layer_index].to_numpy(), + soil_water_potential=data["matric_potential"][top_soil_layer_index].to_numpy(), soil_temp=data["soil_temperature"][top_soil_layer_index].to_numpy(), percent_clay=data["percent_clay"].to_numpy(), mineralisation_rate=data["litter_C_mineralisation_rate"].to_numpy(),