diff --git a/docs/source/api/soil/env_factors.md b/docs/source/api/soil/env_factors.md
new file mode 100644
index 000000000..be4462562
--- /dev/null
+++ b/docs/source/api/soil/env_factors.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.models.soil.env_factors` module
+
+```{eval-rst}
+.. automodule:: virtual_rainforest.models.soil.env_factors
+    :autosummary:
+    :members:
+```
diff --git a/docs/source/index.md b/docs/source/index.md
index efdb0e85c..9b0ce194c 100644
--- a/docs/source/index.md
+++ b/docs/source/index.md
@@ -100,6 +100,7 @@ team.
   Soil Overview <api/soil.md>
   Soil Model <api/soil/soil_model.md>
   Soil Carbon <api/soil/carbon.md>
+  Soil Environmental Factors <api/soil/env_factors.md>
   Soil Constants <api/soil/constants.md>
   Abiotic Simple Overview <api/abiotic_simple.md>
   Abiotic Simple Model <api/abiotic_simple/abiotic_simple_model.md>
diff --git a/docs/source/refs.bib b/docs/source/refs.bib
index bc52a5b6e..d0dd1174e 100644
--- a/docs/source/refs.bib
+++ b/docs/source/refs.bib
@@ -1,3 +1,36 @@
+@article{porporato_hydrologic_2003,
+	title = {Hydrologic controls on soil carbon and nitrogen cycles. {I}. {Modeling} scheme},
+	volume = {26},
+	issn = {03091708},
+	url = {https://linkinghub.elsevier.com/retrieve/pii/S0309170802000945},
+	doi = {10.1016/S0309-1708(02)00094-5},
+	language = {en},
+	number = {1},
+	urldate = {2023-11-14},
+	journal = {Advances in Water Resources},
+	author = {Porporato, A and D’Odorico, P and Laio, F and Rodriguez-Iturbe, I},
+	month = jan,
+	year = {2003},
+	pages = {45--58},
+}
+
+@article{orwin_organic_2011,
+	title = {Organic nutrient uptake by mycorrhizal fungi enhances ecosystem carbon storage: a model-based assessment},
+	volume = {14},
+	issn = {1461023X},
+	shorttitle = {Organic nutrient uptake by mycorrhizal fungi enhances ecosystem carbon storage},
+	url = {https://onlinelibrary.wiley.com/doi/10.1111/j.1461-0248.2011.01611.x},
+	doi = {10.1111/j.1461-0248.2011.01611.x},
+	language = {en},
+	number = {5},
+	urldate = {2023-10-18},
+	journal = {Ecology Letters},
+	author = {Orwin, Kate H. and Kirschbaum, Miko U. F. and St John, Mark G. and Dickie, Ian A.},
+	month = may,
+	year = {2011},
+	pages = {493--502},
+}
+
 
 @book{monteith_light_1969,
 	address = {Madison, Wisconsin, U.S.A.},
@@ -333,22 +366,6 @@ @book{campbell_introduction_1998
 	year = {1998},
 }
 
-@article{mayes_relation_2012,
-	title = {Relation between {Soil} {Order} and {Sorption} of {Dissolved} {Organic} {Carbon} in {Temperate} {Subsoils}},
-	volume = {76},
-	issn = {0361-5995, 1435-0661},
-	url = {https://onlinelibrary.wiley.com/doi/10.2136/sssaj2011.0340},
-	doi = {10.2136/sssaj2011.0340},
-	language = {en},
-	number = {3},
-	urldate = {2023-02-22},
-	journal = {Soil Science Society of America Journal},
-	author = {Mayes, Melanie A. and Heal, Katherine R. and Brandt, Craig C. and Phillips, Jana R. and Jardine, Philip M.},
-	month = may,
-	year = {2012},
-	pages = {1027--1037},
-}
-
 @article{davis_simple_2017,
 	title = {Simple process-led algorithms for simulating habitats ({SPLASH} v.1.0): robust indices of radiation, evapotranspiration and plant-available moisture},
 	volume = {10},
diff --git a/tests/conftest.py b/tests/conftest.py
index f2b432ccf..dac8c7c1a 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -172,12 +172,22 @@ def dummy_carbon_data(layer_roles_fixture):
     """Microbial biomass (carbon) pool (kg C m^-3)"""
     data["soil_c_pool_pom"] = DataArray([0.1, 1.0, 0.7, 0.35], dims=["cell_id"])
     """Particulate organic matter pool (kg C m^-3)"""
+    data["soil_enzyme_pom"] = DataArray(
+        [0.022679, 0.009576, 0.050051, 0.003010], dims=["cell_id"]
+    )
+    """Soil enzyme that breaks down particulate organic matter (kg C m^-3)"""
+    data["soil_enzyme_maom"] = DataArray(
+        [0.0356, 0.0117, 0.02509, 0.00456], dims=["cell_id"]
+    )
+    """Soil enzyme that breaks down mineral associated organic matter (kg C m^-3)"""
     data["pH"] = DataArray([3.0, 7.5, 9.0, 5.7], dims=["cell_id"])
     data["bulk_density"] = DataArray([1350.0, 1800.0, 1000.0, 1500.0], dims=["cell_id"])
-    data["percent_clay"] = DataArray([80.0, 30.0, 10.0, 90.0], dims=["cell_id"])
+    data["clay_fraction"] = DataArray([0.8, 0.3, 0.1, 0.9], dims=["cell_id"])
     data["litter_C_mineralisation_rate"] = DataArray(
         [0.00212106, 0.00106053, 0.00049000, 0.0055], dims=["cell_id"]
     )
+    # Data for combined vertical flow (for entire timestep)
+    data["vertical_flow"] = DataArray([3.0, 15.0, 75.0, 47.7], dims=["cell_id"])
 
     # The layer dependant data has to be handled separately
     data["soil_moisture"] = xr.concat(
@@ -186,7 +196,7 @@ def dummy_carbon_data(layer_roles_fixture):
             # At present the soil model only uses the top soil layer, so this is the
             # only one with real test values in
             DataArray(
-                [[0.472467929, 0.399900047, 0.256053640, 0.153616897]],
+                [[0.9304620050, 0.787549327, 0.504263188, 0.302527807]],
                 dims=["layers", "cell_id"],
             ),
             DataArray(np.full((1, 4), np.nan), dims=["layers", "cell_id"]),
@@ -200,7 +210,6 @@ def dummy_carbon_data(layer_roles_fixture):
             "cell_id": data.grid.cell_id,
         }
     )
-    # TODO - Eventually this should replace the dummy soil moisture entirely
     data["matric_potential"] = xr.concat(
         [
             DataArray(np.full((13, 4), np.nan), dims=["layers", "cell_id"]),
diff --git a/tests/core/test_constants_class.py b/tests/core/test_constants_class.py
index 30b4d12b4..1d3837145 100644
--- a/tests/core/test_constants_class.py
+++ b/tests/core/test_constants_class.py
@@ -29,14 +29,14 @@ class Test(ConstantsDataclass):  # type: ignore [misc]
         pytest.param(
             {},
             does_not_raise(),
-            123.4,
+            0.25,
             (),
             id="defaults_with_no_config",
         ),
         pytest.param(
-            {"placeholder": 432.1},
+            {"depth_of_active_soil_layer": 1.55},
             does_not_raise(),
-            432.1,
+            1.55,
             (),
             id="configured",
         ),
@@ -70,6 +70,6 @@ def test_ConstantsDataclass_from_config(caplog, config, raises, exp_val, exp_log
         constants_instance = CoreConsts.from_config(config)
 
         if isinstance(raises, does_not_raise):
-            assert constants_instance.placeholder == exp_val
+            assert constants_instance.depth_of_active_soil_layer == exp_val
 
         log_check(caplog=caplog, expected_log=exp_log)
diff --git a/tests/core/test_constants_loader.py b/tests/core/test_constants_loader.py
index 8d499cfd1..9764f99fa 100644
--- a/tests/core/test_constants_loader.py
+++ b/tests/core/test_constants_loader.py
@@ -17,16 +17,16 @@
             "core",
             "CoreConsts",
             does_not_raise(),
-            123.4,
+            0.25,
             ((INFO, "Initialised core.CoreConsts from config"),),
             id="default_values",
         ),
         pytest.param(
-            "[core.constants.CoreConsts]\nplaceholder=432.1",
+            "[core.constants.CoreConsts]\ndepth_of_active_soil_layer=1.5",
             "core",
             "CoreConsts",
             does_not_raise(),
-            432.1,
+            1.5,
             ((INFO, "Initialised core.CoreConsts from config"),),
             id="configured_value",
         ),
@@ -96,7 +96,7 @@ def test_load_constants(
             assert isinstance(constants_instance, CoreConsts)
             # The unconfigurable zero_Celsius should take the default value
             assert constants_instance.zero_Celsius == constants.zero_Celsius
-            # Check the placeholder constant has been configured
-            assert constants_instance.placeholder == exp_val
+            # Check the depth_of_active_soil_layer constant has been configured
+            assert constants_instance.depth_of_active_soil_layer == exp_val
 
         log_check(caplog=caplog, expected_log=exp_log)
diff --git a/tests/core/test_data.py b/tests/core/test_data.py
index 59c8423f2..e982af0c3 100644
--- a/tests/core/test_data.py
+++ b/tests/core/test_data.py
@@ -986,6 +986,8 @@ def test_output_current_state(mocker, dummy_carbon_data, time_index):
             "soil_c_pool_lmwc",
             "soil_c_pool_microbe",
             "soil_c_pool_pom",
+            "soil_enzyme_pom",
+            "soil_enzyme_maom",
         ],
         time_index,
     )
diff --git a/tests/models/litter/test_litter_model.py b/tests/models/litter/test_litter_model.py
index b9e296a0a..81bf0d958 100644
--- a/tests/models/litter/test_litter_model.py
+++ b/tests/models/litter/test_litter_model.py
@@ -33,8 +33,8 @@ def litter_model_fixture(dummy_litter_data):
 
 def test_litter_model_initialization(caplog, dummy_litter_data):
     """Test `LitterModel` initialization."""
-
     from virtual_rainforest.core.base_model import BaseModel
+    from virtual_rainforest.core.constants import CoreConsts
     from virtual_rainforest.models.litter.constants import LitterConsts
     from virtual_rainforest.models.litter.litter_model import LitterModel
 
@@ -43,7 +43,8 @@ def test_litter_model_initialization(caplog, dummy_litter_data):
         update_interval=pint.Quantity("1 week"),
         soil_layers=[-0.5, -1.0],
         canopy_layers=10,
-        constants=LitterConsts,
+        model_constants=LitterConsts,
+        core_constants=CoreConsts,
     )
 
     # In cases where it passes then checks that the object has the right properties
@@ -76,7 +77,7 @@ def test_litter_model_initialization(caplog, dummy_litter_data):
 
 def test_litter_model_initialization_no_data(caplog):
     """Test `LitterModel` initialization fails when all data is missing."""
-
+    from virtual_rainforest.core.constants import CoreConsts
     from virtual_rainforest.core.data import Data
     from virtual_rainforest.core.grid import Grid
     from virtual_rainforest.models.litter.constants import LitterConsts
@@ -94,7 +95,8 @@ def test_litter_model_initialization_no_data(caplog):
             update_interval=pint.Quantity("1 week"),
             soil_layers=2,  # FIXME - incorrect soil layer spec in model
             canopy_layers=10,
-            constants=LitterConsts,
+            model_constants=LitterConsts,
+            core_constants=CoreConsts,
         )
 
     # Final check that expected logging entries are produced
@@ -149,7 +151,7 @@ def test_litter_model_initialization_no_data(caplog):
 
 def test_litter_model_initialization_bad_pool_bounds(caplog, dummy_litter_data):
     """Test `LitterModel` initialization fails when litter pools are out of bounds."""
-
+    from virtual_rainforest.core.constants import CoreConsts
     from virtual_rainforest.models.litter.constants import LitterConsts
     from virtual_rainforest.models.litter.litter_model import LitterModel
 
@@ -164,7 +166,8 @@ def test_litter_model_initialization_bad_pool_bounds(caplog, dummy_litter_data):
             update_interval=pint.Quantity("1 week"),
             soil_layers=2,
             canopy_layers=10,
-            constants=LitterConsts,
+            model_constants=LitterConsts,
+            core_constants=CoreConsts,
         )
 
     # Final check that the last log entry is as expected
@@ -177,7 +180,7 @@ def test_litter_model_initialization_bad_pool_bounds(caplog, dummy_litter_data):
 
 def test_litter_model_initialization_bad_lignin_bounds(caplog, dummy_litter_data):
     """Test `LitterModel` initialization fails for lignin proportions not in bounds."""
-
+    from virtual_rainforest.core.constants import CoreConsts
     from virtual_rainforest.models.litter.constants import LitterConsts
     from virtual_rainforest.models.litter.litter_model import LitterModel
 
@@ -187,7 +190,14 @@ def test_litter_model_initialization_bad_lignin_bounds(caplog, dummy_litter_data
         # Put incorrect data in for woody lignin
         litter_data["lignin_woody"] = DataArray([0.5, 0.4, 1.1], dims=["cell_id"])
 
-        LitterModel(litter_data, pint.Quantity("1 week"), 2, 10, constants=LitterConsts)
+        LitterModel(
+            litter_data,
+            pint.Quantity("1 week"),
+            2,
+            10,
+            model_constants=LitterConsts,
+            core_constants=CoreConsts,
+        )
 
     # Final check that expected logging entries are produced
     log_check(
@@ -207,6 +217,7 @@ def test_litter_model_initialization_bad_lignin_bounds(caplog, dummy_litter_data
             does_not_raise(),
             (
                 (INFO, "Initialised litter.LitterConsts from config"),
+                (INFO, "Initialised core.CoreConsts from config"),
                 (
                     INFO,
                     "Information required to initialise the litter model successfully "
@@ -255,6 +266,7 @@ def test_litter_model_initialization_bad_lignin_bounds(caplog, dummy_litter_data
             does_not_raise(),
             (
                 (INFO, "Initialised litter.LitterConsts from config"),
+                (INFO, "Initialised core.CoreConsts from config"),
                 (
                     INFO,
                     "Information required to initialise the litter model successfully "
@@ -327,7 +339,7 @@ def test_generate_litter_model(
             pint.Quantity(config["core"]["timing"]["update_interval"]),
         )
         assert model.update_interval == time_interval
-        assert model.constants.litter_decomp_temp_response == temp_response
+        assert model.model_constants.litter_decomp_temp_response == temp_response
 
     # Final check that expected logging entries are produced
     log_check(caplog, expected_log_entries)
diff --git a/tests/models/litter/test_litter_pools.py b/tests/models/litter/test_litter_pools.py
index 447cc86de..7f1927400 100644
--- a/tests/models/litter/test_litter_pools.py
+++ b/tests/models/litter/test_litter_pools.py
@@ -131,6 +131,7 @@ def test_calculate_change_in_litter_variables(
     dummy_litter_data, surface_layer_index, top_soil_layer_index
 ):
     """Test that litter pool update calculation is correct."""
+    from virtual_rainforest.core.constants import CoreConsts
     from virtual_rainforest.models.litter.litter_pools import (
         calculate_change_in_litter_variables,
     )
@@ -168,7 +169,8 @@ def test_calculate_change_in_litter_variables(
         decomposed_excrement=dummy_litter_data["decomposed_excrement"].to_numpy(),
         decomposed_carcasses=dummy_litter_data["decomposed_carcasses"].to_numpy(),
         update_interval=1.0,
-        constants=LitterConsts,
+        model_constants=LitterConsts,
+        core_constants=CoreConsts,
     )
 
     for name in expected_pools.keys():
@@ -206,6 +208,7 @@ def test_calculate_decay_rates(dummy_litter_data, temp_and_water_factors):
 
 def test_calculate_total_C_mineralised(decay_rates):
     """Test that calculation of total C mineralised is as expected."""
+    from virtual_rainforest.core.constants import CoreConsts
     from virtual_rainforest.models.litter.litter_pools import (
         calculate_total_C_mineralised,
     )
@@ -213,8 +216,7 @@ def test_calculate_total_C_mineralised(decay_rates):
     expected_mineralisation = [0.0212182, 0.0274272, 0.00617274]
 
     actual_mineralisation = calculate_total_C_mineralised(
-        decay_rates=decay_rates,
-        constants=LitterConsts,
+        decay_rates=decay_rates, model_constants=LitterConsts, core_constants=CoreConsts
     )
 
     assert np.allclose(actual_mineralisation, expected_mineralisation)
diff --git a/tests/models/soil/conftest.py b/tests/models/soil/conftest.py
new file mode 100644
index 000000000..d98147bc8
--- /dev/null
+++ b/tests/models/soil/conftest.py
@@ -0,0 +1,51 @@
+"""Collection of fixtures to assist the testing of the soil model."""
+
+import pytest
+
+from virtual_rainforest.models.soil.constants import SoilConsts
+
+
+@pytest.fixture
+def environmental_factors(dummy_carbon_data, top_soil_layer_index):
+    """Environmental factors based on dummy carbon data."""
+    from virtual_rainforest.models.soil.env_factors import (
+        calculate_clay_impact_on_enzyme_saturation,
+        calculate_clay_impact_on_necromass_decay,
+        calculate_pH_suitability,
+        calculate_water_potential_impact_on_microbes,
+    )
+
+    water_factors = calculate_water_potential_impact_on_microbes(
+        water_potential=dummy_carbon_data["matric_potential"][
+            top_soil_layer_index
+        ].to_numpy(),
+        water_potential_halt=SoilConsts.soil_microbe_water_potential_halt,
+        water_potential_opt=SoilConsts.soil_microbe_water_potential_optimum,
+        response_curvature=SoilConsts.microbial_water_response_curvature,
+    )
+
+    pH_factors = calculate_pH_suitability(
+        soil_pH=dummy_carbon_data["pH"].to_numpy(),
+        maximum_pH=SoilConsts.max_pH_microbes,
+        minimum_pH=SoilConsts.min_pH_microbes,
+        lower_optimum_pH=SoilConsts.lowest_optimal_pH_microbes,
+        upper_optimum_pH=SoilConsts.highest_optimal_pH_microbes,
+    )
+
+    clay_saturation_factors = calculate_clay_impact_on_enzyme_saturation(
+        clay_fraction=dummy_carbon_data["clay_fraction"].to_numpy(),
+        base_protection=SoilConsts.base_soil_protection,
+        protection_with_clay=SoilConsts.soil_protection_with_clay,
+    )
+
+    clay_decay_factors = calculate_clay_impact_on_necromass_decay(
+        clay_fraction=dummy_carbon_data["clay_fraction"].to_numpy(),
+        decay_exponent=SoilConsts.clay_necromass_decay_exponent,
+    )
+
+    return {
+        "water": water_factors,
+        "pH": pH_factors,
+        "clay_saturation": clay_saturation_factors,
+        "clay_decay": clay_decay_factors,
+    }
diff --git a/tests/models/soil/test_carbon.py b/tests/models/soil/test_carbon.py
index 8c7da068a..af5772f57 100644
--- a/tests/models/soil/test_carbon.py
+++ b/tests/models/soil/test_carbon.py
@@ -3,72 +3,24 @@
 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 tests.conftest import log_check
 from virtual_rainforest.models.soil.constants import SoilConsts
 
 
-@pytest.fixture
-def moist_scalars(dummy_carbon_data, top_soil_layer_index):
-    """Moisture scalars based on dummy carbon data."""
-    from virtual_rainforest.models.soil.carbon import convert_moisture_to_scalar
-
-    moist_scalars = convert_moisture_to_scalar(
-        np.array([0.5, 0.7, 0.6, 0.2]),
-        SoilConsts.moisture_scalar_coefficient,
-        SoilConsts.moisture_scalar_exponent,
-    )
-
-    return moist_scalars
-
-
-@pytest.fixture
-def water_factors(dummy_carbon_data, top_soil_layer_index):
-    """Water potential factors based on dummy carbon data."""
-    from virtual_rainforest.models.soil.carbon import (
-        calculate_water_potential_impact_on_microbes,
-    )
-
-    moist_scalars = calculate_water_potential_impact_on_microbes(
-        water_potential=dummy_carbon_data["matric_potential"][top_soil_layer_index],
-        water_potential_halt=SoilConsts.soil_microbe_water_potential_halt,
-        water_potential_opt=SoilConsts.soil_microbe_water_potential_optimum,
-        moisture_response_curvature=SoilConsts.moisture_response_curvature,
-    )
-
-    return moist_scalars
-
-
-def test_top_soil_data_extraction(dummy_carbon_data, top_soil_layer_index):
-    """Test that top soil data can be extracted from the data object correctly."""
-
-    top_soil_temps = [35.0, 37.5, 40.0, 25.0]
-    top_soil_water_potentials = [-3.0, -10.0, -250.0, -10000.0]
-
-    assert np.allclose(
-        dummy_carbon_data["soil_temperature"][top_soil_layer_index], top_soil_temps
-    )
-    assert np.allclose(
-        dummy_carbon_data["matric_potential"][top_soil_layer_index],
-        top_soil_water_potentials,
-    )
-
-
 def test_calculate_soil_carbon_updates(dummy_carbon_data, top_soil_layer_index):
     """Test that the two pool update functions work correctly."""
-
+    from virtual_rainforest.core.constants import CoreConsts
     from virtual_rainforest.models.soil.carbon import calculate_soil_carbon_updates
 
     change_in_pools = {
-        "soil_c_pool_lmwc": [0.25847352, 0.59685734, 0.97122869, 0.02112077],
-        "soil_c_pool_maom": [0.10296645, 0.04445693, -0.31401747, 0.00422143],
-        "soil_c_pool_microbe": [-0.16002809, -0.07621711, -0.36452355, -0.01140071],
-        "soil_c_pool_pom": [-0.25377786, -0.58704913, -0.41245236, -0.01566563],
+        "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],
     }
 
     # Make order of pools object
@@ -81,17 +33,24 @@ def test_calculate_soil_carbon_updates(dummy_carbon_data, top_soil_layer_index):
         soil_c_pool_maom=dummy_carbon_data["soil_c_pool_maom"].to_numpy(),
         soil_c_pool_microbe=dummy_carbon_data["soil_c_pool_microbe"].to_numpy(),
         soil_c_pool_pom=dummy_carbon_data["soil_c_pool_pom"].to_numpy(),
+        soil_enzyme_pom=dummy_carbon_data["soil_enzyme_pom"].to_numpy(),
+        soil_enzyme_maom=dummy_carbon_data["soil_enzyme_maom"].to_numpy(),
         pH=dummy_carbon_data["pH"],
         bulk_density=dummy_carbon_data["bulk_density"],
-        soil_moisture=np.array([0.5, 0.7, 0.6, 0.2]),
+        soil_moisture=dummy_carbon_data["soil_moisture"][
+            top_soil_layer_index
+        ].to_numpy(),
         soil_water_potential=dummy_carbon_data["matric_potential"][
             top_soil_layer_index
         ].to_numpy(),
+        # TODO - Change this once average vertical flow is used
+        vertical_flow_rate=dummy_carbon_data["vertical_flow"] / 30.0,
         soil_temp=dummy_carbon_data["soil_temperature"][top_soil_layer_index],
-        percent_clay=dummy_carbon_data["percent_clay"],
+        clay_fraction=dummy_carbon_data["clay_fraction"],
         mineralisation_rate=dummy_carbon_data["litter_C_mineralisation_rate"],
         delta_pools_ordered=pool_order,
-        constants=SoilConsts,
+        model_constants=SoilConsts,
+        core_constants=CoreConsts,
     )
 
     # Check that the updates are correctly calculated. Using a loop here implicitly
@@ -100,241 +59,60 @@ def test_calculate_soil_carbon_updates(dummy_carbon_data, top_soil_layer_index):
         assert np.allclose(delta_pools[i * 4 : (i + 1) * 4], change_in_pools[pool])
 
 
-def test_calculate_mineral_association(dummy_carbon_data, moist_scalars):
-    """Test that mineral_association runs and generates the correct values."""
-
-    from virtual_rainforest.models.soil.carbon import calculate_mineral_association
-
-    output_l_to_m = [-5.78872135e-03, -1.00408341e-02, -5.62807109e-01, 2.60743689e-05]
-
-    # Then calculate mineral association rate
-    lmwc_to_maom = calculate_mineral_association(
-        dummy_carbon_data["soil_c_pool_lmwc"],
-        dummy_carbon_data["soil_c_pool_maom"],
-        dummy_carbon_data["pH"],
-        dummy_carbon_data["bulk_density"],
-        moist_scalars,
-        dummy_carbon_data["percent_clay"],
-        constants=SoilConsts,
-    )
-
-    # Check that expected values are generated
-    assert np.allclose(lmwc_to_maom, output_l_to_m)
-
-
-def test_calculate_equilibrium_maom(dummy_carbon_data):
-    """Test that equilibrium maom calculation works as expected."""
-    from virtual_rainforest.models.soil.carbon import calculate_equilibrium_maom
+def test_determine_microbial_biomass_losses(
+    dummy_carbon_data, top_soil_layer_index, environmental_factors
+):
+    """Check that the determination of microbial biomass losses works correctly."""
+    from virtual_rainforest.models.soil.carbon import determine_microbial_biomass_losses
 
-    Q_max = [2.38520786, 1.98025934, 0.64714262, 2.80537157]
-    output_eqb_maoms = [2.13182275, 0.65105909, 0.36433141, 0.58717765]
+    expected_maintenance = [0.05443078, 0.02298407, 0.12012258, 0.00722288]
+    expected_pom_enzyme = [0.0005443078, 0.0002298407, 0.0012012258, 7.22288e-5]
+    expected_maom_enzyme = [0.0005443078, 0.0002298407, 0.0012012258, 7.22288e-5]
+    expected_decay_to_pom = [0.04631043, 0.01809481, 0.09055279, 0.00621707]
+    expected_decay_to_lmwc = [0.007031729, 0.004429577, 0.027167343, 8.613595e-4]
 
-    equib_maoms = calculate_equilibrium_maom(
-        dummy_carbon_data["pH"],
-        Q_max,
-        dummy_carbon_data["soil_c_pool_lmwc"],
+    losses = determine_microbial_biomass_losses(
+        soil_c_pool_microbe=dummy_carbon_data["soil_c_pool_microbe"],
+        soil_temp=dummy_carbon_data["soil_temperature"][top_soil_layer_index],
+        clay_factor_decay=environmental_factors["clay_decay"],
         constants=SoilConsts,
     )
-    assert np.allclose(equib_maoms, output_eqb_maoms)
 
+    # Check that each rate matches expectation
+    assert np.allclose(losses.maintenance_synthesis, expected_maintenance)
+    assert np.allclose(losses.pom_enzyme_production, expected_pom_enzyme)
+    assert np.allclose(losses.maom_enzyme_production, expected_maom_enzyme)
+    assert np.allclose(losses.necromass_decay_to_lmwc, expected_decay_to_lmwc)
+    assert np.allclose(losses.necromass_decay_to_pom, expected_decay_to_pom)
 
-@pytest.mark.parametrize(
-    "alternative,output_capacities,raises,expected_log_entries",
-    [
-        (
-            None,
-            [2.38520786, 1.98025934, 0.64714262, 2.80537157],
-            does_not_raise(),
-            (),
-        ),
-        (
-            [156.0],
-            [],
-            pytest.raises(ValueError),
-            ((ERROR, "Relative clay content must be expressed as a percentage!"),),
-        ),
-        (
-            [-9.0],
-            [],
-            pytest.raises(ValueError),
-            ((ERROR, "Relative clay content must be expressed as a percentage!"),),
-        ),
-    ],
-)
-def test_calculate_max_sorption_capacity(
-    caplog,
-    dummy_carbon_data,
-    alternative,
-    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:
-        if alternative:
-            max_capacities = calculate_max_sorption_capacity(
-                dummy_carbon_data["bulk_density"],
-                np.array(alternative, dtype=np.float32),
-                SoilConsts.max_sorption_with_clay_slope,
-                SoilConsts.max_sorption_with_clay_intercept,
-            )
-        else:
-            max_capacities = calculate_max_sorption_capacity(
-                dummy_carbon_data["bulk_density"],
-                dummy_carbon_data["percent_clay"],
-                SoilConsts.max_sorption_with_clay_slope,
-                SoilConsts.max_sorption_with_clay_intercept,
-            )
-
-        assert np.allclose(max_capacities, output_capacities)
-
-    log_check(caplog, expected_log_entries)
-
-
-def test_calculate_binding_coefficient(dummy_carbon_data):
-    """Test that Langmuir binding coefficient calculation works as expected."""
-    from virtual_rainforest.models.soil.carbon import calculate_binding_coefficient
-
-    output_coefs = [168.26740611, 24.49063242, 12.88249552, 52.9419581]
-
-    binding_coefs = calculate_binding_coefficient(
-        dummy_carbon_data["pH"],
-        SoilConsts.binding_with_ph_slope,
-        SoilConsts.binding_with_ph_intercept,
+    # Then check that sum of other rates is the same as the overall
+    # maintenance_synthesis rate
+    assert np.allclose(
+        losses.maintenance_synthesis,
+        losses.pom_enzyme_production
+        + losses.maom_enzyme_production
+        + losses.necromass_decay_to_lmwc
+        + losses.necromass_decay_to_pom,
     )
 
-    assert np.allclose(binding_coefs, output_coefs)
 
-
-@pytest.mark.parametrize(
-    "activation_energy,expected_factors",
-    [
-        (30000.0, [2.57153601, 2.82565326, 3.10021393, 1.73629781]),
-        (45000.0, [4.12371761, 4.74983258, 5.45867825, 2.28789625]),
-        (57000.0, [6.01680536, 7.19657491, 8.58309980, 2.85289648]),
-    ],
-)
-def calculate_temperature_effect_on_microbes(
-    dummy_carbon_data, top_soil_layer_index, activation_energy, expected_factors
+def test_calculate_maintenance_biomass_synthesis(
+    dummy_carbon_data, top_soil_layer_index
 ):
-    """Test function to calculate microbial temperature response."""
-    from virtual_rainforest.models.soil.carbon import (
-        calculate_temperature_effect_on_microbes,
-    )
-
-    actual_factors = calculate_temperature_effect_on_microbes(
-        soil_temperature=dummy_carbon_data["soil_temperature"][top_soil_layer_index],
-        activation_energy=activation_energy,
-        reference_temperature=SoilConsts.arrhenius_reference_temp,
-        gas_constant=SoilConsts.universal_gas_constant,
-    )
-
-    assert np.allclose(expected_factors, actual_factors)
-
-
-@pytest.mark.parametrize(
-    "alternative,output_scalars,raises,expected_log_entries",
-    [
-        (None, [0.750035, 0.947787, 0.880671, 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,
-    dummy_carbon_data,
-    alternative,
-    output_scalars,
-    raises,
-    expected_log_entries,
-    top_soil_layer_index,
-):
-    """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:
-        if alternative:
-            moist_scalar = convert_moisture_to_scalar(
-                np.array(alternative, dtype=np.float32),
-                SoilConsts.moisture_scalar_coefficient,
-                SoilConsts.moisture_scalar_exponent,
-            )
-        else:
-            moist_scalar = convert_moisture_to_scalar(
-                np.array([0.5, 0.7, 0.6, 0.2]),
-                SoilConsts.moisture_scalar_coefficient,
-                SoilConsts.moisture_scalar_exponent,
-            )
-
-        assert np.allclose(moist_scalar, output_scalars)
-
-    log_check(caplog, expected_log_entries)
-
-
-def test_calculate_water_potential_impact_on_microbes():
-    """Test the calculation of the impact of soil water on microbial rates."""
+    """Check maintenance respiration cost calculates correctly."""
     from virtual_rainforest.models.soil.carbon import (
-        calculate_water_potential_impact_on_microbes,
-    )
-
-    water_potentials = np.array([-3.0, -10.0, -250.0, -10000.0])
-
-    expected_factor = [1.0, 0.94414168, 0.62176357, 0.07747536]
-
-    actual_factor = calculate_water_potential_impact_on_microbes(
-        water_potential=water_potentials,
-        water_potential_halt=SoilConsts.soil_microbe_water_potential_halt,
-        water_potential_opt=SoilConsts.soil_microbe_water_potential_optimum,
-        moisture_response_curvature=SoilConsts.moisture_response_curvature,
+        calculate_maintenance_biomass_synthesis,
     )
 
-    assert np.allclose(actual_factor, expected_factor)
-
-
-def test_calculate_maintenance_respiration(
-    dummy_carbon_data, moist_scalars, top_soil_layer_index
-):
-    """Check maintenance respiration cost calculates correctly."""
-    from virtual_rainforest.models.soil.carbon import calculate_maintenance_respiration
-
-    expected_resps = [0.05443078, 0.02298407, 0.12012258, 0.00722288]
+    expected_loss = [0.05443078, 0.02298407, 0.12012258, 0.00722288]
 
-    main_resps = calculate_maintenance_respiration(
+    actual_loss = calculate_maintenance_biomass_synthesis(
         soil_c_pool_microbe=dummy_carbon_data["soil_c_pool_microbe"],
         soil_temp=dummy_carbon_data["soil_temperature"][top_soil_layer_index],
         constants=SoilConsts,
     )
 
-    assert np.allclose(main_resps, expected_resps)
-
-
-def test_calculate_necromass_adsorption(dummy_carbon_data, moist_scalars):
-    """Check maintenance respiration cost calculates correctly."""
-    from virtual_rainforest.models.soil.carbon import calculate_necromass_adsorption
-
-    expected_adsorps = [0.10875517, 0.05449776, 0.24878964, 0.00419536]
-
-    actual_adsorps = calculate_necromass_adsorption(
-        dummy_carbon_data["soil_c_pool_microbe"],
-        moist_scalars,
-        SoilConsts.necromass_adsorption_rate,
-    )
-
-    assert np.allclose(actual_adsorps, expected_adsorps)
+    assert np.allclose(actual_loss, expected_loss)
 
 
 def test_calculate_carbon_use_efficiency(dummy_carbon_data, top_soil_layer_index):
@@ -353,99 +131,80 @@ def test_calculate_carbon_use_efficiency(dummy_carbon_data, top_soil_layer_index
     assert np.allclose(actual_cues, expected_cues)
 
 
-def test_calculate_microbial_saturation(dummy_carbon_data):
-    """Check microbial activity saturation calculates correctly."""
-    from virtual_rainforest.models.soil.carbon import calculate_microbial_saturation
-
-    expected_saturated = [0.99876016, 0.99687933, 0.99936324, 0.99285147]
-
-    actual_saturated = calculate_microbial_saturation(
-        dummy_carbon_data["soil_c_pool_microbe"],
-        SoilConsts.half_sat_microbial_activity,
-    )
-
-    assert np.allclose(actual_saturated, expected_saturated)
-
-
-def test_calculate_microbial_pom_mineralisation_saturation(dummy_carbon_data):
-    """Check microbial mineralisation saturation calculates correctly."""
-    from virtual_rainforest.models.soil.carbon import (
-        calculate_microbial_pom_mineralisation_saturation,
-    )
-
-    expected_saturated = [0.99793530, 0.99480968, 0.99893917, 0.98814229]
-
-    actual_saturated = calculate_microbial_pom_mineralisation_saturation(
-        dummy_carbon_data["soil_c_pool_microbe"],
-        SoilConsts.half_sat_microbial_pom_mineralisation,
-    )
-
-    assert np.allclose(actual_saturated, expected_saturated)
-
-
-def test_calculate_pom_decomposition_saturation(dummy_carbon_data):
-    """Check POM decomposition saturation calculates correctly."""
-    from virtual_rainforest.models.soil.carbon import (
-        calculate_pom_decomposition_saturation,
-    )
-
-    expected_saturated = [0.4, 0.86956521, 0.82352941, 0.7]
+@pytest.mark.parametrize(
+    "turnover,expected_decay",
+    [
+        (
+            2.4e-2,
+            [0.000544296, 0.000229824, 0.001201224, 7.224e-5],
+        ),
+        (
+            6.5e-2,
+            [0.001474135, 0.00062244, 0.003253315, 0.00019565],
+        ),
+        (
+            2.4e-3,
+            [5.44296e-5, 2.29824e-5, 0.0001201224, 7.224e-6],
+        ),
+    ],
+)
+def test_calculate_enzyme_turnover(dummy_carbon_data, turnover, expected_decay):
+    """Check that enzyme turnover rates are calculated correctly."""
+    from virtual_rainforest.models.soil.carbon import calculate_enzyme_turnover
 
-    actual_saturated = calculate_pom_decomposition_saturation(
-        dummy_carbon_data["soil_c_pool_pom"], SoilConsts.half_sat_pom_decomposition
+    actual_decay = calculate_enzyme_turnover(
+        enzyme_pool=dummy_carbon_data["soil_enzyme_pom"], turnover_rate=turnover
     )
 
-    assert np.allclose(actual_saturated, expected_saturated)
+    assert np.allclose(actual_decay, expected_decay)
 
 
 def test_calculate_microbial_carbon_uptake(
-    dummy_carbon_data, top_soil_layer_index, water_factors
+    dummy_carbon_data, top_soil_layer_index, environmental_factors
 ):
     """Check microbial carbon uptake calculates correctly."""
     from virtual_rainforest.models.soil.carbon import calculate_microbial_carbon_uptake
 
-    expected_uptake = [3.15786124e-3, 1.26472838e-3, 4.38868085e-3, 1.75270792e-5]
+    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 = calculate_microbial_carbon_uptake(
+    actual_uptake, actual_assimilation = calculate_microbial_carbon_uptake(
         soil_c_pool_lmwc=dummy_carbon_data["soil_c_pool_lmwc"],
         soil_c_pool_microbe=dummy_carbon_data["soil_c_pool_microbe"],
-        water_factor=water_factors,
-        soil_temp=dummy_carbon_data["soil_temperature"][top_soil_layer_index],
+        water_factor=environmental_factors["water"],
+        pH_factor=environmental_factors["pH"],
+        soil_temp=dummy_carbon_data["soil_temperature"][
+            top_soil_layer_index
+        ].to_numpy(),
         constants=SoilConsts,
     )
 
     assert np.allclose(actual_uptake, expected_uptake)
+    assert np.allclose(actual_assimilation, expected_assimilation)
 
 
-def test_calculate_labile_carbon_leaching(dummy_carbon_data, moist_scalars):
-    """Check leaching of labile carbon is calculated correctly."""
-    from virtual_rainforest.models.soil.carbon import calculate_labile_carbon_leaching
-
-    expected_leaching = [5.62526764e-05, 2.84336164e-05, 1.32100695e-04, 1.25860748e-06]
-
-    actual_leaching = calculate_labile_carbon_leaching(
-        dummy_carbon_data["soil_c_pool_lmwc"],
-        moist_scalars,
-        SoilConsts.leaching_rate_labile_carbon,
-    )
-
-    assert np.allclose(actual_leaching, expected_leaching)
-
-
-def test_calculate_pom_decomposition(
-    dummy_carbon_data, top_soil_layer_index, water_factors
+def test_calculate_enzyme_mediated_decomposition(
+    dummy_carbon_data, top_soil_layer_index, environmental_factors
 ):
     """Check that particulate organic matter decomposition is calculated correctly."""
-    from virtual_rainforest.models.soil.carbon import calculate_pom_decomposition
+    from virtual_rainforest.models.soil.carbon import (
+        calculate_enzyme_mediated_decomposition,
+    )
 
-    expected_decomp = [0.25589892, 0.58810966, 0.41294236, 0.02116563]
+    expected_decomp = [3.39844565e-4, 8.91990315e-3, 1.25055119e-2, 4.14247999e-5]
 
-    actual_decomp = calculate_pom_decomposition(
-        dummy_carbon_data["soil_c_pool_pom"],
-        dummy_carbon_data["soil_c_pool_microbe"],
-        water_factor=water_factors,
+    actual_decomp = calculate_enzyme_mediated_decomposition(
+        soil_c_pool=dummy_carbon_data["soil_c_pool_pom"],
+        soil_enzyme=dummy_carbon_data["soil_enzyme_pom"],
+        water_factor=environmental_factors["water"],
+        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)
diff --git a/tests/models/soil/test_env_factors.py b/tests/models/soil/test_env_factors.py
new file mode 100644
index 000000000..87e7a70b4
--- /dev/null
+++ b/tests/models/soil/test_env_factors.py
@@ -0,0 +1,233 @@
+"""Test module for soil.env_factors.py.
+
+This module tests the functions which calculate environmental impacts on soil processes.
+"""
+
+import numpy as np
+import pytest
+
+
+def test_top_soil_data_extraction(dummy_carbon_data, top_soil_layer_index):
+    """Test that top soil data can be extracted from the data object correctly."""
+
+    top_soil_temps = [35.0, 37.5, 40.0, 25.0]
+    top_soil_water_potentials = [-3.0, -10.0, -250.0, -10000.0]
+
+    assert np.allclose(
+        dummy_carbon_data["soil_temperature"][top_soil_layer_index], top_soil_temps
+    )
+    assert np.allclose(
+        dummy_carbon_data["matric_potential"][top_soil_layer_index],
+        top_soil_water_potentials,
+    )
+
+
+def test_calculate_environmental_effect_factors(
+    dummy_carbon_data, top_soil_layer_index
+):
+    """Test that function to calculate all set of environmental factors works."""
+    from virtual_rainforest.models.soil.constants import SoilConsts
+    from virtual_rainforest.models.soil.env_factors import (
+        calculate_environmental_effect_factors,
+    )
+
+    expected_water = [1.0, 0.94414168, 0.62176357, 0.07747536]
+    expected_pH = [0.25, 1.0, 0.428571428, 1.0]
+    expected_clay_sat = [1.782, 1.102, 0.83, 1.918]
+    expected_clay_decay = [0.52729242, 0.78662786, 0.92311634, 0.48675225]
+
+    env_factors = calculate_environmental_effect_factors(
+        soil_water_potential=dummy_carbon_data["matric_potential"][
+            top_soil_layer_index
+        ],
+        pH=dummy_carbon_data["pH"],
+        clay_fraction=dummy_carbon_data["clay_fraction"],
+        constants=SoilConsts,
+    )
+
+    assert np.allclose(env_factors.water, expected_water)
+    assert np.allclose(env_factors.pH, expected_pH)
+    assert np.allclose(env_factors.clay_saturation, expected_clay_sat)
+    assert np.allclose(env_factors.clay_decay, expected_clay_decay)
+
+
+@pytest.mark.parametrize(
+    "activation_energy,expected_factors",
+    [
+        (30000.0, [2.57153601, 2.82565326, 3.10021393, 1.73629781]),
+        (45000.0, [4.12371761, 4.74983258, 5.45867825, 2.28789625]),
+        (57000.0, [6.01680536, 7.19657491, 8.58309980, 2.85289648]),
+    ],
+)
+def calculate_temperature_effect_on_microbes(
+    dummy_carbon_data, top_soil_layer_index, activation_energy, expected_factors
+):
+    """Test function to calculate microbial temperature response."""
+    from virtual_rainforest.models.soil.constants import SoilConsts
+    from virtual_rainforest.models.soil.env_factors import (
+        calculate_temperature_effect_on_microbes,
+    )
+
+    actual_factors = calculate_temperature_effect_on_microbes(
+        soil_temperature=dummy_carbon_data["soil_temperature"][top_soil_layer_index],
+        activation_energy=activation_energy,
+        reference_temperature=SoilConsts.arrhenius_reference_temp,
+        gas_constant=SoilConsts.universal_gas_constant,
+    )
+
+    assert np.allclose(expected_factors, actual_factors)
+
+
+def test_calculate_water_potential_impact_on_microbes(
+    dummy_carbon_data, top_soil_layer_index
+):
+    """Test the calculation of the impact of soil water on microbial rates."""
+    from virtual_rainforest.models.soil.constants import SoilConsts
+    from virtual_rainforest.models.soil.env_factors import (
+        calculate_water_potential_impact_on_microbes,
+    )
+
+    expected_factor = [1.0, 0.94414168, 0.62176357, 0.07747536]
+
+    actual_factor = calculate_water_potential_impact_on_microbes(
+        water_potential=dummy_carbon_data["matric_potential"][top_soil_layer_index],
+        water_potential_halt=SoilConsts.soil_microbe_water_potential_halt,
+        water_potential_opt=SoilConsts.soil_microbe_water_potential_optimum,
+        response_curvature=SoilConsts.microbial_water_response_curvature,
+    )
+
+    assert np.allclose(actual_factor, expected_factor)
+
+
+def test_soil_water_potential_too_high(dummy_carbon_data, top_soil_layer_index):
+    """Test that too high soil water potential results in an error."""
+    from virtual_rainforest.models.soil.constants import SoilConsts
+    from virtual_rainforest.models.soil.env_factors import (
+        calculate_water_potential_impact_on_microbes,
+    )
+
+    water_potentials = np.array([-2.0, -10.0, -250.0, -10000.0])
+
+    with pytest.raises(ValueError):
+        calculate_water_potential_impact_on_microbes(
+            water_potential=water_potentials,
+            water_potential_halt=SoilConsts.soil_microbe_water_potential_halt,
+            water_potential_opt=SoilConsts.soil_microbe_water_potential_optimum,
+            response_curvature=SoilConsts.microbial_water_response_curvature,
+        )
+
+
+def test_calculate_pH_suitability():
+    """Test that calculation of pH suitability is correct."""
+    from virtual_rainforest.models.soil.constants import SoilConsts
+    from virtual_rainforest.models.soil.env_factors import calculate_pH_suitability
+
+    pH_values = np.array([3.0, 7.5, 9.0, 5.7, 2.0, 11.5])
+    expected_inhib = [0.25, 1.0, 0.428571428, 1.0, 0.0, 0.0]
+
+    actual_inhib = calculate_pH_suitability(
+        soil_pH=pH_values,
+        maximum_pH=SoilConsts.max_pH_microbes,
+        minimum_pH=SoilConsts.min_pH_microbes,
+        lower_optimum_pH=SoilConsts.lowest_optimal_pH_microbes,
+        upper_optimum_pH=SoilConsts.highest_optimal_pH_microbes,
+    )
+
+    assert np.allclose(expected_inhib, actual_inhib)
+
+
+@pytest.mark.parametrize(
+    argnames=["params"],
+    argvalues=[
+        pytest.param(
+            {
+                "maximum_pH": 7.0,
+                "minimum_pH": 2.5,
+                "lower_optimum_pH": 4.5,
+                "upper_optimum_pH": 7.5,
+            },
+            id="maximum_pH too low",
+        ),
+        pytest.param(
+            {
+                "maximum_pH": 11.0,
+                "minimum_pH": 2.5,
+                "lower_optimum_pH": 1.5,
+                "upper_optimum_pH": 7.5,
+            },
+            id="lower_optimum_pH too low",
+        ),
+        pytest.param(
+            {
+                "maximum_pH": 11.0,
+                "minimum_pH": 2.5,
+                "lower_optimum_pH": 4.5,
+                "upper_optimum_pH": 3.5,
+            },
+            id="upper_optimum_pH too low",
+        ),
+    ],
+)
+def test_calculate_pH_suitability_errors(params):
+    """Test that calculation of pH suitability generates errors if constants are bad."""
+    from virtual_rainforest.models.soil.env_factors import calculate_pH_suitability
+
+    pH_values = np.array([3.0, 7.5, 9.0, 5.7, 2.0, 11.5])
+
+    with pytest.raises(ValueError):
+        calculate_pH_suitability(soil_pH=pH_values, **params)
+
+
+def test_calculate_clay_impact_on_enzyme_saturation(dummy_carbon_data):
+    """Test calculation of the effect of soil clay fraction on saturation constants."""
+    from virtual_rainforest.models.soil.constants import SoilConsts
+    from virtual_rainforest.models.soil.env_factors import (
+        calculate_clay_impact_on_enzyme_saturation,
+    )
+
+    expected_factor = [1.782, 1.102, 0.83, 1.918]
+
+    actual_factor = calculate_clay_impact_on_enzyme_saturation(
+        clay_fraction=dummy_carbon_data["clay_fraction"],
+        base_protection=SoilConsts.base_soil_protection,
+        protection_with_clay=SoilConsts.soil_protection_with_clay,
+    )
+
+    assert np.allclose(expected_factor, actual_factor)
+
+
+def test_calculate_clay_impact_on_necromass_decay(dummy_carbon_data):
+    """Test calculation of the effect of soil clay fraction on necromass decay."""
+    from virtual_rainforest.models.soil.constants import SoilConsts
+    from virtual_rainforest.models.soil.env_factors import (
+        calculate_clay_impact_on_necromass_decay,
+    )
+
+    expected_factor = [0.52729242, 0.78662786, 0.92311634, 0.48675225]
+
+    actual_factor = calculate_clay_impact_on_necromass_decay(
+        clay_fraction=dummy_carbon_data["clay_fraction"],
+        decay_exponent=SoilConsts.clay_necromass_decay_exponent,
+    )
+
+    assert np.allclose(expected_factor, actual_factor)
+
+
+def test_calculate_leaching_rate(dummy_carbon_data, top_soil_layer_index):
+    """Test calculation of solute leaching rates."""
+    from virtual_rainforest.core.constants import CoreConsts
+    from virtual_rainforest.models.soil.constants import SoilConsts
+    from virtual_rainforest.models.soil.env_factors import calculate_leaching_rate
+
+    expected_rate = [1.07473723e-6, 2.53952130e-6, 9.91551977e-5, 5.25567712e-5]
+    vertical_flow_per_day = np.array([0.1, 0.5, 2.5, 15.9])
+
+    actual_rate = calculate_leaching_rate(
+        solute_density=dummy_carbon_data["soil_c_pool_lmwc"],
+        vertical_flow_rate=vertical_flow_per_day,
+        soil_moisture=dummy_carbon_data["soil_moisture"][top_soil_layer_index],
+        solubility_coefficient=SoilConsts.solubility_coefficient_lmwc,
+        soil_layer_thickness=CoreConsts.depth_of_active_soil_layer,
+    )
+
+    assert np.allclose(expected_rate, actual_rate)
diff --git a/tests/models/soil/test_soil_model.py b/tests/models/soil/test_soil_model.py
index 794d87bbe..c69613990 100644
--- a/tests/models/soil/test_soil_model.py
+++ b/tests/models/soil/test_soil_model.py
@@ -17,13 +17,9 @@
 @pytest.fixture
 def soil_model_fixture(dummy_carbon_data):
     """Create a soil model fixture based on the dummy carbon data."""
-
     from virtual_rainforest.core.config import Config
-    from virtual_rainforest.core.registry import register_module
     from virtual_rainforest.models.soil.soil_model import SoilModel
 
-    # Register the module components to access constants classes
-    register_module("virtual_rainforest.models.abiotic_simple")
     # Build the config object
     config = Config(
         cfg_strings="[core]\n[core.timing]\nupdate_interval = '12 hours'\n[soil]\n"
@@ -36,8 +32,8 @@ def soil_model_fixture(dummy_carbon_data):
 
 def test_soil_model_initialization(caplog, dummy_carbon_data):
     """Test `SoilModel` initialization with good data."""
-
     from virtual_rainforest.core.base_model import BaseModel
+    from virtual_rainforest.core.constants import CoreConsts
     from virtual_rainforest.models.soil.constants import SoilConsts
     from virtual_rainforest.models.soil.soil_model import SoilModel
 
@@ -46,7 +42,8 @@ def test_soil_model_initialization(caplog, dummy_carbon_data):
         pint.Quantity("1 week"),
         [-0.5, -1.0],
         10,
-        constants=SoilConsts,
+        model_constants=SoilConsts,
+        core_constants=CoreConsts,
     )
 
     # In cases where it passes then checks that the object has the right properties
@@ -66,7 +63,7 @@ def test_soil_model_initialization(caplog, dummy_carbon_data):
             (DEBUG, "soil model: required var 'soil_c_pool_pom' checked"),
             (DEBUG, "soil model: required var 'pH' checked"),
             (DEBUG, "soil model: required var 'bulk_density' checked"),
-            (DEBUG, "soil model: required var 'percent_clay' checked"),
+            (DEBUG, "soil model: required var 'clay_fraction' checked"),
         ),
     )
 
@@ -74,6 +71,7 @@ def test_soil_model_initialization(caplog, dummy_carbon_data):
 def test_soil_model_initialization_no_data(caplog, dummy_carbon_data):
     """Test `SoilModel` initialization with no data."""
 
+    from virtual_rainforest.core.constants import CoreConsts
     from virtual_rainforest.core.data import Data
     from virtual_rainforest.core.grid import Grid
     from virtual_rainforest.models.soil.constants import SoilConsts
@@ -90,7 +88,8 @@ def test_soil_model_initialization_no_data(caplog, dummy_carbon_data):
             pint.Quantity("1 week"),
             [-0.5, -1.0],
             10,
-            constants=SoilConsts,
+            model_constants=SoilConsts,
+            core_constants=CoreConsts,
         )
 
     # Final check that expected logging entries are produced
@@ -103,7 +102,7 @@ def test_soil_model_initialization_no_data(caplog, dummy_carbon_data):
             (ERROR, "soil model: init data missing required var 'soil_c_pool_pom'"),
             (ERROR, "soil model: init data missing required var 'pH'"),
             (ERROR, "soil model: init data missing required var 'bulk_density'"),
-            (ERROR, "soil model: init data missing required var 'percent_clay'"),
+            (ERROR, "soil model: init data missing required var 'clay_fraction'"),
             (ERROR, "soil model: error checking required_init_vars, see log."),
         ),
     )
@@ -112,6 +111,7 @@ def test_soil_model_initialization_no_data(caplog, dummy_carbon_data):
 def test_soil_model_initialization_bounds_error(caplog, dummy_carbon_data):
     """Test `SoilModel` initialization."""
 
+    from virtual_rainforest.core.constants import CoreConsts
     from virtual_rainforest.models.soil.constants import SoilConsts
     from virtual_rainforest.models.soil.soil_model import SoilModel
 
@@ -127,7 +127,8 @@ def test_soil_model_initialization_bounds_error(caplog, dummy_carbon_data):
             pint.Quantity("1 week"),
             [-0.5, -1.0],
             10,
-            constants=SoilConsts,
+            model_constants=SoilConsts,
+            core_constants=CoreConsts,
         )
 
     # Final check that expected logging entries are produced
@@ -141,7 +142,7 @@ def test_soil_model_initialization_bounds_error(caplog, dummy_carbon_data):
             (DEBUG, "soil model: required var 'soil_c_pool_pom' checked"),
             (DEBUG, "soil model: required var 'pH' checked"),
             (DEBUG, "soil model: required var 'bulk_density' checked"),
-            (DEBUG, "soil model: required var 'percent_clay' checked"),
+            (DEBUG, "soil model: required var 'clay_fraction' checked"),
             (ERROR, "Initial carbon pools contain at least one negative value!"),
         ),
     )
@@ -153,10 +154,11 @@ def test_soil_model_initialization_bounds_error(caplog, dummy_carbon_data):
         pytest.param(
             "[core]\n[core.timing]\nupdate_interval = '12 hours'\n[soil]\n",
             pint.Quantity("12 hours"),
-            0.2,
+            60.0,
             does_not_raise(),
             (
                 (INFO, "Initialised soil.SoilConsts from config"),
+                (INFO, "Initialised core.CoreConsts from config"),
                 (
                     INFO,
                     "Information required to initialise the soil model successfully "
@@ -168,7 +170,7 @@ def test_soil_model_initialization_bounds_error(caplog, dummy_carbon_data):
                 (DEBUG, "soil model: required var 'soil_c_pool_pom' checked"),
                 (DEBUG, "soil model: required var 'pH' checked"),
                 (DEBUG, "soil model: required var 'bulk_density' checked"),
-                (DEBUG, "soil model: required var 'percent_clay' checked"),
+                (DEBUG, "soil model: required var 'clay_fraction' checked"),
             ),
             id="default_config",
         ),
@@ -180,6 +182,7 @@ def test_soil_model_initialization_bounds_error(caplog, dummy_carbon_data):
             does_not_raise(),
             (
                 (INFO, "Initialised soil.SoilConsts from config"),
+                (INFO, "Initialised core.CoreConsts from config"),
                 (
                     INFO,
                     "Information required to initialise the soil model successfully "
@@ -191,7 +194,7 @@ def test_soil_model_initialization_bounds_error(caplog, dummy_carbon_data):
                 (DEBUG, "soil model: required var 'soil_c_pool_pom' checked"),
                 (DEBUG, "soil model: required var 'pH' checked"),
                 (DEBUG, "soil model: required var 'bulk_density' checked"),
-                (DEBUG, "soil model: required var 'percent_clay' checked"),
+                (DEBUG, "soil model: required var 'clay_fraction' checked"),
             ),
             id="modified_config_correct",
         ),
@@ -240,7 +243,7 @@ def test_generate_soil_model(
             pint.Quantity(config["core"]["timing"]["update_interval"]),
         )
         assert model.update_interval == time_interval
-        assert model.constants.max_decomp_rate_pom == max_decomp
+        assert model.model_constants.max_decomp_rate_pom == max_decomp
 
     # Final check that expected logging entries are produced
     log_check(caplog, expected_log_entries)
@@ -287,17 +290,22 @@ def test_update(mocker, soil_model_fixture, dummy_carbon_data):
             Dataset(
                 data_vars=dict(
                     lmwc=DataArray(
-                        [0.13122264, 0.30432201, 0.50601794, 0.01541629], dims="cell_id"
+                        [0.04823845, 0.02119714, 0.0895863, 0.00528887], dims="cell_id"
                     ),
                     maom=DataArray(
-                        [2.54647825, 1.71347325, 4.32317632, 0.50157362], dims="cell_id"
+                        [2.49936689, 1.70118553, 4.50085129, 0.50000614], dims="cell_id"
                     ),
                     microbe=DataArray(
-                        [5.72558633, 2.27805994, 11.21111949, 0.99495352],
-                        dims="cell_id",
+                        [5.77512315, 2.2899636, 11.24827514, 0.99640928], dims="cell_id"
                     ),
                     pom=DataArray(
-                        [0.02068157, 0.71317448, 0.50004912, 0.34220329], dims="cell_id"
+                        [0.12397575, 1.00508662, 0.7389913, 0.35583206], dims="cell_id"
+                    ),
+                    enzyme_pom=DataArray(
+                        [0.02267842, 0.00957576, 0.05004963, 0.00300993], dims="cell_id"
+                    ),
+                    enzyme_maom=DataArray(
+                        [0.0354453, 0.01167442, 0.02538637, 0.00454144], dims="cell_id"
                     ),
                 )
             ),
@@ -337,6 +345,8 @@ def test_integrate_soil_model(
         assert np.allclose(new_pools["soil_c_pool_maom"], final_pools["maom"])
         assert np.allclose(new_pools["soil_c_pool_microbe"], final_pools["microbe"])
         assert np.allclose(new_pools["soil_c_pool_pom"], final_pools["pom"])
+        assert np.allclose(new_pools["soil_enzyme_pom"], final_pools["enzyme_pom"])
+        assert np.allclose(new_pools["soil_enzyme_maom"], final_pools["enzyme_maom"])
 
     # Check that integrator is called once (and once only)
     if mock_output:
@@ -369,8 +379,9 @@ def test_order_independance(dummy_carbon_data, soil_model_fixture):
         "bulk_density",
         "soil_moisture",
         "matric_potential",
+        "vertical_flow",
         "soil_temperature",
-        "percent_clay",
+        "clay_fraction",
         "litter_C_mineralisation_rate",
     ]
     for not_pool in not_pools:
@@ -380,7 +391,7 @@ def test_order_independance(dummy_carbon_data, soil_model_fixture):
     pool_names = [
         str(name)
         for name in dummy_carbon_data.data.keys()
-        if str(name).startswith("soil_c_pool_")
+        if str(name).startswith("soil_c_pool_") or str(name).startswith("soil_enzyme_")
     ]
 
     # Add pool values from object in reversed order
@@ -404,29 +415,36 @@ def test_order_independance(dummy_carbon_data, soil_model_fixture):
 
 def test_construct_full_soil_model(dummy_carbon_data, top_soil_layer_index):
     """Test that the function that creates the object to integrate exists and works."""
+    from virtual_rainforest.core.constants import CoreConsts
     from virtual_rainforest.models.soil.constants import SoilConsts
     from virtual_rainforest.models.soil.soil_model import construct_full_soil_model
 
-    delta_pools = (
-        [
-            0.25809707,
-            0.59264789,
-            0.56851013,
-            0.02112901,
-            0.0962045,
-            0.02576618,
-            -0.08926844,
-            0.0029497,
-            -0.15288599,
-            -0.05330495,
-            -0.18645947,
-            -0.01013683,
-            -0.25377786,
-            -0.58704913,
-            -0.41245236,
-            -0.01566563,
-        ],
-    )
+    delta_pools = [
+        -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,
+        1.16756409e-2,
+        1.17571917e-8,
+        1.67442231e-8,
+        1.83311362e-9,
+        -1.11675865e-08,
+        -0.00031009,
+        -5.09593e-5,
+        0.0005990658,
+        -3.72112e-5,
+    ]
 
     # make pools
     pools = np.concatenate(
@@ -434,6 +452,7 @@ def test_construct_full_soil_model(dummy_carbon_data, top_soil_layer_index):
             dummy_carbon_data[str(name)].to_numpy()
             for name in dummy_carbon_data.data.keys()
             if str(name).startswith("soil_c_pool_")
+            or str(name).startswith("soil_enzyme_")
         ]
     )
 
@@ -441,17 +460,22 @@ def test_construct_full_soil_model(dummy_carbon_data, top_soil_layer_index):
     delta_pools_ordered = {
         str(name): np.array([])
         for name in dummy_carbon_data.data.keys()
-        if str(name).startswith("soil_c_pool_")
+        if str(name).startswith("soil_c_pool_") or str(name).startswith("soil_enzyme_")
     }
 
+    # Find rate of flow per day
+    vertical_flow_per_day = dummy_carbon_data["vertical_flow"].to_numpy() / 30.0
+
     rate_of_change = construct_full_soil_model(
         0.0,
-        pools,
-        dummy_carbon_data,
-        4,
-        top_soil_layer_index,
+        pools=pools,
+        data=dummy_carbon_data,
+        vertical_flow_per_day=vertical_flow_per_day,
+        no_cells=4,
+        top_soil_layer_index=top_soil_layer_index,
         delta_pools_ordered=delta_pools_ordered,
-        constants=SoilConsts,
+        model_constants=SoilConsts,
+        core_constants=CoreConsts,
     )
 
     assert np.allclose(delta_pools, rate_of_change)
diff --git a/tests/test_main.py b/tests/test_main.py
index 8721f0dee..016b90778 100644
--- a/tests/test_main.py
+++ b/tests/test_main.py
@@ -28,6 +28,7 @@
             (
                 (INFO, "Initialising models: soil"),
                 (INFO, "Initialised soil.SoilConsts from config"),
+                (INFO, "Initialised core.CoreConsts from config"),
                 (
                     INFO,
                     "Information required to initialise the soil model successfully "
@@ -39,7 +40,7 @@
                 (DEBUG, "soil model: required var 'soil_c_pool_pom' checked"),
                 (DEBUG, "soil model: required var 'pH' checked"),
                 (DEBUG, "soil model: required var 'bulk_density' checked"),
-                (DEBUG, "soil model: required var 'percent_clay' checked"),
+                (DEBUG, "soil model: required var 'clay_fraction' checked"),
             ),
             id="valid config",
         ),
@@ -51,6 +52,7 @@
             (
                 (INFO, "Initialising models: soil"),
                 (INFO, "Initialised soil.SoilConsts from config"),
+                (INFO, "Initialised core.CoreConsts from config"),
                 (
                     INFO,
                     "Information required to initialise the soil model successfully "
@@ -69,6 +71,7 @@
             (
                 (INFO, "Initialising models: soil"),
                 (INFO, "Initialised soil.SoilConsts from config"),
+                (INFO, "Initialised core.CoreConsts from config"),
                 (
                     INFO,
                     "Information required to initialise the soil model successfully "
diff --git a/virtual_rainforest/core/constants.py b/virtual_rainforest/core/constants.py
index eb8b550c9..06b6019db 100644
--- a/virtual_rainforest/core/constants.py
+++ b/virtual_rainforest/core/constants.py
@@ -21,5 +21,11 @@ class CoreConsts(ConstantsDataclass):
     zero_Celsius: ClassVar[float] = constants.zero_Celsius
     """Conversion constant from Kelvin to Celsius (°)."""
 
-    placeholder: float = 123.4
-    """A placeholder configurable constant."""
+    depth_of_active_soil_layer: float = 0.25
+    """Depth of the biogeochemically active soil layer [m].
+
+    The soil model considered a homogenous layer in which all significant nutrient
+    processes take place. This is a major assumption of the model. The value is taken
+    from :cite:t:`fatichi_mechanistic_2019`. No empirical source is provided for this
+    value.
+    """
diff --git a/virtual_rainforest/core/module_schema.json b/virtual_rainforest/core/module_schema.json
index 34aa7e8af..557c0cef5 100644
--- a/virtual_rainforest/core/module_schema.json
+++ b/virtual_rainforest/core/module_schema.json
@@ -196,7 +196,7 @@
                      "minItems": 1,
                      "uniqueItems": true,
                      "default": [
-                        -0.5,
+                        -0.25,
                         -1.0
                      ]
                   },
diff --git a/virtual_rainforest/data_variables.toml b/virtual_rainforest/data_variables.toml
index 8d348543e..75e66f4f4 100644
--- a/virtual_rainforest/data_variables.toml
+++ b/virtual_rainforest/data_variables.toml
@@ -477,14 +477,6 @@ initialised_by = "external"
 updated_by = "soil"
 used_by = [ "abiotic", "soil"]
 
-[[variable]]
-name = "percent_clay"
-description = "Clay as a percentage of total soil"
-unit = "%"
-initialised_by = "external"
-updated_by = ""
-used_by = [ "abiotic", "soil"]
-
 [[variable]]
 name = "pH"
 description = "Soil pH"
diff --git a/virtual_rainforest/example_data/dummy_soil_data.nc b/virtual_rainforest/example_data/dummy_soil_data.nc
index 6e8f8a04d..5a839a471 100644
Binary files a/virtual_rainforest/example_data/dummy_soil_data.nc and b/virtual_rainforest/example_data/dummy_soil_data.nc differ
diff --git a/virtual_rainforest/example_data/dummy_total_config.toml b/virtual_rainforest/example_data/dummy_total_config.toml
index 0dddf7cf2..1a8b2828f 100644
--- a/virtual_rainforest/example_data/dummy_total_config.toml
+++ b/virtual_rainforest/example_data/dummy_total_config.toml
@@ -48,7 +48,7 @@ file = "dummy_soil_data.nc"
 var_name = "bulk_density"
 [[core.data.variable]]
 file = "dummy_soil_data.nc"
-var_name = "percent_clay"
+var_name = "clay_fraction"
 [[core.data.variable]]
 file = "dummy_soil_data.nc"
 var_name = "soil_c_pool_lmwc"
@@ -61,6 +61,12 @@ var_name = "soil_c_pool_microbe"
 [[core.data.variable]]
 file = "dummy_soil_data.nc"
 var_name = "soil_c_pool_pom"
+[[core.data.variable]]
+file = "dummy_soil_data.nc"
+var_name = "soil_enzyme_pom"
+[[core.data.variable]]
+file = "dummy_soil_data.nc"
+var_name = "soil_enzyme_maom"
 
 [[core.data.variable]]
 file = "dummy_litter_data.nc"
diff --git a/virtual_rainforest/example_data/soil_dummy.py b/virtual_rainforest/example_data/soil_dummy.py
index 6ba0d2aed..dd5b1d292 100644
--- a/virtual_rainforest/example_data/soil_dummy.py
+++ b/virtual_rainforest/example_data/soil_dummy.py
@@ -6,87 +6,42 @@
 function sensibly for.
 """
 
+import numpy as np
 from xarray import Dataset
 
+# How far the center of each cell is from the origin. This applies to both the x and y
+# direction independently, so cell (0,0) is at the origin, whereas cell (2,3) is 180m
+# from the origin in the x direction and 270m in the y direction.
+cell_spacing = np.arange(0, 721, 90)
+gradient = np.multiply.outer(cell_spacing / 90, cell_spacing / 90)
 
-def generate_pH_values(x: float, y: float) -> float:
-    """Function to generate a reasonable range of pH values.
-
-    We're looking at acidic soils so a range of 3.5-4.5 seems plausible.
-    """
-    return 3.5 + (x * y) / (64)
-
-
-def generate_BD_values(x: float, y: float) -> float:
-    """Function to generate a reasonable range of bulk density values.
-
-    Bulk density can vary quite a lot so a range of 1200-1800 kg m^-3 seems sensible.
-    """
-    return 1200.0 + 600.0 * (x * y) / (64)
-
-
-def generate_clay_values(x: float, y: float) -> float:
-    """Function to generate a reasonable range of clay content values.
-
-    We're considering fairly clayey soils, so look at a range of 27.0-40.0 % clay.
-    """
-    return 27.0 + 13.0 * (x * y) / (64)
-
-
-def generate_lmwc_values(x: float, y: float) -> float:
-    """Function to generate a reasonable range of lmwc values.
-
-    LMWC generally a very small carbon pool, so a range of 0.005-0.01 kg C m^-3 is used.
-    """
-    return 0.005 + 0.005 * (x * y) / (64)
-
-
-def generate_maom_values(x: float, y: float) -> float:
-    """Function to generate a reasonable range of maom values.
-
-    A huge amount of carbon can be locked away as MAOM, so a range of 1.0-3.0 kg C m^-3
-    is used.
-    """
-    return 1.0 + 2.0 * (x * y) / (64)
-
+# Generate a range of plausible values (3.5-4.5) for the soil pH [unitless].
+pH_values = 3.5 + 1.00 * (gradient) / (64)
 
-def generate_microbial_C_values(x: float, y: float) -> float:
-    """Function to generate a reasonable range of microbial C values.
+# Generate a range of plausible values (1200-1800) for the bulk density [kg m^-3].
+bulk_density_values = 1200.0 + 600.0 * (gradient) / (64)
 
-    The carbon locked up as microbial biomass is tiny, so a range of 0.0015-0.005
-    kg C m^-3 is used.
-    """
-    return 0.0015 + 0.0035 * (x * y) / (64)
+# Generate a range of plausible values (0.27-0.40) for the clay fraction [fraction].
+clay_fraction_values = 0.27 + 0.13 * (gradient) / (64)
 
+# Generate a range of plausible values (0.005-0.01) for the lmwc pool [kg C m^-3].
+lmwc_values = 0.005 + 0.005 * (gradient) / (64)
 
-def generate_pom_values(x: float, y: float) -> float:
-    """Function to generate a reasonable range of pom values.
+# Generate a range of plausible values (1.0-3.0) for the maom pool [kg C m^-3].
+maom_values = 1.0 + 2.0 * (gradient) / (64)
 
-    A reasonable amount of carbon is stored as particulate organic matter (POM), so a
-    range of 0.1-1.0 kg C m^-3 is used.
-    """
-    return 0.1 + 0.9 * (x * y) / (64)
+# Generate a range of plausible values (0.0015-0.005) for the microbial C pool
+# [kg C m^-3].
+microbial_C_values = 0.0015 + 0.0035 * (gradient) / (64)
 
+# Generate a range of plausible values (0.1-1.0) for the POM pool [kg C m^-3].
+pom_values = 0.1 + 0.9 * (gradient) / (64)
 
-# Generate range of cell numbers in the a x and y directions. Here we have a 9x9 grid,
-# so cells are numbered from 0 to 8 in each direction.
-x_cell_ids = range(0, 9)
-y_cell_ids = range(0, 9)
+# Generate a range of plausible values (0.01-0.5) for the POM enzyme pool [kg C m^-3].
+pom_enzyme_values = 0.01 + 0.49 * (gradient) / (64)
 
-# Make matrices containing all the relevant values
-pH_values = [[generate_pH_values(x, y) for y in y_cell_ids] for x in x_cell_ids]
-bulk_density_values = [
-    [generate_BD_values(x, y) for y in y_cell_ids] for x in x_cell_ids
-]
-percent_clay_values = [
-    [generate_clay_values(x, y) for y in y_cell_ids] for x in x_cell_ids
-]
-lmwc_values = [[generate_lmwc_values(x, y) for y in y_cell_ids] for x in x_cell_ids]
-maom_values = [[generate_maom_values(x, y) for y in y_cell_ids] for x in x_cell_ids]
-microbial_C_values = [
-    [generate_microbial_C_values(x, y) for y in y_cell_ids] for x in x_cell_ids
-]
-pom_values = [[generate_pom_values(x, y) for y in y_cell_ids] for x in x_cell_ids]
+# Generate a range of plausible values (0.01-0.5) for the MAOM enzyme pool [kg C m^-3].
+maom_enzyme_values = 0.01 + 0.49 * (gradient) / (64)
 
 # How far the center of each cell is from the origin. This applies to both the x and y
 # direction independently, so cell (0,0) is at the origin, whereas cell (2,3) is 180m
@@ -98,11 +53,13 @@ def generate_pom_values(x: float, y: float) -> float:
     data_vars=dict(
         pH=(["x", "y"], pH_values),
         bulk_density=(["x", "y"], bulk_density_values),
-        percent_clay=(["x", "y"], percent_clay_values),
+        clay_fraction=(["x", "y"], clay_fraction_values),
         soil_c_pool_lmwc=(["x", "y"], lmwc_values),
         soil_c_pool_maom=(["x", "y"], maom_values),
         soil_c_pool_microbe=(["x", "y"], microbial_C_values),
         soil_c_pool_pom=(["x", "y"], pom_values),
+        soil_enzyme_pom=(["x", "y"], pom_enzyme_values),
+        soil_enzyme_maom=(["x", "y"], maom_enzyme_values),
     ),
     coords=dict(
         x=(["x"], cell_displacements),
diff --git a/virtual_rainforest/models/litter/constants.py b/virtual_rainforest/models/litter/constants.py
index cb7e3bdc4..b49743c81 100644
--- a/virtual_rainforest/models/litter/constants.py
+++ b/virtual_rainforest/models/litter/constants.py
@@ -151,19 +151,6 @@ class LitterConsts(ConstantsDataclass):
     documentation for :attr:`cue_metabolic` for details.
     """
 
-    depth_of_active_layer: float = 0.25
-    """Depth of the biogeochemically active soil layer [m].
-
-    The soil model considered a homogenous layer in which all significant nutrient
-    processes take place. This is a major assumption of the model. The value is taken
-    from :cite:t:`fatichi_mechanistic_2019`. No empirical source is provided for this
-    value.
-
-    This is really a core constant as it is shared across models. However, the core
-    constants are not yet setup, so this constant is being stored here for the time
-    being.
-    """
-
     lignin_inhibition_factor: float = -5.0
     """Exponential factor expressing the extent that lignin inhibits litter breakdown.
 
diff --git a/virtual_rainforest/models/litter/litter_model.py b/virtual_rainforest/models/litter/litter_model.py
index 626b84cb2..2e4583121 100644
--- a/virtual_rainforest/models/litter/litter_model.py
+++ b/virtual_rainforest/models/litter/litter_model.py
@@ -41,6 +41,7 @@ class instance. If errors crop here when converting the information from the con
 
 from virtual_rainforest.core.base_model import BaseModel
 from virtual_rainforest.core.config import Config
+from virtual_rainforest.core.constants import CoreConsts
 from virtual_rainforest.core.constants_loader import load_constants
 from virtual_rainforest.core.data import Data
 from virtual_rainforest.core.exceptions import InitialisationError
@@ -106,7 +107,8 @@ def __init__(
         update_interval: Quantity,
         soil_layers: list[float],
         canopy_layers: int,
-        constants: LitterConsts,
+        model_constants: LitterConsts,
+        core_constants: CoreConsts,
         **kwargs: Any,
     ):
         super().__init__(data, update_interval, **kwargs)
@@ -150,8 +152,10 @@ def __init__(
             LOGGER.error(to_raise)
             raise to_raise
 
-        self.constants = constants
-        """Set of constants for the litter model"""
+        self.model_constants = model_constants
+        """Set of constants for the litter model."""
+        self.core_constants = core_constants
+        """Set of constants shared across models."""
 
         # create a list of layer roles
         layer_roles = set_layer_roles(canopy_layers, soil_layers)
@@ -187,13 +191,21 @@ def from_config(
         canopy_layers = config["core"]["layers"]["canopy_layers"]
 
         # Load in the relevant constants
-        constants = load_constants(config, "litter", "LitterConsts")
+        model_constants = load_constants(config, "litter", "LitterConsts")
+        core_constants = load_constants(config, "core", "CoreConsts")
 
         LOGGER.info(
             "Information required to initialise the litter model successfully "
             "extracted."
         )
-        return cls(data, update_interval, soil_layers, canopy_layers, constants)
+        return cls(
+            data=data,
+            update_interval=update_interval,
+            soil_layers=soil_layers,
+            canopy_layers=canopy_layers,
+            model_constants=model_constants,
+            core_constants=core_constants,
+        )
 
     def setup(self) -> None:
         """Placeholder function to setup up the litter model."""
@@ -219,7 +231,8 @@ def update(self, time_index: int, **kwargs: Any) -> None:
             water_potential=self.data["matric_potential"][
                 self.top_soil_layer_index
             ].to_numpy(),
-            constants=self.constants,
+            model_constants=self.model_constants,
+            core_constants=self.core_constants,
             update_interval=self.update_interval.to("day").magnitude,
             above_metabolic=self.data["litter_pool_above_metabolic"].to_numpy(),
             above_structural=self.data["litter_pool_above_structural"].to_numpy(),
diff --git a/virtual_rainforest/models/litter/litter_pools.py b/virtual_rainforest/models/litter/litter_pools.py
index 130389431..54f283fb6 100644
--- a/virtual_rainforest/models/litter/litter_pools.py
+++ b/virtual_rainforest/models/litter/litter_pools.py
@@ -21,6 +21,7 @@
 import numpy as np
 from numpy.typing import NDArray
 
+from virtual_rainforest.core.constants import CoreConsts
 from virtual_rainforest.models.litter.constants import LitterConsts
 
 
@@ -39,7 +40,8 @@ def calculate_change_in_litter_variables(
     decomposed_excrement: NDArray[np.float32],
     decomposed_carcasses: NDArray[np.float32],
     update_interval: float,
-    constants: LitterConsts,
+    model_constants: LitterConsts,
+    core_constants: CoreConsts,
 ) -> dict[str, NDArray[np.float32]]:
     """Calculate changes for all the litter variables (pool sizes and chemistries).
 
@@ -64,7 +66,8 @@ def calculate_change_in_litter_variables(
         decomposed_carcasses: Input rate of (partially) decomposed carcass biomass from
             the animal model [kg C m^-2 day^-1]
         update_interval: Interval that the litter pools are being updated for [days]
-        constants: Set of constants for the litter model
+        model_constants: Set of constants for the litter model
+        core_constants: Set of core constants shared between all models
 
     Returns:
         The new value for each of the litter pools, and the total mineralisation rate.
@@ -76,7 +79,7 @@ def calculate_change_in_litter_variables(
         surface_temp=surface_temp,
         topsoil_temp=topsoil_temp,
         water_potential=water_potential,
-        constants=constants,
+        constants=model_constants,
     )
 
     # Calculate the litter pool decay rates
@@ -90,11 +93,13 @@ def calculate_change_in_litter_variables(
         lignin_woody=lignin_woody,
         lignin_below_structural=lignin_below_structural,
         environmental_factors=environmental_factors,
-        constants=constants,
+        constants=model_constants,
     )
 
     # Calculate the total mineralisation of carbon from the litter
-    total_C_mineralisation_rate = calculate_total_C_mineralised(decay_rates, constants)
+    total_C_mineralisation_rate = calculate_total_C_mineralised(
+        decay_rates, model_constants=model_constants, core_constants=core_constants
+    )
 
     # Calculate the updated pool masses
     updated_pools = calculate_updated_pools(
@@ -107,7 +112,7 @@ def calculate_change_in_litter_variables(
         decomposed_carcasses=decomposed_carcasses,
         decay_rates=decay_rates,
         update_interval=update_interval,
-        constants=constants,
+        constants=model_constants,
     )
 
     # Find the changes in the lignin concentrations of the 3 relevant pools
@@ -117,7 +122,7 @@ def calculate_change_in_litter_variables(
         lignin_below_structural=lignin_below_structural,
         updated_pools=updated_pools,
         update_interval=update_interval,
-        constants=constants,
+        constants=model_constants,
     )
 
     # Construct dictionary of data arrays to return
@@ -217,14 +222,17 @@ def calculate_decay_rates(
 
 
 def calculate_total_C_mineralised(
-    decay_rates: dict[str, NDArray[np.float32]], constants: LitterConsts
+    decay_rates: dict[str, NDArray[np.float32]],
+    model_constants: LitterConsts,
+    core_constants: CoreConsts,
 ) -> NDArray[np.float32]:
     """Calculate the total carbon mineralisation rate from all five litter pools.
 
     Args:
         decay_rates: Dictionary containing the rates of decay for all 5 litter pools
             [kg C m^-2 day^-1]
-        constants: Set of constants for the litter model
+        model_constants: Set of constants for the litter model
+        core_constants: Set of core constants shared between all models
 
     Returns:
         Rate of carbon mineralisation from litter into soil [kg C m^-3 day^-1].
@@ -232,22 +240,24 @@ def calculate_total_C_mineralised(
 
     # Calculate mineralisation from each pool
     metabolic_above_mineral = calculate_carbon_mineralised(
-        decay_rates["metabolic_above"], carbon_use_efficiency=constants.cue_metabolic
+        decay_rates["metabolic_above"],
+        carbon_use_efficiency=model_constants.cue_metabolic,
     )
     structural_above_mineral = calculate_carbon_mineralised(
         decay_rates["structural_above"],
-        carbon_use_efficiency=constants.cue_structural_above_ground,
+        carbon_use_efficiency=model_constants.cue_structural_above_ground,
     )
     woody_mineral = calculate_carbon_mineralised(
         decay_rates["woody"],
-        carbon_use_efficiency=constants.cue_woody,
+        carbon_use_efficiency=model_constants.cue_woody,
     )
     metabolic_below_mineral = calculate_carbon_mineralised(
-        decay_rates["metabolic_below"], carbon_use_efficiency=constants.cue_metabolic
+        decay_rates["metabolic_below"],
+        carbon_use_efficiency=model_constants.cue_metabolic,
     )
     structural_below_mineral = calculate_carbon_mineralised(
         decay_rates["structural_below"],
-        carbon_use_efficiency=constants.cue_structural_below_ground,
+        carbon_use_efficiency=model_constants.cue_structural_below_ground,
     )
 
     # Calculate mineralisation rate
@@ -260,7 +270,7 @@ def calculate_total_C_mineralised(
     )
 
     # Convert mineralisation rate into kg m^-3 units (from kg m^-2)
-    return total_C_mineralisation_rate / constants.depth_of_active_layer
+    return total_C_mineralisation_rate / core_constants.depth_of_active_soil_layer
 
 
 def calculate_updated_pools(
diff --git a/virtual_rainforest/models/soil/__init__.py b/virtual_rainforest/models/soil/__init__.py
index 9e5512130..155b4eaad 100644
--- a/virtual_rainforest/models/soil/__init__.py
+++ b/virtual_rainforest/models/soil/__init__.py
@@ -8,6 +8,8 @@
   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.
+* The :mod:`~virtual_rainforest.models.soil.env_factors` provides functions that capture
+  the impact of environmental factors on microbial rates.
 * The :mod:`~virtual_rainforest.models.soil.constants` provides a set of dataclasses
   containing the constants required by the broader soil model.
 """  # noqa: D205, D415
diff --git a/virtual_rainforest/models/soil/carbon.py b/virtual_rainforest/models/soil/carbon.py
index f1ef4d2a4..ea8ec0026 100644
--- a/virtual_rainforest/models/soil/carbon.py
+++ b/virtual_rainforest/models/soil/carbon.py
@@ -1,34 +1,60 @@
 """The ``models.soil.carbon`` module  simulates the soil carbon cycle for the Virtual
-Rainforest. At the moment four pools are modelled, these are low molecular weight carbon
-(LMWC), mineral associated organic matter (MAOM), microbial biomass, and particulate
-organic matter (POM).
+Rainforest. At the moment five pools are modelled, these are low molecular weight carbon
+(LMWC), mineral associated organic matter (MAOM), microbial biomass, particulate organic
+matter (POM), and POM degrading enzymes.
 """  # noqa: D205, D415
 
+from dataclasses import dataclass
+
 import numpy as np
 from numpy.typing import NDArray
-from scipy.constants import convert_temperature, gas_constant
 
-from virtual_rainforest.core.logger import LOGGER
+from virtual_rainforest.core.constants import CoreConsts
 from virtual_rainforest.models.soil.constants import SoilConsts
+from virtual_rainforest.models.soil.env_factors import (
+    calculate_environmental_effect_factors,
+    calculate_leaching_rate,
+    calculate_temperature_effect_on_microbes,
+)
 
-# TODO - Once enzymes are added, temperature dependence of saturation constants should
-# be added.
 
+@dataclass
+class MicrobialBiomassLoss:
+    """A data class to store the various biomass losses from microbial biomass."""
 
+    maintenance_synthesis: NDArray[np.float32]
+    """Rate at which biomass must be synthesised to balance losses [kg C m^-3 day^-1].
+    """
+    pom_enzyme_production: NDArray[np.float32]
+    """Rate at which POM degrading enzymes are produced [kg C m^-3 day^-1]."""
+    maom_enzyme_production: NDArray[np.float32]
+    """Rate at which MAOM degrading enzymes are produced [kg C m^-3 day^-1]."""
+    necromass_decay_to_lmwc: NDArray[np.float32]
+    """Rate at which biomass is lost to the LMWC pool [kg C m^-3 day^-1]."""
+    necromass_decay_to_pom: NDArray[np.float32]
+    """Rate at which biomass is lost to the POM pool [kg C m^-3 day^-1]."""
+
+
+# TODO - This function should probably be shortened, leaving as is for the moment as a
+# sensible split will probably be more obvious once more is added to this function.
 def calculate_soil_carbon_updates(
     soil_c_pool_lmwc: NDArray[np.float32],
     soil_c_pool_maom: NDArray[np.float32],
     soil_c_pool_microbe: NDArray[np.float32],
     soil_c_pool_pom: NDArray[np.float32],
+    soil_enzyme_pom: NDArray[np.float32],
+    soil_enzyme_maom: NDArray[np.float32],
     pH: NDArray[np.float32],
     bulk_density: NDArray[np.float32],
     soil_moisture: NDArray[np.float32],
     soil_water_potential: NDArray[np.float32],
+    vertical_flow_rate: NDArray[np.float32],
     soil_temp: NDArray[np.float32],
-    percent_clay: NDArray[np.float32],
+    clay_fraction: NDArray[np.float32],
     mineralisation_rate: NDArray[np.float32],
     delta_pools_ordered: dict[str, NDArray[np.float32]],
-    constants: SoilConsts,
+    core_constants: CoreConsts,
+    model_constants: SoilConsts,
 ) -> NDArray[np.float32]:
     """Calculate net change for each carbon pool.
 
@@ -41,355 +67,195 @@ def calculate_soil_carbon_updates(
         soil_c_pool_maom: Mineral associated organic matter pool [kg C m^-3]
         soil_c_pool_microbe: Microbial biomass (carbon) pool [kg C m^-3]
         soil_c_pool_pom: Particulate organic matter pool [kg C m^-3]
-        pH: pH values for each soil grid cell
+        soil_enzyme_pom: Amount of enzyme class which breaks down particulate organic
+            matter [kg C m^-3]
+        soil_enzyme_maom: Amount of enzyme class which breaks down mineral associated
+            organic matter [kg C m^-3]
+        pH: pH values for each soil grid cell [unitless]
         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_water_potential: Soil water potential for each grid cell [kPa]
         soil_temp: soil temperature for each soil grid cell [degrees C]
-        percent_clay: Percentage clay for each soil grid cell
+        clay_fraction: The clay fraction for each soil grid cell [unitless]
         mineralisation_rate: Amount of litter mineralised into POM pool [kg C m^-3
             day^-1]
         delta_pools_ordered: Dictionary to store pool changes in the order that pools
             are stored in the initial condition vector.
-        constants: Set of constants for the soil model.
+        core_constants: Set of constants shared between models.
+        model_constants: Set of constants for the soil model.
 
     Returns:
         A vector containing net changes to each pool. Order [lmwc, maom].
     """
 
-    # TODO - At present we have two different factors that capture the impact of soil
-    # water. water_factor is based on soil water potential and applies to the microbial
-    # processes, whereas moist_scalar is based on soil water content and applies to the
-    # chemical processes. In future, all processes will be microbially mediated, meaning
-    # that moist_scalar can be removed.
-
-    # Find the impact of soil water content on chemical soil processes
-    moist_scalar = convert_moisture_to_scalar(
-        soil_moisture,
-        constants.moisture_scalar_coefficient,
-        constants.moisture_scalar_exponent,
-    )
-    # Find the impact of soil water potential on the biochemical soil processes
-    water_factor = calculate_water_potential_impact_on_microbes(
-        water_potential=soil_water_potential,
-        water_potential_halt=constants.soil_microbe_water_potential_halt,
-        water_potential_opt=constants.soil_microbe_water_potential_optimum,
-        moisture_response_curvature=constants.moisture_response_curvature,
-    )
-    # Calculate transfers between pools
-    lmwc_to_maom = calculate_mineral_association(
-        soil_c_pool_lmwc=soil_c_pool_lmwc,
-        soil_c_pool_maom=soil_c_pool_maom,
+    # Find environmental factors which impact biogeochemical soil processes
+    env_factors = calculate_environmental_effect_factors(
+        soil_water_potential=soil_water_potential,
         pH=pH,
-        bulk_density=bulk_density,
-        moisture_scalar=moist_scalar,
-        percent_clay=percent_clay,
-        constants=constants,
+        clay_fraction=clay_fraction,
+        constants=model_constants,
     )
-    microbial_uptake = calculate_microbial_carbon_uptake(
+
+    microbial_uptake, microbial_assimilation = calculate_microbial_carbon_uptake(
         soil_c_pool_lmwc=soil_c_pool_lmwc,
         soil_c_pool_microbe=soil_c_pool_microbe,
-        water_factor=water_factor,
+        water_factor=env_factors.water,
+        pH_factor=env_factors.pH,
         soil_temp=soil_temp,
-        constants=constants,
+        constants=model_constants,
     )
-    microbial_respiration = calculate_maintenance_respiration(
+    biomass_losses = determine_microbial_biomass_losses(
         soil_c_pool_microbe=soil_c_pool_microbe,
         soil_temp=soil_temp,
-        constants=constants,
+        clay_factor_decay=env_factors.clay_decay,
+        constants=model_constants,
     )
-    necromass_adsorption = calculate_necromass_adsorption(
-        soil_c_pool_microbe=soil_c_pool_microbe,
-        moisture_scalar=moist_scalar,
-        necromass_adsorption_rate=constants.necromass_adsorption_rate,
+    pom_enzyme_turnover = calculate_enzyme_turnover(
+        enzyme_pool=soil_enzyme_pom,
+        turnover_rate=model_constants.pom_enzyme_turnover_rate,
     )
-    labile_carbon_leaching = calculate_labile_carbon_leaching(
-        soil_c_pool_lmwc=soil_c_pool_lmwc,
-        moisture_scalar=moist_scalar,
-        leaching_rate=constants.leaching_rate_labile_carbon,
+    maom_enzyme_turnover = calculate_enzyme_turnover(
+        enzyme_pool=soil_enzyme_maom,
+        turnover_rate=model_constants.maom_enzyme_turnover_rate,
     )
-    pom_decomposition_to_lmwc = calculate_pom_decomposition(
-        soil_c_pool_pom=soil_c_pool_pom,
-        soil_c_pool_microbe=soil_c_pool_microbe,
-        water_factor=water_factor,
+    labile_carbon_leaching = calculate_leaching_rate(
+        solute_density=soil_c_pool_lmwc,
+        vertical_flow_rate=vertical_flow_rate,
+        soil_moisture=soil_moisture,
+        solubility_coefficient=model_constants.solubility_coefficient_lmwc,
+        soil_layer_thickness=core_constants.depth_of_active_soil_layer,
+    )
+    pom_decomposition_rate = calculate_enzyme_mediated_decomposition(
+        soil_c_pool=soil_c_pool_pom,
+        soil_enzyme=soil_enzyme_pom,
+        water_factor=env_factors.water,
+        pH_factor=env_factors.pH,
+        clay_factor_saturation=env_factors.clay_saturation,
         soil_temp=soil_temp,
-        constants=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
+    )
+    maom_decomposition_to_lmwc = calculate_enzyme_mediated_decomposition(
+        soil_c_pool=soil_c_pool_maom,
+        soil_enzyme=soil_enzyme_maom,
+        water_factor=env_factors.water,
+        pH_factor=env_factors.pH,
+        clay_factor_saturation=env_factors.clay_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
-        - lmwc_to_maom
+        + 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"] = lmwc_to_maom + necromass_adsorption
+    delta_pools_ordered["soil_c_pool_maom"] = (
+        pom_decomposition_to_maom - maom_decomposition_to_lmwc
+    )
     delta_pools_ordered["soil_c_pool_microbe"] = (
-        microbial_uptake - microbial_respiration - necromass_adsorption
+        microbial_assimilation - biomass_losses.maintenance_synthesis
     )
     delta_pools_ordered["soil_c_pool_pom"] = (
-        mineralisation_rate - pom_decomposition_to_lmwc
+        mineralisation_rate
+        + biomass_losses.necromass_decay_to_pom
+        - pom_decomposition_rate
+    )
+    delta_pools_ordered["soil_enzyme_pom"] = (
+        biomass_losses.pom_enzyme_production - pom_enzyme_turnover
+    )
+    delta_pools_ordered["soil_enzyme_maom"] = (
+        biomass_losses.maom_enzyme_production - maom_enzyme_turnover
     )
 
     # Create output array of pools in desired order
     return np.concatenate(list(delta_pools_ordered.values()))
 
 
-def calculate_mineral_association(
-    soil_c_pool_lmwc: NDArray[np.float32],
-    soil_c_pool_maom: NDArray[np.float32],
-    pH: NDArray[np.float32],
-    bulk_density: NDArray[np.float32],
-    moisture_scalar: NDArray[np.float32],
-    percent_clay: NDArray[np.float32],
-    constants: SoilConsts,
-) -> NDArray[np.float32]:
-    """Calculates net rate of LMWC association with soil minerals.
-
-    Following :cite:t:`abramoff_millennial_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:
-        soil_c_pool_lmwc: Low molecular weight carbon pool [kg C m^-3]
-        soil_c_pool_maom: Mineral associated organic matter pool [kg C m^-3]
-        pH: pH values for each soil grid cell
-        bulk_density: bulk density values for each soil grid cell [kg m^-3]
-        moisture_scalar: A scalar capturing the impact of soil moisture and on process
-            rates [unitless]
-        percent_clay: Percentage clay for each soil grid cell
-        constants: Set of constants for the soil model.
-
-    Returns:
-        The net flux from LMWC to MAOM [kg C m^-3 day^-1]
-    """
-
-    # Calculate maximum sorption
-    Q_max = calculate_max_sorption_capacity(
-        bulk_density,
-        percent_clay,
-        constants.max_sorption_with_clay_slope,
-        constants.max_sorption_with_clay_intercept,
-    )
-    equib_maom = calculate_equilibrium_maom(pH, Q_max, soil_c_pool_lmwc, constants)
-
-    return moisture_scalar * soil_c_pool_lmwc * (equib_maom - soil_c_pool_maom) / Q_max
-
-
-def calculate_max_sorption_capacity(
-    bulk_density: NDArray[np.float32],
-    percent_clay: NDArray[np.float32],
-    max_sorption_with_clay_slope: float,
-    max_sorption_with_clay_intercept: float,
-) -> 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 :cite:t:`mayes_relation_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
-        max_sorption_with_clay_slope: Slope of relationship between clay content and
-            maximum organic matter sorption [(% clay)^-1]
-        max_sorption_with_clay_intercept: Intercept of relationship between clay content
-            and maximum organic matter sorption [log(kg C kg soil ^-1)]
-
-    Returns:
-        Maximum sorption capacity [kg C 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 ** (
-        max_sorption_with_clay_slope * np.log10(percent_clay)
-        + max_sorption_with_clay_intercept
-    )
-    return Q_max
-
-
-def calculate_equilibrium_maom(
-    pH: NDArray[np.float32],
-    Q_max: NDArray[np.float32],
-    lmwc: NDArray[np.float32],
+def determine_microbial_biomass_losses(
+    soil_c_pool_microbe: NDArray[np.float32],
+    soil_temp: NDArray[np.float32],
+    clay_factor_decay: NDArray[np.float32],
     constants: SoilConsts,
-) -> NDArray[np.float32]:
-    """Calculate equilibrium MAOM concentration based on Langmuir coefficients.
+) -> MicrobialBiomassLoss:
+    """Calculate all of the losses from the microbial biomass pool.
 
-    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.
+    Microbes need to synthesis new biomass at a certain rate just to maintain their
+    current biomass. This function calculates this overall rate and the various losses
+    that contribute to this rate. The main sources of this loss are the external
+    excretion of enzymes, cell death, and protein degradation.
 
     Args:
-        pH: pH values for each soil grid cell
-        Q_max: Maximum sorption capacities [kg C m^-3]
-        lmwc: Low molecular weight carbon pool [kg C m^-3]
+        soil_c_pool_microbe: Microbial biomass (carbon) pool [kg C m^-3]
+        soil_temp: soil temperature for each soil grid cell [degrees C]
+        clay_factor_decay: A factor capturing the impact of soil clay fraction on
+            necromass decay destination [unitless]
         constants: Set of constants for the soil model.
 
     Returns:
-        Equilibrium concentration of MAOM [kg C m^-3]
-    """
-
-    binding_coefficient = calculate_binding_coefficient(
-        pH, constants.binding_with_ph_slope, constants.binding_with_ph_intercept
-    )
-    return (binding_coefficient * Q_max * lmwc) / (1 + lmwc * binding_coefficient)
-
-
-def calculate_binding_coefficient(
-    pH: NDArray[np.float32],
-    binding_with_ph_slope: float,
-    binding_with_ph_intercept: float,
-) -> NDArray[np.float32]:
-    """Calculate Langmuir binding coefficient based on pH.
-
-    This specific expression and its parameters are drawn from
-    :cite:t:`mayes_relation_2012`.
-
-    Args:
-        pH: pH values for each soil grid cell
-        binding_with_ph_slope: Slope of relationship between pH and binding coefficient
-            [pH^-1]
-        binding_with_ph_intercept: Intercept of relationship between pH and binding
-            coefficient [log(m^3 kg^-1)]
-
-    Returns:
-        Langmuir binding coefficients for mineral association of labile carbon [m^3
-        kg^-1]
-    """
-
-    return 10.0 ** (binding_with_ph_slope * pH + binding_with_ph_intercept)
-
-
-def calculate_temperature_effect_on_microbes(
-    soil_temperature: NDArray[np.float32],
-    activation_energy: float,
-    reference_temperature: float,
-) -> NDArray[np.float32]:
-    """Calculate the effect that temperature has on microbial metabolic rates.
-
-    This uses a standard Arrhenius equation to calculate the impact of temperature.
-
-    This function takes temperatures in Celsius as inputs and converts them into Kelvin
-    for use in the Arrhenius equation. TODO - review this after we have decided how to
-    handle these conversions in general.
-
-    Args:
-        soil_temperature: The temperature of the soil [C]
-        activation_energy: Energy of activation [J mol^-1]
-        soil_temperature: The reference temperature of the Arrhenius equation [C]
-
-    Returns:
-        A multiplicative factor capturing the effect of temperature on microbial rates
+        A dataclass containing all the losses from the microbial biomass pool.
     """
 
-    # Convert the temperatures to Kelvin
-    soil_temp_in_kelvin = convert_temperature(
-        soil_temperature, old_scale="Celsius", new_scale="Kelvin"
-    )
-    ref_temp_in_kelvin = convert_temperature(
-        reference_temperature, old_scale="Celsius", new_scale="Kelvin"
-    )
-
-    return np.exp(
-        (-activation_energy / gas_constant)
-        * ((1 / (soil_temp_in_kelvin)) - (1 / (ref_temp_in_kelvin)))
+    # Calculate the rate of maintenance synthesis
+    maintenance_synthesis = calculate_maintenance_biomass_synthesis(
+        soil_c_pool_microbe=soil_c_pool_microbe,
+        soil_temp=soil_temp,
+        constants=constants,
     )
 
-
-def calculate_water_potential_impact_on_microbes(
-    water_potential: NDArray[np.float32],
-    water_potential_halt: float,
-    water_potential_opt: float,
-    moisture_response_curvature: float,
-) -> NDArray[np.float32]:
-    """Calculate the effect that soil water potential has on microbial rates.
-
-    This function only returns valid output for soil water potentials that are less than
-    the optimal water potential.
-
-    Args:
-        water_potential: Soil water potential [kPa]
-        water_potential_halt: Water potential at which all microbial activity stops
-            [kPa]
-        water_potential_opt: Optimal water potential for microbial activity [kPa]
-        moisture_response_curvature: Parameter controlling the curvature of the moisture
-            response function [unitless]
-
-    Returns:
-        A multiplicative factor capturing the impact of moisture on soil microbe rates
-        decomposition [unitless]
-    """
-
-    # If the water potential is greater than the optimal then the function produces NaNs
-    # so the simulation should be interrupted
-    if any(water_potential > water_potential_opt):
-        err = ValueError("Water potential greater than minimum value")
-        LOGGER.critical(err)
-        raise err
-
-    # Calculate how much moisture suppresses microbial activity
-    supression = (
-        (np.log10(-water_potential) - np.log10(-water_potential_opt))
-        / (np.log10(-water_potential_halt) - np.log10(-water_potential_opt))
-    ) ** moisture_response_curvature
-
-    return 1 - supression
-
-
-def convert_moisture_to_scalar(
-    soil_moisture: NDArray[np.float32],
-    moisture_scalar_coefficient: float,
-    moisture_scalar_exponent: float,
-) -> NDArray[np.float32]:
-    """Convert soil moisture into a factor to multiply rates by.
-
-    This form is used in :cite:t:`abramoff_millennial_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]
-        moisture_scalar_coefficient: [unit less]
-        moisture_scalar_exponent: [(Relative water content)^-1]
-
-    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
-        + moisture_scalar_coefficient
-        * np.exp(-moisture_scalar_exponent * soil_moisture)
+    # Calculation the production of each enzyme class
+    pom_enzyme_production = constants.maintenance_pom_enzyme * maintenance_synthesis
+    maom_enzyme_production = constants.maintenance_maom_enzyme * maintenance_synthesis
+
+    # Remaining maintenance synthesis is used to replace degraded proteins and cells
+    replacement_synthesis = (
+        1 - constants.maintenance_pom_enzyme - constants.maintenance_maom_enzyme
+    ) * maintenance_synthesis
+
+    # TODO - This split will change when a necromass pool is introduced
+    # Calculate fraction of necromass that decays to LMWC
+    necromass_proportion_to_lmwc = constants.necromass_to_lmwc * clay_factor_decay
+    # These proteins and cells that are replaced decay into either the POM or LMWC pool
+    necromass_to_lmwc = necromass_proportion_to_lmwc * replacement_synthesis
+    necromass_to_pom = (1 - necromass_proportion_to_lmwc) * replacement_synthesis
+
+    return MicrobialBiomassLoss(
+        maintenance_synthesis=maintenance_synthesis,
+        pom_enzyme_production=pom_enzyme_production,
+        maom_enzyme_production=maom_enzyme_production,
+        necromass_decay_to_lmwc=necromass_to_lmwc,
+        necromass_decay_to_pom=necromass_to_pom,
     )
 
 
-def calculate_maintenance_respiration(
+def calculate_maintenance_biomass_synthesis(
     soil_c_pool_microbe: NDArray[np.float32],
     soil_temp: NDArray[np.float32],
     constants: SoilConsts,
 ) -> NDArray[np.float32]:
-    """Calculate the maintenance respiration of the microbial pool.
+    """Calculate microbial biomass synthesis rate required to offset losses.
+
+    In order for a microbial population to not decline it must synthesise enough new
+    biomass to offset losses. These losses mostly come from cell death and protein
+    decay, but also include loses due to extracellular enzyme excretion.
 
     Args:
         soil_c_pool_microbe: Microbial biomass (carbon) pool [kg C m^-3]
@@ -397,7 +263,8 @@ def calculate_maintenance_respiration(
         constants: Set of constants for the soil model.
 
     Returns:
-        Total maintenance respiration for all microbial biomass
+        The rate of microbial biomass loss that must be matched to maintain a steady
+        population [kg C m^-3 day^-1]
     """
 
     temp_factor = calculate_temperature_effect_on_microbes(
@@ -409,27 +276,6 @@ def calculate_maintenance_respiration(
     return constants.microbial_turnover_rate * temp_factor * soil_c_pool_microbe
 
 
-def calculate_necromass_adsorption(
-    soil_c_pool_microbe: NDArray[np.float32],
-    moisture_scalar: NDArray[np.float32],
-    necromass_adsorption_rate: float,
-) -> NDArray[np.float32]:
-    """Calculate adsorption of microbial necromass to soil minerals.
-
-    Args:
-        soil_c_pool_microbe: Microbial biomass (carbon) pool [kg C m^-3]
-        moisture_scalar: A scalar capturing the impact of soil moisture on process rates
-            [unitless]
-        necromass_adsorption_rate: Rate at which necromass is adsorbed by soil minerals
-            [day^-1]
-
-    Returns:
-        Adsorption of microbial biomass to mineral associated organic matter (MAOM)
-    """
-
-    return necromass_adsorption_rate * moisture_scalar * soil_c_pool_microbe
-
-
 def calculate_carbon_use_efficiency(
     soil_temp: NDArray[np.float32],
     reference_cue: float,
@@ -454,196 +300,157 @@ def calculate_carbon_use_efficiency(
     return reference_cue - cue_with_temperature * (soil_temp - cue_reference_temp)
 
 
-def calculate_microbial_saturation(
-    soil_c_pool_microbe: NDArray[np.float32],
-    half_sat_microbial_activity: float,
-) -> NDArray[np.float32]:
-    """Calculate microbial activity saturation.
-
-    This ensures that microbial activity (per unit biomass) drops as biomass density
-    increases. This is adopted from Abramoff et al. It feels like an assumption that
-    should be revised as the Virtual Rainforest develops.
-
-    Args:
-        soil_c_pool_microbe: Microbial biomass (carbon) pool [kg C m^-3]
-        half_sat_microbial_activity: Half saturation constant for microbial activity
-
-    Returns:
-        A rescaling of microbial biomass that takes into account activity saturation
-        with increasing biomass density
-    """
-
-    return soil_c_pool_microbe / (soil_c_pool_microbe + half_sat_microbial_activity)
-
-
-def calculate_microbial_pom_mineralisation_saturation(
-    soil_c_pool_microbe: NDArray[np.float32],
-    half_sat_microbial_mineralisation: float,
+def calculate_enzyme_turnover(
+    enzyme_pool: NDArray[np.float32], turnover_rate: float
 ) -> NDArray[np.float32]:
-    """Calculate microbial POM mineralisation saturation (with increasing biomass).
-
-    This ensures that microbial mineralisation of POM (per unit biomass) drops as
-    biomass density increases. This is adopted from Abramoff et al. This function is
-    very similar to the
-    :func:`~virtual_rainforest.models.soil.carbon.calculate_microbial_saturation`
-    function. They could in theory be reworked into a single function, but it doesn't
-    seem worth the effort as we do not anticipate using biomass saturation functions
-    beyond the first model draft.
+    """Calculate the turnover rate of a specific enzyme class.
 
     Args:
-        soil_c_pool_microbe: Microbial biomass (carbon) pool [kg C m^-3]
-        half_sat_microbial_mineralisation: Half saturation constant for microbial
-            mineralisation of POM
+        enzyme_pool: The pool size for the enzyme class in question [kg C m^-3]
+        turnover_rate: The rate at which enzymes in the pool turnover [day^-1]
 
     Returns:
-        A rescaling of microbial biomass that takes into account POM mineralisation rate
-        saturation with increasing biomass density
+        The rate at which enzymes are lost from the pool [kg C m^-3 day^-1]
     """
 
-    return soil_c_pool_microbe / (
-        soil_c_pool_microbe + half_sat_microbial_mineralisation
-    )
-
-
-def calculate_pom_decomposition_saturation(
-    soil_c_pool_pom: NDArray[np.float32],
-    half_sat_pom_decomposition: float,
-) -> NDArray[np.float32]:
-    """Calculate particulate organic matter (POM) decomposition saturation.
-
-    This ensures that decomposition of POM to low molecular weight carbon (LMWC)
-    saturates with increasing POM. This effect arises from the saturation of enzymes
-    with increasing substrate.
-
-    Args:
-        soil_c_pool_pom: Particulate organic matter (carbon) pool [kg C m^-3]
-        half_sat_pom_decomposition: Half saturation constant for POM decomposition
-
-    Returns:
-        The saturation of the decomposition process
-    """
-
-    return soil_c_pool_pom / (soil_c_pool_pom + half_sat_pom_decomposition)
+    return turnover_rate * enzyme_pool
 
 
 def calculate_microbial_carbon_uptake(
     soil_c_pool_lmwc: NDArray[np.float32],
     soil_c_pool_microbe: NDArray[np.float32],
     water_factor: NDArray[np.float32],
+    pH_factor: NDArray[np.float32],
     soil_temp: NDArray[np.float32],
     constants: SoilConsts,
-) -> NDArray[np.float32]:
-    """Calculate amount of labile carbon taken up by microbes.
+) -> tuple[NDArray[np.float32], NDArray[np.float32]]:
+    """Calculate uptake and assimilation of labile carbon by microbes.
+
+    This function starts by calculating the impact that environmental factors have on
+    the rate and saturation constants for microbial uptake. These constants are then
+    used to calculate the rate of uptake of labile carbon. Carbon use efficiency is then
+    calculated and used to find how much of this carbon ends up assimilated as biomass
+    (rather than respired).
 
     Args:
         soil_c_pool_lmwc: Low molecular weight carbon pool [kg C m^-3]
         soil_c_pool_microbe: Microbial biomass (carbon) pool [kg C m^-3]
         water_factor: A factor capturing the impact of soil water potential on microbial
             rates [unitless]
+        pH_factor: A factor capturing the impact of soil pH on microbial rates
+            [unitless]
         soil_temp: soil temperature for each soil grid cell [degrees C]
         constants: Set of constants for the soil model.
 
     Returns:
-        Uptake of low molecular weight carbon (LMWC) by the soil microbial biomass.
+        A tuple containing the uptake rate of low molecular weight carbon (LMWC) by the
+        soil microbial biomass, and the rate at which this causes microbial biomass to
+        increase.
     """
 
-    # Calculate carbon use efficiency and microbial saturation
+    # Calculate carbon use efficiency
     carbon_use_efficency = calculate_carbon_use_efficiency(
         soil_temp,
         constants.reference_cue,
         constants.cue_reference_temp,
         constants.cue_with_temperature,
     )
-    microbial_saturation = calculate_microbial_saturation(
-        soil_c_pool_microbe, constants.half_sat_microbial_activity
-    )
-    temp_factor = calculate_temperature_effect_on_microbes(
+    # Calculate impact of temperature on the rate and saturation constants
+    temp_factor_rate = calculate_temperature_effect_on_microbes(
         soil_temperature=soil_temp,
         activation_energy=constants.activation_energy_microbial_uptake,
         reference_temperature=constants.arrhenius_reference_temp,
     )
+    temp_factor_saturation = calculate_temperature_effect_on_microbes(
+        soil_temperature=soil_temp,
+        activation_energy=constants.activation_energy_labile_C_saturation,
+        reference_temperature=constants.arrhenius_reference_temp,
+    )
+    # Then use to calculate rate constant and saturation constant (which also change
+    # with other environmental conditions)
+    rate_constant = (
+        constants.max_uptake_rate_labile_C * temp_factor_rate * water_factor * pH_factor
+    )
+    saturation_constant = constants.half_sat_labile_C_uptake * temp_factor_saturation
+
+    # Calculate both the rate of carbon uptake, and the rate at which this carbon is
+    # assimilated into microbial biomass.
+    uptake_rate = rate_constant * np.divide(
+        (soil_c_pool_lmwc * soil_c_pool_microbe),
+        (soil_c_pool_lmwc + saturation_constant),
+    )
+    assimilation_rate = uptake_rate * carbon_use_efficency
 
     # TODO - the quantities calculated above can be used to calculate the carbon
     # respired instead of being uptaken. This isn't currently of interest, but will be
     # in future
 
-    return (
-        constants.max_uptake_rate_labile_C
-        * water_factor
-        * temp_factor
-        * soil_c_pool_lmwc
-        * microbial_saturation
-        * carbon_use_efficency
-    )
-
-
-def calculate_labile_carbon_leaching(
-    soil_c_pool_lmwc: NDArray[np.float32],
-    moisture_scalar: NDArray[np.float32],
-    leaching_rate: float,
-) -> NDArray[np.float32]:
-    """Calculate rate at which labile carbon is leached.
-
-    This is adopted from Abramoff et al. We definitely need to give more thought to how
-    we model leaching.
-
-    Args:
-        soil_c_pool_lmwc: Low molecular weight carbon pool [kg C m^-3]
-        moisture_scalar: A scalar capturing the impact of soil moisture on process rates
-            [unitless]
-        leaching_rate: The rate at which labile carbon leaches from the soil [day^-1]
-
-    Returns:
-        The amount of labile carbon leached
-    """
-
-    return leaching_rate * moisture_scalar * soil_c_pool_lmwc
+    return uptake_rate, assimilation_rate
 
 
-def calculate_pom_decomposition(
-    soil_c_pool_pom: NDArray[np.float32],
-    soil_c_pool_microbe: NDArray[np.float32],
+def calculate_enzyme_mediated_decomposition(
+    soil_c_pool: NDArray[np.float32],
+    soil_enzyme: NDArray[np.float32],
     water_factor: NDArray[np.float32],
+    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 decomposition of particulate organic matter into labile carbon (LMWC).
+    """Calculate rate of a enzyme mediated decomposition process.
 
-    This is adopted from Abramoff et al. We definitely want to change this down the line
-    to something that uses enzymes explicitly.
+    This function calculates various environmental factors that effect enzyme activity,
+    then uses these to find environmental adjusted rate and saturation constants. These
+    are then used to find the decomposition rate of the pool in question.
 
     Args:
-        soil_c_pool_pom: Particulate organic matter pool [kg C m^-3]
-        soil_c_pool_microbe: Microbial biomass (carbon) pool [kg C m^-3]
+        soil_c_pool: Size of organic matter pool [kg C m^-3]
+        soil_enzyme: Amount of enzyme class which breaks down the organic matter pool in
+            question [kg C m^-3]
         water_factor: A factor capturing the impact of soil water potential on microbial
             rates [unitless]
+        pH_factor: A factor capturing the impact of soil pH on microbial rates
+            [unitless]
+        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 amount of particulate organic matter (POM) decomposed into labile carbon
-            (LMWC)
+        The rate of decomposition of the organic matter pool in question [kg C m^-3
+        day^-1]
     """
 
-    # Calculate the two relevant saturations
-    saturation_with_biomass = calculate_microbial_pom_mineralisation_saturation(
-        soil_c_pool_microbe, constants.half_sat_microbial_pom_mineralisation
+    # Calculate the factors which impact the rate and saturation constants
+    temp_factor_rate = calculate_temperature_effect_on_microbes(
+        soil_temperature=soil_temp,
+        activation_energy=activation_energy_rate,
+        reference_temperature=reference_temp,
     )
-    saturation_with_pom = calculate_pom_decomposition_saturation(
-        soil_c_pool_pom, constants.half_sat_pom_decomposition
+    temp_factor_saturation = calculate_temperature_effect_on_microbes(
+        soil_temperature=soil_temp,
+        activation_energy=activation_energy_sat,
+        reference_temperature=reference_temp,
     )
 
-    # Calculate the impact of temperature on the rate
-    temp_factor = calculate_temperature_effect_on_microbes(
-        soil_temperature=soil_temp,
-        activation_energy=constants.activation_energy_pom_decomp,
-        reference_temperature=constants.arrhenius_reference_temp,
+    # Calculate the adjusted rate and saturation constants
+    rate_constant = max_decomp_rate * temp_factor_rate * water_factor * pH_factor
+    saturation_constant = (
+        half_saturation * temp_factor_saturation * clay_factor_saturation
     )
 
     return (
-        constants.max_decomp_rate_pom
-        * saturation_with_pom
-        * saturation_with_biomass
-        * water_factor
-        * temp_factor
+        rate_constant * soil_enzyme * soil_c_pool / (saturation_constant + soil_c_pool)
     )
diff --git a/virtual_rainforest/models/soil/constants.py b/virtual_rainforest/models/soil/constants.py
index a4858ec01..9e7d5e3c1 100644
--- a/virtual_rainforest/models/soil/constants.py
+++ b/virtual_rainforest/models/soil/constants.py
@@ -7,53 +7,12 @@
 
 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)
 class SoilConsts(ConstantsDataclass):
-    """Dataclass to store all constants for the `soil` model.
-
-    All constants are taken from :cite:t:`abramoff_millennial_2018` unless otherwise
-    stated.
-    """
-
-    binding_with_ph_slope: float = -0.186
-    """Change in the binding affinity of soil mineral with pH.
-
-    Units of [pH^-1]. From linear regression :cite:p:`mayes_relation_2012`."""
-
-    binding_with_ph_intercept: float = -0.216 + 3.0
-    """Binding affinity of soil minerals at zero pH.
-
-    Unit of [log(m^3 kg^-1)]. n.b. +3 converts from mg^-1 to kg^-1 and L to m^3. From
-    linear regression :cite:p:`mayes_relation_2012`.
-    """
-
-    max_sorption_with_clay_slope: float = 0.483
-    """Change in the maximum size of the MAOM pool with increasing clay content.
-
-    Units of [(% clay)^-1]. From linear regression :cite:p:`mayes_relation_2012`.
-    """
-
-    max_sorption_with_clay_intercept: float = 2.328 - 6.0
-    """Maximum size of the MAOM pool at zero clay content.
-
-    Unit of [log(kg C kg soil ^-1)]. n.b. -6 converts from mg to kg. From linear
-    regression :cite:p:`mayes_relation_2012`.
-    """
-
-    moisture_scalar_coefficient: float = 30.0
-    """Used in :cite:t:`abramoff_millennial_2018`, can't find original source.
-
-    Value at zero relative water content (RWC) [unit less].
-    """
-
-    moisture_scalar_exponent: float = 9.0
-    """Used in :cite:t:`abramoff_millennial_2018`, can't find original source.
-
-    Units of [(RWC)^-1]
-    """
+    """Dataclass to store all constants for the `soil` model."""
 
     reference_cue: float = 0.6
     """Carbon use efficiency of community at the reference temperature [no units].
@@ -73,27 +32,6 @@ class SoilConsts(ConstantsDataclass):
     Default value taken from :cite:t:`abramoff_millennial_2018`.
     """
 
-    necromass_adsorption_rate: float = 0.025
-    """Rate at which necromass is adsorbed by soil minerals [day^-1].
-
-    Taken from :cite:t:`abramoff_millennial_2018`, where it was obtained by calibration.
-    """
-
-    half_sat_microbial_activity: float = 0.0072
-    """Half saturation constant for microbial activity (with increasing biomass).
-
-    Units of [kg C m^-2].
-    """
-
-    half_sat_microbial_pom_mineralisation: float = 0.012
-    """Half saturation constant for microbial POM mineralisation [kg C m^-2]."""
-
-    leaching_rate_labile_carbon: float = 1.5e-3
-    """Leaching rate for labile carbon (lmwc) [day^-1]."""
-
-    half_sat_pom_decomposition: float = 0.150
-    """Half saturation constant for POM decomposition to LMWC [kg C m^-2]."""
-
     soil_microbe_water_potential_optimum: float = -3.0
     """The water potential at which soil microbial rates are maximised [kPa].
 
@@ -106,10 +44,10 @@ class SoilConsts(ConstantsDataclass):
     Value is taken from :cite:t`moyano_responses_2013`.
     """
 
-    moisture_response_curvature: float = 1.47
-    """Curvature of the soil microbial moisture response function [unitless].
+    microbial_water_response_curvature: float = 1.47
+    """Curvature of function for response of soil microbial rates to water potential.
 
-    Value is taken from :cite:t`moyano_responses_2013`.
+    [unitless]. Value is taken from :cite:t`moyano_responses_2013`.
     """
 
     arrhenius_reference_temp: float = 12.0
@@ -131,7 +69,7 @@ class SoilConsts(ConstantsDataclass):
     once fungi are added.
     """
 
-    activation_energy_microbial_uptake = 47000
+    activation_energy_microbial_uptake: float = 47000
     """Activation energy for microbial uptake of low molecular weight carbon [J K^-1].
 
     Value taken from :cite:t:`wang_development_2013`. The maximum labile carbon uptake
@@ -139,26 +77,84 @@ class SoilConsts(ConstantsDataclass):
     :attr:`max_uptake_rate_labile_C`.
     """
 
-    # TODO - Add another set of constants once we start tracking lignin
-    max_decomp_rate_pom: float = 0.2
+    half_sat_labile_C_uptake: float = 0.364
+    """Half saturation constant for microbial uptake of labile carbon (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_labile_C_saturation`.
+    """
+
+    activation_energy_labile_C_saturation: float = 30000
+    """Activation energy for labile C uptake saturation constant [J K^-1].
+
+    Taken from :cite:t:`wang_development_2013`.
+    """
+
+    half_sat_pom_decomposition: float = 70.0
+    """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]. 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
+    """Activation energy for POM decomposition saturation constant [J K^-1].
+
+    Taken from :cite:t:`wang_development_2013`.
+    """
+
+    max_decomp_rate_pom: float = 60.0
     """Maximum rate for particulate organic matter break down (at reference temp).
 
     Units of [day^-1]. The reference temperature is given by
     :attr:`arrhenius_reference_temp`, and the corresponding activation energy is given
-    by :attr:`activation_energy_pom_decomp`.
-
-    TODO - Once enzymes are included this should take the value of 200.0 [day^-1].
+    by :attr:`activation_energy_pom_decomp_rate`.
 
     TODO - Source of this constant is not completely clear, investigate this further
     once lignin chemistry is added.
     """
 
-    activation_energy_pom_decomp = 37000
+    activation_energy_pom_decomp_rate: float = 37000
     """Activation energy for decomposition of particulate organic matter [J K^-1].
 
     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].
@@ -177,3 +173,103 @@ class SoilConsts(ConstantsDataclass):
     Value taken from :cite:t:`wang_development_2013`. The microbial turnover rate that
     this activation energy corresponds to is given by :attr:`microbial_turnover_rate`.
     """
+
+    # TODO - At some point I need to split these enzyme constants into fungi and
+    # bacteria specific constants
+    pom_enzyme_turnover_rate: float = 2.4e-2
+    """Turnover rate for POM degrading enzymes [day^-1].
+
+    Value taken from :cite:t:`wang_development_2013`.
+    """
+
+    maom_enzyme_turnover_rate: float = 2.4e-2
+    """Turnover rate for MAOM degrading enzymes [day^-1].
+
+    Value taken from :cite:t:`wang_development_2013`.
+    """
+
+    maintenance_pom_enzyme: float = 1e-2
+    """Fraction of maintenance synthesis used to produce POM degrading enzymes.
+
+    [unitless]. Value taken from :cite:t:`wang_development_2013`.
+    """
+
+    maintenance_maom_enzyme: float = 1e-2
+    """Fraction of maintenance synthesis used to produce MAOM degrading enzymes.
+
+    [unitless]. Value taken from :cite:t:`wang_development_2013`.
+    """
+
+    necromass_to_lmwc: float = 0.25
+    """Proportion of necromass that flows to LMWC rather than POM [unitless].
+
+    Value taken from :cite:t:`wang_development_2013`.
+    """
+
+    # TODO - The 4 constants below should take different values for fungi and bacteria,
+    # once that separation is implemented.
+    min_pH_microbes: float = 2.5
+    """Soil pH below which microbial activity is completely inhibited [unitless].
+
+    This value cannot be larger than :attr:`lowest_optimal_pH_microbes`. The default
+    value was obtained by averaging the fungi and bacteria specific values given in
+    :cite:t:`orwin_organic_2011`.
+    """
+
+    lowest_optimal_pH_microbes: float = 4.5
+    """Soil pH above which microbial activity is not inhibited at all [unitless].
+
+    This value cannot be smaller than :attr:`min_pH_microbes` or larger than
+    :attr:`highest_optimal_pH_microbes`. The default value was obtained by averaging the
+    fungi and bacteria specific values given in :cite:t:`orwin_organic_2011`.
+    """
+
+    highest_optimal_pH_microbes: float = 7.5
+    """Soil pH below which microbial activity is not inhibited at all [unitless].
+
+    This value cannot be smaller than :attr:`lowest_optimal_pH_microbes` or larger than
+    :attr:`max_pH_microbes`. The default value was obtained by averaging the fungi
+    and bacteria specific values given in :cite:t:`orwin_organic_2011`.
+    """
+
+    max_pH_microbes: float = 11.0
+    """Soil pH above which microbial activity is completely inhibited [unitless].
+
+    This value cannot be smaller than :attr:`highest_optimal_pH_microbes`. The default
+    value was obtained by averaging the fungi and bacteria specific values given in
+    :cite:t:`orwin_organic_2011`.
+    """
+
+    base_soil_protection: float = 0.694
+    """Basal change in saturation constants due to soil structure [unitless]
+
+    This value is multiplicative and is taken from :cite:t:`fatichi_mechanistic_2019`.
+    """
+
+    soil_protection_with_clay: float = 1.36
+    """Rate at which soil protection of carbon increases with clay content [unitless].
+
+    This protection contributes multiplicatively to the effective saturation constant.
+    The value of this constant is taken from :cite:t:`fatichi_mechanistic_2019`.
+    """
+
+    clay_necromass_decay_exponent: float = -0.8
+    """Change in proportion of necromass which decays with increasing soil clay content.
+
+    [unitless]. The function this is used in is an exponential, and the sign should be
+    negative so increases in clay leads to a lower proportion of necromass decaying to
+    LMWC. The value of this constant is taken from :cite:t:`fatichi_mechanistic_2019`.
+    """
+
+    pom_decomposition_fraction_lmwc: float = 0.5
+    """Fraction of decomposed POM that becomes LMWC rather than MAOM [unitless].
+
+    Value taken from :cite:t:`wang_development_2013`.
+    """
+
+    solubility_coefficient_lmwc: float = 0.05
+    """Solubility coefficient for low molecular weight organic carbon [unitless].
+
+    Value taken from :cite:t:`fatichi_mechanistic_2019`, where it is estimated in quite
+    a loose manner.
+    """
diff --git a/virtual_rainforest/models/soil/env_factors.py b/virtual_rainforest/models/soil/env_factors.py
new file mode 100644
index 000000000..778f895d2
--- /dev/null
+++ b/virtual_rainforest/models/soil/env_factors.py
@@ -0,0 +1,306 @@
+"""The ``models.soil.env_factors`` module contains functions that are used to
+capture the impact that environmental factors have on microbial rates. These include
+temperature, soil water potential, pH and soil texture.
+"""  # noqa: D205, D415
+
+from dataclasses import dataclass
+
+import numpy as np
+from numpy.typing import NDArray
+from scipy.constants import convert_temperature, gas_constant
+
+from virtual_rainforest.core.logger import LOGGER
+from virtual_rainforest.models.soil.constants import SoilConsts
+
+
+@dataclass
+class EnvironmentalEffectFactors:
+    """The various factors through which the environment effects soil cycling rates."""
+
+    water: NDArray[np.float32]
+    """Impact of soil water potential on enzymatic rates [unitless]."""
+    pH: NDArray[np.float32]
+    """Impact of soil pH on enzymatic rates [unitless]."""
+    clay_saturation: NDArray[np.float32]
+    """Impact of soil clay fraction on enzyme saturation constants [unitless]."""
+    clay_decay: NDArray[np.float32]
+    """Impact of soil clay fraction on necromass decay destination [unitless]."""
+
+
+def calculate_environmental_effect_factors(
+    soil_water_potential: NDArray[np.float32],
+    pH: NDArray[np.float32],
+    clay_fraction: NDArray[np.float32],
+    constants: SoilConsts,
+) -> EnvironmentalEffectFactors:
+    """Calculate the effects that the environment has on relevant biogeochemical rates.
+
+    For each environmental effect a multiplicative factor is calculated, and all of them
+    are returned in a single object for use elsewhere in the soil model.
+
+    Args:
+        soil_water_potential: Soil water potential for each grid cell [kPa]
+        pH: pH values for each soil grid cell [unitless]
+        clay_fraction: The clay fraction for each soil grid cell [unitless]
+        constants: Set of constants for the soil model
+
+    Returns:
+        An object containing four environmental factors, one for the effect of water
+        potential on enzyme rates, one for the effect of pH on enzyme rates, one for the
+        effect of clay fraction on enzyme saturation, and one for the effect of clay on
+        necromass decay destination.
+    """
+
+    # Calculate the impact that each environment variable has on the relevant
+    # biogeochemical soil processes
+    water_factor = calculate_water_potential_impact_on_microbes(
+        water_potential=soil_water_potential,
+        water_potential_halt=constants.soil_microbe_water_potential_halt,
+        water_potential_opt=constants.soil_microbe_water_potential_optimum,
+        response_curvature=constants.microbial_water_response_curvature,
+    )
+    pH_factor = calculate_pH_suitability(
+        soil_pH=pH,
+        maximum_pH=constants.max_pH_microbes,
+        minimum_pH=constants.min_pH_microbes,
+        lower_optimum_pH=constants.lowest_optimal_pH_microbes,
+        upper_optimum_pH=constants.highest_optimal_pH_microbes,
+    )
+    clay_factor_saturation = calculate_clay_impact_on_enzyme_saturation(
+        clay_fraction=clay_fraction,
+        base_protection=constants.base_soil_protection,
+        protection_with_clay=constants.soil_protection_with_clay,
+    )
+    clay_factor_decay = calculate_clay_impact_on_necromass_decay(
+        clay_fraction=clay_fraction,
+        decay_exponent=constants.clay_necromass_decay_exponent,
+    )
+
+    # Combine all factors into a single EnvironmentalFactors object
+    return EnvironmentalEffectFactors(
+        water=water_factor,
+        pH=pH_factor,
+        clay_saturation=clay_factor_saturation,
+        clay_decay=clay_factor_decay,
+    )
+
+
+def calculate_temperature_effect_on_microbes(
+    soil_temperature: NDArray[np.float32],
+    activation_energy: float,
+    reference_temperature: float,
+) -> NDArray[np.float32]:
+    """Calculate the effect that temperature has on microbial metabolic rates.
+
+    This uses a standard Arrhenius equation to calculate the impact of temperature.
+
+    This function takes temperatures in Celsius as inputs and converts them into Kelvin
+    for use in the Arrhenius equation. TODO - review this after we have decided how to
+    handle these conversions in general.
+
+    Args:
+        soil_temperature: The temperature of the soil [C]
+        activation_energy: Energy of activation [J mol^-1]
+        soil_temperature: The reference temperature of the Arrhenius equation [C]
+
+    Returns:
+        A multiplicative factor capturing the effect of temperature on microbial rates
+    """
+
+    # Convert the temperatures to Kelvin
+    soil_temp_in_kelvin = convert_temperature(
+        soil_temperature, old_scale="Celsius", new_scale="Kelvin"
+    )
+    ref_temp_in_kelvin = convert_temperature(
+        reference_temperature, old_scale="Celsius", new_scale="Kelvin"
+    )
+
+    return np.exp(
+        (-activation_energy / gas_constant)
+        * ((1 / (soil_temp_in_kelvin)) - (1 / (ref_temp_in_kelvin)))
+    )
+
+
+def calculate_water_potential_impact_on_microbes(
+    water_potential: NDArray[np.float32],
+    water_potential_halt: float,
+    water_potential_opt: float,
+    response_curvature: float,
+) -> NDArray[np.float32]:
+    """Calculate the effect that soil water potential has on microbial rates.
+
+    This function only returns valid output for soil water potentials that are less than
+    the optimal water potential.
+
+    Args:
+        water_potential: Soil water potential [kPa]
+        water_potential_halt: Water potential at which all microbial activity stops
+            [kPa]
+        water_potential_opt: Optimal water potential for microbial activity [kPa]
+        response_curvature: Parameter controlling the curvature of function that
+            captures the response of microbial rates to water potential [unitless]
+
+    Returns:
+        A multiplicative factor capturing the impact of moisture on soil microbe rates
+        decomposition [unitless]
+    """
+
+    # If the water potential is greater than the optimal then the function produces NaNs
+    # so the simulation should be interrupted
+    if np.any(water_potential > water_potential_opt):
+        err = ValueError("Water potential greater than minimum value")
+        LOGGER.critical(err)
+        raise err
+
+    # Calculate how much moisture suppresses microbial activity
+    supression = (
+        (np.log10(-water_potential) - np.log10(-water_potential_opt))
+        / (np.log10(-water_potential_halt) - np.log10(-water_potential_opt))
+    ) ** response_curvature
+
+    return 1 - supression
+
+
+def calculate_pH_suitability(
+    soil_pH: NDArray[np.float32],
+    maximum_pH: float,
+    minimum_pH: float,
+    upper_optimum_pH: float,
+    lower_optimum_pH: float,
+) -> NDArray[np.float32]:
+    """Calculate the suitability of the soil pH for microbial activity.
+
+    This function is taken from :cite:t:`orwin_organic_2011`. pH values within the
+    optimal range are assumed to cause no microbial inhibition, and pH values above a
+    certain value or below a certain value are assumed to cause total inhibition. Linear
+    declines then occur between the edges of the optimal range and the zone of total
+    inhibition.
+
+    Args:
+        soil_pH: The pH of the soil [unitless]
+        maximum_pH: pH above which microbial rates are completely inhibited [unitless]
+        minimum_pH: pH below which microbial rates are completely inhibited [unitless]
+        upper_optimum_pH: pH above which suitability declines [unitless]
+        lower_optimum_pH: pH below which suitability declines [unitless]
+
+    Returns:
+        A multiplicative factor capturing the effect of pH on microbial rates
+    """
+
+    # TODO - This check is necessary to prevent nonsensical output being generated,
+    # however it could be done when constants are loaded, rather than for every function
+    # call
+    if (
+        maximum_pH <= upper_optimum_pH
+        or upper_optimum_pH <= lower_optimum_pH
+        or lower_optimum_pH <= minimum_pH
+    ):
+        to_raise = ValueError("At least one pH threshold has an invalid value!")
+        LOGGER.error(to_raise)
+        raise to_raise
+
+    pH_factors = np.full(len(soil_pH), np.nan)
+
+    # zero below minimum or above maximum pH
+    pH_factors[soil_pH < minimum_pH] = 0
+    pH_factors[soil_pH > maximum_pH] = 0
+    # and one between the two thresholds
+    pH_factors[(lower_optimum_pH <= soil_pH) & (soil_pH <= upper_optimum_pH)] = 1
+
+    # Find points that lie between optimal region and maximum/minimum
+    between_opt_and_min = (minimum_pH <= soil_pH) & (soil_pH < lower_optimum_pH)
+    between_opt_and_max = (upper_optimum_pH < soil_pH) & (soil_pH <= maximum_pH)
+
+    # Linear increase from minimum pH value to lower threshold
+    pH_factors[between_opt_and_min] = (soil_pH[between_opt_and_min] - minimum_pH) / (
+        lower_optimum_pH - minimum_pH
+    )
+    # Linear decrease from the upper threshold to maximum pH
+    pH_factors[between_opt_and_max] = (
+        soil_pH[between_opt_and_max] - upper_optimum_pH
+    ) / (maximum_pH - upper_optimum_pH)
+
+    return pH_factors
+
+
+def calculate_clay_impact_on_enzyme_saturation(
+    clay_fraction: NDArray[np.float32],
+    base_protection: float,
+    protection_with_clay: float,
+) -> NDArray[np.float32]:
+    """Calculate the impact that the soil clay fraction has on enzyme saturation.
+
+    This factor impacts enzyme saturation constants, based on the assumption that finely
+    textured soils will restrict enzyme access to available C substrates (which protects
+    them). Its form is taken from :cite:t:`fatichi_mechanistic_2019`.
+
+    Args:
+        clay_fraction: The fraction of the soil which is clay [unitless]
+        base_protection: The protection that a soil with no clay provides [unitless]
+        protection_with_clay: The rate at which protection increases with increasing
+           clay [unitless]
+
+    Returns:
+        A multiplicative factor capturing how much the protection due to soil structure
+        changes the effective saturation constant by [unitless]
+    """
+
+    return base_protection + protection_with_clay * clay_fraction
+
+
+def calculate_clay_impact_on_necromass_decay(
+    clay_fraction: NDArray[np.float32], decay_exponent: float
+) -> NDArray[np.float32]:
+    """Calculate the impact that soil clay has on necromass decay to LMWC.
+
+    Necromass which doesn't breakdown fully gets added to the POM pool instead.
+
+    Args:
+        clay_fraction: The fraction of the soil which is clay [unitless]
+        sorption_exponent: Controls the impact that differences in soil clay content
+            have on the proportion of necromass that decays to LMWC [unitless]
+
+    Returns:
+        A multiplicative factor capturing the impact that soil clay has on the
+        proportion of necromass decay which sorbs to form POM [unitless]
+    """
+
+    return np.exp(decay_exponent * clay_fraction)
+
+
+def calculate_leaching_rate(
+    solute_density: NDArray[np.float32],
+    vertical_flow_rate: NDArray[np.float32],
+    soil_moisture: NDArray[np.float32],
+    solubility_coefficient: float,
+    soil_layer_thickness: float,
+) -> NDArray[np.float32]:
+    """Calculate leaching rate for a given solute based on flow rate.
+
+    This functional form is adapted from :cite:t:`porporato_hydrologic_2003`. The amount
+    of solute that is expected to be found in dissolved form is calculated by
+    multiplying the solute density by its solubility coefficient. This is then
+    multiplied by the frequency with which the water column is completely replaced, i.e.
+    the ratio of vertical flow rate to water column height.
+
+    Args:
+        solute_density: The density of the solute in the soil [kg solute m^-3]
+        vertical_flow_rate: Rate of flow downwards through the soil [mm day^-1]
+        soil_moisture: Volumetric relative water content of the soil [unitless]
+        solubility_coefficient: The solubility coefficient of the solute in question
+            [unitless]
+        soil_layer_thickness: Thickness of the biogeochemically active soil layer [m]
+
+    Returns:
+        The rate at which the solute in question is leached [kg solute m^-3 day^-1]
+    """
+
+    # Vertical flow rate has to be converted into m day^-1
+    vert_flow_meters = vertical_flow_rate / 1e3
+
+    return (
+        solubility_coefficient
+        * solute_density
+        * vert_flow_meters
+        / (soil_moisture * soil_layer_thickness)
+    )
diff --git a/virtual_rainforest/models/soil/soil_model.py b/virtual_rainforest/models/soil/soil_model.py
index 3e9139a5f..3852656ae 100644
--- a/virtual_rainforest/models/soil/soil_model.py
+++ b/virtual_rainforest/models/soil/soil_model.py
@@ -28,6 +28,7 @@
 
 from virtual_rainforest.core.base_model import BaseModel
 from virtual_rainforest.core.config import Config
+from virtual_rainforest.core.constants import CoreConsts
 from virtual_rainforest.core.constants_loader import load_constants
 from virtual_rainforest.core.data import Data
 from virtual_rainforest.core.exceptions import InitialisationError
@@ -70,7 +71,7 @@ class SoilModel(BaseModel):
         ("soil_c_pool_pom", ("spatial",)),
         ("pH", ("spatial",)),
         ("bulk_density", ("spatial",)),
-        ("percent_clay", ("spatial",)),
+        ("clay_fraction", ("spatial",)),
     )
     """Required initialisation variables for the soil model.
 
@@ -82,6 +83,8 @@ class SoilModel(BaseModel):
         "soil_c_pool_lmwc",
         "soil_c_pool_microbe",
         "soil_c_pool_pom",
+        "soil_enzyme_pom",
+        "soil_enzyme_maom",
     )
     """Variables updated by the soil model."""
 
@@ -91,7 +94,8 @@ def __init__(
         update_interval: Quantity,
         soil_layers: list[float],
         canopy_layers: int,
-        constants: SoilConsts,
+        model_constants: SoilConsts,
+        core_constants: CoreConsts,
         **kwargs: Any,
     ):
         super().__init__(data, update_interval, **kwargs)
@@ -118,8 +122,10 @@ def __init__(
         """The layer in the data object representing the first soil layer."""
         # TODO - At the moment the soil model only cares about the very top layer. As
         # both the soil and abiotic models get more complex this might well change.
-        self.constants = constants
-        """Set of constants for the soil model"""
+        self.model_constants = model_constants
+        """Set of constants for the soil model."""
+        self.core_constants = core_constants
+        """Set of constants shared between models."""
 
     @classmethod
     def from_config(
@@ -142,13 +148,21 @@ def from_config(
         canopy_layers = config["core"]["layers"]["canopy_layers"]
 
         # Load in the relevant constants
-        constants = load_constants(config, "soil", "SoilConsts")
+        soil_constants = load_constants(config, "soil", "SoilConsts")
+        core_constants = load_constants(config, "core", "CoreConsts")
 
         LOGGER.info(
             "Information required to initialise the soil model successfully "
             "extracted."
         )
-        return cls(data, update_interval, soil_layers, canopy_layers, constants)
+        return cls(
+            data=data,
+            update_interval=update_interval,
+            soil_layers=soil_layers,
+            canopy_layers=canopy_layers,
+            model_constants=soil_constants,
+            core_constants=core_constants,
+        )
 
     def setup(self) -> None:
         """Placeholder function to setup up the soil model."""
@@ -201,19 +215,25 @@ def integrate(self) -> dict[str, DataArray]:
         # Construct vector of initial values y0
         y0 = np.concatenate(
             [
-                self.data[str(name)].to_numpy()
-                for name in self.data.data.keys()
-                if str(name).startswith("soil_c_pool_")
+                self.data[name].to_numpy()
+                for name in map(str, self.data.data.keys())
+                if name.startswith("soil_c_pool_") or name.startswith("soil_enzyme_")
             ]
         )
 
         # Find and store order of pools
         delta_pools_ordered = {
-            str(name): np.array([])
-            for name in self.data.data.keys()
-            if str(name).startswith("soil_c_pool_")
+            name: np.array([])
+            for name in map(str, self.data.data.keys())
+            if name.startswith("soil_c_pool_") or name.startswith("soil_enzyme_")
         }
 
+        # TODO - This is a work around as (to the best of my knowledge) the hydrology
+        # model only currently gives total flow (over the entire time step) rather than
+        # flow rates
+        # Convert vertical flow into per day units
+        vertical_flow_per_day = self.data["vertical_flow"].to_numpy() / update_time
+
         # Carry out simulation
         output = solve_ivp(
             construct_full_soil_model,
@@ -221,10 +241,12 @@ def integrate(self) -> dict[str, DataArray]:
             y0,
             args=(
                 self.data,
+                vertical_flow_per_day,
                 no_cells,
                 self.top_soil_layer_index,
                 delta_pools_ordered,
-                self.constants,
+                self.model_constants,
+                self.core_constants,
             ),
         )
 
@@ -252,10 +274,13 @@ def construct_full_soil_model(
     t: float,
     pools: NDArray[np.float32],
     data: Data,
+    # TODO - Remove this as an argument once vertical flow is averaged per day
+    vertical_flow_per_day: NDArray[np.float32],
     no_cells: int,
     top_soil_layer_index: int,
     delta_pools_ordered: dict[str, NDArray[np.float32]],
-    constants: SoilConsts,
+    model_constants: SoilConsts,
+    core_constants: CoreConsts,
 ) -> NDArray[np.float32]:
     """Function that constructs the full soil model in a solve_ivp friendly form.
 
@@ -265,11 +290,13 @@ def construct_full_soil_model(
             integrated.
         pools: An array containing all soil pools in a single vector
         data: The data object, used to populate the arguments i.e. pH and bulk density
+        vertical_flow_per_day: Rate of vertical water flow through soil [mm day^-1]
         no_cells: Number of grid cells the integration is being performed over
         top_soil_layer_index: Index for layer in data object representing top soil layer
         delta_pools_ordered: Dictionary to store pool changes in the order that pools
             are stored in the initial condition vector.
-        constants: Set of constants for the soil model.
+        model_constants: Set of constants for the soil model.
+        core_constants: Set of constants shared between models.
 
     Returns:
         The rate of change for each soil pool
@@ -289,11 +316,13 @@ def construct_full_soil_model(
         bulk_density=data["bulk_density"].to_numpy(),
         soil_moisture=data["soil_moisture"][top_soil_layer_index].to_numpy(),
         soil_water_potential=data["matric_potential"][top_soil_layer_index].to_numpy(),
+        vertical_flow_rate=vertical_flow_per_day,
         soil_temp=data["soil_temperature"][top_soil_layer_index].to_numpy(),
-        percent_clay=data["percent_clay"].to_numpy(),
+        clay_fraction=data["clay_fraction"].to_numpy(),
         mineralisation_rate=data["litter_C_mineralisation_rate"].to_numpy(),
         delta_pools_ordered=delta_pools_ordered,
-        constants=constants,
+        model_constants=model_constants,
+        core_constants=core_constants,
         **soil_pools,
     )