Minimal carbon pool model

.github/workflows/ci.yml

@@ -9,7 +9,7 @@ jobs:
       - uses: actions/checkout@v3
       - uses: actions/setup-python@v4
         with:
-          python-version: "3.9"
+          python-version: "3.9.16"
       - uses: pre-commit/[email protected]
 
   test:

tests/test_main.py

@@ -19,7 +19,7 @@
     select_models,
     vr_run,
 )
-from virtual_rainforest.soil.model import SoilModel
+from virtual_rainforest.models.soil.model import SoilModel
 
 from .conftest import log_check
 

tests/test_models.py

@@ -11,7 +11,7 @@
 from numpy import datetime64, timedelta64
 
 from virtual_rainforest.core.model import BaseModel, InitialisationError
-from virtual_rainforest.soil.model import SoilModel
+from virtual_rainforest.models.soil.model import SoilModel
 
 from .conftest import log_check
 

tests/test_soil_carbon.py

@@ -0,0 +1,129 @@
+"""Test module for soil.carbon.py.
+
+This module tests the functionality of the soil carbon module
+"""
+
+from contextlib import nullcontext as does_not_raise
+from logging import CRITICAL
+
+import numpy as np
+import pytest
+
+from virtual_rainforest.core.model import InitialisationError
+from virtual_rainforest.models.soil.carbon import SoilCarbonPools
+
+from .conftest import log_check
+
+
+@pytest.mark.parametrize(
+    "maom,lmwc,raises,expected_log_entries",
+    [
+        (
+            np.array([23.0, 12.0], dtype=np.float32),
+            np.array([98.0, 7.0], dtype=np.float32),
+            does_not_raise(),
+            (),
+        ),
+        (
+            np.array([23.0, 12.0], dtype=np.float32),
+            np.array([98.0], dtype=np.float32),
+            pytest.raises(InitialisationError),
+            (
+                (
+                    CRITICAL,
+                    "Dimension mismatch for initial carbon pools!",
+                ),
+            ),
+        ),
+        (
+            np.array([23.0, 12.0], dtype=np.float32),
+            np.array([98.0, -24.0], dtype=np.float32),
+            pytest.raises(InitialisationError),
+            (
+                (
+                    CRITICAL,
+                    "Initial carbon pools contain at least one negative value!",
+                ),
+            ),
+        ),
+    ],
+)
+def test_soil_carbon_class(caplog, maom, lmwc, raises, expected_log_entries):
+    """Test SoilCarbon class initialisation."""
+
+    # Check that initialisation fails (or doesn't) as expected
+    with raises:
+        soil_carbon = SoilCarbonPools(maom, lmwc)
+
+        assert (soil_carbon.maom == maom).all()
+        assert (soil_carbon.lmwc == lmwc).all()
+
+    log_check(caplog, expected_log_entries)
+
+
+def test_update_pools():
+    """Test that update_pools runs and generates the correct values."""
+
+    # Initialise soil carbon class
+    maom = np.array([23.0, 23.0], dtype=np.float32)
+    lmwc = np.array([98.0, 55.0], dtype=np.float32)
+    soil_carbon = SoilCarbonPools(maom, lmwc)
+
+    # Define all the required variables to run function
+    pH = np.array([7.0, 7.0], dtype=np.float32)
+    bulk_density = np.array([1350, 1350], dtype=np.float32)
+    percent_clay = np.array([50.0, 50.0], dtype=np.float32)
+    soil_moisture = np.array([0.5, 0.5], dtype=np.float32)
+    soil_temp = np.array([35.0, 35.0], dtype=np.float32)
+    dt = 2.0 / 24.0
+
+    soil_carbon.update_pools(
+        pH, bulk_density, soil_moisture, soil_temp, percent_clay, dt
+    )
+
+    # Check that pools are correctly incremented
+    assert np.allclose(soil_carbon.maom, np.array([28.8263, 25.7322]))
+    assert np.allclose(soil_carbon.lmwc, np.array([92.1736, 52.2677]))
+
+
+def test_mineral_association():
+    """Test that mineral_association runs and generates the correct values."""
+
+    # Initialise soil carbon class
+    maom = np.array([23.0, 23.0], dtype=np.float32)
+    lmwc = np.array([98.0, 55.0], dtype=np.float32)
+    soil_carbon = SoilCarbonPools(maom, lmwc)
+
+    # Define all the required variables to run function
+    pH = np.array([7.0, 7.0], dtype=np.float32)
+    bulk_density = np.array([1350, 1350], dtype=np.float32)
+    percent_clay = np.array([50.0, 50.0], dtype=np.float32)
+    soil_moisture = np.array([0.5, 0.5], dtype=np.float32)
+    soil_temp = np.array([35.0, 35.0], dtype=np.float32)
+
+    lmwc_to_maom = soil_carbon.mineral_association(
+        pH, bulk_density, soil_moisture, soil_temp, percent_clay
+    )
+
+    # Check that expected values are generated
+    assert np.allclose(lmwc_to_maom, np.array([69.9158, 32.7868]))
+
+
+def test_convert_temperature_to_scalar():
+    """Test that scalar_temperature runs and generates the correct value."""
+    from virtual_rainforest.models.soil.carbon import convert_temperature_to_scalar
+
+    soil_temperature = np.array([35.0, 37.5], dtype=np.float32)
+    temp_scalar = convert_temperature_to_scalar(soil_temperature)
+
+    assert np.allclose(temp_scalar, np.array([1.27113, 1.27196]))
+
+
+def test_convert_moisture_to_scalar():
+    """Test that scalar_moisture runs and generates the correct value."""
+    from virtual_rainforest.models.soil.carbon import convert_moisture_to_scalar
+
+    soil_moisture = np.array([0.5, 0.7], dtype=np.float32)
+    moist_scalar = convert_moisture_to_scalar(soil_moisture)
+
+    assert np.allclose(moist_scalar, np.array([0.750035, 0.947787]))

virtual_rainforest/__init__.py

@@ -2,10 +2,10 @@
 
 # Import all module schema here to ensure that they are added to the registry
 from virtual_rainforest.core import schema  # noqa
-from virtual_rainforest.plants import schema  # noqa
-from virtual_rainforest.soil import schema  # noqa
+from virtual_rainforest.models.plants import schema  # noqa
+from virtual_rainforest.models.soil import schema  # noqa
 
 # Import models here so that they also end up in the registry
-from virtual_rainforest.soil.model import SoilModel  # noqa
+from virtual_rainforest.models.soil.model import SoilModel  # noqa
 
 __version__ = importlib.metadata.version("virtual_rainforest")

virtual_rainforest/models/soil/carbon.py

@@ -0,0 +1,201 @@
+"""The `soil.carbon` module.
+
+This module simulates the radiation soil carbon cycle for the Virtual Rainforest. At the
+moment only two pools are modelled, these are low molecular weight carbon (LMWC) and
+mineral associated organic matter (MAOM). More pools and their interactions will be
+added at a later date.
+"""
+
+import numpy as np
+from numpy.typing import NDArray
+
+from virtual_rainforest.core.logger import log_and_raise
+from virtual_rainforest.core.model import InitialisationError
+
+# from core.constants import CONSTANTS as C
+# but for meanwhile define all the constants needed here
+BINDING_WITH_PH = {
+    "slope": -0.186,  # unitless
+    "intercept": -0.216,  # unitless
+}
+"""From linear regression (Mayes et al. (2012))."""
+MAX_SORPTION_WITH_CLAY = {
+    "slope": 0.483,  # unitless
+    "intercept": 2.328,  # unitless
+}
+"""From linear regression (Mayes et al. (2012))."""
+MOISTURE_SCALAR = {
+    "coefficient": 30.0,  # unitless
+    "exponent": 9.0,  # unitless
+}
+"""Used in Abramoff et al. (2018), but can't trace it back to anything more concrete."""
+TEMP_SCALAR = np.array([15.4, 11.75, 29.7, 0.031, 30.0])
+"""Used in Abramoff et al. (2018), but can't trace it back to anything more concrete.
+
+Values:
+    t_1: (degrees C)
+    t_2: (unit unclear)
+    t_3: (unit unclear)
+    t_4: (unit unclear)
+    ref_temp: (degrees C)
+"""
+
+
+class SoilCarbonPools:
+    """Class containing the full set of soil carbon pools.
+
+    At the moment, only two pools are included. Functions exist for the transfer of
+    carbon between these pools, but not with either the yet to be implemented soil
+    carbon pools, other pools in the soil module, or other modules.
+    """
+
+    def __init__(self, maom: NDArray[np.float32], lmwc: NDArray[np.float32]) -> None:
+        """Initialise set of carbon pools."""
+
+        # Check that arrays are of equal size and shape
+        if maom.shape != lmwc.shape:
+            log_and_raise(
+                "Dimension mismatch for initial carbon pools!",
+                InitialisationError,
+            )
+
+        # Check that negative initial values are not given
+        if np.any(maom < 0) or np.any(lmwc < 0):
+            log_and_raise(
+                "Initial carbon pools contain at least one negative value!",
+                InitialisationError,
+            )
+
+        self.maom = maom
+        """Mineral associated organic matter pool"""
+        self.lmwc = lmwc
+        """Low molecular weight carbon pool"""
+
+    def update_pools(
+        self,
+        pH: NDArray[np.float32],
+        bulk_density: NDArray[np.float32],
+        soil_moisture: NDArray[np.float32],
+        soil_temp: NDArray[np.float32],
+        percent_clay: NDArray[np.float32],
+        dt: float,
+    ) -> None:
+        """Update all soil carbon pools.
+
+        This function calls lower level functions which calculate the transfers between
+        pools. When all transfers have been calculated the net transfer is used to
+        update the soil pools.
+
+        Args:
+            pH: pH values for each soil grid cell
+            bulk_density: bulk density values for each soil grid cell
+            soil_moisture: soil moisture for each soil grid cell
+            soil_temp: soil temperature for each soil grid cell
+            percent_clay: Percentage clay for each soil grid cell
+            dt: time step (days)
+        """
+        # TODO - Add interactions which involve the three missing carbon pools
+
+        lmwc_to_maom = self.mineral_association(
+            pH, bulk_density, soil_moisture, soil_temp, percent_clay
+        )
+
+        # Once changes are determined update all pools
+        self.lmwc -= lmwc_to_maom * dt
+        self.maom += lmwc_to_maom * dt
+
+    def mineral_association(
+        self,
+        pH: NDArray[np.float32],
+        bulk_density: NDArray[np.float32],
+        soil_moisture: NDArray[np.float32],
+        soil_temp: NDArray[np.float32],
+        percent_clay: NDArray[np.float32],
+    ) -> NDArray[np.float32]:
+        """Calculates net rate of LMWC association with soil minerals.
+
+        Following Abramoff et al. (2018), mineral adsorption of carbon is controlled by
+        a Langmuir saturation function. At present, binding affinity and Q_max are
+        recalculated on every function called based on pH, bulk density and clay
+        content. Once a decision has been reached as to how fast pH and bulk density
+        will change (if at all), this calculation may need to be moved elsewhere.
+
+        Args:
+            pH: pH values for each soil grid cell
+            bulk_density: bulk density values for each soil grid cell
+            soil_moisture: soil moisture for each soil grid cell
+            soil_temp: soil temperature for each soil grid cell
+            percent_clay: Percentage clay for each soil grid cell
+
+        Returns:
+            lmwc_to_maom: The net flux from LMWC to MAOM
+        """
+
+        # This expression is drawn from (Mayes et al. (2012))
+        binding_affinity = 10.0 ** (
+            BINDING_WITH_PH["slope"] * pH + BINDING_WITH_PH["intercept"]
+        )
+
+        # This expression is also drawn from Mayes et al. (2012)
+        # Original paper also depends on Fe concentration, but we are ignoring this for
+        # now
+        Q_max = bulk_density * 10 ** (
+            MAX_SORPTION_WITH_CLAY["slope"] * np.log10(percent_clay)
+            + MAX_SORPTION_WITH_CLAY["intercept"]
+        )
+
+        # Using the above calculate the equilibrium MAOM pool
+        equib_maom = (binding_affinity * Q_max * self.lmwc) / (
+            1 + self.lmwc * binding_affinity
+        )
+
+        # Find scalar factors that multiple rates
+        temp_scalar = convert_temperature_to_scalar(soil_temp)
+        moist_scalar = convert_moisture_to_scalar(soil_moisture)
+
+        return temp_scalar * moist_scalar * self.lmwc * (equib_maom - self.maom) / Q_max
+
+
+def convert_temperature_to_scalar(
+    soil_temp: NDArray[np.float32],
+) -> NDArray[np.float32]:
+    """Convert soil temperature into a factor to multiply rates by.
+
+    This form is used in Abramoff et al. (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.
+
+    Args:
+       soil_temp: soil temperature for each soil grid cell
+    """
+
+    t_1, t_2, t_3, t_4, ref_temp = TEMP_SCALAR
+
+    # This expression is drawn from Abramoff et al. (2018)
+    numerator = t_2 + (t_3 / np.pi) * np.arctan(np.pi * (soil_temp - t_1))
+    denominator = t_2 + (t_3 / np.pi) * np.arctan(np.pi * t_4 * (ref_temp - t_1))
+
+    return np.divide(numerator, denominator)
+
+
+def convert_moisture_to_scalar(
+    soil_moisture: NDArray[np.float32],
+) -> NDArray[np.float32]:
+    """Convert soil moisture into a factor to multiply rates by.
+
+    This form is used in Abramoff et al. (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.
+
+    Args:
+        soil_moisture: soil moisture for each soil grid cell
+    """
+
+    # This expression is drawn from Abramoff et al. (2018)
+    return 1 / (
+        1
+        + MOISTURE_SCALAR["coefficient"]
+        * np.exp(-MOISTURE_SCALAR["exponent"] * soil_moisture)
+    )