diff --git a/docs/source/api/abiotic/wind.md b/docs/source/api/abiotic/wind.md
new file mode 100644
index 000000000..aae9c0aad
--- /dev/null
+++ b/docs/source/api/abiotic/wind.md
@@ -0,0 +1,24 @@
+---
+jupytext:
+  cell_metadata_filter: -all
+  formats: md:myst
+  main_language: python
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.13.8
+kernelspec:
+  display_name: vr_python3
+  language: python
+  name: vr_python3
+---
+
+#  API for the {mod}`~virtual_rainforest.models.abiotic.wind` module
+
+```{eval-rst}
+.. automodule:: virtual_rainforest.models.abiotic.wind
+    :autosummary:
+    :members:
+    :special-members: __init__
+```
diff --git a/docs/source/index.md b/docs/source/index.md
index c9d7f9578..65ad75613 100644
--- a/docs/source/index.md
+++ b/docs/source/index.md
@@ -108,6 +108,7 @@ team.
   Abiotic Mechanistic Model <api/abiotic/abiotic_model.md>
   Abiotic Mechanistic Constants <api/abiotic/abiotic_constants.md>
   Abiotic Mechanistic Tools <api/abiotic/abiotic_tools.md>
+  Abiotic Mechanistic Wind <api/abiotic/wind.md>
   Hydrology Overview <api/hydrology.md>
   Hydrology Model <api/hydrology/hydrology_model.md>
   Hydrology Above-ground <api/hydrology/above_ground.md>
diff --git a/docs/source/refs.bib b/docs/source/refs.bib
index d6a469727..01640a5b1 100644
--- a/docs/source/refs.bib
+++ b/docs/source/refs.bib
@@ -1,4 +1,14 @@
 
+@article{yasuda_turbulent_1988,
+	title = {Turbulent diffusivity and diurnal variations in the atmospheric boundary layer.},
+	volume = {43},
+	doi = {https://doi.org/10.1007/BF00128403},
+	journal = {Boundary-layer meteorology},
+	author = {Yasuda, N.},
+	year = {1988},
+	pages = {209--221},
+}
+
 @article{henderson-sellers_new_1984,
 	title = {A new formula for latent heat of vaporization of water as a function of temperature},
 	volume = {110},
diff --git a/tests/conftest.py b/tests/conftest.py
index f2b432ccf..65278ff65 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -319,6 +319,9 @@ def dummy_climate_data(layer_roles_fixture):
         np.full((3, 3), 30),
         dims=["cell_id", "time_index"],
     )
+    data["wind_speed_ref"] = DataArray(
+        [[0, 5, 10], [0, 5, 10], [0, 5, 10]], dims=["time_index", "cell_id"]
+    )
     data["mean_annual_temperature"] = DataArray(
         np.full((3), 20),
         dims=["cell_id"],
@@ -339,6 +342,18 @@ def dummy_climate_data(layer_roles_fixture):
         np.full((3, 3), 400),
         dims=["cell_id", "time_index"],
     )
+    pressure = np.repeat(a=[96.0, np.nan], repeats=[13, 2])
+    data["atmospheric_pressure"] = DataArray(
+        np.broadcast_to(pressure, (3, 15)).T,
+        dims=["layers", "cell_id"],
+        coords={
+            "layers": np.arange(15),
+            "layer_roles": ("layers", layer_roles_fixture),
+            "cell_id": data.grid.cell_id,
+        },
+        name="atmospheric_pressure",
+    )
+
     evapotranspiration = np.repeat(a=[np.nan, 20.0, np.nan], repeats=[1, 3, 11])
     data["evapotranspiration"] = DataArray(
         np.broadcast_to(evapotranspiration, (3, 15)).T,
@@ -472,5 +487,8 @@ def dummy_climate_data(layer_roles_fixture):
         np.full((2, 3), 450),
         dims=("groundwater_layers", "cell_id"),
     )
+    data["canopy_height"] = DataArray([32, 32, 32], dims=["cell_id"])
+    data["sensible_heat_flux_topofcanopy"] = DataArray([100, 50, 10], dims=["cell_id"])
+    data["friction_velocity"] = DataArray([12, 5, 2], dims=["cell_id"])
 
     return data
diff --git a/tests/models/abiotic/test_abiotic_model.py b/tests/models/abiotic/test_abiotic_model.py
index bea3b7b2d..3d3bf21e0 100644
--- a/tests/models/abiotic/test_abiotic_model.py
+++ b/tests/models/abiotic/test_abiotic_model.py
@@ -1,88 +1,247 @@
 """Test module for abiotic.abiotic_model.py."""
 
 from contextlib import nullcontext as does_not_raise
-from logging import DEBUG, ERROR, INFO
+from logging import CRITICAL, DEBUG, ERROR, INFO
 
+import numpy as np
 import pint
 import pytest
+from xarray import DataArray
 
 from tests.conftest import log_check
 from virtual_rainforest.core.exceptions import ConfigurationError
 
 
-@pytest.mark.parametrize(
-    "raises,expected_log_entries",
-    [
-        (
-            does_not_raise(),
-            (
-                (
-                    DEBUG,
-                    "abiotic model: required var 'air_temperature_ref' checked",
-                ),
-            ),
-        ),
-    ],
-)
 def test_abiotic_model_initialization(
     caplog,
     dummy_climate_data,
-    raises,
-    expected_log_entries,
     layer_roles_fixture,
 ):
     """Test `AbioticModel` initialization."""
     from virtual_rainforest.core.base_model import BaseModel
+    from virtual_rainforest.core.constants import CoreConsts
     from virtual_rainforest.models.abiotic.abiotic_model import AbioticModel
     from virtual_rainforest.models.abiotic.constants import AbioticConsts
 
-    with raises:
-        # Initialize model
-        model = AbioticModel(
-            dummy_climate_data,
+    # Initialize model
+    model = AbioticModel(
+        dummy_climate_data,
+        pint.Quantity("1 day"),
+        soil_layers=[-0.5, -1.0],
+        canopy_layers=10,
+        constants=AbioticConsts(),
+        core_constants=CoreConsts(),
+    )
+
+    # In cases where it passes then checks that the object has the right properties
+    assert isinstance(model, BaseModel)
+    assert model.model_name == "abiotic"
+    assert str(model) == "A abiotic model instance"
+    assert repr(model) == "AbioticModel(update_interval = 1 day)"
+    assert model.layer_roles == layer_roles_fixture
+
+    # Final check that expected logging entries are produced
+    log_check(
+        caplog,
+        expected_log=(
+            (DEBUG, "abiotic model: required var 'air_temperature' checked"),
+            (DEBUG, "abiotic model: required var 'canopy_height' checked"),
+            (DEBUG, "abiotic model: required var 'layer_heights' checked"),
+            (DEBUG, "abiotic model: required var 'leaf_area_index' checked"),
+            (DEBUG, "abiotic model: required var 'atmospheric_pressure' checked"),
+            (
+                DEBUG,
+                (
+                    "abiotic model: required var 'sensible_heat_flux_topofcanopy'"
+                    " checked"
+                ),
+            ),
+            (DEBUG, "abiotic model: required var 'wind_speed_ref' checked"),
+        ),
+    )
+
+
+def test_abiotic_model_initialization_no_data(caplog, dummy_climate_data):
+    """Test `AbioticModel` 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.abiotic.abiotic_model import AbioticModel
+    from virtual_rainforest.models.abiotic.constants import AbioticConsts
+
+    with pytest.raises(ValueError):
+        # Make four cell grid
+        grid = Grid(cell_nx=4, cell_ny=1)
+        empty_data = Data(grid)
+
+        # Try and initialise model with empty data object
+        _ = AbioticModel(
+            empty_data,
             pint.Quantity("1 day"),
             soil_layers=[-0.5, -1.0],
             canopy_layers=10,
-            constants=AbioticConsts(),
+            constants=AbioticConsts,
+            core_constants=CoreConsts,
         )
 
-        # In cases where it passes then checks that the object has the right properties
-        assert isinstance(model, BaseModel)
-        assert model.model_name == "abiotic"
-        assert repr(model) == "AbioticModel(update_interval = 1 day)"
-        assert model.layer_roles == layer_roles_fixture
-
     # Final check that expected logging entries are produced
-    log_check(caplog, expected_log_entries)
+    log_check(
+        caplog,
+        expected_log=(
+            (ERROR, "abiotic model: init data missing required var 'air_temperature'"),
+            (ERROR, "abiotic model: init data missing required var 'canopy_height'"),
+            (ERROR, "abiotic model: init data missing required var 'layer_heights'"),
+            (ERROR, "abiotic model: init data missing required var 'leaf_area_index'"),
+            (
+                ERROR,
+                (
+                    "abiotic model: init data missing required var"
+                    " 'atmospheric_pressure'"
+                ),
+            ),
+            (
+                ERROR,
+                (
+                    "abiotic model: init data missing required var"
+                    " 'sensible_heat_flux_topofcanopy'"
+                ),
+            ),
+            (ERROR, "abiotic model: init data missing required var 'wind_speed_ref'"),
+            (ERROR, "abiotic model: error checking required_init_vars, see log."),
+        ),
+    )
 
 
 @pytest.mark.parametrize(
-    "cfg_string,time_interval,raises,expected_log_entries",
+    "cfg_string, time_interval, drag_coeff, raises, expected_log_entries",
     [
         pytest.param(
-            "[core]\n[core.timing]\nupdate_interval = '1 day'\n[abiotic]\n",
-            pint.Quantity("1 day"),
+            "[core]\n[core.timing]\nupdate_interval = '12 hours'\n[abiotic]\n",
+            pint.Quantity("12 hours"),
+            0.2,
             does_not_raise(),
             (
                 (INFO, "Initialised abiotic.AbioticConsts from config"),
+                (INFO, "Initialised core.CoreConsts from config"),
                 (
                     INFO,
-                    "Information required to initialise the abiotic model "
-                    "successfully extracted.",
+                    "Information required to initialise the abiotic model successfully "
+                    "extracted.",
+                ),
+                (DEBUG, "abiotic model: required var 'air_temperature' checked"),
+                (DEBUG, "abiotic model: required var 'canopy_height' checked"),
+                (DEBUG, "abiotic model: required var 'layer_heights' checked"),
+                (DEBUG, "abiotic model: required var 'leaf_area_index' checked"),
+                (
+                    DEBUG,
+                    "abiotic model: required var 'atmospheric_pressure' checked",
                 ),
                 (
                     DEBUG,
-                    "abiotic model: required var 'air_temperature_ref' checked",
+                    (
+                        "abiotic model: required var 'sensible_heat_flux_topofcanopy'"
+                        " checked"
+                    ),
                 ),
+                (DEBUG, "abiotic model: required var 'wind_speed_ref' checked"),
             ),
-            id="initialises correctly",
+            id="default_config",
         ),
+        pytest.param(
+            "[core]\n[core.timing]\nupdate_interval = '12 hours'\n"
+            "[abiotic.constants.AbioticConsts]\ndrag_coefficient = 0.05\n",
+            pint.Quantity("12 hours"),
+            0.05,
+            does_not_raise(),
+            (
+                (INFO, "Initialised abiotic.AbioticConsts from config"),
+                (INFO, "Initialised core.CoreConsts from config"),
+                (
+                    INFO,
+                    "Information required to initialise the abiotic model successfully "
+                    "extracted.",
+                ),
+                (DEBUG, "abiotic model: required var 'air_temperature' checked"),
+                (DEBUG, "abiotic model: required var 'canopy_height' checked"),
+                (DEBUG, "abiotic model: required var 'layer_heights' checked"),
+                (DEBUG, "abiotic model: required var 'leaf_area_index' checked"),
+                (
+                    DEBUG,
+                    "abiotic model: required var 'atmospheric_pressure' checked",
+                ),
+                (
+                    DEBUG,
+                    (
+                        "abiotic model: required var 'sensible_heat_flux_topofcanopy'"
+                        " checked"
+                    ),
+                ),
+                (DEBUG, "abiotic model: required var 'wind_speed_ref' checked"),
+            ),
+            id="modified_config_correct",
+        ),
+        pytest.param(
+            "[core]\n[core.timing]\nupdate_interval = '12 hours'\n"
+            "[abiotic.constants.AbioticConsts]\ndrag_coefficients = 0.05\n",
+            None,
+            None,
+            pytest.raises(ConfigurationError),
+            (
+                (ERROR, "Unknown names supplied for AbioticConsts: drag_coefficients"),
+                (INFO, "Valid names are: "),
+                (CRITICAL, "Could not initialise abiotic.AbioticConsts from config"),
+            ),
+            id="modified_config_incorrect",
+        ),
+    ],
+)
+def test_generate_abiotic_model(
+    caplog,
+    dummy_climate_data,
+    cfg_string,
+    time_interval,
+    drag_coeff,
+    raises,
+    expected_log_entries,
+):
+    """Test that the function to initialise the abiotic model behaves as expected."""
+
+    from virtual_rainforest.core.config import Config
+    from virtual_rainforest.core.registry import register_module
+    from virtual_rainforest.models.abiotic.abiotic_model import AbioticModel
+
+    # Register the module components to access constants classes
+    register_module("virtual_rainforest.models.abiotic")
+
+    # Build the config object
+    config = Config(cfg_strings=cfg_string)
+    caplog.clear()
+
+    # Check whether model is initialised (or not) as expected
+    with raises:
+        model = AbioticModel.from_config(
+            dummy_climate_data,
+            config,
+            pint.Quantity(config["core"]["timing"]["update_interval"]),
+        )
+        assert model.update_interval == time_interval
+        assert model.constants.drag_coefficient == drag_coeff
+
+    # Final check that expected logging entries are produced
+    log_check(caplog, expected_log_entries)
+
+
+@pytest.mark.parametrize(
+    "cfg_string, time_interval, raises, expected_log_entries",
+    [
         pytest.param(
             "[core]\n[core.timing]\nupdate_interval = '1 month'\n[abiotic]\n",
             pint.Quantity("1 month"),
             pytest.raises(ConfigurationError),
             (
                 (INFO, "Initialised abiotic.AbioticConsts from config"),
+                (INFO, "Initialised core.CoreConsts from config"),
                 (
                     INFO,
                     "Information required to initialise the abiotic model "
@@ -94,7 +253,7 @@ def test_abiotic_model_initialization(
         ),
     ],
 )
-def test_generate_abiotic_model(
+def test_generate_abiotic_model_bounds_error(
     caplog,
     dummy_climate_data,
     cfg_string,
@@ -128,3 +287,78 @@ def test_generate_abiotic_model(
 
     # Final check that expected logging entries are produced
     log_check(caplog, expected_log_entries)
+
+
+@pytest.mark.parametrize(
+    "cfg_string,time_interval",
+    [
+        pytest.param(
+            "[core]\n[core.timing]\nupdate_interval = '1 day'\n[abiotic]\n",
+            pint.Quantity("1 day"),
+            id="updates correctly",
+        ),
+    ],
+)
+def test_update_abiotic_model(
+    dummy_climate_data, layer_roles_fixture, cfg_string, time_interval
+):
+    """Test that update() returns expected output in data object."""
+
+    from virtual_rainforest.core.config import Config
+    from virtual_rainforest.core.registry import register_module
+    from virtual_rainforest.models.abiotic.abiotic_model import AbioticModel
+
+    # Register the module components to access constants classes
+    register_module("virtual_rainforest.models.abiotic")
+
+    # Build the config object
+    config = Config(cfg_strings=cfg_string)
+
+    # initialise model
+    model = AbioticModel.from_config(
+        dummy_climate_data,
+        config,
+        pint.Quantity(config["core"]["timing"]["update_interval"]),
+    )
+
+    model.setup()
+    model.update(time_index=0)
+
+    friction_velocity_exp = np.array(
+        [
+            [0.0, 0.818637, 1.638679],
+            [0.0, 0.81887, 1.638726],
+            [0.0, 0.820036, 1.638959],
+            [0.0, 0.821194, 1.639192],
+            [0.0, 0.822174, 1.63939],
+            [0.0, 0.822336, 1.639422],
+        ]
+    )
+    wind_speed_exp = np.array(
+        [
+            [0.55, 5.536364, 11.07365],
+            [0.54557, 5.491774, 10.984462],
+            [0.523951, 5.274152, 10.549181],
+            [0.503188, 5.065153, 10.13115],
+            [0.486188, 4.89403, 9.788873],
+            [0.483444, 4.866404, 9.733618],
+        ]
+    )
+
+    wind_above_exp = np.array([0.55, 5.536364, 11.07365])
+
+    np.testing.assert_allclose(
+        model.data["wind_speed_above_canopy"], wind_above_exp, rtol=1e-3, atol=1e-3
+    )
+    np.testing.assert_allclose(
+        model.data["friction_velocity"],
+        DataArray(np.concatenate((friction_velocity_exp, np.full((9, 3), np.nan)))),
+        rtol=1e-3,
+        atol=1e-3,
+    )
+    np.testing.assert_allclose(
+        model.data["wind_speed_canopy"],
+        DataArray(np.concatenate((wind_speed_exp, np.full((9, 3), np.nan)))),
+        rtol=1e-3,
+        atol=1e-3,
+    )
diff --git a/tests/models/abiotic/test_abiotic_tools.py b/tests/models/abiotic/test_abiotic_tools.py
index e6008099a..f45d03bab 100644
--- a/tests/models/abiotic/test_abiotic_tools.py
+++ b/tests/models/abiotic/test_abiotic_tools.py
@@ -29,6 +29,8 @@ def test_calculate_molar_density_air():
                 [39.681006, 39.681006, 39.681006],
             ]
         ),
+        rtol=1e-5,
+        atol=1e-5,
     )
 
 
@@ -47,7 +49,7 @@ def test_calculate_specific_heat_air():
     )
 
     exp_result = np.array(
-        [[29.2075, 29.207, 29.207], [29.202, 29.202, 29.202], [29.2, 29.2, 29.2]]
+        [[29.2075, 29.207, 29.207], [29.202, 29.202, 29.202], [29.2, 29.2, 29.2]],
     )
 
     np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3)
@@ -74,4 +76,9 @@ def test_calculate_latent_heat_vaporisation():
         ]
     )
 
-    np.testing.assert_allclose(result, exp_result)
+    np.testing.assert_allclose(
+        result,
+        exp_result,
+        rtol=1e-5,
+        atol=1e-5,
+    )
diff --git a/tests/models/abiotic/test_wind.py b/tests/models/abiotic/test_wind.py
new file mode 100644
index 000000000..e356f02f9
--- /dev/null
+++ b/tests/models/abiotic/test_wind.py
@@ -0,0 +1,381 @@
+"""Test module for wind.py."""
+
+
+import numpy as np
+
+from virtual_rainforest.core.constants import CoreConsts
+from virtual_rainforest.models.abiotic.constants import AbioticConsts
+
+
+def test_calculate_zero_plane_displacement(dummy_climate_data):
+    """Test if calculated correctly and set to zero without vegetation."""
+
+    from virtual_rainforest.models.abiotic.wind import calculate_zero_plane_displacement
+
+    result = calculate_zero_plane_displacement(
+        canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(),
+        leaf_area_index=np.array([0, 3, 7]),
+        zero_plane_scaling_parameter=7.5,
+    )
+
+    np.testing.assert_allclose(result, np.array([0.0, 25.312559, 27.58673]))
+
+
+def test_calculate_roughness_length_momentum(dummy_climate_data):
+    """Test roughness length governing momentum transfer."""
+
+    from virtual_rainforest.models.abiotic.wind import (
+        calculate_roughness_length_momentum,
+    )
+
+    result = calculate_roughness_length_momentum(
+        canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(),
+        leaf_area_index=np.array([0, 3, 7]),
+        zero_plane_displacement=np.array([0.0, 25.312559, 27.58673]),
+        substrate_surface_drag_coefficient=0.003,
+        roughness_element_drag_coefficient=0.3,
+        roughness_sublayer_depth_parameter=0.193,
+        max_ratio_wind_to_friction_velocity=0.3,
+        min_roughness_length=0.05,
+        von_karman_constant=CoreConsts.von_karmans_constant,
+    )
+
+    np.testing.assert_allclose(
+        result, np.array([0.017, 1.4533, 0.9591]), rtol=1e-3, atol=1e-3
+    )
+
+
+def test_calculate_diabatic_correction_above(dummy_climate_data):
+    """Test diabatic correction factors for heat and momentum."""
+
+    from virtual_rainforest.models.abiotic.wind import (
+        calculate_diabatic_correction_above,
+    )
+
+    air_temperature = (
+        dummy_climate_data.data["air_temperature"]
+        .where(dummy_climate_data.data["air_temperature"].layer_roles != "soil")
+        .dropna(dim="layers")
+    )
+    result = calculate_diabatic_correction_above(
+        molar_density_air=np.repeat(28.96, 3),
+        specific_heat_air=np.repeat(1, 3),
+        temperature=air_temperature.to_numpy(),
+        sensible_heat_flux=(
+            dummy_climate_data.data["sensible_heat_flux_topofcanopy"].to_numpy()
+        ),
+        friction_velocity=dummy_climate_data.data["friction_velocity"].to_numpy(),
+        wind_heights=np.array([1, 15, 50]),
+        zero_plane_displacement=np.array([0.0, 25.312559, 27.58673]),
+        celsius_to_kelvin=CoreConsts.zero_Celsius,
+        von_karmans_constant=CoreConsts.von_karmans_constant,
+        yasuda_stability_parameter1=AbioticConsts.yasuda_stability_parameter1,
+        yasuda_stability_parameter2=AbioticConsts.yasuda_stability_parameter2,
+        yasuda_stability_parameter3=AbioticConsts.yasuda_stability_parameter3,
+        diabatic_heat_momentum_ratio=AbioticConsts.diabatic_heat_momentum_ratio,
+    )
+
+    exp_result_h = np.array(
+        [
+            [0.003044, -0.036571, 0.042159],
+            [0.003046, -0.03659, 0.042182],
+            [0.003056, -0.036704, 0.042322],
+            [0.003073, -0.036902, 0.042564],
+            [0.00312, -0.037455, 0.043242],
+            [0.003191, -0.038274, 0.044247],
+        ]
+    )
+    exp_result_m = np.array(
+        [
+            [0.001827, -0.021943, 0.025296],
+            [0.001828, -0.021954, 0.025309],
+            [0.001833, -0.022023, 0.025393],
+            [0.001844, -0.022141, 0.025539],
+            [0.001872, -0.022473, 0.025945],
+            [0.001914, -0.022964, 0.026548],
+        ]
+    )
+    np.testing.assert_allclose(result["psi_h"], exp_result_h, rtol=1e-4, atol=1e-4)
+    np.testing.assert_allclose(result["psi_m"], exp_result_m, rtol=1e-4, atol=1e-4)
+
+
+def test_calculate_mean_mixing_length(dummy_climate_data):
+    """Test mixing length with and without vegetation."""
+
+    from virtual_rainforest.models.abiotic.wind import calculate_mean_mixing_length
+
+    result = calculate_mean_mixing_length(
+        canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(),
+        zero_plane_displacement=np.array([0.0, 25.312559, 27.58673]),
+        roughness_length_momentum=np.array([0.017, 1.4533, 0.9591]),
+        mixing_length_factor=AbioticConsts.mixing_length_factor,
+    )
+
+    np.testing.assert_allclose(
+        result, np.array([1.35804, 1.401984, 0.925228]), rtol=1e-4, atol=1e-4
+    )
+
+
+def test_generate_relative_turbulence_intensity(dummy_climate_data):
+    """Test relative turbulence intensity."""
+
+    from virtual_rainforest.models.abiotic.wind import (
+        generate_relative_turbulence_intensity,
+    )
+
+    layer_heights = (
+        dummy_climate_data.data["layer_heights"]
+        .where(dummy_climate_data.data["layer_heights"].layer_roles != "soil")
+        .dropna(dim="layers")
+    )
+    result_T = generate_relative_turbulence_intensity(
+        layer_heights=layer_heights.to_numpy(),
+        min_relative_turbulence_intensity=0.36,
+        max_relative_turbulence_intensity=0.9,
+        increasing_with_height=True,
+    )
+
+    exp_result_T = np.array(
+        [
+            [17.64, 17.64, 17.64],
+            [16.56, 16.56, 16.56],
+            [11.16, 11.16, 11.166],
+            [5.76, 5.76, 5.76],
+            [1.17, 1.17, 1.17],
+            [0.414, 0.414, 0.414],
+        ]
+    )
+    result_F = generate_relative_turbulence_intensity(
+        layer_heights=layer_heights.to_numpy(),
+        min_relative_turbulence_intensity=0.36,
+        max_relative_turbulence_intensity=0.9,
+        increasing_with_height=False,
+    )
+
+    exp_result_F = np.array(
+        [
+            [-16.38, -16.38, -16.38],
+            [-15.3, -15.3, -15.3],
+            [-9.9, -9.9, -9.9],
+            [-4.5, -4.5, -4.5],
+            [0.09, 0.09, 0.09],
+            [0.846, 0.846, 0.846],
+        ]
+    )
+    np.testing.assert_allclose(result_T, exp_result_T, rtol=1e-3, atol=1e-3)
+    np.testing.assert_allclose(result_F, exp_result_F, rtol=1e-3, atol=1e-3)
+
+
+def test_calculate_wind_attenuation_coefficient(dummy_climate_data):
+    """Test wind attenuation coefficient with and without vegetation."""
+
+    from virtual_rainforest.models.abiotic.wind import (
+        calculate_wind_attenuation_coefficient,
+    )
+
+    leaf_area_index = (
+        dummy_climate_data.data["leaf_area_index"]
+        .where(dummy_climate_data.data["leaf_area_index"].layer_roles != "soil")
+        .dropna(dim="layers")
+    )
+    relative_turbulence_intensity = np.array(
+        [
+            [17.64, 17.64, 17.64],
+            [16.56, 16.56, 16.56],
+            [11.16, 11.16, 11.166],
+            [5.76, 5.76, 5.76],
+            [1.17, 1.17, 1.17],
+            [0.414, 0.414, 0.414],
+        ]
+    )
+    result = calculate_wind_attenuation_coefficient(
+        canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(),
+        leaf_area_index=leaf_area_index,
+        mean_mixing_length=np.array([1.35804, 1.401984, 0.925228]),
+        drag_coefficient=AbioticConsts.drag_coefficient,
+        relative_turbulence_intensity=relative_turbulence_intensity,
+    )
+
+    exp_result = np.array(
+        [
+            [0.133579, 0.129392, 0.196066],
+            [0.142291, 0.137831, 0.208853],
+            [0.211141, 0.204523, 0.309744],
+        ]
+    )
+    np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3)
+
+
+def test_wind_log_profile(dummy_climate_data):
+    """Test log wind profile."""
+
+    from virtual_rainforest.models.abiotic.wind import wind_log_profile
+
+    layer_heights = (
+        dummy_climate_data.data["layer_heights"]
+        .where(dummy_climate_data.data["layer_heights"].layer_roles != "soil")
+        .dropna(dim="layers")
+    )
+    diab_correction = np.array([0.003044, -0.036571, 0.042159])
+    result = wind_log_profile(
+        height=layer_heights.to_numpy(),
+        zeroplane_displacement=np.array([0.0, 25.312559, 27.58673]),
+        roughness_length_momentum=np.array([0.017, 1.4533, 0.9591]),
+        diabatic_correction_momentum=diab_correction,
+    )
+
+    exp_result = np.array(
+        [
+            [7.543322, 1.489823, 1.568535],
+            [7.478785, 1.13446, 0.964925],
+            [7.07333, 1.0, 1.0],
+            [6.3802, 1.0, 1.0],
+            [4.483127, 1.0, 1.0],
+            [1.775148, 1.0, 1.0],
+        ]
+    )
+
+    np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3)
+
+
+def test_calculate_friction_velocity(dummy_climate_data):
+    """Calculate friction velocity."""
+
+    from virtual_rainforest.models.abiotic.wind import calculate_fricition_velocity
+
+    result = calculate_fricition_velocity(
+        wind_speed_ref=(
+            dummy_climate_data.data["wind_speed_ref"].isel(time_index=0).to_numpy()
+        ),
+        reference_height=(dummy_climate_data.data["canopy_height"] + 10).to_numpy(),
+        zeroplane_displacement=np.array([0.0, 25.312559, 27.58673]),
+        roughness_length_momentum=np.array([0.017, 1.4533, 0.9591]),
+        diabatic_correction_momentum=np.array([-0.1, 0.0, 0.1]),
+        von_karmans_constant=CoreConsts.von_karmans_constant,
+    )
+    exp_result = np.array([0.0, 0.819397, 1.423534])
+    np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3)
+
+
+def test_calculate_wind_above_canopy(dummy_climate_data):
+    """Wind speed above canopy."""
+
+    from virtual_rainforest.models.abiotic.wind import calculate_wind_above_canopy
+
+    result = calculate_wind_above_canopy(
+        friction_velocity=np.array([0.0, 0.819397, 1.423534]),
+        wind_height_above=(dummy_climate_data.data["canopy_height"] + 15).to_numpy(),
+        zeroplane_displacement=np.array([0.0, 25.312559, 27.58673]),
+        roughness_length_momentum=np.array([0.017, 1.4533, 0.9591]),
+        diabatic_correction_momentum=np.array([0.003, 0.026, 0.013]),
+        von_karmans_constant=CoreConsts.von_karmans_constant,
+        min_wind_speed_above_canopy=0.55,
+    )
+
+    exp_result = np.array([0.55, 5.590124, 10.750233])
+    np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3)
+
+
+def test_calculate_wind_canopy(dummy_climate_data):
+    """Test below canopy wind profile."""
+
+    from virtual_rainforest.models.abiotic.wind import calculate_wind_canopy
+
+    layer_heights = (
+        dummy_climate_data.data["layer_heights"]
+        .where(dummy_climate_data.data["layer_heights"].layer_roles != "soil")
+        .dropna(dim="layers")
+    )
+
+    result = calculate_wind_canopy(
+        top_of_canopy_wind_speed=np.array([0.55, 5.590124, 10.750233]),
+        wind_layer_heights=layer_heights.to_numpy(),
+        canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(),
+        attenuation_coefficient=np.array([0.133579, 0.129392, 0.196066]),
+        min_windspeed_below_canopy=0.001,
+    )
+
+    exp_result = np.array(
+        [
+            [0.55, 5.590124, 10.750233],
+            [0.545427, 5.545099, 10.619302],
+            [0.523128, 5.325355, 9.988183],
+            [0.50174, 5.11432, 9.394572],
+            [0.48425, 4.941529, 8.917824],
+            [0.481428, 4.913634, 8.841656],
+        ]
+    )
+    np.testing.assert_allclose(result, exp_result, rtol=1e-3, atol=1e-3)
+
+
+def test_calculate_wind_profile(dummy_climate_data):
+    """Test full update of wind profile."""
+
+    from virtual_rainforest.models.abiotic.wind import calculate_wind_profile
+
+    wind_layer_heights = (
+        dummy_climate_data.data["layer_heights"]
+        .where(dummy_climate_data.data["layer_heights"].layer_roles != "soil")
+        .dropna(dim="layers")
+    )
+    leaf_area_index = (
+        dummy_climate_data.data["leaf_area_index"]
+        .where(dummy_climate_data.data["leaf_area_index"].layer_roles != "soil")
+        .dropna(dim="layers")
+    )
+    air_temperature = (
+        dummy_climate_data.data["air_temperature"]
+        .where(dummy_climate_data.data["air_temperature"].layer_roles != "soil")
+        .dropna(dim="layers")
+    )
+
+    wind_update = calculate_wind_profile(
+        canopy_height=dummy_climate_data.data["canopy_height"].to_numpy(),
+        wind_height_above=(dummy_climate_data.data["canopy_height"] + 15).to_numpy(),
+        wind_layer_heights=wind_layer_heights.to_numpy(),
+        leaf_area_index=leaf_area_index.to_numpy(),
+        air_temperature=air_temperature.to_numpy(),
+        atmospheric_pressure=np.array([96, 96, 96]),
+        sensible_heat_flux_topofcanopy=np.array([100, 50, 10]),
+        wind_speed_ref=(
+            dummy_climate_data.data["wind_speed_ref"].isel(time_index=0).to_numpy()
+        ),
+        wind_reference_height=(
+            dummy_climate_data.data["canopy_height"] + 10
+        ).to_numpy(),
+        abiotic_constants=AbioticConsts,
+        core_constants=CoreConsts,
+    )
+
+    friction_velocity_exp = np.array(
+        [
+            [0.0, 0.818637, 1.638679],
+            [0.0, 0.81887, 1.638726],
+            [0.0, 0.820036, 1.638959],
+            [0.0, 0.821194, 1.639192],
+            [0.0, 0.822174, 1.63939],
+            [0.0, 0.822336, 1.639422],
+        ]
+    )
+    wind_speed_exp = np.array(
+        [
+            [0.55, 5.536364, 11.07365],
+            [0.54557, 5.491774, 10.984462],
+            [0.523951, 5.274152, 10.549181],
+            [0.503188, 5.065153, 10.13115],
+            [0.486188, 4.89403, 9.788873],
+            [0.483444, 4.866404, 9.733618],
+        ]
+    )
+
+    wind_above_exp = np.array([0.55, 5.536364, 11.07365])
+
+    np.testing.assert_allclose(
+        wind_update["wind_speed_above_canopy"], wind_above_exp, rtol=1e-3, atol=1e-3
+    )
+    np.testing.assert_allclose(
+        wind_update["friction_velocity"], friction_velocity_exp, rtol=1e-3, atol=1e-3
+    )
+    np.testing.assert_allclose(
+        wind_update["wind_speed_canopy"], wind_speed_exp, rtol=1e-3, atol=1e-3
+    )
diff --git a/tests/models/animals/conftest.py b/tests/models/animals/conftest.py
index da188ca31..0b2039bf4 100644
--- a/tests/models/animals/conftest.py
+++ b/tests/models/animals/conftest.py
@@ -77,31 +77,49 @@ def plant_climate_data_instance(layer_roles_fixture):
 
 
 @pytest.fixture
-def functional_group_list_instance(shared_datadir):
+def constants_instance():
+    """Fixture for an instance of animal constants."""
+    from virtual_rainforest.models.animals.constants import AnimalConsts
+
+    return AnimalConsts()
+
+
+@pytest.fixture
+def functional_group_list_instance(shared_datadir, constants_instance):
     """Fixture for an animal functional group used in tests."""
     from virtual_rainforest.models.animals.functional_group import (
         import_functional_groups,
     )
 
     file = shared_datadir / "example_functional_group_import.csv"
-    fg_list = import_functional_groups(file)
+    fg_list = import_functional_groups(file, constants_instance)
 
     return fg_list
 
 
 @pytest.fixture
-def animal_model_instance(data_instance, functional_group_list_instance):
+def animal_model_instance(
+    data_instance, functional_group_list_instance, constants_instance
+):
     """Fixture for an animal model object used in tests."""
     from pint import Quantity
 
     from virtual_rainforest.models.animals.animal_model import AnimalModel
 
-    return AnimalModel(data_instance, Quantity("1 day"), functional_group_list_instance)
+    return AnimalModel(
+        data_instance,
+        Quantity("1 day"),
+        functional_group_list_instance,
+        constants_instance,
+    )
 
 
 @pytest.fixture
 def animal_community_instance(
-    functional_group_list_instance, animal_model_instance, plant_data_instance
+    functional_group_list_instance,
+    animal_model_instance,
+    plant_data_instance,
+    constants_instance,
 ):
     """Fixture for an animal community used in tests."""
     from virtual_rainforest.models.animals.animal_communities import AnimalCommunity
@@ -112,28 +130,31 @@ def animal_community_instance(
         community_key=4,
         neighbouring_keys=[1, 3, 5, 7],
         get_destination=animal_model_instance.get_community_by_key,
+        constants=constants_instance,
     )
 
 
 @pytest.fixture
-def herbivore_functional_group_instance(shared_datadir):
+def herbivore_functional_group_instance(shared_datadir, constants_instance):
     """Fixture for an animal functional group used in tests."""
     from virtual_rainforest.models.animals.functional_group import (
         import_functional_groups,
     )
 
     file = shared_datadir / "example_functional_group_import.csv"
-    fg_list = import_functional_groups(file)
+    fg_list = import_functional_groups(file, constants_instance)
 
     return fg_list[3]
 
 
 @pytest.fixture
-def herbivore_cohort_instance(herbivore_functional_group_instance):
+def herbivore_cohort_instance(herbivore_functional_group_instance, constants_instance):
     """Fixture for an animal cohort used in tests."""
     from virtual_rainforest.models.animals.animal_cohorts import AnimalCohort
 
-    return AnimalCohort(herbivore_functional_group_instance, 10000.0, 1, 10)
+    return AnimalCohort(
+        herbivore_functional_group_instance, 10000.0, 1, 10, constants_instance
+    )
 
 
 @pytest.fixture
@@ -145,8 +166,10 @@ def excrement_instance():
 
 
 @pytest.fixture
-def plant_instance(plant_data_instance):
+def plant_instance(plant_data_instance, constants_instance):
     """Fixture for a plant community used in tests."""
     from virtual_rainforest.models.animals.plant_resources import PlantResources
 
-    return PlantResources(data=plant_data_instance, cell_id=4)
+    return PlantResources(
+        data=plant_data_instance, cell_id=4, constants=constants_instance
+    )
diff --git a/tests/models/animals/test_animal_cohorts.py b/tests/models/animals/test_animal_cohorts.py
index 343c1091b..6fbae3c54 100644
--- a/tests/models/animals/test_animal_cohorts.py
+++ b/tests/models/animals/test_animal_cohorts.py
@@ -5,53 +5,59 @@
 
 
 @pytest.fixture
-def predator_functional_group_instance(shared_datadir):
+def predator_functional_group_instance(shared_datadir, constants_instance):
     """Fixture for an animal functional group used in tests."""
     from virtual_rainforest.models.animals.functional_group import (
         import_functional_groups,
     )
 
     file = shared_datadir / "example_functional_group_import.csv"
-    fg_list = import_functional_groups(file)
+    fg_list = import_functional_groups(file, constants_instance)
 
     return fg_list[2]
 
 
 @pytest.fixture
-def predator_cohort_instance(predator_functional_group_instance):
+def predator_cohort_instance(predator_functional_group_instance, constants_instance):
     """Fixture for an animal cohort used in tests."""
     from virtual_rainforest.models.animals.animal_cohorts import AnimalCohort
 
-    return AnimalCohort(predator_functional_group_instance, 10000.0, 1, 10)
+    return AnimalCohort(
+        predator_functional_group_instance, 10000.0, 1, 10, constants_instance
+    )
 
 
 @pytest.fixture
-def ectotherm_functional_group_instance(shared_datadir):
+def ectotherm_functional_group_instance(shared_datadir, constants_instance):
     """Fixture for an animal functional group used in tests."""
     from virtual_rainforest.models.animals.functional_group import (
         import_functional_groups,
     )
 
     file = shared_datadir / "example_functional_group_import.csv"
-    fg_list = import_functional_groups(file)
+    fg_list = import_functional_groups(file, constants_instance)
 
     return fg_list[5]
 
 
 @pytest.fixture
-def ectotherm_cohort_instance(ectotherm_functional_group_instance):
+def ectotherm_cohort_instance(ectotherm_functional_group_instance, constants_instance):
     """Fixture for an animal cohort used in tests."""
     from virtual_rainforest.models.animals.animal_cohorts import AnimalCohort
 
-    return AnimalCohort(ectotherm_functional_group_instance, 100.0, 1, 10)
+    return AnimalCohort(
+        ectotherm_functional_group_instance, 100.0, 1, 10, constants_instance
+    )
 
 
 @pytest.fixture
-def prey_cohort_instance(herbivore_functional_group_instance):
+def prey_cohort_instance(herbivore_functional_group_instance, constants_instance):
     """Fixture for an animal cohort used in tests."""
     from virtual_rainforest.models.animals.animal_cohorts import AnimalCohort
 
-    return AnimalCohort(herbivore_functional_group_instance, 100.0, 1, 10)
+    return AnimalCohort(
+        herbivore_functional_group_instance, 100.0, 1, 10, constants_instance
+    )
 
 
 @pytest.fixture
@@ -85,6 +91,7 @@ def test_invalid_animal_cohort_initialization(
         age,
         individuals,
         error_type,
+        constants_instance,
     ):
         """Test for invalid inputs during AnimalCohort initialization."""
         from virtual_rainforest.models.animals.animal_cohorts import AnimalCohort
@@ -95,6 +102,7 @@ def test_invalid_animal_cohort_initialization(
                 mass,
                 age,
                 individuals,
+                constants_instance,
             )
 
     @pytest.mark.parametrize(
@@ -346,38 +354,46 @@ def test_eat(self, herbivore_cohort_instance, mocker):
         )
         assert herbivore_cohort_instance.mass_current == 0
 
-    def test_is_below_mass_threshold(self, herbivore_cohort_instance):
+    def test_is_below_mass_threshold(
+        self, herbivore_cohort_instance, constants_instance
+    ):
         """Test the can_reproduce method of AnimalCohort."""
-        from virtual_rainforest.models.animals.constants import BIRTH_MASS_THRESHOLD
 
         # TODO: test other mass thresholds
         # 1. Test when stored_energy is exactly equal to the threshold
         herbivore_cohort_instance.mass_current = (
-            herbivore_cohort_instance.functional_group.adult_mass * BIRTH_MASS_THRESHOLD
+            herbivore_cohort_instance.functional_group.adult_mass
+            * constants_instance.birth_mass_threshold
         )
         assert not herbivore_cohort_instance.is_below_mass_threshold(
-            BIRTH_MASS_THRESHOLD
+            constants_instance.birth_mass_threshold
         )
 
         # 2. Test when stored_energy is just below the threshold
         herbivore_cohort_instance.mass_current = (
-            herbivore_cohort_instance.functional_group.adult_mass * BIRTH_MASS_THRESHOLD
+            herbivore_cohort_instance.functional_group.adult_mass
+            * constants_instance.birth_mass_threshold
             - 0.01
         )
-        assert herbivore_cohort_instance.is_below_mass_threshold(BIRTH_MASS_THRESHOLD)
+        assert herbivore_cohort_instance.is_below_mass_threshold(
+            constants_instance.birth_mass_threshold
+        )
 
         # 3. Test when stored_energy is above the threshold
         herbivore_cohort_instance.mass_current = (
-            herbivore_cohort_instance.functional_group.adult_mass * BIRTH_MASS_THRESHOLD
+            herbivore_cohort_instance.functional_group.adult_mass
+            * constants_instance.birth_mass_threshold
             + 0.01
         )
         assert not herbivore_cohort_instance.is_below_mass_threshold(
-            BIRTH_MASS_THRESHOLD
+            constants_instance.birth_mass_threshold
         )
 
         # 4. Test with stored_energy set to 0
         herbivore_cohort_instance.mass_current = 0.0
-        assert herbivore_cohort_instance.is_below_mass_threshold(BIRTH_MASS_THRESHOLD)
+        assert herbivore_cohort_instance.is_below_mass_threshold(
+            constants_instance.birth_mass_threshold
+        )
 
     @pytest.mark.parametrize(
         "initial_individuals, number_days, mortality_prob",
diff --git a/tests/models/animals/test_animal_communities.py b/tests/models/animals/test_animal_communities.py
index 209f9e28f..7a033afd4 100644
--- a/tests/models/animals/test_animal_communities.py
+++ b/tests/models/animals/test_animal_communities.py
@@ -5,7 +5,10 @@
 
 @pytest.fixture
 def animal_community_destination_instance(
-    functional_group_list_instance, animal_model_instance, plant_data_instance
+    functional_group_list_instance,
+    animal_model_instance,
+    plant_data_instance,
+    constants_instance,
 ):
     """Fixture for an animal community used in tests."""
     from virtual_rainforest.models.animals.animal_communities import AnimalCommunity
@@ -16,29 +19,34 @@ def animal_community_destination_instance(
         community_key=4,
         neighbouring_keys=[1, 3, 5, 7],
         get_destination=animal_model_instance.get_community_by_key,
+        constants=constants_instance,
     )
 
 
 @pytest.fixture
-def functional_group_instance(shared_datadir):
+def functional_group_instance(shared_datadir, constants_instance):
     """Fixture for an animal functional group used in tests."""
     from virtual_rainforest.models.animals.functional_group import (
         import_functional_groups,
     )
 
     file = shared_datadir / "example_functional_group_import.csv"
-    fg_list = import_functional_groups(file)
+    fg_list = import_functional_groups(file, constants_instance)
 
     return fg_list[3]
 
 
 @pytest.fixture
-def animal_cohort_instance(functional_group_instance):
+def animal_cohort_instance(functional_group_instance, constants_instance):
     """Fixture for an animal cohort used in tests."""
     from virtual_rainforest.models.animals.animal_cohorts import AnimalCohort
 
     return AnimalCohort(
-        functional_group_instance, functional_group_instance.adult_mass, 1.0, 10
+        functional_group_instance,
+        functional_group_instance.adult_mass,
+        1.0,
+        10,
+        constants_instance,
     )
 
 
@@ -164,9 +172,10 @@ def test_die_cohort(self, animal_cohort_instance, animal_community_instance):
             not in animal_community_instance.animal_cohorts["herbivorous_mammal"]
         )
 
-    def test_birth(self, animal_community_instance, animal_cohort_instance):
+    def test_birth(
+        self, animal_community_instance, animal_cohort_instance, constants_instance
+    ):
         """Test the birth method in AnimalCommunity."""
-        from virtual_rainforest.models.animals.constants import BIRTH_MASS_THRESHOLD
 
         # Setup initial conditions
         parent_cohort_name = animal_cohort_instance.name
@@ -179,7 +188,8 @@ def test_birth(self, animal_community_instance, animal_cohort_instance):
 
         # Set the reproductive mass of the parent cohort to ensure it can reproduce
         required_mass_for_birth = (
-            animal_cohort_instance.functional_group.adult_mass * BIRTH_MASS_THRESHOLD
+            animal_cohort_instance.functional_group.adult_mass
+            * constants_instance.birth_mass_threshold
             - animal_cohort_instance.functional_group.adult_mass
         )
 
@@ -198,13 +208,11 @@ def test_birth(self, animal_community_instance, animal_cohort_instance):
         # 2. Check that the reproductive mass of the parent cohort is reduced to 0
         assert animal_cohort_instance.reproductive_mass == 0
 
-    def test_birth_community(self, animal_community_instance):
+    def test_birth_community(self, animal_community_instance, constants_instance):
         """Test the thresholding behavior of birth_community."""
 
         from itertools import chain
 
-        from virtual_rainforest.models.animals.constants import BIRTH_MASS_THRESHOLD
-
         # Preparation: populate the community
         animal_community_instance.populate_community()
 
@@ -216,7 +224,8 @@ def test_birth_community(self, animal_community_instance):
 
         # Set mass to just below the threshold
         threshold_mass = (
-            initial_cohort.functional_group.adult_mass * BIRTH_MASS_THRESHOLD
+            initial_cohort.functional_group.adult_mass
+            * constants_instance.birth_mass_threshold
             - initial_cohort.functional_group.adult_mass
         )
 
diff --git a/tests/models/animals/test_animal_model.py b/tests/models/animals/test_animal_model.py
index 9b5bb2469..d0a18ae63 100644
--- a/tests/models/animals/test_animal_model.py
+++ b/tests/models/animals/test_animal_model.py
@@ -11,7 +11,10 @@
 
 
 def test_animal_model_initialization(
-    caplog, plant_climate_data_instance, functional_group_list_instance
+    caplog,
+    plant_climate_data_instance,
+    functional_group_list_instance,
+    constants_instance,
 ):
     """Test `AnimalModel` initialization."""
     from virtual_rainforest.core.base_model import BaseModel
@@ -22,6 +25,7 @@ def test_animal_model_initialization(
         plant_climate_data_instance,
         pint.Quantity("1 week"),
         functional_group_list_instance,
+        constants=constants_instance,
     )
 
     # In cases where it passes then checks that the object has the right properties
@@ -108,6 +112,7 @@ def test_generate_animal_model(
     time_interval,
     raises,
     expected_log_entries,
+    constants_instance,
 ):
     """Test that the function to initialise the animal model behaves as expected."""
     from virtual_rainforest.models.animals.animal_model import AnimalModel
@@ -115,9 +120,9 @@ def test_generate_animal_model(
     # Check whether model is initialised (or not) as expected
     with raises:
         model = AnimalModel.from_config(
-            plant_climate_data_instance,
-            config,
-            pint.Quantity(config["core"]["timing"]["update_interval"]),
+            data=plant_climate_data_instance,
+            config=config,
+            update_interval=pint.Quantity(config["core"]["timing"]["update_interval"]),
         )
         assert model.update_interval == time_interval
         # Run the update step (once this does something should check output)
@@ -152,7 +157,7 @@ def test_get_community_by_key(animal_model_instance):
 
 
 def test_update_method_sequence(
-    plant_climate_data_instance, functional_group_list_instance
+    plant_climate_data_instance, functional_group_list_instance, constants_instance
 ):
     """Test update to ensure it runs the community methods in order.
 
@@ -166,6 +171,7 @@ def test_update_method_sequence(
         plant_climate_data_instance,
         pint.Quantity("1 week"),
         functional_group_list_instance,
+        constants=constants_instance,
     )
 
     # Mock all the methods that are supposed to be called by update
diff --git a/tests/models/animals/test_functional_group.py b/tests/models/animals/test_functional_group.py
index c89bb1794..79cf752ad 100644
--- a/tests/models/animals/test_functional_group.py
+++ b/tests/models/animals/test_functional_group.py
@@ -99,10 +99,17 @@ def test_initialization(
             MetabolicType,
             TaxaType,
         )
+        from virtual_rainforest.models.animals.constants import AnimalConsts
         from virtual_rainforest.models.animals.functional_group import FunctionalGroup
 
         func_group = FunctionalGroup(
-            name, taxa, diet, metabolic_type, birth_mass, adult_mass
+            name,
+            taxa,
+            diet,
+            metabolic_type,
+            birth_mass,
+            adult_mass,
+            constants=AnimalConsts(),
         )
         assert func_group.name == name
         assert func_group.taxa == TaxaType(taxa)
@@ -133,13 +140,14 @@ def test_import_functional_groups(
         MetabolicType,
         TaxaType,
     )
+    from virtual_rainforest.models.animals.constants import AnimalConsts
     from virtual_rainforest.models.animals.functional_group import (
         FunctionalGroup,
         import_functional_groups,
     )
 
     file = shared_datadir / "example_functional_group_import.csv"
-    fg_list = import_functional_groups(file)
+    fg_list = import_functional_groups(file, constants=AnimalConsts())
     assert len(fg_list) == 6
     assert isinstance(fg_list[index], FunctionalGroup)
     assert fg_list[index].name == name
diff --git a/virtual_rainforest/core/constants.py b/virtual_rainforest/core/constants.py
index 440455ffa..ce992056c 100644
--- a/virtual_rainforest/core/constants.py
+++ b/virtual_rainforest/core/constants.py
@@ -18,12 +18,12 @@
 class CoreConsts(ConstantsDataclass):
     """Core constants for use across the Virtual Rainforest modules."""
 
-    zero_Celsius: ClassVar[float] = constants.zero_Celsius
-    """Conversion constant from Kelvin to Celsius (°)."""
-
     placeholder: float = 123.4
     """A placeholder configurable constant."""
 
+    zero_Celsius: ClassVar[float] = constants.zero_Celsius
+    """Conversion constant from Kelvin to Celsius (°)."""
+
     standard_pressure: float = constants.atmosphere / 1000
     """Standard atmospheric pressure, [kPa]"""
 
diff --git a/virtual_rainforest/models/abiotic/__init__.py b/virtual_rainforest/models/abiotic/__init__.py
index f5e12b00d..b8795fc16 100644
--- a/virtual_rainforest/models/abiotic/__init__.py
+++ b/virtual_rainforest/models/abiotic/__init__.py
@@ -15,6 +15,11 @@
 * The :mod:`~virtual_rainforest.models.abiotic.abiotic_tools` submodule contains a set
   of general functions that are shared across submodules in the
   :mod:`~virtual_rainforest.models.abiotic` model.
+
+* The :mod:`~virtual_rainforest.models.abiotic.wind` submodule calculates the
+  above- and within-canopy wind profiles for the Virtual Rainforest. These profiles will
+  determine the exchange of heat, water, and :math:`CO_{2}` between soil and atmosphere
+  below the canopy as well as the exchange with the atmsophere above the canopy.
 """  # noqa: D205, D415
 
 from virtual_rainforest.models.abiotic.abiotic_model import AbioticModel  # noqa: F401
diff --git a/virtual_rainforest/models/abiotic/abiotic_model.py b/virtual_rainforest/models/abiotic/abiotic_model.py
index 347ab46e4..8e34dfbc7 100644
--- a/virtual_rainforest/models/abiotic/abiotic_model.py
+++ b/virtual_rainforest/models/abiotic/abiotic_model.py
@@ -12,20 +12,26 @@ class as a child of the :class:`~virtual_rainforest.core.base_model.BaseModel` c
 config dictionary to the required types they are caught and then logged, and at the end
 of the unpacking an error is thrown. This error should be caught and handled by
 downstream functions so that all model configuration failures can be reported as one.
+
+TODO temperatures in Kelvin
 """  # noqa: D205, D415
 
 from __future__ import annotations
 
 from typing import Any
 
+import numpy as np
 from pint import Quantity
+from xarray import DataArray
 
 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.logger import LOGGER
 from virtual_rainforest.core.utils import set_layer_roles
+from virtual_rainforest.models.abiotic import wind
 from virtual_rainforest.models.abiotic.constants import AbioticConsts
 
 
@@ -46,7 +52,15 @@ class AbioticModel(BaseModel):
     """Shortest time scale that abiotic model can sensibly capture."""
     upper_bound_on_time_scale = "1 day"
     """Longest time scale that abiotic model can sensibly capture."""
-    required_init_vars = (("air_temperature_ref", ("spatial",)),)
+    required_init_vars = (
+        ("air_temperature", ("spatial",)),
+        ("canopy_height", ("spatial",)),
+        ("layer_heights", ("spatial",)),
+        ("leaf_area_index", ("spatial",)),
+        ("atmospheric_pressure", ("spatial",)),
+        ("sensible_heat_flux_topofcanopy", ("spatial",)),
+        ("wind_speed_ref", ("spatial",)),
+    )
     """The required variables and axes for the abiotic model"""
     vars_updated = ()
     """Variables updated by the abiotic model"""
@@ -58,6 +72,7 @@ def __init__(
         soil_layers: list[float],
         canopy_layers: int,
         constants: AbioticConsts,
+        core_constants: CoreConsts,
         **kwargs: Any,
     ):
         super().__init__(data, update_interval, **kwargs)
@@ -72,7 +87,9 @@ def __init__(
         self.update_interval
         """The time interval between model updates."""
         self.constants = constants
-        """Set of constants for the abiotic model"""
+        """Set of constants for the abiotic model."""
+        self.core_constants = core_constants
+        """Set of universal constants that are used across all models."""
 
     @classmethod
     def from_config(
@@ -96,6 +113,7 @@ def from_config(
 
         # Load in the relevant constants
         constants = load_constants(config, "abiotic", "AbioticConsts")
+        core_constants = load_constants(config, "core", "CoreConsts")
 
         LOGGER.info(
             "Information required to initialise the abiotic model successfully "
@@ -107,6 +125,7 @@ def from_config(
             soil_layers,
             canopy_layers,
             constants,
+            core_constants,
         )
 
     def setup(self) -> None:
@@ -122,5 +141,62 @@ def update(self, time_index: int, **kwargs: Any) -> None:
             time_index: The index of the current time step in the data object.
         """
 
+        wind_update_inputs: dict[str, DataArray] = {}
+
+        for var in ["layer_heights", "leaf_area_index", "air_temperature"]:
+            selection = (
+                self.data[var]
+                .where(self.data[var].layer_roles != "soil")
+                .dropna(dim="layers")
+            )
+            wind_update_inputs[var] = selection
+
+        wind_update = wind.calculate_wind_profile(
+            canopy_height=self.data["canopy_height"].to_numpy(),
+            wind_height_above=(self.data["canopy_height"] + 15).to_numpy(),  # TODO
+            wind_layer_heights=wind_update_inputs["layer_heights"].to_numpy(),
+            leaf_area_index=wind_update_inputs["leaf_area_index"].to_numpy(),
+            air_temperature=wind_update_inputs["air_temperature"].to_numpy(),
+            atmospheric_pressure=self.data["atmospheric_pressure"].to_numpy()[0],
+            sensible_heat_flux_topofcanopy=(
+                self.data["sensible_heat_flux_topofcanopy"].to_numpy()
+            ),
+            wind_speed_ref=(
+                self.data["wind_speed_ref"].isel(time_index=time_index).to_numpy()
+            ),
+            wind_reference_height=(self.data["canopy_height"] + 10).to_numpy(),
+            abiotic_constants=self.constants,
+            core_constants=self.core_constants,
+        )
+
+        wind_output = {}
+        wind_speed_above_canopy = DataArray(
+            wind_update["wind_speed_above_canopy"],
+            dims="cell_id",
+            coords={"cell_id": self.data.grid.cell_id},
+        )
+        wind_output["wind_speed_above_canopy"] = wind_speed_above_canopy
+
+        for var in ["wind_speed_canopy", "friction_velocity"]:
+            var_out = DataArray(
+                np.concatenate(
+                    (
+                        wind_update[var],
+                        np.full(
+                            (
+                                len(self.layer_roles) - len(wind_update[var]),
+                                self.data.grid.n_cells,
+                            ),
+                            np.nan,
+                        ),
+                    ),
+                ),
+                dims=self.data["layer_heights"].dims,
+                coords=self.data["layer_heights"].coords,
+            )
+            wind_output[var] = var_out
+
+        self.data.add_from_dict(output_dict=wind_output)
+
     def cleanup(self) -> None:
         """Placeholder function for abiotic model cleanup."""
diff --git a/virtual_rainforest/models/abiotic/abiotic_tools.py b/virtual_rainforest/models/abiotic/abiotic_tools.py
index 758fac76d..62aae0a2a 100644
--- a/virtual_rainforest/models/abiotic/abiotic_tools.py
+++ b/virtual_rainforest/models/abiotic/abiotic_tools.py
@@ -2,6 +2,7 @@
 are shared across submodules in the :mod:`~virtual_rainforest.models.abiotic` model.
 
 TODO cross-check with pyrealm for duplication/ different implementation
+TODO change temperatures to Kelvin
 """  # noqa: D205, D415
 
 import numpy as np
diff --git a/virtual_rainforest/models/abiotic/constants.py b/virtual_rainforest/models/abiotic/constants.py
index ed67860e7..a06bc3ce8 100644
--- a/virtual_rainforest/models/abiotic/constants.py
+++ b/virtual_rainforest/models/abiotic/constants.py
@@ -35,3 +35,131 @@ class AbioticConsts(ConstantsDataclass):
     Implementation after :cite:t:`maclean_microclimc_2021`, value is taken from
     :cite:t:`henderson-sellers_new_1984`.
     """
+
+    zero_plane_scaling_parameter: float = 7.5
+    """Control parameter for scaling zero displacement to height, dimensionless.
+
+    Implementation after :cite:t:`maclean_microclimc_2021`, value is taken from
+    :cite:t:`raupach_simplified_1994`."""
+
+    substrate_surface_drag_coefficient: float = 0.003
+    """Substrate-surface drag coefficient, dimensionless.
+
+    The substrate-surface drag coefficient represents the resistance encountered by an
+    object moving on or through a surface and varies based on the nature of the surface
+    and the object's properties. Here, it affects how wind speed is altered by a surface
+    . Implementation and value from :cite:t:`maclean_microclimc_2021`."""
+
+    roughness_element_drag_coefficient: float = 0.3
+    """Roughness-element drag coefficient, dimensionless.
+
+    The roughness-element drag coefficient refers to the dimensionless coefficient used
+    to quantify the drag force exerted by individual roughness elements (such as
+    buildings, trees, or surface irregularities) on airflow, influencing the overall
+    aerodynamic characteristics of a surface within the atmospheric boundary layer.
+    Implementation and value from :cite:t:`maclean_microclimc_2021`."""
+
+    roughness_sublayer_depth_parameter: float = 0.193
+    """Parameter characterizes the roughness sublayer depth.
+
+    The roughness sublayer depth refers to the layer near the surface where the effects
+    of surface roughness significantly influence airflow, turbulence, and momentum
+    transfer, typically extending up to about 10% of the height of the roughness
+    elements or features on the surface. This layer is characterized by intense
+    turbulence and rapid velocity changes due to surface irregularities.
+    Implentation and value taken from :cite:p:`maclean_microclimc_2021`."""
+
+    max_ratio_wind_to_friction_velocity: float = 0.3
+    """Maximum ratio of wind velocity to friction velocity, dimensionless.
+
+    Implementation and value from :cite:t:`maclean_microclimc_2021`."""
+
+    drag_coefficient: float = 0.2
+    """Drag coefficient, dimensionless.
+
+    The drag coefficient is a dimensionless quantity that characterizes the drag or
+    resistance experienced by an object moving through a fluid (here the atmosphere) and
+    is defined as the ratio of the drag force on the object to the dynamic pressure of
+    the fluid flow and the reference area of the object.
+    Implementation and value from :cite:t:`maclean_microclimc_2021`."""
+
+    relative_turbulence_intensity: float = 0.5
+    """Relative turbulence intensity, dimensionless.
+
+    The relative turbulence intensity is a proportionality factor that relates the mean
+    eddy velocity is assumed to the local wind speed below the canopy. Implementation
+    and value from :cite:t:`maclean_microclimc_2021`."""
+
+    diabatic_correction_factor_below: float = 1
+    """Diabatic correction factor below canopy, dimensionless.
+
+    The diabatic correction factor is a scaling adjustment used to compensate for the
+    effects of vertical heat transfer or thermal non-adiabaticity on atmospheric
+    variables or processes, particularly when estimating or interpreting measurements
+    across different heights or conditions. This factor is used to adjust wind profiles
+    below the canopy. Implementation and value from :cite:t:`maclean_microclimc_2021`.
+    """
+
+    mixing_length_factor: float = 0.32
+    """Factor in calculation of mixing length, dimensionless.
+
+    Implementation and value from :cite:t:`maclean_microclimc_2021`."""
+
+    min_relative_turbulence_intensity: float = 0.36
+    """Minimum relative turbulence intensity, dimensionless.
+
+    See :attr:`relative_turbulence_intensity`.
+    The default value is taken from Shaw et al (1974) Agricultural Meteorology, 13:
+    419-425. TODO this is not representative of a rainforest environment and needs to be
+    adjusted.
+    """
+
+    max_relative_turbulence_intensity: float = 0.9
+    """Maximum relative turbulence intensity, dimensionless.
+
+    See :attr:`relative_turbulence_intensity`.
+    The default value from Shaw et al (1974) Agricultural Meteorology, 13: 419-425.
+    TODO this is not representative of a rainforest environment and needs to be
+    adjusted."""
+
+    min_wind_speed_above_canopy: float = 0.55
+    """Minimum wind speed above the canopy, [m/s].
+
+    Implementation and value from :cite:t:`maclean_microclimc_2021`."""
+
+    min_windspeed_below_canopy: float = 0.001
+    """Minimum wind speed below the canopy or in absence of vegetation, [m/s]."""
+
+    min_roughness_length: float = 0.05
+    """Minimum roughness length, [m].
+
+    The minimum roughness length represents the lowest height at which the surface
+    roughness significantly affects the wind flow over a particular terrain or surface.
+    Implementation and value from :cite:t:`maclean_microclimc_2021`."""
+
+    yasuda_stability_parameter1: float = 6
+    """Parameter to approximate diabatic correction factors for heat and momentum.
+
+    Dimenionless parameter, implementation after :cite:t:`maclean_microclimc_2021` and
+    values taken from :cite:t:`yasuda_turbulent_1988`."""
+
+    yasuda_stability_parameter2: float = 2
+    """Parameter to approximate diabatic correction factors for heat and momentum.
+
+    Dimenionless parameter, implementation after :cite:t:`maclean_microclimc_2021` and
+    values taken from :cite:t:`yasuda_turbulent_1988`."""
+
+    yasuda_stability_parameter3: float = 16
+    """Parameter to approximate diabatic correction factors for heat and momentum.
+
+    Dimenionless parameter, implementation after :cite:t:`maclean_microclimc_2021` and
+    values taken from :cite:t:`yasuda_turbulent_1988`."""
+
+    diabatic_heat_momentum_ratio: float = 0.6
+    """Factor that relates diabatic correction factors for heat and momentum.
+
+    Dimenionless parameter, implementation after :cite:t:`maclean_microclimc_2021` and
+    values taken from :cite:t:`yasuda_turbulent_1988`."""
+
+    turbulence_sign: bool = True
+    """Flag indicating if turbulence increases or decreases with height."""
diff --git a/virtual_rainforest/models/abiotic/wind.py b/virtual_rainforest/models/abiotic/wind.py
new file mode 100644
index 000000000..42106f8ff
--- /dev/null
+++ b/virtual_rainforest/models/abiotic/wind.py
@@ -0,0 +1,653 @@
+r"""The wind module calculates the above- and within-canopy wind profile for the
+Virtual Rainforest. The wind profile determines the exchange of heat, water, and
+:math:`CO_{2}` between soil and atmosphere below the canopy as well as the exchange with
+the atmosphere above the canopy.
+
+TODO: add sanity checks, errors and logging
+TODO replace leaf area index by plant area index when we have more info about vertical
+distribution of leaf and woody parts
+TODO change temperatures to Kelvin
+"""  # noqa: D205, D415
+
+from typing import Union
+
+import numpy as np
+from numpy.typing import NDArray
+
+from virtual_rainforest.core.constants import CoreConsts
+from virtual_rainforest.models.abiotic.abiotic_tools import (
+    calculate_molar_density_air,
+    calculate_specific_heat_air,
+)
+from virtual_rainforest.models.abiotic.constants import AbioticConsts
+
+
+def calculate_zero_plane_displacement(
+    canopy_height: NDArray[np.float32],
+    leaf_area_index: NDArray[np.float32],
+    zero_plane_scaling_parameter: float,
+) -> NDArray[np.float32]:
+    """Calculate zero plane displacement height, [m].
+
+    The zero plane displacement height of a vegetated surface is the height at which the
+    wind speed would go to zero if the logarithmic wind profile was maintained from the
+    outer flow all the way down to the surface (that is, in the absence of the
+    vegetation when the value is set to zero). Implementation after
+    :cite:t:`maclean_microclimc_2021`.
+
+    Args:
+        canopy_height: Canopy height, [m]
+        leaf_area_index: Total leaf area index, [m m-1]
+        zero_plane_scaling_parameter: Control parameter for scaling d/h, dimensionless
+            :cite:p:`raupach_simplified_1994`
+
+    Returns:
+        zero place displacement height, [m]
+    """
+
+    # Select grid cells where vegetation is present
+    displacement = np.where(leaf_area_index > 0, leaf_area_index, np.nan)
+
+    # Calculate zero displacement height
+    scale_displacement = np.sqrt(zero_plane_scaling_parameter * displacement)
+    zero_place_displacement = (
+        (1 - (1 - np.exp(-scale_displacement)) / scale_displacement) * canopy_height,
+    )
+
+    # No displacement in absence of vegetation
+    return np.nan_to_num(zero_place_displacement, nan=0.0).squeeze()
+
+
+def calculate_roughness_length_momentum(
+    canopy_height: NDArray[np.float32],
+    leaf_area_index: NDArray[np.float32],
+    zero_plane_displacement: NDArray[np.float32],
+    substrate_surface_drag_coefficient: float,
+    roughness_element_drag_coefficient: float,
+    roughness_sublayer_depth_parameter: float,
+    max_ratio_wind_to_friction_velocity: float,
+    min_roughness_length: float,
+    von_karman_constant: float,
+) -> NDArray[np.float32]:
+    """Calculate roughness length governing momentum transfer, [m].
+
+    Roughness length is defined as the height at which the mean velocity is zero due to
+    substrate roughness. Real surfaces such as the ground or vegetation are not smooth
+    and often have varying degrees of roughness. Roughness length accounts for that
+    effect. Implementation after :cite:t:`maclean_microclimc_2021`.
+
+    Args:
+        canopy_height: Canopy height, [m]
+        leaf_area_index: Total leaf area index, [m m-1]
+        zero_plane_displacement: Height above ground within the canopy where the wind
+            profile extrapolates to zero, [m]
+        substrate_surface_drag_coefficient: Substrate-surface drag coefficient,
+            dimensionless
+        roughness_element_drag_coefficient: Roughness-element drag coefficient
+        roughness_sublayer_depth_parameter: Parameter that characterizes the roughness
+            sublayer depth, dimensionless
+        max_ratio_wind_to_friction_velocity: Maximum ratio of wind velocity to friction
+            velocity, dimensionless
+        min_roughness_length: Minimum roughness length, [m]
+        von_karman_constant: Von Karman's constant, dimensionless constant describing
+            the logarithmic velocity profile of a turbulent fluid near a no-slip
+            boundary.
+
+    Returns:
+        momentum roughness length, [m]
+    """
+
+    # calculate ratio of wind velocity to friction velocity
+    ratio_wind_to_friction_velocity = np.sqrt(
+        substrate_surface_drag_coefficient
+        + (roughness_element_drag_coefficient * leaf_area_index) / 2
+    )
+
+    # if the ratio of wind velocity to friction velocity is larger than the set maximum,
+    # set the value to set maximum
+    set_maximum_ratio = np.where(
+        ratio_wind_to_friction_velocity > max_ratio_wind_to_friction_velocity,
+        max_ratio_wind_to_friction_velocity,
+        ratio_wind_to_friction_velocity,
+    )
+
+    # calculate initial roughness length
+    initial_roughness_length = (canopy_height - zero_plane_displacement) * np.exp(
+        -von_karman_constant * (1 / set_maximum_ratio)
+        - roughness_sublayer_depth_parameter
+    )
+
+    # if roughness smaller than the substrate surface drag coefficient, set to value to
+    # the substrate surface drag coefficient
+    roughness_length = np.where(
+        initial_roughness_length < substrate_surface_drag_coefficient,
+        substrate_surface_drag_coefficient,
+        initial_roughness_length,
+    )
+
+    return np.where(roughness_length <= 0, min_roughness_length, roughness_length)
+
+
+def calculate_diabatic_correction_above(
+    molar_density_air: Union[float, NDArray[np.float32]],
+    specific_heat_air: Union[float, NDArray[np.float32]],
+    temperature: NDArray[np.float32],
+    sensible_heat_flux: NDArray[np.float32],
+    friction_velocity: NDArray[np.float32],
+    wind_heights: NDArray[np.float32],
+    zero_plane_displacement: NDArray[np.float32],
+    celsius_to_kelvin: float,
+    von_karmans_constant: float,
+    yasuda_stability_parameter1: float,
+    yasuda_stability_parameter2: float,
+    yasuda_stability_parameter3: float,
+    diabatic_heat_momentum_ratio: float,
+) -> dict[str, NDArray[np.float32]]:
+    r"""Calculate the diabatic correction factors for momentum and heat above canopy.
+
+    Diabatic correction factor for heat and momentum are used to adjust wind profiles
+    for surface heating and cooling :cite:p:`maclean_microclimc_2021`. When the surface
+    is strongly heated, the diabatic correction factor for momemtum :math:`psi_{m}`
+    becomes negative and drops to values of around -1.5. In contrast, when the surface
+    is much cooler than the air above it, it increases to values around 4.
+
+    Args:
+        molar_density_air: molar density of air, [mol m-3]
+        specific_heat_air: specific heat of air, [J mol-1 K-1]
+        temperature: 2m temperature, [C]
+        sensible_heat_flux: Sensible heat flux from canopy to atmosphere above,
+            [W m-2], # TODO: could be the top entry of the general sensible heat flux
+        friction_velocity: Friction velocity, [m s-1]
+        wind_heights: Vector of heights for which wind speed is calculated, [m]
+        zero_plane_displacement: Height above ground within the canopy where the wind
+            profile extrapolates to zero, [m]
+        celsius_to_kelvin: Factor to convert temperature in Celsius to absolute
+            temperature in Kelvin
+        von_karmans_constant: Von Karman's constant, dimensionless constant describing
+            the logarithmic velocity profile of a turbulent fluid near a no-slip
+            boundary.
+        yasuda_stability_parameter1: Parameter to approximate diabatic correction
+            factors for heat and momentum after :cite:t:`yasuda_turbulent_1988`
+        yasuda_stability_parameter2: Parameter to approximate diabatic correction
+            factors for heat and momentum after :cite:t:`yasuda_turbulent_1988`
+        yasuda_stability_parameter3: Parameter to approximate diabatic correction
+            factors for heat and momentum after :cite:t:`yasuda_turbulent_1988`
+        diabatic_heat_momentum_ratio: Factor that relates diabatic correction
+            factors for heat and momentum after :cite:t:`yasuda_turbulent_1988`
+
+    Returns:
+        diabatic correction factors for heat (psi_h) and momentum (psi_m) transfer
+    """
+
+    # calculate atmospheric stability
+    stability = (
+        von_karmans_constant
+        * (wind_heights - zero_plane_displacement)
+        * sensible_heat_flux
+    ) / (
+        molar_density_air
+        * specific_heat_air
+        * (temperature + celsius_to_kelvin)
+        * friction_velocity
+    )
+
+    stable_condition = yasuda_stability_parameter1 * np.log(1 - stability)
+    unstable_condition = -yasuda_stability_parameter2 * np.log(
+        (1 + np.sqrt(1 - yasuda_stability_parameter3 * stability)) / 2
+    )
+
+    diabatic_correction_heat = np.where(
+        sensible_heat_flux < 0, stable_condition, unstable_condition
+    )
+
+    diabatic_correction_momentum = np.where(
+        sensible_heat_flux < 0,
+        diabatic_correction_heat,
+        diabatic_heat_momentum_ratio * diabatic_correction_heat,
+    )
+
+    return {"psi_m": diabatic_correction_momentum, "psi_h": diabatic_correction_heat}
+
+
+def calculate_mean_mixing_length(
+    canopy_height: NDArray[np.float32],
+    zero_plane_displacement: NDArray[np.float32],
+    roughness_length_momentum: NDArray[np.float32],
+    mixing_length_factor: float,
+) -> NDArray[np.float32]:
+    """Calculate mixing length for canopy air transport, [m].
+
+    The mean mixing length is used to calculate turbulent air transport inside vegetated
+    canopies. It is made equivalent to the above canopy value at the canopy surface. In
+    absence of vegetation, it is set to zero. Implementation after
+    :cite:t:`maclean_microclimc_2021`.
+
+    Args:
+        canopy_height: Canopy height, [m]
+        zero_plane_displacement: Height above ground within the canopy where the wind
+            profile extrapolates to zero, [m]
+        roughness_length_momentum: Momentum roughness length, [m]
+        mixing_length_factor: Factor in calculation of mean mixing length, dimensionless
+
+    Returns:
+        mixing length for canopy air transport, [m]
+    """
+
+    mean_mixing_length = (
+        mixing_length_factor * (canopy_height - zero_plane_displacement)
+    ) / np.log((canopy_height - zero_plane_displacement) / roughness_length_momentum)
+
+    return np.nan_to_num(mean_mixing_length, nan=0)
+
+
+def generate_relative_turbulence_intensity(
+    layer_heights: NDArray[np.float32],
+    min_relative_turbulence_intensity: float,
+    max_relative_turbulence_intensity: float,
+    increasing_with_height: bool,
+) -> NDArray[np.float32]:
+    """Generate relative turbulence intensity profile, dimensionless.
+
+    At the moment, default values are for a maize crop Shaw et al (1974)
+    Agricultural Meteorology, 13: 419-425. TODO adjust default to environment
+
+    Args:
+        layer_heights: heights of above ground layers, [m]
+        min_relative_turbulence_intensity: minimum relative turbulence intensity,
+            dimensionless
+        max_relative_turbulence_intensity: maximum relative turbulence intensity,
+            dimensionless
+        increasing_with_height: increasing logical indicating whether turbulence
+            intensity increases (TRUE) or decreases (FALSE) with height
+
+    Returns:
+        relative turbulence intensity for each node, dimensionless
+    """
+
+    if increasing_with_height:
+        turbulence_intensity = (
+            min_relative_turbulence_intensity
+            + (max_relative_turbulence_intensity - min_relative_turbulence_intensity)
+            * layer_heights
+        )
+    else:
+        turbulence_intensity = (
+            max_relative_turbulence_intensity
+            - (max_relative_turbulence_intensity - min_relative_turbulence_intensity)
+            * layer_heights
+        )
+    return turbulence_intensity
+
+
+def calculate_wind_attenuation_coefficient(
+    canopy_height: NDArray[np.float32],
+    leaf_area_index: NDArray[np.float32],
+    mean_mixing_length: NDArray[np.float32],
+    drag_coefficient: float,
+    relative_turbulence_intensity: NDArray[np.float32],
+) -> NDArray[np.float32]:
+    """Calculate wind attenuation coefficient, dimensionless.
+
+    The wind attenuation coefficient describes how wind is slowed down by the presence
+    of vegetation. In absence of vegetation, the coefficient is set to zero.
+    Implementation after :cite:t:`maclean_microclimc_2021`.
+
+    Args:
+        canopy_height: Canopy height, [m]
+        leaf_area_index: Total leaf area index, [m m-1]
+        mean_mixing_length: Mixing length for canopy air transport, [m]
+        drag_coefficient: Drag coefficient, dimensionless
+        relative_turbulence_intensity: Relative turbulence intensity, dimensionless
+
+    Returns:
+        wind attenuation coefficient, dimensionless
+    """
+
+    intermediate_coefficient = (
+        (drag_coefficient * leaf_area_index * canopy_height)
+        / (
+            2
+            * mean_mixing_length
+            * relative_turbulence_intensity[0 : len(leaf_area_index)]
+        ),
+    )
+
+    return np.nan_to_num(intermediate_coefficient, nan=0).squeeze()
+
+
+def wind_log_profile(
+    height: Union[float, NDArray[np.float32]],
+    zeroplane_displacement: Union[float, NDArray[np.float32]],
+    roughness_length_momentum: Union[float, NDArray[np.float32]],
+    diabatic_correction_momentum: Union[float, NDArray[np.float32]],
+) -> NDArray[np.float32]:
+    """Calculate logarithmic wind profile.
+
+    Args:
+        height: Array of heights for which wind speed is calculated, [m]
+        zeroplane_displacement: Height above ground within the canopy where the wind
+            profile extrapolates to zero, [m]
+        roughness_length_momentum: Momentum roughness length, [m]
+        diabatic_correction_momentum: Diabatic correction factor for momentum
+
+    Returns:
+        logarithmic wind profile
+    """
+
+    wind_profile = (
+        np.log((height - zeroplane_displacement) / roughness_length_momentum)
+        + diabatic_correction_momentum,
+    )
+
+    return np.nan_to_num(wind_profile, nan=1).squeeze()
+
+
+def calculate_fricition_velocity(
+    wind_speed_ref: NDArray[np.float32],
+    reference_height: Union[float, NDArray[np.float32]],
+    zeroplane_displacement: NDArray[np.float32],
+    roughness_length_momentum: NDArray[np.float32],
+    diabatic_correction_momentum: Union[float, NDArray[np.float32]],
+    von_karmans_constant: float,
+) -> NDArray[np.float32]:
+    """Calculate friction velocity from wind speed at reference height, [m s-1].
+
+    Args:
+        wind_speed_ref: Wind speed at reference height, [m s-1]
+        reference_height: Height of wind measurement, [m]
+        zero_plane_displacement: Height above ground within the canopy where the wind
+            profile extrapolates to zero, [m]
+        roughness_length_momentum: Momentum roughness length, [m]
+        diabatic_correction_momentum: Diabatic correction factor for momentum as
+            returned by
+            :func:`~virtual_rainforest.models.abiotic.wind.calculate_diabatic_correction_above`
+        von_karmans_constant: Von Karman's constant, dimensionless constant describing
+            the logarithmic velocity profile of a turbulent fluid near a no-slip
+            boundary.
+
+    Returns:
+        friction velocity
+    """
+
+    wind_profile_reference = wind_log_profile(
+        height=reference_height,
+        zeroplane_displacement=zeroplane_displacement,
+        roughness_length_momentum=roughness_length_momentum,
+        diabatic_correction_momentum=diabatic_correction_momentum,
+    )
+
+    return von_karmans_constant * (wind_speed_ref / wind_profile_reference)
+
+
+def calculate_wind_above_canopy(
+    friction_velocity: NDArray[np.float32],
+    wind_height_above: NDArray[np.float32],
+    zeroplane_displacement: NDArray[np.float32],
+    roughness_length_momentum: NDArray[np.float32],
+    diabatic_correction_momentum: NDArray[np.float32],
+    von_karmans_constant: float,
+    min_wind_speed_above_canopy: float,
+) -> NDArray[np.float32]:
+    """Calculate wind speed above canopy from wind speed at reference height, [m s-1].
+
+    Wind speed above the canopy dictates heat and vapor exchange between the canopy
+    and the air above it, and therefore ultimately determines temperature and vapor
+    profiles.
+    The wind profile above canopy typically follows a logarithmic height profile, which
+    extrapolates to zero roughly two thirds of the way to the top of the canopy. The
+    profile itself is thus dependent on the height of the canopy, but also on the
+    roughness of the vegetation layer, which causes wind shear. We follow the
+    implementation by :cite:t:`campbell_introduction_1998` as described in
+    :cite:t:`maclean_microclimc_2021`.
+
+    Args:
+        friction_velocity: friction velocity, [m s-1]
+        wind_height_above: Height above canopy for which wind speed is required, [m]
+        zero_plane_displacement: Height above ground within the canopy where the wind
+            profile extrapolates to zero, [m]
+        roughness_length_momentum: Momentum roughness length, [m]
+        diabatic_correction_momentum: Diabatic correction factor for momentum as
+            returned by
+            :func:`~virtual_rainforest.models.abiotic.wind.calculate_diabatic_correction_above`
+        von_karmans_constant: Von Karman's constant, dimensionless constant describing
+            the logarithmic velocity profile of a turbulent fluid near a no-slip
+            boundary.
+
+    Returns:
+        wind speed at required heights above canopy, [m s-1]
+    """
+
+    wind_profile_above = wind_log_profile(
+        height=wind_height_above,
+        zeroplane_displacement=zeroplane_displacement,
+        roughness_length_momentum=roughness_length_momentum,
+        diabatic_correction_momentum=diabatic_correction_momentum,
+    )
+    wind_profile = (friction_velocity / von_karmans_constant) * wind_profile_above
+
+    return np.where(
+        wind_profile < min_wind_speed_above_canopy,
+        min_wind_speed_above_canopy,
+        wind_profile,
+    )
+
+
+def calculate_wind_canopy(
+    top_of_canopy_wind_speed: NDArray[np.float32],
+    wind_layer_heights: NDArray[np.float32],
+    canopy_height: NDArray[np.float32],
+    attenuation_coefficient: NDArray[np.float32],
+    min_windspeed_below_canopy: float,
+) -> NDArray[np.float32]:
+    """Calculate wind speed in a multi-layer canopy, [m s-1].
+
+    This function can be extended to account for edge distance effects.
+
+    Args:
+        top_of_canopy_wind_speed: Wind speed at top of canopy layer, [m s-1]
+        wind_layer_heights: Heights of canopy layer nodes, [m]
+        canopy_height: Height to top of canopy layer, [m]
+        attenuation_coefficient: Mean attenuation coefficient based on the profile
+            calculated by
+            :func:`~virtual_rainforest.models.abiotic.wind.calculate_wind_attenuation_coefficient`
+        min_windspeed_below_canopy: Minimum wind speed below the canopy or in absence of
+            vegetation, [m/s]. This value is set to avoid dividion by zero.
+
+    Returns:
+        wind speed at height of canopy node, [m s-1]
+    """
+
+    return np.nan_to_num(
+        top_of_canopy_wind_speed
+        * np.exp(attenuation_coefficient * ((wind_layer_heights / canopy_height) - 1)),
+        nan=min_windspeed_below_canopy,
+    ).squeeze()
+
+
+def calculate_wind_profile(
+    canopy_height: NDArray[np.float32],
+    wind_height_above: NDArray[np.float32],
+    wind_layer_heights: NDArray[np.float32],
+    leaf_area_index: NDArray[np.float32],
+    air_temperature: NDArray[np.float32],
+    atmospheric_pressure: NDArray[np.float32],
+    sensible_heat_flux_topofcanopy: NDArray[np.float32],
+    wind_speed_ref: NDArray[np.float32],
+    wind_reference_height: Union[float, NDArray[np.float32]],
+    abiotic_constants: AbioticConsts,
+    core_constants: CoreConsts,
+) -> dict[str, NDArray[np.float32]]:
+    r"""Calculate wind speed above and below the canopy, [m s-1].
+
+    The wind profile above the canopy is described as follows (based on
+    :cite:p:`campbell_introduction_1998` as implemented in
+    :cite:t:`maclean_microclimc_2021`):
+
+    :math:`u_z = \frac{u^{*}}{0.4} ln \frac{z-d}{z_M} + \phi_M`
+
+    where :math:`u_z` is wind speed at height :math:`z` above the canopy, :math:`d` is
+    the height above ground within the canopy where the wind profile extrapolates to
+    zero, :math:`z_m` the roughness length for momentum, :math:`\phi_M` is a diabatic
+    correction for momentum and :math:`u^{*}` is the friction velocity, which gives the
+    wind speed at height :math:`d + z_m`.
+
+    The wind profile below canopy is derived as follows:
+
+    :math:`u_z = u_h exp(a(\frac{z}{h} - 1))`
+
+    where :math:`u_z` is wind speed at height :math:`z` within the canopy, :math:`u_h`
+    is wind speed at the top of the canopy at height :math:`h`, and :math:`a` is a wind
+    attenuation coefficient given by :math:`a = 2 l_m i_w`, where :math:`c_d` is a drag
+    coefficient that varies with leaf inclination and shape, :math:`i_w` is a
+    coefficient describing relative turbulence intensity and :math:`l_m` is the mean
+    mixing length, equivalent to the free space between the leaves and stems. For
+    details, see :cite:t:`maclean_microclimc_2021`.
+
+    Args:
+        canopy_height: Canopy height, [m]
+        wind_height_above: Height above canopy for which wind speed is required, [m]
+        wind_layer_heights: Heights of canopy layer nodes, [m]
+        leaf_area_index: Leaf area index, [m m-1]
+        air_temperature: Air temperature, [C]
+        atmospheric_pressure: Atmospheric pressure, [kPa]
+        sensible_heat_flux_topofcanopy: Sensible heat flux from the top of the canopy to
+            the atmosphere, [W m-2],
+        wind_speed_ref: Wind speed at reference height, [m s-1]
+        wind_reference_height: Reference height for wind measurement, [m]
+        abiotic_constants: Specific constants for the abiotic model
+        core_constants: Universal constants shared across all models
+
+    Returns:
+        dictionnary that contains wind speed above the canopy, [m s-1], wind speed
+        within and below the canopy, [m s-1], and friction velocity, [m s-1]
+    """
+
+    output = {}
+
+    # TODO adjust wind to 2m above canopy?
+
+    molar_density_air = calculate_molar_density_air(
+        temperature=air_temperature,
+        atmospheric_pressure=atmospheric_pressure,
+        standard_mole=core_constants.standard_mole,
+        standard_pressure=core_constants.standard_pressure,
+        celsius_to_kelvin=core_constants.zero_Celsius,
+    )
+
+    specific_heat_air = calculate_specific_heat_air(
+        temperature=air_temperature,
+        molar_heat_capacity_air=core_constants.molar_heat_capacity_air,
+        specific_heat_equ_factor_1=abiotic_constants.specific_heat_equ_factor_1,
+        specific_heat_equ_factor_2=abiotic_constants.specific_heat_equ_factor_2,
+    )
+
+    leaf_area_index_sum = leaf_area_index.sum(axis=1)
+
+    zero_plane_displacement = calculate_zero_plane_displacement(
+        canopy_height=canopy_height,
+        leaf_area_index=leaf_area_index_sum,
+        zero_plane_scaling_parameter=abiotic_constants.zero_plane_scaling_parameter,
+    )
+
+    roughness_length_momentum = calculate_roughness_length_momentum(
+        canopy_height=canopy_height,
+        leaf_area_index=leaf_area_index_sum,
+        zero_plane_displacement=zero_plane_displacement,
+        substrate_surface_drag_coefficient=(
+            AbioticConsts.substrate_surface_drag_coefficient
+        ),
+        roughness_element_drag_coefficient=(
+            AbioticConsts.roughness_element_drag_coefficient
+        ),
+        roughness_sublayer_depth_parameter=(
+            AbioticConsts.roughness_sublayer_depth_parameter
+        ),
+        max_ratio_wind_to_friction_velocity=(
+            AbioticConsts.max_ratio_wind_to_friction_velocity
+        ),
+        min_roughness_length=AbioticConsts.min_roughness_length,
+        von_karman_constant=CoreConsts.von_karmans_constant,
+    )
+
+    friction_velocity_uncorrected = calculate_fricition_velocity(
+        wind_speed_ref=wind_speed_ref,
+        reference_height=wind_reference_height,
+        zeroplane_displacement=zero_plane_displacement,
+        roughness_length_momentum=roughness_length_momentum,
+        diabatic_correction_momentum=0.0,
+        von_karmans_constant=core_constants.von_karmans_constant,
+    )
+
+    diabatic_correction_above = calculate_diabatic_correction_above(
+        molar_density_air=molar_density_air,
+        specific_heat_air=specific_heat_air,
+        temperature=air_temperature,
+        sensible_heat_flux=sensible_heat_flux_topofcanopy,
+        friction_velocity=friction_velocity_uncorrected,
+        wind_heights=wind_layer_heights,
+        zero_plane_displacement=zero_plane_displacement,
+        celsius_to_kelvin=CoreConsts.zero_Celsius,
+        von_karmans_constant=CoreConsts.von_karmans_constant,
+        yasuda_stability_parameter1=AbioticConsts.yasuda_stability_parameter1,
+        yasuda_stability_parameter2=AbioticConsts.yasuda_stability_parameter2,
+        yasuda_stability_parameter3=AbioticConsts.yasuda_stability_parameter3,
+        diabatic_heat_momentum_ratio=AbioticConsts.diabatic_heat_momentum_ratio,
+    )
+
+    friction_velocity = calculate_fricition_velocity(
+        wind_speed_ref=wind_speed_ref,
+        reference_height=wind_reference_height,
+        zeroplane_displacement=zero_plane_displacement,
+        roughness_length_momentum=roughness_length_momentum,
+        diabatic_correction_momentum=diabatic_correction_above["psi_m"],
+        von_karmans_constant=core_constants.von_karmans_constant,
+    )
+    output["friction_velocity"] = friction_velocity
+
+    mean_mixing_length = calculate_mean_mixing_length(
+        canopy_height=canopy_height,
+        zero_plane_displacement=zero_plane_displacement,
+        roughness_length_momentum=roughness_length_momentum,
+        mixing_length_factor=abiotic_constants.mixing_length_factor,
+    )
+
+    relative_turbulence_intensity = generate_relative_turbulence_intensity(
+        layer_heights=wind_layer_heights,
+        min_relative_turbulence_intensity=(
+            abiotic_constants.min_relative_turbulence_intensity
+        ),
+        max_relative_turbulence_intensity=(
+            abiotic_constants.max_relative_turbulence_intensity
+        ),
+        increasing_with_height=AbioticConsts.turbulence_sign,
+    )
+
+    attennuation_coefficient = calculate_wind_attenuation_coefficient(
+        canopy_height=canopy_height,
+        leaf_area_index=leaf_area_index,
+        mean_mixing_length=mean_mixing_length,
+        drag_coefficient=abiotic_constants.drag_coefficient,
+        relative_turbulence_intensity=relative_turbulence_intensity,
+    )
+
+    wind_speed_above_canopy = calculate_wind_above_canopy(
+        friction_velocity=friction_velocity[0],
+        wind_height_above=wind_height_above,
+        zeroplane_displacement=zero_plane_displacement,
+        roughness_length_momentum=roughness_length_momentum,
+        diabatic_correction_momentum=diabatic_correction_above["psi_m"][0],
+        von_karmans_constant=core_constants.von_karmans_constant,
+        min_wind_speed_above_canopy=abiotic_constants.min_wind_speed_above_canopy,
+    )
+    output["wind_speed_above_canopy"] = wind_speed_above_canopy
+
+    wind_speed_canopy = calculate_wind_canopy(
+        top_of_canopy_wind_speed=wind_speed_above_canopy,
+        wind_layer_heights=wind_layer_heights,
+        canopy_height=canopy_height,
+        attenuation_coefficient=attennuation_coefficient[0],
+        min_windspeed_below_canopy=abiotic_constants.min_windspeed_below_canopy,
+    )
+    output["wind_speed_canopy"] = wind_speed_canopy
+
+    return output
diff --git a/virtual_rainforest/models/abiotic_simple/abiotic_simple_model.py b/virtual_rainforest/models/abiotic_simple/abiotic_simple_model.py
index 1bd888263..9a944aa73 100644
--- a/virtual_rainforest/models/abiotic_simple/abiotic_simple_model.py
+++ b/virtual_rainforest/models/abiotic_simple/abiotic_simple_model.py
@@ -13,6 +13,8 @@ class as a child of the :class:`~virtual_rainforest.core.base_model.BaseModel` c
 config dictionary to the required types they are caught and then logged, and at the end
 of the unpacking an error is thrown. This error should be caught and handled by
 downstream functions so that all model configuration failures can be reported as one.
+
+TODO update temperatures to Kelvin
 """  # noqa: D205, D415
 
 from __future__ import annotations
diff --git a/virtual_rainforest/models/abiotic_simple/microclimate.py b/virtual_rainforest/models/abiotic_simple/microclimate.py
index b045d3d44..765344034 100644
--- a/virtual_rainforest/models/abiotic_simple/microclimate.py
+++ b/virtual_rainforest/models/abiotic_simple/microclimate.py
@@ -8,6 +8,8 @@
 1 m depth which equals the mean annual temperature.
 The module also provides a constant vertical profile of atmospheric pressure and
 :math:`\ce{CO2}`.
+
+TODO change tenperatures to Kelvin
 """  # noqa: D205, D415
 
 import numpy as np
diff --git a/virtual_rainforest/models/animals/animal_cohorts.py b/virtual_rainforest/models/animals/animal_cohorts.py
index 375d5e770..33facb677 100644
--- a/virtual_rainforest/models/animals/animal_cohorts.py
+++ b/virtual_rainforest/models/animals/animal_cohorts.py
@@ -16,11 +16,7 @@
 
 from virtual_rainforest.core.logger import LOGGER
 from virtual_rainforest.models.animals.animal_traits import DietType
-from virtual_rainforest.models.animals.constants import (
-    DECAY_FRACTION_CARCASSES,
-    DECAY_FRACTION_EXCREMENT,
-    FLOW_TO_REPRODUCTIVE_MASS_THRESHOLD,
-)
+from virtual_rainforest.models.animals.constants import AnimalConsts
 from virtual_rainforest.models.animals.decay import CarcassPool
 from virtual_rainforest.models.animals.functional_group import FunctionalGroup
 from virtual_rainforest.models.animals.protocols import Consumer, DecayPool, Resource
@@ -42,6 +38,7 @@ def __init__(
         mass: float,
         age: float,
         individuals: int,
+        constants: AnimalConsts = AnimalConsts(),
     ) -> None:
         if age < 0:
             raise ValueError("Age must be a positive number.")
@@ -63,6 +60,8 @@ def __init__(
         """The age of the animal cohort [days]."""
         self.individuals = individuals
         """The number of individuals in this cohort."""
+        self.constants = constants
+        """Animal constants."""
         self.damuth_density: int = damuths_law(
             self.functional_group.adult_mass, self.functional_group.damuths_law_terms
         )
@@ -92,9 +91,9 @@ def __init__(
 
         # TODO - In future this should be parameterised using a constants dataclass, but
         # this hasn't yet been implemented for the animal model
-        self.decay_fraction_excrement: float = DECAY_FRACTION_EXCREMENT
+        self.decay_fraction_excrement: float = self.constants.decay_fraction_excrement
         """The fraction of excrement which decays before it gets consumed."""
-        self.decay_fraction_carcasses: float = DECAY_FRACTION_CARCASSES
+        self.decay_fraction_carcasses: float = self.constants.decay_fraction_carcasses
         """The fraction of carcass biomass which decays before it gets consumed."""
 
     def metabolize(self, temperature: float, dt: timedelta64) -> None:
@@ -326,7 +325,9 @@ def eat(self, food: Resource, pool: DecayPool) -> float:
         # get the per-individual energetic gain from the bulk value
         mass_consumed = food.get_eaten(self, pool) / self.individuals
 
-        if self.is_below_mass_threshold(FLOW_TO_REPRODUCTIVE_MASS_THRESHOLD):
+        if self.is_below_mass_threshold(
+            self.constants.flow_to_reproductive_mass_threshold
+        ):
             # if current mass equals or exceeds standard adult mass, gains to repro mass
             self.mass_current += mass_consumed
         else:
diff --git a/virtual_rainforest/models/animals/animal_communities.py b/virtual_rainforest/models/animals/animal_communities.py
index 94f2e9579..8c6b78abc 100644
--- a/virtual_rainforest/models/animals/animal_communities.py
+++ b/virtual_rainforest/models/animals/animal_communities.py
@@ -18,10 +18,7 @@
 from virtual_rainforest.core.data import Data
 from virtual_rainforest.core.logger import LOGGER
 from virtual_rainforest.models.animals.animal_cohorts import AnimalCohort
-from virtual_rainforest.models.animals.constants import (
-    BIRTH_MASS_THRESHOLD,
-    DISPERSAL_MASS_THRESHOLD,
-)
+from virtual_rainforest.models.animals.constants import AnimalConsts
 from virtual_rainforest.models.animals.decay import CarcassPool, ExcrementPool
 from virtual_rainforest.models.animals.functional_group import FunctionalGroup
 from virtual_rainforest.models.animals.plant_resources import PlantResources
@@ -47,6 +44,7 @@ def __init__(
         community_key: int,
         neighbouring_keys: list[int],
         get_destination: Callable[[int], "AnimalCommunity"],
+        constants: AnimalConsts = AnimalConsts(),
     ) -> None:
         # The constructor of the AnimalCommunity class.
         self.data = data
@@ -59,6 +57,8 @@ def __init__(
         """List of integer keys of neighbouring communities."""
         self.get_destination = get_destination
         """Callable get_destination from AnimalModel."""
+        self.constants = constants
+        """Animal constants."""
 
         self.animal_cohorts: dict[str, list[AnimalCohort]] = {
             k.name: [] for k in self.functional_groups
@@ -96,7 +96,11 @@ def populate_community(self) -> None:
             )
 
             cohort = AnimalCohort(
-                functional_group, functional_group.adult_mass, 0.0, individuals
+                functional_group,
+                functional_group.adult_mass,
+                0.0,
+                individuals,
+                self.constants,
             )
             self.animal_cohorts[functional_group.name].append(cohort)
 
@@ -121,7 +125,7 @@ def migrate(self, migrant: AnimalCohort, destination: AnimalCommunity) -> None:
     def migrate_community(self) -> None:
         """This handles migrating all cohorts in a community."""
         for cohort in self.all_animal_cohorts:
-            if cohort.is_below_mass_threshold(DISPERSAL_MASS_THRESHOLD):
+            if cohort.is_below_mass_threshold(self.constants.dispersal_mass_threshold):
                 # Random walk destination from the neighbouring keys
                 destination_key = choice(self.neighbouring_keys)
                 destination = self.get_destination(destination_key)
@@ -177,6 +181,7 @@ def birth(self, parent_cohort: AnimalCohort) -> None:
                 parent_cohort.functional_group.birth_mass,
                 0.0,
                 number_offspring,
+                self.constants,
             )
         )
 
@@ -185,7 +190,7 @@ def birth_community(self) -> None:
 
         # reproduction occurs for cohorts with sufficient reproductive mass
         for cohort in self.all_animal_cohorts:
-            if not cohort.is_below_mass_threshold(BIRTH_MASS_THRESHOLD):
+            if not cohort.is_below_mass_threshold(self.constants.birth_mass_threshold):
                 self.birth(cohort)
 
     def forage_community(self) -> None:
@@ -203,7 +208,9 @@ def forage_community(self) -> None:
         """
         # Generate the plant resources for foraging.
         plant_community: PlantResources = PlantResources(
-            data=self.data, cell_id=self.community_key
+            data=self.data,
+            cell_id=self.community_key,
+            constants=self.constants,
         )
 
         plant_list = [plant_community]
diff --git a/virtual_rainforest/models/animals/animal_model.py b/virtual_rainforest/models/animals/animal_model.py
index e6de497fa..9aef7bf0e 100644
--- a/virtual_rainforest/models/animals/animal_model.py
+++ b/virtual_rainforest/models/animals/animal_model.py
@@ -31,6 +31,7 @@
 from virtual_rainforest.core.data import Data
 from virtual_rainforest.core.logger import LOGGER
 from virtual_rainforest.models.animals.animal_communities import AnimalCommunity
+from virtual_rainforest.models.animals.constants import AnimalConsts
 from virtual_rainforest.models.animals.functional_group import FunctionalGroup
 
 
@@ -69,6 +70,7 @@ def __init__(
         data: Data,
         update_interval: Quantity,
         functional_groups: list[FunctionalGroup],
+        constants: AnimalConsts = AnimalConsts(),
         **kwargs: Any,
     ):
         super().__init__(data, update_interval, **kwargs)
@@ -82,6 +84,8 @@ def __init__(
 
         self.communities: dict[int, AnimalCommunity] = {}
         """Set empty dict for populating with communities."""
+        self.constants = constants
+        """Animal constants."""
         self._initialize_communities(functional_groups)
         """Create the dictionary of animal communities and populate each community with
         animal cohorts."""
@@ -126,6 +130,7 @@ def _initialize_communities(self, functional_groups: list[FunctionalGroup]) -> N
                 community_key=k,
                 neighbouring_keys=list(self.data.grid.neighbours[k]),
                 get_destination=self.get_community_by_key,
+                constants=self.constants,
             )
             for k in self.data.grid.cell_id
         }
@@ -155,7 +160,10 @@ def from_config(
 
         functional_groups_raw = config["animals"]["functional_groups"]
 
-        animal_functional_groups = [FunctionalGroup(**k) for k in functional_groups_raw]
+        animal_functional_groups = [
+            FunctionalGroup(**k, constants=AnimalConsts())
+            for k in functional_groups_raw
+        ]
 
         LOGGER.info(
             "Information required to initialise the animal model successfully "
diff --git a/virtual_rainforest/models/animals/constants.py b/virtual_rainforest/models/animals/constants.py
index dbd156c4c..1b599ea90 100644
--- a/virtual_rainforest/models/animals/constants.py
+++ b/virtual_rainforest/models/animals/constants.py
@@ -11,6 +11,7 @@
 
 from dataclasses import dataclass, field
 
+from virtual_rainforest.core.constants_class import ConstantsDataclass
 from virtual_rainforest.models.animals.animal_traits import (
     DietType,
     MetabolicType,
@@ -19,7 +20,7 @@
 
 
 @dataclass(frozen=True)
-class AnimalConsts:
+class AnimalConsts(ConstantsDataclass):
     """Dataclass to store all constants related to metabolic rates.
 
     TODO: The entire constants fille will be reworked in this style after the energy to
@@ -41,7 +42,99 @@ class AnimalConsts:
         }
     )
 
+    damuths_law_terms: dict[TaxaType, dict[DietType, tuple[float, float]]] = field(
+        default_factory=lambda: {
+            TaxaType.MAMMAL: {
+                DietType.HERBIVORE: (-0.75, 4.23),
+                DietType.CARNIVORE: (-0.75, 1.00),
+            },
+            TaxaType.BIRD: {
+                DietType.HERBIVORE: (-0.75, 5.00),
+                DietType.CARNIVORE: (-0.75, 2.00),
+            },
+            TaxaType.INSECT: {
+                DietType.HERBIVORE: (-0.75, 5.00),
+                DietType.CARNIVORE: (-0.75, 2.00),
+            },
+        }
+    )
 
+    fat_mass_terms: dict[TaxaType, tuple[float, float]] = field(
+        default_factory=lambda: {
+            TaxaType.MAMMAL: (1.19, 0.02),  # Scaling of mammalian herbivore fat mass
+            TaxaType.BIRD: (1.19, 0.05),  # Toy Values
+            TaxaType.INSECT: (1.19, 0.05),  # Toy Values
+        }
+    )
+
+    muscle_mass_terms: dict[TaxaType, tuple[float, float]] = field(
+        default_factory=lambda: {
+            TaxaType.MAMMAL: (1.0, 0.38),  # Scaling of mammalian herbivore muscle mass
+            TaxaType.BIRD: (1.0, 0.40),  # Toy Values
+            TaxaType.INSECT: (1.0, 0.40),  # Toy Values
+        }
+    )
+
+    intake_rate_terms: dict[TaxaType, tuple[float, float]] = field(
+        default_factory=lambda: {
+            TaxaType.MAMMAL: (0.71, 0.63),  # Mammalian maximum intake rate
+            TaxaType.BIRD: (0.7, 0.50),  # Toy Values
+            TaxaType.INSECT: (0.7, 0.50),  # Toy Values
+        }
+    )
+
+    energy_density: dict[str, float] = field(
+        default_factory=lambda: {
+            "meat": 7000.0,  # Energy of mammal meat [J/g]
+            "plant": 18200000.0,  # Energy of plant food [J/g]
+        }
+    )
+
+    conversion_efficiency: dict[DietType, float] = field(
+        default_factory=lambda: {
+            DietType.HERBIVORE: 0.1,  # Toy value
+            DietType.CARNIVORE: 0.25,  # Toy value
+        }
+    )
+
+    mechanical_efficiency: dict[DietType, float] = field(
+        default_factory=lambda: {
+            DietType.HERBIVORE: 0.9,  # Toy value
+            DietType.CARNIVORE: 0.8,  # Toy value
+        }
+    )
+
+    prey_mass_scaling_terms: dict[
+        MetabolicType, dict[TaxaType, tuple[float, float]]
+    ] = field(
+        default_factory=lambda: {
+            MetabolicType.ENDOTHERMIC: {
+                TaxaType.MAMMAL: (1.0, 1.0),  # Toy values
+                TaxaType.BIRD: (1.0, 1.0),  # Toy values
+            },
+            MetabolicType.ECTOTHERMIC: {TaxaType.INSECT: (1.0, 1.0)},  # Toy values
+        }
+    )
+
+    longevity_scaling_terms: dict[TaxaType, tuple[float, float]] = field(
+        default_factory=lambda: {
+            TaxaType.MAMMAL: (0.25, 0.02),  # Toy values
+            TaxaType.BIRD: (0.25, 0.05),  # Toy values
+            TaxaType.INSECT: (0.25, 0.05),  # Toy values
+        }
+    )
+
+    birth_mass_threshold: float = 1.5  # Threshold for reproduction
+    flow_to_reproductive_mass_threshold: float = (
+        1.0  # Threshold of trophic flow to reproductive mass
+    )
+    dispersal_mass_threshold: float = 0.75  # Threshold for dispersal
+    energy_percentile_threshold: float = 0.5  # Threshold for initiating migration
+    decay_fraction_excrement: float = 0.5  # Decay fraction for excrement
+    decay_fraction_carcasses: float = 0.2  # Decay fraction for carcasses
+
+
+"""
 METABOLIC_RATE_TERMS: dict[MetabolicType, dict[str, tuple[float, float]]] = {
     # Parameters from Madingley, mass based metabolic rates
     MetabolicType.ENDOTHERMIC: {
@@ -161,7 +254,7 @@ class AnimalConsts:
 DISPERSAL_MASS_THRESHOLD: float = 0.75  # Toy value for thesholding dispersal.
 
 ENERGY_PERCENTILE_THRESHOLD: float = 0.5  # Toy value for initiating migration
-
+"""
 DECAY_FRACTION_EXCREMENT: float = 0.5
 """Fraction of excrement that is assumed to decay rather than be consumed [unitless].
 
@@ -176,3 +269,6 @@ class AnimalConsts:
 [unitless]. TODO - The number given here is very much made up, see
 :attr:`DECAY_FRACTION_EXCREMENT` for details of how this should be changed in future.
 """
+BOLTZMANN_CONSTANT: float = 8.617333262145e-5  # Boltzmann constant [eV/K]
+
+TEMPERATURE: float = 37.0  # Toy temperature for setting up metabolism [C].
diff --git a/virtual_rainforest/models/animals/functional_group.py b/virtual_rainforest/models/animals/functional_group.py
index dc92adb62..1d6cd5a09 100644
--- a/virtual_rainforest/models/animals/functional_group.py
+++ b/virtual_rainforest/models/animals/functional_group.py
@@ -10,17 +10,7 @@
     MetabolicType,
     TaxaType,
 )
-from virtual_rainforest.models.animals.constants import (
-    CONVERSION_EFFICIENCY,
-    DAMUTHS_LAW_TERMS,
-    FAT_MASS_TERMS,
-    INTAKE_RATE_TERMS,
-    LONGEVITY_SCALING_TERMS,
-    MECHANICAL_EFFICIENCY,
-    METABOLIC_RATE_TERMS,
-    MUSCLE_MASS_TERMS,
-    PREY_MASS_SCALING_TERMS,
-)
+from virtual_rainforest.models.animals.constants import AnimalConsts
 
 
 class FunctionalGroup:
@@ -44,6 +34,7 @@ def __init__(
         metabolic_type: str,
         birth_mass: float,
         adult_mass: float,
+        constants: AnimalConsts = AnimalConsts(),
     ) -> None:
         """The constructor for the FunctionalGroup class."""
 
@@ -59,27 +50,35 @@ def __init__(
         """The mass of the functional group at birth."""
         self.adult_mass = adult_mass
         """The mass of the functional group at adulthood."""
-        self.metabolic_rate_terms = METABOLIC_RATE_TERMS[self.metabolic_type]
+        self.constants = constants
+        """Animal constants."""
+        self.metabolic_rate_terms = self.constants.metabolic_rate_terms[
+            self.metabolic_type
+        ]
         """The coefficient and exponent of metabolic rate."""
-        self.damuths_law_terms = DAMUTHS_LAW_TERMS[self.taxa][self.diet]
+        self.damuths_law_terms = self.constants.damuths_law_terms[self.taxa][self.diet]
         """The coefficient and exponent of damuth's law for population density."""
-        self.muscle_mass_terms = MUSCLE_MASS_TERMS[self.taxa]
+        self.muscle_mass_terms = self.constants.muscle_mass_terms[self.taxa]
         """The coefficient and exponent of muscle mass allometry."""
-        self.fat_mass_terms = FAT_MASS_TERMS[self.taxa]
+        self.fat_mass_terms = self.constants.fat_mass_terms[self.taxa]
         """The coefficient and exponent of fat mass allometry."""
-        self.intake_rate_terms = INTAKE_RATE_TERMS[self.taxa]
+        self.intake_rate_terms = self.constants.intake_rate_terms[self.taxa]
         """The coefficient and exponent of intake allometry."""
-        self.conversion_efficiency = CONVERSION_EFFICIENCY[self.diet]
+        self.conversion_efficiency = self.constants.conversion_efficiency[self.diet]
         """The conversion efficiency of the functional group based on diet."""
-        self.mechanical_efficiency = MECHANICAL_EFFICIENCY[self.diet]
+        self.mechanical_efficiency = self.constants.mechanical_efficiency[self.diet]
         """The mechanical transfer efficiency of a functional group based on diet."""
-        self.prey_scaling = PREY_MASS_SCALING_TERMS[self.metabolic_type][self.taxa]
+        self.prey_scaling = self.constants.prey_mass_scaling_terms[self.metabolic_type][
+            self.taxa
+        ]
         """The predator-prey mass ratio scaling relationship."""
-        self.longevity_scaling = LONGEVITY_SCALING_TERMS[self.taxa]
+        self.longevity_scaling = self.constants.longevity_scaling_terms[self.taxa]
         """The coefficient and exponent for lifespan allometry."""
 
 
-def import_functional_groups(fg_csv_file: str) -> list[FunctionalGroup]:
+def import_functional_groups(
+    fg_csv_file: str, constants: AnimalConsts
+) -> list[FunctionalGroup]:
     """The function to import pre-defined functional groups.
 
     This function is a first-pass of how we might import pre-defined functional groups.
@@ -117,6 +116,7 @@ def import_functional_groups(fg_csv_file: str) -> list[FunctionalGroup]:
             row.metabolic_type,
             row.birth_mass,
             row.adult_mass,
+            constants=constants,
         )
         for row in fg.itertuples()
     ]
diff --git a/virtual_rainforest/models/animals/plant_resources.py b/virtual_rainforest/models/animals/plant_resources.py
index 4d38af6d7..b355137e6 100644
--- a/virtual_rainforest/models/animals/plant_resources.py
+++ b/virtual_rainforest/models/animals/plant_resources.py
@@ -5,7 +5,7 @@
 from __future__ import annotations
 
 from virtual_rainforest.core.data import Data
-from virtual_rainforest.models.animals.constants import ENERGY_DENSITY
+from virtual_rainforest.models.animals.constants import AnimalConsts
 from virtual_rainforest.models.animals.protocols import Consumer, DecayPool
 
 
@@ -26,7 +26,7 @@ class PlantResources:
         cell_id: The cell id for the plant community to expose.
     """
 
-    def __init__(self, data: Data, cell_id: int) -> None:
+    def __init__(self, data: Data, cell_id: int, constants: AnimalConsts) -> None:
         # Store the data and extract the appropriate plant data
         self.data = data
         """A reference to the core data object."""
@@ -34,11 +34,12 @@ def __init__(self, data: Data, cell_id: int) -> None:
             data["layer_leaf_mass"].sel(cell_id=cell_id).sum(dim="layers").item()
         )
         """The mass of the plant leaf mass [kg]."""
-
+        self.constants = constants
+        """The animals constants."""
         # Calculate energy availability
         # TODO - this needs to be handed back to the plants model, which will define PFT
         #        specific conversions to different resources.
-        self.energy_density: float = ENERGY_DENSITY["plant"]
+        self.energy_density: float = self.constants.energy_density["plant"]
         """The energy (J) in a kg of plant [currently set to toy value of Alfalfa]."""
         self.energy_max: float = self.mass_current * self.energy_density
         """The maximum amount of energy that the cohort can have [J] [Alfalfa]."""
diff --git a/virtual_rainforest/models/hydrology/above_ground.py b/virtual_rainforest/models/hydrology/above_ground.py
index bb0dedc41..37ec45c85 100644
--- a/virtual_rainforest/models/hydrology/above_ground.py
+++ b/virtual_rainforest/models/hydrology/above_ground.py
@@ -2,6 +2,8 @@
 processes for the Virtual Rainforest. At the moment, this includes rain water
 interception by the canopy, soil evaporation, and functions related to surface
 runoff, bypass flow, and river discharge.
+
+TODO change temperatures to Kelvin
 """  # noqa: D205, D415
 
 from math import sqrt
@@ -30,7 +32,7 @@ def calculate_soil_evaporation(
     heat_transfer_coefficient: float,
     extinction_coefficient_global_radiation: float,
 ) -> NDArray[np.float32]:
-    r"""Calculate soil evaporation based classical bulk aerodynamic formulation.
+    r"""Calculate soil evaporation based on classical bulk aerodynamic formulation.
 
     This function uses the so-called 'alpha' method to estimate the evaporative flux
     :cite:p:`mahfouf_comparative_1991`.
@@ -51,7 +53,7 @@ def calculate_soil_evaporation(
 
     :math:`E_{act} = E_{g} * exp(-\kappa_{gb}*LAI)`
 
-    where :math:`kappa_{gb}` is the extinction coefficient for global radiation, and
+    where :math:`\kappa_{gb}` is the extinction coefficient for global radiation, and
     :math:`LAI` is the total leaf area index.
 
     Args:
diff --git a/virtual_rainforest/models/hydrology/hydrology_model.py b/virtual_rainforest/models/hydrology/hydrology_model.py
index 14c1b4ca2..b443adb8d 100644
--- a/virtual_rainforest/models/hydrology/hydrology_model.py
+++ b/virtual_rainforest/models/hydrology/hydrology_model.py
@@ -13,6 +13,8 @@ class as a child of the :class:`~virtual_rainforest.core.base_model.BaseModel` c
 config dictionary to the required types they are caught and then logged, and at the end
 of the unpacking an error is thrown. This error should be caught and handled by
 downstream functions so that all model configuration failures can be reported as one.
+
+TODO change temperature to Kelvin
 """  # noqa: D205, D415
 
 from __future__ import annotations