Added step to break MAOM down into LMWC

tests/models/soil/test_carbon.py

@@ -15,9 +15,9 @@ 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": [-3.70969075e-2, -2.2791314e-2, -0.15100693, 7.2342438e-4],
-        "soil_c_pool_maom": [1.6992235e-4, 4.45995156e-3, 6.252756e-3, 2.0712399e-5],
-        "soil_c_pool_microbe": [-0.03828773, -0.01245439, -0.06446385, -0.00711458],
+        "soil_c_pool_lmwc": [-0.00371115, 0.00278502, -0.01849181, 0.00089995],
+        "soil_c_pool_maom": [-1.28996257e-3, 2.35822401e-3, 1.5570399e-3, 1.2082886e-5],
+        "soil_c_pool_microbe": [-0.04978105, -0.02020101, -0.10280967, -0.00719517],
         "soil_c_pool_pom": [0.04809165, 0.01023544, 0.07853728, 0.01167564],
         "soil_enzyme_pom": [1.18e-8, 1.67e-8, 1.8e-9, -1.12e-8],
         "soil_enzyme_maom": [-0.00031009, -5.09593e-5, 0.0005990658, -3.72112e-5],
@@ -165,8 +165,8 @@ def test_calculate_microbial_carbon_uptake(
     """Check microbial carbon uptake calculates correctly."""
     from virtual_rainforest.models.soil.carbon import calculate_microbial_carbon_uptake
 
-    expected_uptake = [0.04484178, 0.03190812, 0.18552910, 0.00022563]
-    expected_assimilation = [0.01614304, 0.01052968, 0.05565873, 1.08303e-4]
+    expected_uptake = [1.29159055e-2, 8.43352433e-3, 5.77096991e-2, 5.77363558e-5]
+    expected_assimilation = [4.64972597e-3, 2.78306303e-3, 1.73129097e-2, 2.77134508e-5]
 
     actual_uptake, actual_assimilation = calculate_microbial_carbon_uptake(
         soil_c_pool_lmwc=dummy_carbon_data["soil_c_pool_lmwc"],
@@ -200,7 +200,11 @@ def test_calculate_enzyme_mediated_decomposition(
         pH_factor=environmental_factors["pH"],
         clay_factor_saturation=environmental_factors["clay_saturation"],
         soil_temp=dummy_carbon_data["soil_temperature"][top_soil_layer_index],
-        constants=SoilConsts,
+        reference_temp=SoilConsts.arrhenius_reference_temp,
+        max_decomp_rate=SoilConsts.max_decomp_rate_pom,
+        activation_energy_rate=SoilConsts.activation_energy_pom_decomp_rate,
+        half_saturation=SoilConsts.half_sat_pom_decomposition,
+        activation_energy_sat=SoilConsts.activation_energy_pom_decomp_saturation,
     )
 
     assert np.allclose(actual_decomp, expected_decomp)

tests/models/soil/test_soil_model.py

@@ -294,22 +294,22 @@ def test_update(mocker, soil_model_fixture, dummy_carbon_data):
             Dataset(
                 data_vars=dict(
                     lmwc=DataArray(
-                        [0.03460932, 0.01198232, 0.04607224, 0.00520078], dims="cell_id"
+                        [0.04823845, 0.02119714, 0.0895863, 0.00528887], dims="cell_id"
                     ),
                     maom=DataArray(
-                        [2.50009515, 1.70223563, 4.50321333, 0.50001044], dims="cell_id"
+                        [2.49936689, 1.70118553, 4.50085129, 0.50000614], dims="cell_id"
                     ),
                     microbe=DataArray(
-                        [5.779758, 2.2926573, 11.26080985, 0.99645009], dims="cell_id"
+                        [5.77512315, 2.2899636, 11.24827514, 0.99640928], dims="cell_id"
                     ),
                     pom=DataArray(
-                        [0.12398572, 1.00509262, 0.73902017, 0.35583213], dims="cell_id"
+                        [0.12397575, 1.00508662, 0.7389913, 0.35583206], dims="cell_id"
                     ),
                     enzyme_pom=DataArray(
-                        [0.02267854, 0.00957584, 0.05005016, 0.00300993], dims="cell_id"
+                        [0.02267842, 0.00957576, 0.05004963, 0.00300993], dims="cell_id"
                     ),
                     enzyme_maom=DataArray(
-                        [0.03544542, 0.0116745, 0.02538689, 0.00454144], dims="cell_id"
+                        [0.0354453, 0.01167442, 0.02538637, 0.00454144], dims="cell_id"
                     ),
                 )
             ),
@@ -424,18 +424,18 @@ 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 = [
-        -3.70969075e-2,
-        -2.2791314e-2,
-        -0.15100693,
-        7.2342438e-4,
-        1.6992235e-4,
-        4.45995156e-3,
-        6.252756e-3,
-        2.0712399e-5,
-        -3.82877348e-2,
-        -1.24543917e-2,
-        -6.44638512e-2,
-        -7.11457978e-3,
+        -0.00371115,
+        0.00278502,
+        -0.01849181,
+        0.00089995,
+        -1.28996257e-3,
+        2.35822401e-3,
+        1.5570399e-3,
+        1.2082886e-5,
+        -0.04978105,
+        -0.02020101,
+        -0.10280967,
+        -0.00719517,
         4.80916464e-2,
         1.02354410e-2,
         7.85372753e-2,

virtual_rainforest/models/soil/carbon.py

@@ -38,7 +38,9 @@ class MicrobialBiomassLoss:
     """Rate at which biomass is lost to the POM pool [kg C m^-3 day^-1]."""
 
 
-# TODO - This function should probably be shortened
+# TODO - This function should probably be shortened. First step is group the calculation
+# of environmental factors. Otherwise there's probably stuff that I'm missing, which I
+# should tackle later
 def calculate_soil_carbon_updates(
     soil_c_pool_lmwc: NDArray[np.float32],
     soil_c_pool_maom: NDArray[np.float32],
@@ -150,26 +152,45 @@ def calculate_soil_carbon_updates(
         pH_factor=pH_factor,
         clay_factor_saturation=clay_factor_saturation,
         soil_temp=soil_temp,
-        constants=model_constants,
+        reference_temp=model_constants.arrhenius_reference_temp,
+        max_decomp_rate=model_constants.max_decomp_rate_pom,
+        activation_energy_rate=model_constants.activation_energy_pom_decomp_rate,
+        half_saturation=model_constants.half_sat_pom_decomposition,
+        activation_energy_sat=model_constants.activation_energy_pom_decomp_saturation,
     )
-
+    # Calculate how pom decomposition is split between lmwc and maom pools
     pom_decomposition_to_lmwc = (
         pom_decomposition_rate * model_constants.pom_decomposition_fraction_lmwc
     )
     pom_decomposition_to_maom = pom_decomposition_rate * (
         1 - model_constants.pom_decomposition_fraction_lmwc
     )
-    # TODO - Use enzyme decomp equation to calculate lmwc to maom transfer
+    maom_decomposition_to_lmwc = calculate_enzyme_mediated_decomposition(
+        soil_c_pool=soil_c_pool_maom,
+        soil_enzyme=soil_enzyme_maom,
+        water_factor=water_factor,
+        pH_factor=pH_factor,
+        clay_factor_saturation=clay_factor_saturation,
+        soil_temp=soil_temp,
+        reference_temp=model_constants.arrhenius_reference_temp,
+        max_decomp_rate=model_constants.max_decomp_rate_maom,
+        activation_energy_rate=model_constants.activation_energy_maom_decomp_rate,
+        half_saturation=model_constants.half_sat_maom_decomposition,
+        activation_energy_sat=model_constants.activation_energy_maom_decomp_saturation,
+    )
 
     # Determine net changes to the pools
     delta_pools_ordered["soil_c_pool_lmwc"] = (
         pom_decomposition_to_lmwc
         + biomass_losses.necromass_decay_to_lmwc
         + pom_enzyme_turnover
+        + maom_decomposition_to_lmwc
         - microbial_uptake
         - labile_carbon_leaching
     )
-    delta_pools_ordered["soil_c_pool_maom"] = pom_decomposition_to_maom
+    delta_pools_ordered["soil_c_pool_maom"] = (
+        pom_decomposition_to_maom - maom_decomposition_to_lmwc
+    )
     delta_pools_ordered["soil_c_pool_microbe"] = (
         microbial_assimilation - biomass_losses.maintenance_synthesis
     )
@@ -394,7 +415,11 @@ def calculate_enzyme_mediated_decomposition(
     pH_factor: NDArray[np.float32],
     clay_factor_saturation: NDArray[np.float32],
     soil_temp: NDArray[np.float32],
-    constants: SoilConsts,
+    reference_temp: float,
+    max_decomp_rate: float,
+    activation_energy_rate: float,
+    half_saturation: float,
+    activation_energy_sat: float,
 ) -> NDArray[np.float32]:
     """Calculate rate of a enzyme mediated decomposition process.
 
@@ -413,7 +438,15 @@ def calculate_enzyme_mediated_decomposition(
         clay_factor_saturation: A factor capturing the impact of soil clay fraction on
             enzyme saturation constants [unitless]
         soil_temp: soil temperature for each soil grid cell [degrees C]
-        constants: Set of constants for the soil model.
+        reference_temp: The reference temperature that enzyme rates were determined
+            relative to [degrees C]
+        max_decomp_rate: The maximum rate of substrate decomposition (at the reference
+            temperature) [day^-1]
+        activation_energy_rate: Activation energy for maximum decomposition rate
+            [J K^-1]
+        half_saturation: Half saturation constant for decomposition (at the reference
+            temperature) [kg C m^-3]
+        activation_energy_sat: Activation energy for decomposition saturation [J K^-1]
 
     Returns:
         The rate of decomposition of the organic matter pool in question [kg C m^-3
@@ -423,23 +456,19 @@ def calculate_enzyme_mediated_decomposition(
     # Calculate the factors which impact the rate and saturation constants
     temp_factor_rate = calculate_temperature_effect_on_microbes(
         soil_temperature=soil_temp,
-        activation_energy=constants.activation_energy_pom_decomp_rate,
-        reference_temperature=constants.arrhenius_reference_temp,
+        activation_energy=activation_energy_rate,
+        reference_temperature=reference_temp,
     )
     temp_factor_saturation = calculate_temperature_effect_on_microbes(
         soil_temperature=soil_temp,
-        activation_energy=constants.activation_energy_pom_decomp_saturation,
-        reference_temperature=constants.arrhenius_reference_temp,
+        activation_energy=activation_energy_sat,
+        reference_temperature=reference_temp,
     )
 
     # Calculate the adjusted rate and saturation constants
-    rate_constant = (
-        constants.max_decomp_rate_pom * temp_factor_rate * water_factor * pH_factor
-    )
+    rate_constant = max_decomp_rate * temp_factor_rate * water_factor * pH_factor
     saturation_constant = (
-        constants.half_sat_pom_decomposition
-        * temp_factor_saturation
-        * clay_factor_saturation
+        half_saturation * temp_factor_saturation * clay_factor_saturation
     )
 
     return (

virtual_rainforest/models/soil/constants.py

@@ -7,7 +7,7 @@
 
 from virtual_rainforest.core.constants_class import ConstantsDataclass
 
-# TODO - Need to figure out a sensible area to volume conversion
+# TODO - Once lignin is tracked a large number of constants will have to be duplicated
 
 
 @dataclass(frozen=True)
@@ -77,14 +77,14 @@ class SoilConsts(ConstantsDataclass):
     :attr:`max_uptake_rate_labile_C`.
     """
 
-    half_sat_labile_C_uptake: float = 0.091
+    half_sat_labile_C_uptake: float = 0.364
     """Half saturation constant for microbial uptake of labile carbon (LMWC).
 
-    [kg C m^-2]. This was calculated from the value provided in
-    :cite:t:`wang_development_2013` assuming an average bulk density of 1400 [kg m^-3],
-    and a topsoil depth of 0.25 [m]. The reference temperature is given by
-    :attr:`arrhenius_reference_temp`, and the corresponding activation energy is given
-    by :attr:`activation_energy_labile_C_saturation`.
+    [kg C m^-3]. This was calculated from the value provided in
+    :cite:t:`wang_development_2013` assuming an average bulk density of 1400 [kg m^-3].
+    The reference temperature is given by :attr:`arrhenius_reference_temp`, and the
+    corresponding activation energy is given by
+    :attr:`activation_energy_labile_C_saturation`.
     """
 
     activation_energy_labile_C_saturation: float = 30000
@@ -93,15 +93,13 @@ class SoilConsts(ConstantsDataclass):
     Taken from :cite:t:`wang_development_2013`.
     """
 
-    # TODO - Add another set of constants once we start tracking lignin
     half_sat_pom_decomposition: float = 70.0
-    """Half saturation constant for POM decomposition to LMWC [kg C m^-2].
+    """Half saturation constant for POM decomposition to LMWC [kg C m^-3].
 
     This was calculated from the value provided in :cite:t:`wang_development_2013`
-    assuming an average bulk density of 1400 [kg m^-3], and a topsoil depth of 0.25 [m].
-    The reference temperature is given by :attr:`arrhenius_reference_temp`, and the
-    corresponding activation energy is given by
-    :attr:`activation_energy_pom_decomp_saturation`.
+    assuming an average bulk density of 1400 [kg m^-3]. The reference temperature is
+    given by :attr:`arrhenius_reference_temp`, and the corresponding activation energy
+    is given by :attr:`activation_energy_pom_decomp_saturation`.
     """
 
     activation_energy_pom_decomp_saturation: float = 30000
@@ -110,7 +108,6 @@ class SoilConsts(ConstantsDataclass):
     Taken from :cite:t:`wang_development_2013`.
     """
 
-    # TODO - Add another set of constants once we start tracking lignin
     max_decomp_rate_pom: float = 60.0
     """Maximum rate for particulate organic matter break down (at reference temp).
 
@@ -128,6 +125,36 @@ class SoilConsts(ConstantsDataclass):
     Taken from :cite:t:`wang_development_2013`.
     """
 
+    half_sat_maom_decomposition: float = 350.0
+    """Half saturation constant for MAOM decomposition to LMWC [kg C m^-3].
+
+    This was calculated from the value provided in :cite:t:`wang_development_2013`
+    assuming an average bulk density of 1400 [kg m^-3]. The reference temperature is
+    given by :attr:`arrhenius_reference_temp`, and the corresponding activation energy
+    is given by :attr:`activation_energy_maom_decomp_saturation`.
+    """
+
+    activation_energy_maom_decomp_saturation: float = 30000
+    """Activation energy for MAOM decomposition saturation constant [J K^-1].
+
+    Taken from :cite:t:`wang_development_2013`.
+    """
+
+    max_decomp_rate_maom: float = 24.0
+    """Maximum rate for mineral associated organic matter decomposition enzyme.
+
+    Units of [day^-1]. The rate is for a reference temperature which is given by
+    :attr:`arrhenius_reference_temp`, and the corresponding activation energy is given
+    by :attr:`activation_energy_maom_decomp_rate`. The value is taken from
+    :cite:t:`wang_development_2013`.
+    """
+
+    activation_energy_maom_decomp_rate: float = 47000
+    """Activation energy for decomposition of mineral associated organic matter.
+
+    Units of [J K^-1]. Taken from :cite:t:`wang_development_2013`.
+    """
+
     # 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].