Merge pull request #724 from ImperialCollegeLondon/708-add-functional…

…ity-for-nitrogen-fixation

Add nitrogen fixation to the soil model

docs/source/refs.bib

@@ -920,3 +920,37 @@ @article{xu-ri_terrestrial_2008
 	year = {2008},
 	pages = {1745--1764},
 }
+
+@article{brzostek_modeling_2014,
+	title = {Modeling the carbon cost of plant nitrogen acquisition: {Mycorrhizal} trade-offs and multipath resistance uptake improve predictions of retranslocation: {Carbon} cost of mycorrhizae},
+	volume = {119},
+	issn = {21698953},
+	shorttitle = {Modeling the carbon cost of plant nitrogen acquisition},
+	url = {http://doi.wiley.com/10.1002/2014JG002660},
+	doi = {10.1002/2014JG002660},
+	language = {en},
+	number = {8},
+	urldate = {2021-10-27},
+	journal = {Journal of Geophysical Research: Biogeosciences},
+	author = {Brzostek, Edward R. and Fisher, Joshua B. and Phillips, Richard P.},
+	month = aug,
+	year = {2014},
+	pages = {1684--1697},
+}
+
+@article{lin_modelling_2000,
+	title = {Modelling a global biogeochemical nitrogen cycle in terrestrial ecosystems},
+	volume = {135},
+	copyright = {https://www.elsevier.com/tdm/userlicense/1.0/},
+	issn = {03043800},
+	url = {https://linkinghub.elsevier.com/retrieve/pii/S0304380000003720},
+	doi = {10.1016/S0304-3800(00)00372-0},
+	language = {en},
+	number = {1},
+	urldate = {2025-02-04},
+	journal = {Ecological Modelling},
+	author = {Lin, Bin-Le and Sakoda, Akiyoshi and Shibasaki, Ryosuke and Goto, Naohiro and Suzuki, Motoyuki},
+	month = nov,
+	year = {2000},
+	pages = {89--110},
+}

tests/conftest.py

@@ -281,6 +281,7 @@ def dummy_carbon_data(fixture_core_components):
         "litter_N_mineralisation_rate": [3.5351e-5, 7.0702e-5, 0.000183, 1.63333e-5],
         "litter_P_mineralisation_rate": [7.32e-6, 1.41404e-6, 2.82808e-6, 6.53332e-7],
         "vertical_flow": [0.1, 0.5, 2.5, 1.59],
+        "nitrogen_fixation_carbon_supply": [0.01, 0.25, 0.0075, 0.0047],
     }
 
     for var_name, var_values in data_values.items():

tests/models/soil/test_env_factors.py

@@ -315,6 +315,62 @@ def test_denitrification_temperature_factor_bad_temp(
     assert np.allclose(expected_factor, actual_factor)
 
 
+def test_calculate_symbiotic_nitrogen_fixation_carbon_cost(
+    dummy_carbon_data, fixture_core_components
+):
+    """Test calculation of symbiotic nitrogen fixation cost."""
+    from virtual_ecosystem.models.soil.constants import SoilConsts
+    from virtual_ecosystem.models.soil.env_factors import (
+        calculate_symbiotic_nitrogen_fixation_carbon_cost,
+    )
+
+    expected_cost = [45.8029373, 49.464657, 52.689498, 36.073368]
+
+    actual_cost = calculate_symbiotic_nitrogen_fixation_carbon_cost(
+        soil_temp=dummy_carbon_data["soil_temperature"][
+            fixture_core_components.layer_structure.index_topsoil_scalar
+        ],
+        cost_at_zero_celsius=SoilConsts.nitrogen_fixation_cost_zero_celcius,
+        infinite_temp_cost_offset=SoilConsts.nitrogen_fixation_cost_infinite_temp_offset,
+        thermal_sensitivity=SoilConsts.nitrogen_fixation_cost_thermal_sensitivity,
+        cost_equality_temp=SoilConsts.nitrogen_fixation_cost_equality_temperature,
+    )
+
+    assert np.allclose(expected_cost, actual_cost)
+
+
+def test_calculate_symbiotic_nitrogen_fixation_carbon_cost_bad_temp(
+    dummy_carbon_data, fixture_core_components
+):
+    """Check calculation of nitrogen fixation cost handles bad temperature values."""
+    from virtual_ecosystem.models.soil.constants import SoilConsts
+    from virtual_ecosystem.models.soil.env_factors import (
+        calculate_symbiotic_nitrogen_fixation_carbon_cost,
+    )
+
+    soil_temp = dummy_carbon_data["soil_temperature"][
+        fixture_core_components.layer_structure.index_topsoil_scalar
+    ]
+
+    # Modify some of the soil temps to be below the minimum
+    soil_temp[1] = -23.3
+    soil_temp[3] = -200.0
+
+    expected_cost = [45.8029373, np.inf, 52.689498, np.inf]
+
+    actual_cost = calculate_symbiotic_nitrogen_fixation_carbon_cost(
+        soil_temp=dummy_carbon_data["soil_temperature"][
+            fixture_core_components.layer_structure.index_topsoil_scalar
+        ],
+        cost_at_zero_celsius=SoilConsts.nitrogen_fixation_cost_zero_celcius,
+        infinite_temp_cost_offset=SoilConsts.nitrogen_fixation_cost_infinite_temp_offset,
+        thermal_sensitivity=SoilConsts.nitrogen_fixation_cost_thermal_sensitivity,
+        cost_equality_temp=SoilConsts.nitrogen_fixation_cost_equality_temperature,
+    )
+
+    assert np.allclose(expected_cost, actual_cost)
+
+
 def test_calculate_leaching_rate(dummy_carbon_data, fixture_core_components):
     """Test calculation of solute leaching rates."""
     from virtual_ecosystem.models.soil.constants import SoilConsts

tests/models/soil/test_pools.py

@@ -53,7 +53,7 @@ def test_calculate_all_pool_updates(dummy_carbon_data, fixture_core_components):
         "soil_n_pool_particulate": [1.102338e-5, 6.422491e-5, 0.000131687, 1.461799e-5],
         "soil_n_pool_necromass": [0.00786114, -0.01209909, 0.00432363, -0.00891218],
         "soil_n_pool_maom": [0.00148604, 0.01179891, 0.01365197, 0.0077315],
-        "soil_n_pool_ammonium": [-6.657325e-6, -0.000415127, -0.00021181, -9.40141e-5],
+        "soil_n_pool_ammonium": [0.000952008, 0.019913667, 0.000505414, 0.000455603],
         "soil_n_pool_nitrate": [-0.000293386, -1.292735e-5, -3.576543e-5, -0.000255954],
         "soil_p_pool_dop": [0.000194453, 7.1014337e-5, 0.0001851685, 0.0001017010],
         "soil_p_pool_particulate": [7.22218e-6, -1.13464e-6, 7.86083e-7, 5.85634364e-7],
@@ -649,6 +649,86 @@ def test_calculate_rate_of_denitrification(dummy_carbon_data, fixture_core_compo
     assert np.allclose(actual_rate, expected_rate)
 
 
+def test_calculate_symbiotic_nitrogen_fixation(
+    dummy_carbon_data, fixture_core_components
+):
+    """Check calculation of the rate of symbiotic nitrogen fixation is correct."""
+    from virtual_ecosystem.core.constants import CoreConsts
+    from virtual_ecosystem.models.soil.constants import SoilConsts
+    from virtual_ecosystem.models.soil.pools import (
+        calculate_symbiotic_nitrogen_fixation,
+    )
+
+    expected_fixation = [0.000873306, 0.02021645, 0.00056937, 0.00052116]
+
+    actual_fixation = calculate_symbiotic_nitrogen_fixation(
+        carbon_supply=dummy_carbon_data["nitrogen_fixation_carbon_supply"],
+        soil_temp=dummy_carbon_data["soil_temperature"][
+            fixture_core_components.layer_structure.index_topsoil_scalar
+        ],
+        active_depth=CoreConsts.max_depth_of_microbial_activity,
+        constants=SoilConsts,
+    )
+
+    assert np.allclose(actual_fixation, expected_fixation)
+
+
+def test_calculate_symbiotic_nitrogen_fixation_negative_temps(
+    dummy_carbon_data, fixture_core_components
+):
+    """Check symbiotic nitrogen fixation functions handles negative temperatures."""
+    from virtual_ecosystem.core.constants import CoreConsts
+    from virtual_ecosystem.models.soil.constants import SoilConsts
+    from virtual_ecosystem.models.soil.pools import (
+        calculate_symbiotic_nitrogen_fixation,
+    )
+
+    # Modify some of the soil temps to be below the minimum
+    soil_temp = dummy_carbon_data["soil_temperature"][
+        fixture_core_components.layer_structure.index_topsoil_scalar
+    ]
+    soil_temp[1] = -23.3
+    soil_temp[3] = -200.0
+
+    expected_fixation = [0.000873306, 0.0, 0.00056937, 0.0]
+
+    actual_fixation = calculate_symbiotic_nitrogen_fixation(
+        carbon_supply=dummy_carbon_data["nitrogen_fixation_carbon_supply"],
+        soil_temp=dummy_carbon_data["soil_temperature"][
+            fixture_core_components.layer_structure.index_topsoil_scalar
+        ],
+        active_depth=CoreConsts.max_depth_of_microbial_activity,
+        constants=SoilConsts,
+    )
+
+    assert np.allclose(actual_fixation, expected_fixation)
+
+
+def test_calculate_free_living_nitrogen_fixation(
+    dummy_carbon_data, fixture_core_components
+):
+    """Check calculation of the rate of free-living nitrogen fixation is correct."""
+    from virtual_ecosystem.core.constants import CoreConsts
+    from virtual_ecosystem.models.soil.constants import SoilConsts
+    from virtual_ecosystem.models.soil.pools import (
+        calculate_free_living_nitrogen_fixation,
+    )
+
+    expected_fixation = [8.535774e-5, 0.0001123371, 0.0001478439, 2.845258e-5]
+
+    actual_fixation = calculate_free_living_nitrogen_fixation(
+        soil_temp=dummy_carbon_data["soil_temperature"][
+            fixture_core_components.layer_structure.index_topsoil_scalar
+        ],
+        fixation_at_reference=SoilConsts.free_living_N_fixation_reference_rate,
+        reference_temperature=SoilConsts.free_living_N_fixation_reference_temp,
+        q10_nitrogen_fixation=SoilConsts.free_living_N_fixation_q10_coefficent,
+        active_depth=CoreConsts.max_depth_of_microbial_activity,
+    )
+
+    assert np.allclose(actual_fixation, expected_fixation)
+
+
 def test_calculate_net_formation_of_secondary_P(dummy_carbon_data):
     """Test that calculation of the net formation of secondary P is correct."""
     from virtual_ecosystem.models.soil.pools import (

tests/models/soil/test_soil_model.py

@@ -333,11 +333,11 @@ def test_update(mocker, fixture_soil_model, dummy_carbon_data):
                         [0.86671423, 0.48576345, 0.33406677, 0.09935391], dims="cell_id"
                     ),
                     soil_n_pool_ammonium=DataArray(
-                        [6.47546289e-5, 4.93023264e-3, 1.07788527e-4, 5.11474885e-3],
+                        [0.00053642, 0.01499882, 0.00044842, 0.00538707],
                         dims="cell_id",
                     ),
                     soil_n_pool_nitrate=DataArray(
-                        [0.00188967, 0.00375189, 0.00029596, 0.01290329], dims="cell_id"
+                        [0.00189682, 0.0038413, 0.00031329, 0.01290568], dims="cell_id"
                     ),
                     soil_p_pool_dop=DataArray(
                         [1.68559250e-4, 9.03050817e-5, 3.15038568e-4, 1.66029558e-4],
@@ -444,6 +444,7 @@ def test_order_independance(
         "litter_C_mineralisation_rate",
         "litter_N_mineralisation_rate",
         "litter_P_mineralisation_rate",
+        "nitrogen_fixation_carbon_supply",
     ]
     for not_pool in not_pools:
         new_data[not_pool] = dummy_carbon_data[not_pool]
@@ -527,10 +528,10 @@ def test_construct_full_soil_model(dummy_carbon_data, fixture_core_components):
         0.01179891,
         0.01365197,
         0.0077315,
-        -6.657325e-6,
-        -0.000415127,
-        -0.00021181,
-        -9.40141e-5,
+        0.000952008,
+        0.019913667,
+        0.000505414,
+        0.000455603,
         -0.000293386,
         -1.292735e-5,
         -3.576543e-5,

virtual_ecosystem/data_variables.toml

@@ -1041,6 +1041,13 @@ name = "root_turnover_c_p_ratio"
 unit = "-"
 variable_type = "float"
 
+[[variable]]
+axis = ["spatial"]
+description = "Rate at which plants are supplying carbon to their soil-living nitrogen-fixing symbiotic partners"
+name = "nitrogen_fixation_carbon_supply"
+unit = "kg C m^-2 day^-1"
+variable_type = "float"
+
 [[variable]]
 axis = [] # as_yet_undefined_cohort_setup_axis
 description = "Cell ID of plant cohorts"

virtual_ecosystem/models/plants/plants_model.py

@@ -73,6 +73,7 @@ class PlantsModel(
         "leaf_turnover_c_p_ratio",
         "plant_reproductive_tissue_turnover_c_p_ratio",
         "root_turnover_c_p_ratio",
+        "nitrogen_fixation_carbon_supply",
     ),
     vars_populated_by_first_update=(
         "evapotranspiration",
@@ -92,6 +93,7 @@ class PlantsModel(
         "leaf_turnover_c_p_ratio",
         "plant_reproductive_tissue_turnover_c_p_ratio",
         "root_turnover_c_p_ratio",
+        "nitrogen_fixation_carbon_supply",
     ),
 ):
     """A class defining the plants model.
@@ -467,7 +469,8 @@ def calculate_turnover(self) -> None:
         This function calculates the turnover rate for each plant biomass pool (wood,
         leaves, roots, and reproductive tissues). As well as this the lignin
         concentration, carbon nitrogen ratio and carbon phosphorus ratio of each
-        turnover flow is calculated.
+        turnover flow is calculated. It also returns the rate at which plants supply
+        carbon to their nitrogen fixing symbionts in the soil.
 
         Warning:
             At present, this function literally just returns constant values for each of
@@ -507,3 +510,7 @@ def calculate_turnover(self) -> None:
         self.data["root_turnover_c_p_ratio"] = xr.full_like(
             self.data["elevation"], 656.7
         )
+
+        self.data["nitrogen_fixation_carbon_supply"] = xr.full_like(
+            self.data["elevation"], 0.01
+        )

virtual_ecosystem/models/soil/constants.py

@@ -517,6 +517,57 @@ class SoilConsts(ConstantsDataclass):
     :attr:`denitrification_minimum_temperature`!
     """
 
+    nitrogen_fixation_cost_zero_celcius: float = 59.19651970522086
+    """Cost (in carbon) that plants pay to their symbiotic partners at zero Celsius.
+    
+    Units of [kg C kg N^-1]. This is cost per unit of nitrogen received, and will be
+    higher than the symbiotic partners actually spend to fix the nitrogen. Value is
+    obtained from :cite:t:`brzostek_modeling_2014`.
+    """
+
+    nitrogen_fixation_cost_infinite_temp_offset: float = -0.8034802947791453
+    """Difference in nitrogen fixation cost between zero Celsius and infinite limit.
+    
+    Units of [kg C kg N^-1]. This limit of infinite temperature is not biologically
+    meaningful and is instead just a way of characterising the form of the empirical
+    function. A negative value means that the cost in the infinite temperature limit is
+    higher than at zero Celsius. Value is obtained from
+    :cite:t:`brzostek_modeling_2014`.
+    """
+
+    nitrogen_fixation_cost_thermal_sensitivity: float = 0.27
+    """Sensitivity of symbiotic nitrogen fixation cost to changes in temperature.
+    
+    Units of [C^-1]. Value is obtained from :cite:t:`brzostek_modeling_2014`.
+    """
+
+    nitrogen_fixation_cost_equality_temperature: float = 50.28
+    """Positive temperature at which nitrogen fixation cost is the same at zero Celsius.
+    
+    Units of [C]. Value is obtained from :cite:t:`brzostek_modeling_2014`.
+    """
+
+    free_living_N_fixation_reference_rate: float = 15.0 * 1e-4 / 365.25
+    """Rate at which free living microbes fix nitrogen (at the reference temperature).
+    
+    Units of [kg N m^-2 day^-1]. Value specific to tropical forests, and is taken from
+    :cite:t:`lin_modelling_2000` (with the units adjusted). Should not be changed
+    independently from :attr:`free_living_N_fixation_reference_temp`.
+    """
+
+    free_living_N_fixation_reference_temp: float = 293.15
+    """Temperature reference rate of free-living nitrogen fixation was measured at.
+
+    Units of [K]. Value taken from :cite:t:`lin_modelling_2000`. Should not be changed
+    independently from :attr:`free_living_N_fixation_reference_rate`.
+    """
+
+    free_living_N_fixation_q10_coefficent: float = 3.0
+    """Q10 coefficient for free-living fixation of nitrogen [unitless].
+
+    Value taken from :cite:t:`lin_modelling_2000`.
+    """
+
     primary_phosphorus_breakdown_rate: float = 1.0 / 4.38e6
     """Rate constant for breakdown of primary phosphorus to labile phosphorus [day^-1].