diff --git a/docs/source/api/soil.md b/docs/source/api/soil.md
index 562aef3df..d58988a84 100644
--- a/docs/source/api/soil.md
+++ b/docs/source/api/soil.md
@@ -16,13 +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 {mod}`~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:
-```
+* 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.
diff --git a/docs/source/api/soil/carbon.md b/docs/source/api/soil/carbon.md
new file mode 100644
index 000000000..90077fae1
--- /dev/null
+++ b/docs/source/api/soil/carbon.md
@@ -0,0 +1,23 @@
+---
+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
+---
+
+# API documentation for the {mod}`~virtual_rainforest.soil.carbon` module
+
+```{eval-rst}
+.. automodule:: virtual_rainforest.models.soil.carbon
+    :autosummary:
+    :members:
+```
diff --git a/docs/source/api/soil/model.md b/docs/source/api/soil/model.md
new file mode 100644
index 000000000..429221000
--- /dev/null
+++ b/docs/source/api/soil/model.md
@@ -0,0 +1,23 @@
+---
+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
+---
+
+# API documentation for the {mod}`~virtual_rainforest.soil.model` module
+
+```{eval-rst}
+.. automodule:: virtual_rainforest.models.soil.model
+    :autosummary:
+    :members:
+```
diff --git a/docs/source/index.md b/docs/source/index.md
index 2faea5c9a..f8710513a 100644
--- a/docs/source/index.md
+++ b/docs/source/index.md
@@ -86,8 +86,9 @@ team.
   File readers <api/core/readers.md>
   Core axes <api/core/axes.md>
   Base Model <api/core/model.md>
-  Soil Model <api/soil.md>
-
+  Soil Overview <api/soil.md>
+  Soil Model <api/soil/model.md>
+  Soil Carbon <api/soil/carbon.md>
 ```
 
 ```{eval-rst}
@@ -110,6 +111,7 @@ team.
   Docstring Style <development/documentation/docstring_style.md>
   API Generation <development/documentation/api_generation.md>
   Core Design <development/design/core.md>
+  Adding New Models <development/defining_new_models.md>
 ```
 
 ```{eval-rst}
diff --git a/docs/source/refs.bib b/docs/source/refs.bib
index b933ca2a8..127241c9a 100644
--- a/docs/source/refs.bib
+++ b/docs/source/refs.bib
@@ -183,47 +183,26 @@ @article{Metcalfe2015
 }
 
 @article{cheng_wei_2021,
-  author       = {Cheng, Wei and
-                  Dan, Li and
-                  Deng, Xiangzheng and
-                  Feng, Jinming and
-                  Wang, Yongli and
-                  Peng, Jing and
-                  Tian, Jing and
-                  Qi, Wei and
-                  Liu, Zhu and
-                  Zheng, Xinqi and
-                  Zhou, Demin and
-                  Jiang, Sijian and
-                  Zhao, Haipeng and
-                  Wang, Xiaoyu},
-  title        = {{Global monthly distributions of atmospheric CO2 
-                   concentrations under the historical and future
-                   scenarios}},
-  month        = jun,
-  year         = 2021,
-  note         = {{The data records include 1 file Network Common 
-                   Data Form (NetCDF) format for CO2 distributions in
-                   historical period named
-                   CO2\_1deg\_month\_1850-2013.nc, and 8 files NetCDF
-                   format with the naming convention
-                   CO2\_SSP{XYY}\_2015\_2150.nc, where X and YY are the
-                   shared socioeconomic pathway and radiative forcing
-                   level at 2100, respectively, for CO2 distributions
-                   in the future scenarios. Each NetCDF file includes
-                   3 dimensions: time (month of the year expressed as
-                   days since the first day of 1850, n = 1968 and
-                   1632 for the historical and the future,
-                   respectively); latitude (Degrees North of the
-                   equator [cell centres], n = 180); longitude
-                   (Degrees East of the Prime Meridian [cell
-                   centres], n = 360). Each NetCDF file contains a
-                   monthly variable representing mole fraction of
-                   carbon dioxide in air (variable name: value in the
-                   historical file and CO2 in the future scenario
-                   files) with the unit ppm and the 1º × 1º
-                   resolution.}},
-  publisher    = {Zenodo},
-  doi          = {10.5281/zenodo.5021361},
-  url          = {https://doi.org/10.5281/zenodo.5021361}
+  author  = {Cheng, Wei and
+             Dan, Li and
+             Deng, Xiangzheng and
+             Feng, Jinming and
+             Wang, Yongli and
+             Peng, Jing and
+             Tian, Jing and
+             Qi, Wei and
+             Liu, Zhu and
+             Zheng, Xinqi and
+             Zhou, Demin and
+             Jiang, Sijian and
+             Zhao, Haipeng and
+             Wang, Xiaoyu},
+  title   = {{Global monthly distributions of atmospheric CO2 
+             concentrations under the historical and future
+             scenarios [Data set]}},
+  journal = {Zenodo},
+  month   = jun,
+  year    = 2021,
+  doi     = {10.5281/zenodo.5021361},
+  url     = {https://doi.org/10.5281/zenodo.5021361}
 }
\ No newline at end of file
diff --git a/tests/test_main.py b/tests/test_main.py
index 78d8ae1b6..deb843516 100644
--- a/tests/test_main.py
+++ b/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
 
diff --git a/tests/test_models.py b/tests/test_models.py
index 081f21619..57fad558d 100644
--- a/tests/test_models.py
+++ b/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
 
diff --git a/tests/test_soil_carbon.py b/tests/test_soil_carbon.py
new file mode 100644
index 000000000..f916dc76e
--- /dev/null
+++ b/tests/test_soil_carbon.py
@@ -0,0 +1,325 @@
+"""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 ERROR
+
+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),
+            (
+                (
+                    ERROR,
+                    "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),
+            (
+                (
+                    ERROR,
+                    "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)
+
+
+@pytest.mark.parametrize(
+    argnames=[
+        "maom",
+        "lmwc",
+        "pH",
+        "bulk_density",
+        "moistures",
+        "temperatures",
+        "percent_clay",
+        "dt",
+        "end_maom",
+        "end_lmwc",
+    ],
+    argvalues=[
+        (
+            [2.5, 1.7],
+            [0.05, 0.02],
+            [3.0, 7.5],
+            [1350.0, 1800.0],
+            [0.5, 0.7],
+            [35.0, 37.5],
+            [80.0, 30.0],
+            [2.0 / 24.0, 1.0 / 24.0],
+            [2.500033, 1.70000049],
+            [0.0499668, 0.0199995],
+        ),
+        (
+            [4.5],
+            [0.1],
+            [9.0],
+            [1000.0],
+            [0.6],
+            [40.0],
+            [10.0],
+            [0.5],
+            [4.500071],
+            [0.0999282],
+        ),
+        (
+            [0.5],
+            [0.005],
+            [5.7],
+            [1500.0],
+            [0.2],
+            [25.0],
+            [90.0],
+            [1.0 / 30.0],
+            [0.500000],
+            [0.00499999],
+        ),
+    ],
+)
+def test_update_pools(
+    maom,
+    lmwc,
+    pH,
+    bulk_density,
+    moistures,
+    temperatures,
+    percent_clay,
+    dt,
+    end_maom,
+    end_lmwc,
+):
+    """Test that update_pools runs and generates the correct values."""
+
+    # Initialise soil carbon class
+    soil_carbon = SoilCarbonPools(
+        maom=np.array(maom, dtype=np.float32), lmwc=np.array(lmwc, dtype=np.float32)
+    )
+
+    soil_pH = np.array(pH, dtype=np.float32)
+    soil_BD = np.array(bulk_density, dtype=np.float32)
+    soil_moisture = np.array(moistures, dtype=np.float32)
+    soil_temp = np.array(temperatures, dtype=np.float32)
+    soil_clay = np.array(percent_clay, dtype=np.float32)
+
+    soil_carbon.update_pools(soil_pH, soil_BD, soil_moisture, soil_temp, soil_clay, dt)
+
+    # Check that pools are correctly incremented
+    assert np.allclose(soil_carbon.maom, np.array(end_maom))
+    assert np.allclose(soil_carbon.lmwc, np.array(end_lmwc))
+
+
+@pytest.mark.parametrize(
+    "maom,lmwc,pH,bulk_density,moistures,temperatures,percent_clay,output_l_to_m",
+    [
+        (
+            [2.5, 1.7],
+            [0.05, 0.02],
+            [3.0, 7.5],
+            [1350.0, 1800.0],
+            [0.5, 0.7],
+            [35.0, 37.5],
+            [80.0, 30.0],
+            [0.000397665, 1.178336e-5],
+        ),
+        ([4.5], [0.1], [9.0], [1000.0], [0.6], [40.0], [10.0], [0.0001434178]),
+        ([0.5], [0.005], [5.7], [1500.0], [0.2], [25.0], [90.0], [2.80359e-7]),
+    ],
+)
+def test_mineral_association(
+    maom, lmwc, pH, bulk_density, moistures, temperatures, percent_clay, output_l_to_m
+):
+    """Test that mineral_association runs and generates the correct values."""
+
+    # Initialise soil carbon class
+    soil_carbon = SoilCarbonPools(
+        maom=np.array(maom, dtype=np.float32), lmwc=np.array(lmwc, dtype=np.float32)
+    )
+
+    soil_pH = np.array(pH, dtype=np.float32)
+    soil_BD = np.array(bulk_density, dtype=np.float32)
+    soil_moisture = np.array(moistures, dtype=np.float32)
+    soil_temp = np.array(temperatures, dtype=np.float32)
+    soil_clay = np.array(percent_clay, dtype=np.float32)
+
+    # Then calculate mineral association rate
+    lmwc_to_maom = soil_carbon.mineral_association(
+        soil_pH, soil_BD, soil_moisture, soil_temp, soil_clay
+    )
+
+    # Check that expected values are generated
+    assert np.allclose(lmwc_to_maom, output_l_to_m)
+
+
+@pytest.mark.parametrize(
+    "pH,Q_max,lmwc,output_eqb_maoms",
+    [
+        (
+            [3.0, 7.5],
+            [2.385207e6, 1.980259e6],
+            [0.05, 0.02],
+            [19900.19, 969.4813],
+        ),
+        ([9.0], [647142.61], [0.1], [832.6088]),
+        ([5.7], [2.805371e6], [0.005], [742.4128]),
+    ],
+)
+def test_calculate_equilibrium_maom(pH, Q_max, lmwc, output_eqb_maoms):
+    """Test that equilibrium maom calculation works as expected."""
+    from virtual_rainforest.models.soil.carbon import calculate_equilibrium_maom
+
+    soil_pH = np.array(pH, 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_pH, 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),
+            ((ERROR, "Relative clay content must be expressed as a percentage!"),),
+        ),
+        (
+            [1500.0],
+            [-9.0],
+            [],
+            pytest.raises(ValueError),
+            ((ERROR, "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))
+
+
+@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(temperatures, dtype=np.float32)
+    temp_scalar = convert_temperature_to_scalar(soil_temperature)
+
+    assert np.allclose(temp_scalar, np.array(output_scalars))
+
+
+@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),
+            ((ERROR, "Relative water content cannot go below zero or above one!"),),
+        ),
+        (
+            [2.7],
+            [],
+            pytest.raises(ValueError),
+            ((ERROR, "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
+
+    # 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 np.allclose(moist_scalar, np.array(output_scalars))
+
+    log_check(caplog, expected_log_entries)
diff --git a/virtual_rainforest/__init__.py b/virtual_rainforest/__init__.py
index 65317d397..f312c26d8 100644
--- a/virtual_rainforest/__init__.py
+++ b/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")
diff --git a/virtual_rainforest/core/config.py b/virtual_rainforest/core/config.py
index 73b7019f3..886016082 100644
--- a/virtual_rainforest/core/config.py
+++ b/virtual_rainforest/core/config.py
@@ -1,7 +1,7 @@
 """The :mod:`core.config` module is used to read in the various configuration files,
 validate their contents, and then configure a ready to run instance of the virtual
 rainforest model. The basic details of how this system is used can be found
-:ref:`here<virtual_rainforest/core/config:the configuration module>`.
+:doc:`here </virtual_rainforest/core/config>`.
 
 When a new module is defined a ``JSON`` file should be written, which includes the
 expected configuration tags, their expected types, and any constraints on their values
@@ -518,7 +518,7 @@ def validate_with_defaults(
     This function also adds default values into the configuration dictionary where it is
     appropriate.
 
-     Args:
+    Args:
         config_dict: The complete configuration settings for the particular model
             instance
         comb_schema: Combined schema for all modules that are being configured
diff --git a/virtual_rainforest/plants/__init__.py b/virtual_rainforest/models/plants/__init__.py
similarity index 100%
rename from virtual_rainforest/plants/__init__.py
rename to virtual_rainforest/models/plants/__init__.py
diff --git a/virtual_rainforest/plants/plants_schema.json b/virtual_rainforest/models/plants/plants_schema.json
similarity index 100%
rename from virtual_rainforest/plants/plants_schema.json
rename to virtual_rainforest/models/plants/plants_schema.json
diff --git a/virtual_rainforest/soil/__init__.py b/virtual_rainforest/models/soil/__init__.py
similarity index 100%
rename from virtual_rainforest/soil/__init__.py
rename to virtual_rainforest/models/soil/__init__.py
diff --git a/virtual_rainforest/models/soil/carbon.py b/virtual_rainforest/models/soil/carbon.py
new file mode 100644
index 000000000..3492ee47b
--- /dev/null
+++ b/virtual_rainforest/models/soil/carbon.py
@@ -0,0 +1,255 @@
+"""The `models.soil.carbon` module  simulates the 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.
+"""  # noqa: D205, D415
+
+import numpy as np
+from numpy.typing import NDArray
+
+from virtual_rainforest.core.logger import LOGGER
+from virtual_rainforest.core.model import InitialisationError
+from virtual_rainforest.models.soil.constants import (
+    BindingWithPH,
+    MaxSorptionWithClay,
+    MoistureScalar,
+    TempScalar,
+)
+
+# TODO - I'm basically certain that the paper I've taken this model structure from has
+# not used units consistently (in particular the BINDING_WITH_PH). Down the line I need
+# to track down a reliable parameterisation for this section.
+
+
+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:
+            to_raise = InitialisationError(
+                "Dimension mismatch for initial carbon pools!",
+            )
+            LOGGER.error(to_raise)
+            raise to_raise
+
+        # Check that negative initial values are not given
+        if np.any(maom < 0) or np.any(lmwc < 0):
+            to_raise = InitialisationError(
+                "Initial carbon pools contain at least one negative value!"
+            )
+            LOGGER.error(to_raise)
+            raise to_raise
+
+        self.maom = maom
+        """Mineral associated organic matter pool (kg C m^-3)"""
+        self.lmwc = lmwc
+        """Low molecular weight carbon pool (kg C m^-3)"""
+
+    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 (kg m^-3)
+            soil_moisture: relative water content for each soil grid cell (unitless)
+            soil_temp: soil temperature for each soil grid cell (degrees C)
+            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 (kg m^-3)
+            soil_moisture: relative water content for each soil grid cell (unitless)
+            soil_temp: soil temperature for each soil grid cell (degrees C)
+            percent_clay: Percentage clay for each soil grid cell
+
+        Returns:
+            The net flux from LMWC to MAOM (kg C m^-3 day^-1)
+        """
+
+        # Calculate
+        Q_max = calculate_max_sorption_capacity(bulk_density, percent_clay)
+        equib_maom = calculate_equilibrium_maom(pH, Q_max, self.lmwc)
+
+        # 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 calculate_max_sorption_capacity(
+    bulk_density: NDArray[np.float32],
+    percent_clay: NDArray[np.float32],
+    coef: MaxSorptionWithClay = MaxSorptionWithClay(),
+) -> NDArray[np.float32]:
+    """Calculate maximum sorption capacity based on bulk density and clay content.
+
+    The maximum sorption capacity is the maximum amount of mineral associated organic
+    matter that can exist per unit volume. This expression and its parameters are also
+    drawn from Mayes et al. (2012). In that paper max sorption also depends on Fe
+    content, but we are ignoring this for now.
+
+    Args:
+        bulk_density: bulk density values for each soil grid cell (kg m^-3)
+        percent_clay: Percentage clay for each soil grid cell
+
+    Returns:
+        Maximum sorption capacity (kg m^-3)
+    """
+
+    # Check that negative initial values are not given
+    if np.any(percent_clay > 100.0) or np.any(percent_clay < 0.0):
+        to_raise = ValueError(
+            "Relative clay content must be expressed as a percentage!"
+        )
+        LOGGER.error(to_raise)
+        raise to_raise
+
+    Q_max = bulk_density * 10 ** (coef.slope * np.log10(percent_clay) + coef.intercept)
+    return Q_max
+
+
+def calculate_equilibrium_maom(
+    pH: NDArray[np.float32],
+    Q_max: NDArray[np.float32],
+    lmwc: NDArray[np.float32],
+) -> NDArray[np.float32]:
+    """Calculate equilibrium MAOM concentration based on Langmuir coefficients.
+
+    Equilibrium concentration of mineral associated organic matter (MAOM) is calculated
+    by this function under the assumption that the concentration of low molecular weight
+    carbon (LMWC) is fixed.
+
+    Args:
+        pH: pH values for each soil grid cell Q_max: Maximum sorption capacities (kg
+        m^-3)
+        lmwc: Low molecular weight carbon pool (kg C m^-3)
+
+    Returns:
+        Equilibrium concentration of MAOM (kg C m^-3)
+    """
+
+    binding_coefficient = calculate_binding_coefficient(pH)
+    return (binding_coefficient * Q_max * lmwc) / (1 + lmwc * binding_coefficient)
+
+
+def calculate_binding_coefficient(
+    pH: NDArray[np.float32], coef: BindingWithPH = BindingWithPH()
+) -> NDArray[np.float32]:
+    """Calculate Langmuir binding coefficient based on pH.
+
+    This specific expression and its parameters are drawn from (Mayes et al. (2012)).
+
+    Args:
+        pH: pH values for each soil grid cell
+
+    Returns:
+        Langmuir binding coefficients for mineral association of labile carbon (m^3
+        kg^-1)
+    """
+
+    return 10.0 ** (coef.slope * pH + coef.intercept)
+
+
+def convert_temperature_to_scalar(
+    soil_temp: NDArray[np.float32], coef: TempScalar = TempScalar()
+) -> 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 (degrees C)
+
+    Returns:
+        A scalar that captures the impact of soil temperature on process rates
+    """
+
+    # This expression is drawn from Abramoff et al. (2018)
+    numerator = coef.t_2 + (coef.t_3 / np.pi) * np.arctan(
+        np.pi * (soil_temp - coef.t_1)
+    )
+    denominator = coef.t_2 + (coef.t_3 / np.pi) * np.arctan(
+        np.pi * coef.t_4 * (coef.ref_temp - coef.t_1)
+    )
+
+    return np.divide(numerator, denominator)
+
+
+def convert_moisture_to_scalar(
+    soil_moisture: NDArray[np.float32],
+    coef: MoistureScalar = MoistureScalar(),
+) -> 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: relative water content for each soil grid cell (unitless)
+
+    Returns:
+        A scalar that captures the impact of soil moisture on process rates
+    """
+
+    if np.any(soil_moisture > 1.0) or np.any(soil_moisture < 0.0):
+        to_raise = ValueError(
+            "Relative water content cannot go below zero or above one!"
+        )
+        LOGGER.error(to_raise)
+        raise to_raise
+
+    # This expression is drawn from Abramoff et al. (2018)
+    return 1 / (1 + coef.coefficient * np.exp(-coef.exponent * soil_moisture))
diff --git a/virtual_rainforest/models/soil/constants.py b/virtual_rainforest/models/soil/constants.py
new file mode 100644
index 000000000..2cd27fb84
--- /dev/null
+++ b/virtual_rainforest/models/soil/constants.py
@@ -0,0 +1,53 @@
+"""The `models.soil.carbon` module  simulates the 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.
+"""  # noqa: D205, D415
+
+from dataclasses import dataclass
+
+
+@dataclass
+class BindingWithPH:
+    """From linear regression (Mayes et al. (2012))."""
+
+    slope: float = -0.186
+    """Units of pH^-1."""
+    intercept: float = -0.216
+    """Unit less."""
+
+
+@dataclass
+class MaxSorptionWithClay:
+    """From linear regression (Mayes et al. (2012))."""
+
+    slope: float = 0.483
+    """Units of (% clay)^-1."""
+    intercept: float = 2.328
+    """Unit less."""
+
+
+@dataclass
+class MoistureScalar:
+    """Used in Abramoff et al. (2018), but can't trace it back to anything concrete."""
+
+    coefficient: float = 30.0
+    """Value at zero relative water content (RWC) [unit less]."""
+    exponent: float = 9.0
+    """Units of (RWC)^-1"""
+
+
+@dataclass
+class TempScalar:
+    """Used in Abramoff et al. (2018), but can't trace it back to anything concrete."""
+
+    t_1: float = 15.4
+    """Unclear exactly what this parameter is [degrees C]"""
+    t_2: float = 11.75
+    """Unclear exactly what this parameter is [units unclear]"""
+    t_3: float = 29.7
+    """Unclear exactly what this parameter is [units unclear]"""
+    t_4: float = 0.031
+    """Unclear exactly what this parameter is [units unclear]"""
+    ref_temp: float = 30.0
+    """Reference temperature [degrees C]"""
diff --git a/virtual_rainforest/soil/model.py b/virtual_rainforest/models/soil/model.py
similarity index 100%
rename from virtual_rainforest/soil/model.py
rename to virtual_rainforest/models/soil/model.py
diff --git a/virtual_rainforest/soil/soil_schema.json b/virtual_rainforest/models/soil/soil_schema.json
similarity index 100%
rename from virtual_rainforest/soil/soil_schema.json
rename to virtual_rainforest/models/soil/soil_schema.json