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.16"
+          python-version: "3.9"
       - uses: pre-commit/[email protected]
 
   test:

docs/source/api/soil.md

@@ -16,14 +16,13 @@ kernelspec:
 
 # API reference for `soil` modules
 
-This page contains the detailed documentation of the functions and classes in the
-`virtual_rainforest.soil` modules.
+The {mod}`~virtual_rainforest.models.soil` module is one of the component models of the
+Virtual Rainforest. It is comprised of a number of submodules.
 
-## The `virtual_rainforest.soil.model` module
+Each of the soil sub-modules has its own API reference page:
 
-```{eval-rst}
-.. automodule:: virtual_rainforest.soil.model
-    :autosummary:
-    :members:
-    :special-members: __init__
-```
+* The {mod}`~virtual_rainforest.models.soil.model` submodule instantiates the SoilModel
+  class which consolidates the functionality of the soil module into a single class,
+  which the high level functions of the Virtual Rainforest can then make use of.
+* The {mod}`~virtual_rainforest.models.soil.carbon` provides a model for the soil carbon
+  cycle.

docs/source/api/soil/carbon.md

@@ -0,0 +1,24 @@
+---
+jupytext:
+  cell_metadata_filter: -all
+  formats: md:myst
+  main_language: python
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.13.8
+kernelspec:
+  display_name: vr_python3
+  language: python
+  name: vr_python3
+---
+
+<!-- markdownlint-disable MD041 -->
+
+```{eval-rst}
+.. automodule:: virtual_rainforest.models.soil.carbon
+    :autosummary:
+    :members:
+    :special-members: __init__
+```

docs/source/api/soil/model.md

@@ -0,0 +1,24 @@
+---
+jupytext:
+  cell_metadata_filter: -all
+  formats: md:myst
+  main_language: python
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.13.8
+kernelspec:
+  display_name: vr_python3
+  language: python
+  name: vr_python3
+---
+
+<!-- markdownlint-disable MD041 -->
+
+```{eval-rst}
+.. automodule:: virtual_rainforest.models.soil.model
+    :autosummary:
+    :members:
+    :special-members: __init__
+```

docs/source/index.md

@@ -85,7 +85,9 @@ team.
   File readers <api/core/readers.md>
   Core axes <api/core/axes.md>
   Base Model <api/core/model.md>
-  Soil <api/soil.md>
+  Soil Overview <api/soil.md>
+  Soil Model <api/soil/model.md>
+  Soil Carbon <api/soil/carbon.md>
 ```
 
 ```{eval-rst}

tests/test_soil_carbon.py

@@ -5,17 +5,12 @@
 
 from contextlib import nullcontext as does_not_raise
 from logging import CRITICAL
-from math import isclose
 
 import numpy as np
 import pytest
 
 from virtual_rainforest.core.model import InitialisationError
-from virtual_rainforest.models.soil.carbon import (
-    SoilCarbonPools,
-    convert_moisture_to_scalar,
-    convert_temperature_to_scalar,
-)
+from virtual_rainforest.models.soil.carbon import SoilCarbonPools
 
 from .conftest import log_check
 
@@ -66,72 +61,233 @@ def test_soil_carbon_class(caplog, maom, lmwc, raises, expected_log_entries):
     log_check(caplog, expected_log_entries)
 
 
-def test_update_pools():
+@pytest.mark.parametrize(
+    "maom,lmwc,lmwc_to_maom,dt,end_maom,end_lmwc",
+    [
+        (
+            [2.5, 1.7],
+            [0.05, 0.02],
+            [0.000397665, 1.178336e-5],
+            [2.0 / 24.0, 1.0 / 24.0],
+            [2.500033, 1.70000049],
+            [0.0499668, 0.0199995],
+        ),
+        ([4.5], [0.1], [0.0001434178], [0.5], [4.500071], [0.0999282]),
+        ([0.5], [0.005], [2.80359e-7], [1.0 / 30.0], [0.500000], [0.00499999]),
+    ],
+)
+def test_update_pools(mocker, maom, lmwc, lmwc_to_maom, dt, end_maom, end_lmwc):
     """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
+    soil_carbon = SoilCarbonPools(
+        maom=np.array(maom, dtype=np.float32), lmwc=np.array(lmwc, dtype=np.float32)
     )
 
+    # Mock required values
+    mock_l_to_m = mocker.patch(
+        "virtual_rainforest.models.soil.carbon.SoilCarbonPools.mineral_association"
+    )
+    mock_l_to_m.return_value = np.array(lmwc_to_maom, dtype=np.float32)
+
+    soil_carbon.update_pools([], [], [], [], [], dt)
+
     # Check that pools are correctly incremented
-    assert isclose(soil_carbon.maom[0].item(), 28.82632255)
-    assert isclose(soil_carbon.lmwc[0].item(), 92.17367553)
-    assert isclose(soil_carbon.maom[1].item(), 25.73223495)
-    assert isclose(soil_carbon.lmwc[1].item(), 52.26776504)
+    assert np.allclose(soil_carbon.maom, np.array(end_maom))
+    assert np.allclose(soil_carbon.lmwc, np.array(end_lmwc))
 
 
-def test_mineral_association():
+@pytest.mark.parametrize(
+    "maom,lmwc,temp_scalar,moist_scalar,equib_maom,Q_max,output_l_to_m",
+    [
+        (
+            [2.5, 1.7],
+            [0.05, 0.02],
+            [1.27113, 1.27196],
+            [0.750035, 0.947787],
+            [19900.19, 969.4813],
+            [2.385207e6, 1.980259e6],
+            [0.000397665, 1.178336e-5],
+        ),
+        ([4.5], [0.1], [1.27263], [0.880671], [832.6088], [647142.61], [0.0001434178]),
+        ([0.5], [0.005], [1.26344], [0.167814], [742.4128], [2.805371e6], [2.80359e-7]),
+    ],
+)
+def test_mineral_association(
+    mocker, maom, lmwc, temp_scalar, moist_scalar, equib_maom, Q_max, output_l_to_m
+):
     """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
+    soil_carbon = SoilCarbonPools(
+        maom=np.array(maom, dtype=np.float32), lmwc=np.array(lmwc, dtype=np.float32)
+    )
+
+    # Mock required values
+    mock_t_scalar = mocker.patch(
+        "virtual_rainforest.models.soil.carbon.convert_temperature_to_scalar"
+    )
+    mock_t_scalar.return_value = np.array(temp_scalar, dtype=np.float32)
+    mock_m_scalar = mocker.patch(
+        "virtual_rainforest.models.soil.carbon.convert_moisture_to_scalar"
     )
+    mock_m_scalar.return_value = np.array(moist_scalar, dtype=np.float32)
+    mock_equib_maom = mocker.patch(
+        "virtual_rainforest.models.soil.carbon.calculate_equilibrium_maom"
+    )
+    mock_equib_maom.return_value = np.array(equib_maom, dtype=np.float32)
+    mock_Q_max = mocker.patch(
+        "virtual_rainforest.models.soil.carbon.calculate_max_sorption_capacity"
+    )
+    mock_Q_max.return_value = np.array(Q_max, dtype=np.float32)
+
+    # Then calculate mineral association rate
+    lmwc_to_maom = soil_carbon.mineral_association([], [], [], [], [])
 
     # Check that expected values are generated
-    assert isclose(lmwc_to_maom[0].item(), 69.9158630)
-    assert isclose(lmwc_to_maom[1].item(), 32.78682708)
+    assert np.allclose(lmwc_to_maom, output_l_to_m)
+
+
+@pytest.mark.parametrize(
+    "binding_coefficients,Q_max,lmwc,output_eqb_maoms",
+    [
+        (
+            [0.16826738, 0.02449064],
+            [2.385207e6, 1.980259e6],
+            [0.05, 0.02],
+            [19900.19, 969.4813],
+        ),
+        ([0.0128825], [647142.61], [0.1], [832.6088]),
+        ([0.05294197], [2.805371e6], [0.005], [742.4128]),
+    ],
+)
+def test_calculate_equilibrium_maom(
+    mocker, binding_coefficients, Q_max, lmwc, output_eqb_maoms
+):
+    """Test that equilibrium maom calculation works as expected."""
+    from virtual_rainforest.models.soil.carbon import calculate_equilibrium_maom
+
+    # Configure the mock to return specific binding coefficients
+    mock_binding = mocker.patch(
+        "virtual_rainforest.models.soil.carbon.calculate_binding_coefficient"
+    )
+    mock_binding.return_value = np.array(binding_coefficients, dtype=np.float32)
+
+    soil_Q_max = np.array(Q_max, dtype=np.float32)
+    soil_lmwc = np.array(lmwc, dtype=np.float32)
+
+    equib_maoms = calculate_equilibrium_maom([], soil_Q_max, soil_lmwc)
+    assert np.allclose(equib_maoms, np.array(output_eqb_maoms))
+
+
+@pytest.mark.parametrize(
+    "bulk_density,percent_clay,output_capacities,raises,expected_log_entries",
+    [
+        (
+            [1350.0, 1800.0],
+            [80.0, 30.0],
+            [2.385207e6, 1.980259e6],
+            does_not_raise(),
+            (),
+        ),
+        ([1000.0], [10.0], [647142.61], does_not_raise(), ()),
+        ([1500.0], [90.0], [2.805371e6], does_not_raise(), ()),
+        (
+            [1500.0],
+            [156.0],
+            [],
+            pytest.raises(ValueError),
+            ((CRITICAL, "Relative clay content must be expressed as a percentage!"),),
+        ),
+        (
+            [1500.0],
+            [-9.0],
+            [],
+            pytest.raises(ValueError),
+            ((CRITICAL, "Relative clay content must be expressed as a percentage!"),),
+        ),
+    ],
+)
+def test_calculate_max_sorption_capacity(
+    caplog, bulk_density, percent_clay, output_capacities, raises, expected_log_entries
+):
+    """Test that max sorption capacity calculation works as expected."""
+    from virtual_rainforest.models.soil.carbon import calculate_max_sorption_capacity
+
+    # Check that initialisation fails (or doesn't) as expected
+    with raises:
+        soil_BD = np.array(bulk_density, dtype=np.float32)
+        soil_clay = np.array(percent_clay, dtype=np.float32)
+        max_capacities = calculate_max_sorption_capacity(soil_BD, soil_clay)
+
+        assert np.allclose(max_capacities, np.array(output_capacities))
+
+    log_check(caplog, expected_log_entries)
+
+
+@pytest.mark.parametrize(
+    "pH,output_coefs",
+    [
+        ([3.0, 7.5], [0.16826738, 0.02449064]),
+        ([9.0], [0.0128825]),
+        ([5.7], [0.05294197]),
+    ],
+)
+def test_calculate_binding_coefficient(pH, output_coefs):
+    """Test that Langmuir binding coefficient calculation works as expected."""
+    from virtual_rainforest.models.soil.carbon import calculate_binding_coefficient
+
+    soil_pH = np.array(pH, dtype=np.float32)
+    binding_coefs = calculate_binding_coefficient(soil_pH)
+
+    assert np.allclose(binding_coefs, np.array(output_coefs))
 
 
-def test_convert_temperature_to_scalar():
+@pytest.mark.parametrize(
+    "temperatures,output_scalars",
+    [([35.0, 37.5], [1.27113, 1.27196]), ([40.0], [1.27263]), ([25.0], [1.26344])],
+)
+def test_convert_temperature_to_scalar(temperatures, output_scalars):
     """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)
+    soil_temperature = np.array(temperatures, dtype=np.float32)
     temp_scalar = convert_temperature_to_scalar(soil_temperature)
 
-    assert isclose(temp_scalar[0].item(), 1.271131634)
-    assert isclose(temp_scalar[1].item(), 1.271966338)
+    assert np.allclose(temp_scalar, np.array(output_scalars))
 
 
-def test_convert_moisture_to_scalar():
+@pytest.mark.parametrize(
+    "moistures,output_scalars,raises,expected_log_entries",
+    [
+        ([0.5, 0.7], [0.750035, 0.947787], does_not_raise(), ()),
+        ([0.6], [0.880671], does_not_raise(), ()),
+        ([0.2], [0.167814], does_not_raise(), ()),
+        (
+            [-0.2],
+            [],
+            pytest.raises(ValueError),
+            ((CRITICAL, "Relative water content cannot go below zero or above one!"),),
+        ),
+        (
+            [2.7],
+            [],
+            pytest.raises(ValueError),
+            ((CRITICAL, "Relative water content cannot go below zero or above one!"),),
+        ),
+    ],
+)
+def test_convert_moisture_to_scalar(
+    caplog, moistures, output_scalars, raises, expected_log_entries
+):
     """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)
+    # Check that initialisation fails (or doesn't) as expected
+    with raises:
+        soil_moisture = np.array(moistures, dtype=np.float32)
+        moist_scalar = convert_moisture_to_scalar(soil_moisture)
 
-    assert isclose(moist_scalar[0].item(), 0.750035703)
-    assert isclose(moist_scalar[1].item(), 0.947787225)
+        assert np.allclose(moist_scalar, np.array(output_scalars))
+
+    log_check(caplog, expected_log_entries)