diff --git a/tests/models/abiotic/test_abiotic_model.py b/tests/models/abiotic/test_abiotic_model.py
index fd0facac5..fd252fe91 100644
--- a/tests/models/abiotic/test_abiotic_model.py
+++ b/tests/models/abiotic/test_abiotic_model.py
@@ -369,40 +369,49 @@ def test_setup_abiotic_model(
 
         model.update(time_index=0)
 
-    expected_soil_temp1 = lyr_strct.from_template()
-    expected_soil_temp1[lyr_strct.index_all_soil] = np.array(
-        [
-            [22.3935, 24.064858, 25.343753, 25.343753],
-            [20.040154, 20.068193, 20.089648, 20.089648],
-        ],
+        # Check that values fall within a reasonable expected range
+    soil_temps = model.data["soil_temperature"].isel(layers=lyr_strct.index_all_soil)
+
+    # To test with varying canopy layers, need to mask
+    canopy_mask = ~np.isnan(
+        dummy_climate_data_varying_canopy["canopy_temperature"].isel(
+            layers=lyr_strct.index_filled_canopy
+        )
     )
-    expected_soil_moist = lyr_strct.from_template()
-    expected_soil_moist[lyr_strct.index_all_soil] = np.array([5.0, 500])[:, None]
-    xr.testing.assert_allclose(
-        model.data["soil_temperature"], expected_soil_temp1, rtol=1e-3
+    atm_mask = ~np.isnan(
+        dummy_climate_data_varying_canopy["air_temperature"].isel(
+            layers=lyr_strct.index_filled_atmosphere
+        )
     )
-    xr.testing.assert_allclose(model.data["soil_moisture"], expected_soil_moist)
-
-    exp_air_temp = lyr_strct.from_template()
 
-    exp_air_temp[lyr_strct.index_filled_atmosphere] = np.array(
-        [
-            [30.000105, 30.000091, 30.000047, 30.000047],
-            [29.925696, 29.938505, 29.970583, 29.970583],
-            [29.42155, 29.613321, np.nan, np.nan],
-            [28.558176, np.nan, np.nan, np.nan],
-            [23.158518, 26.130326, 29.416104, 29.416104],
-        ]
+    canopy_temp_result = model.data["canopy_temperature"].isel(
+        layers=lyr_strct.index_filled_canopy
     )
-    xr.testing.assert_allclose(model.data["air_temperature"], exp_air_temp)
-
-    exp_canopytemp = lyr_strct.from_template()
-    exp_canopytemp[lyr_strct.index_filled_canopy] = np.array(
-        [
-            [28.489317, 31.736854, 31.631768, 31.631768],
-            [29.414784, 29.609833, np.nan, np.nan],
-            [28.551829, np.nan, np.nan, np.nan],
-        ]
+    air_temp_result = model.data["air_temperature"].isel(
+        layers=lyr_strct.index_filled_atmosphere
+    )
+    rel_hum_result = model.data["relative_humidity"].isel(
+        layers=lyr_strct.index_filled_atmosphere
     )
 
-    xr.testing.assert_allclose(model.data["canopy_temperature"], exp_canopytemp)
+    # Use the mask as a DataArray for .where()
+    valid_values_can_temp = canopy_temp_result.where(canopy_mask)
+    valid_values_air_temp = air_temp_result.where(atm_mask)
+    valid_values_rel_hum = rel_hum_result.where(atm_mask)
+
+    # Now drop the NaNs (i.e., masked values)
+    valid_values_can_temp_clean = valid_values_can_temp.dropna(dim="layers", how="any")
+    valid_values_air_temp_clean = valid_values_air_temp.dropna(dim="layers", how="any")
+    valid_values_rel_hum_clean = valid_values_rel_hum.dropna(dim="layers", how="any")
+
+    # Now do the test
+    assert ((soil_temps >= 18.0) & (soil_temps <= 28.0)).all()
+    assert (
+        (valid_values_can_temp_clean >= 15.0) & (valid_values_can_temp_clean <= 40.0)
+    ).all()
+    assert (
+        (valid_values_air_temp_clean >= 15.0) & (valid_values_air_temp_clean <= 40.0)
+    ).all()
+    assert (
+        (valid_values_rel_hum_clean >= 0.0) & (valid_values_rel_hum_clean <= 100.0)
+    ).all()
diff --git a/tests/models/abiotic/test_abiotic_tools.py b/tests/models/abiotic/test_abiotic_tools.py
index 11d924f29..466c084ed 100644
--- a/tests/models/abiotic/test_abiotic_tools.py
+++ b/tests/models/abiotic/test_abiotic_tools.py
@@ -255,3 +255,35 @@ def test_compute_layer_thickness_for_varying_canopy(
     )
 
     np.testing.assert_allclose(result, exp_result)
+
+
+def test_calculate_specific_humidity(
+    dummy_climate_data_varying_canopy, fixture_core_components
+):
+    """Test specific humidity."""
+
+    from virtual_ecosystem.models.abiotic.abiotic_tools import (
+        calculate_specific_humidity,
+    )
+
+    data = dummy_climate_data_varying_canopy
+    atm_index = fixture_core_components.layer_structure.index_filled_atmosphere
+
+    result = calculate_specific_humidity(
+        air_temperature=data["air_temperature"][atm_index].to_numpy(),
+        relative_humidity=data["relative_humidity"][atm_index].to_numpy(),
+        atmospheric_pressure=data["atmospheric_pressure"][atm_index].to_numpy(),
+        molecular_weight_ratio_water_to_dry_air=0.622,
+        pyrealm_const=PyrealmConst(),
+    )
+
+    exp_result = np.array(
+        [
+            [0.025064, 0.025064, 0.025064, 0.025064],
+            [0.024934, 0.024934, 0.024934, 0.024934],
+            [0.02412, 0.02412, np.nan, np.nan],
+            [0.02274, np.nan, np.nan, np.nan],
+            [0.01638, 0.01638, 0.01638, 0.01638],
+        ]
+    )
+    np.testing.assert_allclose(result, exp_result, rtol=1e-4, atol=1e-4)
diff --git a/tests/models/abiotic/test_energy_balance.py b/tests/models/abiotic/test_energy_balance.py
index 60c203a5f..39fc79b20 100644
--- a/tests/models/abiotic/test_energy_balance.py
+++ b/tests/models/abiotic/test_energy_balance.py
@@ -415,12 +415,12 @@ def test_update_humidity_vpd(
         density_air=data["density_air"][atm_index].to_numpy(),
         mixing_coefficient=mixing_coefficient,
         ventilation_rate=ventilation_rate,
-        wind_speed=data["wind_speed_ref"].isel(time_index=0).to_numpy(),
         molecular_weight_ratio_water_to_dry_air=(
             CoreConsts.molecular_weight_ratio_water_to_dry_air
         ),
         dry_air_factor=1 - CoreConsts.molecular_weight_ratio_water_to_dry_air,
         cell_area=fixture_core_components.grid.cell_area,
+        limits=(0, 60),
         time_interval=time_interval,
     )
 
diff --git a/tests/models/abiotic/test_microclimate.py b/tests/models/abiotic/test_microclimate.py
index 02b610150..5643b1fc8 100644
--- a/tests/models/abiotic/test_microclimate.py
+++ b/tests/models/abiotic/test_microclimate.py
@@ -5,9 +5,9 @@
 
 from virtual_ecosystem.core.constants import CoreConsts
 from virtual_ecosystem.models.abiotic.constants import AbioticConsts
+from virtual_ecosystem.models.abiotic_simple.constants import AbioticSimpleBounds
 
 
-# Test integration (TODO add structural and value range check)
 def test_run_microclimate(dummy_climate_data_varying_canopy, fixture_core_components):
     """Test microclimate function."""
 
@@ -25,6 +25,7 @@ def test_run_microclimate(dummy_climate_data_varying_canopy, fixture_core_compon
         abiotic_constants=AbioticConsts(),
         core_constants=CoreConsts(),
         pyrealm_const=PyrealmConst(),
+        abiotic_bounds=AbioticSimpleBounds(),
     )
 
     for var in [
@@ -35,73 +36,52 @@ def test_run_microclimate(dummy_climate_data_varying_canopy, fixture_core_compon
     ]:
         assert var in result
 
-    exp_soiltemp = lyr_str.from_template()
-    exp_soiltemp[lyr_str.index_all_soil] = np.array(
-        [
-            [21.095, 21.053, 20.627, 20.627],
-            [20.018, 20.017, 20.010, 20.010],
-        ]
-    )
-    np.testing.assert_allclose(
-        result["soil_temperature"][lyr_str.index_all_soil],
-        exp_soiltemp[lyr_str.index_all_soil],
-        rtol=1e-02,
-        atol=1e-02,
-    )
+    # Check that values fall within a reasonable expected range
+    soil_temps = result["soil_temperature"].isel(layers=lyr_str.index_all_soil)
 
-    exp_cantemp = lyr_str.from_template()
-    exp_cantemp[lyr_str.index_filled_canopy] = np.array(
-        [
-            [27.223404, 27.296955, 27.407225, 27.407225],
-            [28.871276, 28.871276, np.nan, np.nan],
-            [27.206485, np.nan, np.nan, np.nan],
-        ]
+    # To test with varying canopy layers, need to mask
+    canopy_mask = ~np.isnan(
+        dummy_climate_data_varying_canopy["canopy_temperature"].isel(
+            layers=lyr_str.index_filled_canopy
+        )
     )
-    np.testing.assert_allclose(
-        result["canopy_temperature"][lyr_str.index_filled_canopy],
-        exp_cantemp[lyr_str.index_filled_canopy],
-        rtol=1e-02,
-        atol=1e-02,
+    atm_mask = ~np.isnan(
+        dummy_climate_data_varying_canopy["air_temperature"].isel(
+            layers=lyr_str.index_filled_atmosphere
+        )
     )
 
-    # Check that all air temperature values fall within a reasonable expected range
+    canopy_temp_result = result["canopy_temperature"].isel(
+        layers=lyr_str.index_filled_canopy
+    )
     air_temp_result = result["air_temperature"].isel(
         layers=lyr_str.index_filled_atmosphere
     )
-
-    mask = ~np.isnan(
-        dummy_climate_data_varying_canopy["air_temperature"].isel(
-            layers=lyr_str.index_filled_atmosphere
-        )
+    rel_hum_result = result["relative_humidity"].isel(
+        layers=lyr_str.index_filled_atmosphere
     )
 
     # Use the mask as a DataArray for .where()
-    valid_values = air_temp_result.where(mask)
+    valid_values_can_temp = canopy_temp_result.where(canopy_mask)
+    valid_values_air_temp = air_temp_result.where(atm_mask)
+    valid_values_rel_hum = rel_hum_result.where(atm_mask)
 
     # Now drop the NaNs (i.e., masked values)
-    valid_values_clean = valid_values.dropna(
-        dim="layers", how="any"
-    )  # or "all" if you're doing 2D
+    valid_values_can_temp_clean = valid_values_can_temp.dropna(dim="layers", how="any")
+    valid_values_air_temp_clean = valid_values_air_temp.dropna(dim="layers", how="any")
+    valid_values_rel_hum_clean = valid_values_rel_hum.dropna(dim="layers", how="any")
 
     # Now do the test
-    assert ((valid_values_clean >= 16.0) & (valid_values_clean <= 36.0)).all()
-
-    exp_relhum = lyr_str.from_template()
-    exp_relhum[lyr_str.index_filled_atmosphere] = np.array(
-        [
-            [100, 100, 100, 100],
-            [100, 100, 100, 100],
-            [100, np.nan, np.nan, np.nan],
-            [100, np.nan, np.nan, np.nan],
-            [100, 100, 100, 100],
-        ]
-    )
-    np.testing.assert_allclose(
-        result["relative_humidity"],
-        exp_relhum,
-        rtol=1e-02,
-        atol=1e-02,
-    )
+    assert ((soil_temps >= 18.0) & (soil_temps <= 25.0)).all()
+    assert (
+        (valid_values_can_temp_clean >= 15.0) & (valid_values_can_temp_clean <= 40.0)
+    ).all()
+    assert (
+        (valid_values_air_temp_clean >= 15.0) & (valid_values_air_temp_clean <= 40.0)
+    ).all()
+    assert (
+        (valid_values_rel_hum_clean >= 0.0) & (valid_values_rel_hum_clean <= 100.0)
+    ).all()
 
 
 def test_run_microclimate_subdaily(
@@ -116,6 +96,14 @@ def test_run_microclimate_subdaily(
     )
 
     lyr_str = fixture_core_components.layer_structure
+
+    # Adjust input data to more realistsic values over this time scale
+    dummy_climate_data_varying_canopy["canopy_evaporation"] = (
+        dummy_climate_data_varying_canopy["canopy_evaporation"] / 86400 * 3600
+    )
+    dummy_climate_data_varying_canopy["transpiration"] = (
+        dummy_climate_data_varying_canopy["transpiration"] / 86400 * 3600
+    )
     result = run_microclimate(
         data=dummy_climate_data_varying_canopy,
         time_index=0,
@@ -125,30 +113,63 @@ def test_run_microclimate_subdaily(
         abiotic_constants=AbioticConsts(),
         core_constants=CoreConsts(),
         pyrealm_const=PyrealmConst(),
+        abiotic_bounds=AbioticSimpleBounds(),
     )
 
-    # Check that all air temperature values fall within a reasonable expected range
-    # Check that all air temperature values fall within a reasonable expected range
-    air_temp_result = result["air_temperature"].isel(
-        layers=lyr_str.index_filled_atmosphere
-    )
+    for var in [
+        "canopy_temperature",
+        "air_temperature",
+        "sensible_heat_flux",
+        "latent_heat_flux",
+    ]:
+        assert var in result
+
+    # Check that values fall within a reasonable expected range
+    soil_temps = result["soil_temperature"].isel(layers=lyr_str.index_all_soil)
 
-    mask = ~np.isnan(
+    # To test with varying canopy layers, need to mask
+    canopy_mask = ~np.isnan(
+        dummy_climate_data_varying_canopy["canopy_temperature"].isel(
+            layers=lyr_str.index_filled_canopy
+        )
+    )
+    atm_mask = ~np.isnan(
         dummy_climate_data_varying_canopy["air_temperature"].isel(
             layers=lyr_str.index_filled_atmosphere
         )
     )
 
+    canopy_temp_result = result["canopy_temperature"].isel(
+        layers=lyr_str.index_filled_canopy
+    )
+    air_temp_result = result["air_temperature"].isel(
+        layers=lyr_str.index_filled_atmosphere
+    )
+    rel_hum_result = result["relative_humidity"].isel(
+        layers=lyr_str.index_filled_atmosphere
+    )
+
     # Use the mask as a DataArray for .where()
-    valid_values = air_temp_result.where(mask)
+    valid_values_can_temp = canopy_temp_result.where(canopy_mask)
+    valid_values_air_temp = air_temp_result.where(atm_mask)
+    valid_values_rel_hum = rel_hum_result.where(atm_mask)
 
     # Now drop the NaNs (i.e., masked values)
-    valid_values_clean = valid_values.dropna(
-        dim="layers", how="any"
-    )  # or "all" if you're doing 2D
+    valid_values_can_temp_clean = valid_values_can_temp.dropna(dim="layers", how="any")
+    valid_values_air_temp_clean = valid_values_air_temp.dropna(dim="layers", how="any")
+    valid_values_rel_hum_clean = valid_values_rel_hum.dropna(dim="layers", how="any")
 
     # Now do the test
-    assert ((valid_values_clean >= 16.0) & (valid_values_clean <= 36.0)).all()
+    assert ((soil_temps >= 18.0) & (soil_temps <= 25.0)).all()
+    assert (
+        (valid_values_can_temp_clean >= 15.0) & (valid_values_can_temp_clean <= 40.0)
+    ).all()
+    assert (
+        (valid_values_air_temp_clean >= 15.0) & (valid_values_air_temp_clean <= 40.0)
+    ).all()
+    assert (
+        (valid_values_rel_hum_clean >= 0.0) & (valid_values_rel_hum_clean <= 100.0)
+    ).all()
 
 
 def test_run_microclimate_minutes(
@@ -161,6 +182,14 @@ def test_run_microclimate_minutes(
     )
 
     lyr_str = fixture_core_components.layer_structure
+    # Adjust input data to more realistsic values over this time scale
+    dummy_climate_data_varying_canopy["canopy_evaporation"] = (
+        dummy_climate_data_varying_canopy["canopy_evaporation"] / 86400 * 60
+    )
+    dummy_climate_data_varying_canopy["transpiration"] = (
+        dummy_climate_data_varying_canopy["transpiration"] / 86400 * 60
+    )
+
     result = run_microclimate(
         data=dummy_climate_data_varying_canopy,
         time_index=0,
@@ -170,27 +199,52 @@ def test_run_microclimate_minutes(
         abiotic_constants=AbioticConsts(),
         core_constants=CoreConsts(),
         pyrealm_const=PyrealmConst(),
+        abiotic_bounds=AbioticSimpleBounds(),
     )
 
-    # Check that all air temperature values fall within a reasonable expected range
-    # Check that all air temperature values fall within a reasonable expected range
-    air_temp_result = result["air_temperature"].isel(
-        layers=lyr_str.index_filled_atmosphere
-    )
+    # Check that values fall within a reasonable expected range
+    soil_temps = result["soil_temperature"].isel(layers=lyr_str.index_all_soil)
 
-    mask = ~np.isnan(
+    # To test with varying canopy layers, need to mask
+    canopy_mask = ~np.isnan(
+        dummy_climate_data_varying_canopy["canopy_temperature"].isel(
+            layers=lyr_str.index_filled_canopy
+        )
+    )
+    atm_mask = ~np.isnan(
         dummy_climate_data_varying_canopy["air_temperature"].isel(
             layers=lyr_str.index_filled_atmosphere
         )
     )
 
+    canopy_temp_result = result["canopy_temperature"].isel(
+        layers=lyr_str.index_filled_canopy
+    )
+    air_temp_result = result["air_temperature"].isel(
+        layers=lyr_str.index_filled_atmosphere
+    )
+    rel_hum_result = result["relative_humidity"].isel(
+        layers=lyr_str.index_filled_atmosphere
+    )
+
     # Use the mask as a DataArray for .where()
-    valid_values = air_temp_result.where(mask)
+    valid_values_can_temp = canopy_temp_result.where(canopy_mask)
+    valid_values_air_temp = air_temp_result.where(atm_mask)
+    valid_values_rel_hum = rel_hum_result.where(atm_mask)
 
     # Now drop the NaNs (i.e., masked values)
-    valid_values_clean = valid_values.dropna(
-        dim="layers", how="any"
-    )  # or "all" if you're doing 2D
+    valid_values_can_temp_clean = valid_values_can_temp.dropna(dim="layers", how="any")
+    valid_values_air_temp_clean = valid_values_air_temp.dropna(dim="layers", how="any")
+    valid_values_rel_hum_clean = valid_values_rel_hum.dropna(dim="layers", how="any")
 
     # Now do the test
-    assert ((valid_values_clean >= 16.0) & (valid_values_clean <= 36.0)).all()
+    assert ((soil_temps >= 18.0) & (soil_temps <= 25.0)).all()
+    assert (
+        (valid_values_can_temp_clean >= 15.0) & (valid_values_can_temp_clean <= 40.0)
+    ).all()
+    assert (
+        (valid_values_air_temp_clean >= 15.0) & (valid_values_air_temp_clean <= 40.0)
+    ).all()
+    assert (
+        (valid_values_rel_hum_clean >= 0.0) & (valid_values_rel_hum_clean <= 100.0)
+    ).all()
diff --git a/tests/models/abiotic/test_wind.py b/tests/models/abiotic/test_wind.py
index fc3cbfc44..98add30d0 100644
--- a/tests/models/abiotic/test_wind.py
+++ b/tests/models/abiotic/test_wind.py
@@ -1,6 +1,7 @@
 """Test module for abiotic.wind.py."""
 
 import numpy as np
+from numpy.testing import assert_allclose
 
 
 def test_calculate_zero_plane_displacement(dummy_climate_data_varying_canopy):
@@ -16,7 +17,7 @@ def test_calculate_zero_plane_displacement(dummy_climate_data_varying_canopy):
         zero_plane_scaling_parameter=7.5,
     )
 
-    np.testing.assert_allclose(result, np.array([0.0, 0.0, 25.86256, 25.86256]))
+    assert_allclose(result, np.array([0.0, 0.0, 25.86256, 25.86256]))
 
 
 def test_calculate_roughness_length_momentum(dummy_climate_data_varying_canopy):
@@ -38,7 +39,7 @@ def test_calculate_roughness_length_momentum(dummy_climate_data_varying_canopy):
         min_roughness_length=0.01,
     )
 
-    np.testing.assert_allclose(
+    assert_allclose(
         result, np.array([0.01, 0.01666, 0.524479, 0.524479]), rtol=1e-3, atol=1e-3
     )
 
@@ -72,7 +73,7 @@ def test_calculate_wind_profile(
         ]
     )
 
-    np.testing.assert_allclose(result, exp_wind, rtol=1e-3, atol=1e-3)
+    assert_allclose(result, exp_wind, rtol=1e-3, atol=1e-3)
 
 
 def test_calculate_friction_velocity(dummy_climate_data_varying_canopy):
@@ -92,7 +93,7 @@ def test_calculate_friction_velocity(dummy_climate_data_varying_canopy):
         von_karman_constant=0.4,
     )
     exp_friction_velocity = np.array([0.080945, 0.085658, 0.099079, 0.099079])
-    np.testing.assert_allclose(result, exp_friction_velocity, rtol=1e-3, atol=1e-3)
+    assert_allclose(result, exp_friction_velocity, rtol=1e-3, atol=1e-3)
 
 
 def test_calculate_ventilation_rate_scalar():
@@ -122,7 +123,7 @@ def test_calculate_ventilation_rate_array():
     expected = np.array([5.0e-02, 1.0e-03, 1.0e03])
 
     result = calculate_ventilation_rate(ra, h)
-    np.testing.assert_allclose(result, expected)
+    assert_allclose(result, expected)
 
 
 def test_calculate_ventilation_rate_zero_denominator():
@@ -158,12 +159,58 @@ def test_calculate_mixing_coefficients():
 
     assert result.shape == layer_midpoints.shape
     assert np.all(result >= 0)
-    np.testing.assert_allclose(result, expected, rtol=1e-6)
+    assert_allclose(result, expected, rtol=1e-6)
+
+
+def test_clamp_variable_within_limits():
+    """Test clamping of variable within limits."""
+
+    from virtual_ecosystem.models.abiotic.wind import clamp_variable_within_limits
+
+    # Set limits
+    limits = (4, 6)
+
+    # Cells with mix of issues: undershoot and overshoot, can be absorbed by a single
+    # cell vs need to spread further up into canopy, empty layers and complete, and mix
+    # of under and overshoot. Two cases explicitly test the edge case where undershoot
+    # and overshoot cannot be absorbed within the canopy, leading to out of limit values
+    # in the top layer.
+
+    n = np.nan
+    variable = np.array(
+        [
+            [5, 5, 5, 5, 5, 4, 5, 5, 5, 5, 5, 5],
+            [5, 5, 5, 5, 5, 4, 5, 5, 5, 5, 5, 8],
+            [5, 5, 5, 5, 3, 4, 5, 5, 5, 7, 5, 2],
+            [5, n, 5, n, n, 4, 4, 5, n, n, 7, 7],
+            [5, 5, 3, 3, 5, 0, 2, 7, 7, 7, 9, 3],
+        ]
+    )
+
+    variable_expected = np.array(
+        [
+            [5, 5, 5, 5, 5, 0, 5, 5, 5, 6, 7, 5],
+            [5, 5, 5, 5, 4, 4, 4, 5, 5, 6, 6, 6],
+            [5, 5, 5, 4, 4, 4, 4, 5, 6, 6, 6, 4],
+            [5, n, 4, n, n, 4, 4, 6, n, n, 6, 6],
+            [5, 5, 4, 4, 5, 4, 4, 6, 6, 6, 6, 4],
+        ]
+    )
+
+    clamped_variable = clamp_variable_within_limits(variable=variable, limits=limits)
+
+    assert_allclose(clamped_variable, variable_expected)
+
+    # sanity checks: the column sums should be maintained.
+    assert_allclose(variable.sum(axis=0), clamped_variable.sum(axis=0))
 
 
 def test_mix_and_ventilate(dummy_climate_data_varying_canopy, fixture_core_components):
-    """Test mixing and ventilation."""
+    """Test mixing and ventilation within bounds."""
 
+    from virtual_ecosystem.models.abiotic.abiotic_tools import (
+        compute_layer_thickness_for_varying_canopy,
+    )
     from virtual_ecosystem.models.abiotic.wind import (
         mix_and_ventilate,
     )
@@ -173,38 +220,50 @@ def test_mix_and_ventilate(dummy_climate_data_varying_canopy, fixture_core_compo
     atm_index = lystr.index_filled_atmosphere
 
     heights = data["layer_heights"][atm_index].to_numpy()
-    heights_with_base = np.vstack([heights, np.zeros(heights.shape[1])])
-    above_ground_layer_thickness = -np.diff(heights_with_base, axis=0)
+    above_ground_layer_thickness = compute_layer_thickness_for_varying_canopy(
+        heights=heights
+    )
 
     mixing_coefficient = np.array(
         [
             [0.001, 0.001, 0.001, 0.001],
-            [0.005, 0.005, np.nan, np.nan],
+            [0.005, 0.005, 0.005, np.nan],
             [0.01, 0.01, np.nan, np.nan],
             [0.001, np.nan, np.nan, np.nan],
             [0.012, 0.012, 0.012, 0.012],
         ]
     )
-    ventilation_rate = np.array([0.01, 0.01, 0.01, 0.01])
+    ventilation_rate = np.array([0.001, 0.001, 0.001, 0.001])
+
+    input_humidity = np.array(
+        [
+            [95.0, 95.0, 95.0, 95.0],
+            [100.0, 100.0, 100.0, 100.0],
+            [100.0, 100.0, 100.0, 100.0],
+            [90.0, 90.0, 90.0, 90.0],
+            [100.0, 100.0, 100.0, 100.0],
+        ],
+    )
 
     exp_result = np.array(
         [
-            [77.700816, 77.700816, 77.700816, 77.700816],
-            [90.341644, 90.341644, 90.341644, 90.341644],
-            [92.233778, np.nan, np.nan, np.nan],
-            [96.503298, np.nan, np.nan, np.nan],
+            [94.82, 94.82, 94.979866, 94.979866],
+            [100.0, 100.0, 100.0, 100.0],
+            [99.64, 98.909118, np.nan, np.nan],
+            [98.08080808, np.nan, np.nan, np.nan],
             [100.0, 100.0, 100.0, 100.0],
         ]
     )
 
     result = mix_and_ventilate(
-        input_variable=data["relative_humidity"][atm_index].to_numpy(),
+        input_variable=input_humidity,
         layer_thickness=above_ground_layer_thickness,
         mixing_coefficient=mixing_coefficient,
         ventilation_rate=ventilation_rate,
+        limits=(0, 100),
         time_interval=3600.0,
     )
-    np.testing.assert_allclose(result, exp_result, rtol=1e-6, atol=1e-6)
+    assert_allclose(result, exp_result, rtol=1e-6, atol=1e-6)
 
 
 def test_advect_from_toplayer():
@@ -231,7 +290,7 @@ def test_advect_from_toplayer():
         time_interval=time_interval,
     )
 
-    np.testing.assert_allclose(result, expected_specific_humidity)
+    assert_allclose(result, expected_specific_humidity)
 
 
 def test_calculate_aerodynamic_resistance(
@@ -261,4 +320,4 @@ def test_calculate_aerodynamic_resistance(
         wind_speed=np.array([1.0, 2.0, 0.5, 0.01]),
         von_karman_constant=0.4,
     )
-    np.testing.assert_allclose(result, exp_ra, rtol=1e-3, atol=1e-3)
+    assert_allclose(result, exp_ra, rtol=1e-3, atol=1e-3)
diff --git a/virtual_ecosystem/models/abiotic/abiotic_model.py b/virtual_ecosystem/models/abiotic/abiotic_model.py
index aa7f30de5..ed916c150 100644
--- a/virtual_ecosystem/models/abiotic/abiotic_model.py
+++ b/virtual_ecosystem/models/abiotic/abiotic_model.py
@@ -238,6 +238,7 @@ def _update(self, time_index: int, **kwargs: Any) -> None:
             abiotic_constants=self.model_constants,
             core_constants=self.core_constants,
             pyrealm_const=PyrealmConst,
+            abiotic_bounds=AbioticSimpleBounds(),
         )
 
         self.data.add_from_dict(output_dict=update_dict)
diff --git a/virtual_ecosystem/models/abiotic/abiotic_tools.py b/virtual_ecosystem/models/abiotic/abiotic_tools.py
index d6dd3dae3..598b59c9d 100644
--- a/virtual_ecosystem/models/abiotic/abiotic_tools.py
+++ b/virtual_ecosystem/models/abiotic/abiotic_tools.py
@@ -226,3 +226,43 @@ def compute_layer_thickness_for_varying_canopy(
                 thickness[row, col] = current - 0.0
 
     return thickness
+
+
+def calculate_specific_humidity(
+    air_temperature: NDArray[np.floating],
+    relative_humidity: NDArray[np.floating],
+    atmospheric_pressure: NDArray[np.floating],
+    molecular_weight_ratio_water_to_dry_air: float,
+    pyrealm_const: PyrealmConst,
+) -> NDArray[np.floating]:
+    """Calculate specific humidity.
+
+    Args:
+        air_temperature: Air temperature, [C]
+        relative_humidity: Relative humidity, [%]
+        atmospheric_pressure: Atmospheric pressure, [kPa]
+        molecular_weight_ratio_water_to_dry_air: The ratio of the molar mass of water
+            vapour to the molar mass of dry air
+        pyrealm_const: Pyrealm constants
+
+    Returns:
+        Specific humidity (kg/kg)
+    """
+    # Saturation vapor pressure
+    saturation_vapour_pressure = calc_vp_sat(
+        ta=air_temperature,
+        core_const=pyrealm_const,
+    )
+
+    # Actual vapor pressure (hPa)
+    actual_vapour_pressure = (relative_humidity / 100.0) * saturation_vapour_pressure
+
+    # Specific humidity formula
+    specific_humidity = (
+        molecular_weight_ratio_water_to_dry_air * actual_vapour_pressure
+    ) / (
+        atmospheric_pressure
+        - ((1 - molecular_weight_ratio_water_to_dry_air) * actual_vapour_pressure)
+    )
+
+    return specific_humidity
diff --git a/virtual_ecosystem/models/abiotic/energy_balance.py b/virtual_ecosystem/models/abiotic/energy_balance.py
index 8e6bba426..7cf1bed31 100644
--- a/virtual_ecosystem/models/abiotic/energy_balance.py
+++ b/virtual_ecosystem/models/abiotic/energy_balance.py
@@ -326,10 +326,10 @@ def calculate_energy_balance_residual(
     # Energy balance residual, [W m-2]
     energy_balance_residual = (
         absorbed_radiation_canopy
-        - np.abs(longwave_emission_canopy)
-        - np.abs(sensible_heat_flux_canopy)
-        - np.abs(latent_heat_flux_canopy)
-        # - np.abs(absorption_par)
+        - longwave_emission_canopy
+        + sensible_heat_flux_canopy
+        + latent_heat_flux_canopy
+        # - absorption_par
     )
 
     if return_fluxes:
@@ -469,7 +469,9 @@ def residual_func(canopy_temp_scalar):
                     aerodynamic_resistance=np.array(
                         [[aerodynamic_resistance[i, j]]], dtype=np.float64
                     ),
-                    latent_heat_vapourisation=latent_heat_vapourisation,
+                    latent_heat_vapourisation=np.array(
+                        [[latent_heat_vapourisation[i, j]]], dtype=np.float64
+                    ),
                     leaf_emissivity=emissivity_leaf,
                     stefan_boltzmann_constant=stefan_boltzmann_constant,
                     zero_Celsius=zero_Celsius,
@@ -569,11 +571,11 @@ def update_humidity_vpd(
     atmospheric_pressure: NDArray[np.floating],
     density_air: NDArray[np.floating],
     mixing_coefficient: NDArray[np.floating],
-    ventilation_rate: float | NDArray[np.floating],
-    wind_speed: NDArray[np.floating],
+    ventilation_rate: NDArray[np.floating],
     molecular_weight_ratio_water_to_dry_air: float,
     dry_air_factor: float,
     cell_area: float,
+    limits: tuple[float, float],
     time_interval: float,
 ) -> dict[str, NDArray[np.floating]]:
     """Update specific humidity and vapour pressure deficit for a multilayer canopy.
@@ -591,11 +593,11 @@ def update_humidity_vpd(
         density_air: Density of air, [kg m-3]
         mixing_coefficient: Turbulent mixing coefficient, [m2 s-1]
         ventilation_rate: Ventilation rate, [s-1]
-        wind_speed: Horizontal wind speed above canopy, [m s-1]
         molecular_weight_ratio_water_to_dry_air: Molecular weight ratio of water to dry
             air, dimensionless
         dry_air_factor: Complement of water_to_air_mass_ratio, accounting for dry air
         cell_area: Grid cell area, [m2]
+        limits: Realistic bounds of specific humidity
         time_interval: Time interval, [s]
 
     Returns:
@@ -630,19 +632,20 @@ def update_humidity_vpd(
         layer_thickness=layer_thickness,
         mixing_coefficient=mixing_coefficient,
         ventilation_rate=ventilation_rate,
+        limits=limits,
         time_interval=time_interval,
     )
 
-    # Advection
-    specific_humidity_advected = wind.advect_water_from_toplayer(
-        specific_humidity=specific_humidity_updated[0],
-        layer_thickness=layer_thickness[0],
-        density_air=density_air[0],
-        wind_speed=wind_speed,
-        characteristic_length=np.sqrt(cell_area),
-        time_interval=time_interval,
-    )
-    specific_humidity_updated[0] = specific_humidity_advected
+    # NOTE Advection not implemented as everything is removed with time interval > 1h
+    # specific_humidity_advected = wind.advect_water_from_toplayer(
+    #     specific_humidity=specific_humidity_updated[0],
+    #     layer_thickness=layer_thickness[0],
+    #     density_air=density_air[0],
+    #     wind_speed=wind_speed,
+    #     characteristic_length=np.sqrt(cell_area),
+    #     time_interval=time_interval,
+    # )
+    # specific_humidity_updated[0] = specific_humidity_advected
 
     # Vapour pressure [kPa]
     vapour_pressure_updated = (specific_humidity_updated * atmospheric_pressure) / (
@@ -651,7 +654,7 @@ def update_humidity_vpd(
     )
 
     # Ensure vapor pressure doesn't exceed the saturated vapor pressure
-    # NOTE we need to make sure that we do not loose water here
+    # TODO we need to make sure that we do not loose water here
     vapour_pressure_updated = np.minimum(
         vapour_pressure_updated, saturated_vapour_pressure
     )
diff --git a/virtual_ecosystem/models/abiotic/microclimate.py b/virtual_ecosystem/models/abiotic/microclimate.py
index f1470aec5..323c0c6dc 100644
--- a/virtual_ecosystem/models/abiotic/microclimate.py
+++ b/virtual_ecosystem/models/abiotic/microclimate.py
@@ -13,6 +13,7 @@
 from virtual_ecosystem.core.logger import LOGGER
 from virtual_ecosystem.models.abiotic import abiotic_tools, energy_balance, wind
 from virtual_ecosystem.models.abiotic.constants import AbioticConsts
+from virtual_ecosystem.models.abiotic_simple.constants import AbioticSimpleBounds
 
 
 def run_microclimate(
@@ -24,6 +25,7 @@ def run_microclimate(
     abiotic_constants: AbioticConsts,
     core_constants: CoreConsts,
     pyrealm_const: PyrealmConst,
+    abiotic_bounds: AbioticSimpleBounds,
 ) -> dict[str, DataArray]:
     """Run microclimate model.
 
@@ -43,6 +45,7 @@ def run_microclimate(
         abiotic_constants: Set of constants for abiotic model
         core_constants: Set of constants that are shared across all models
         pyrealm_const: Set of constants from pyrealm
+        abiotic_bounds: Bounds for vertical mixing of atmospheric variables
 
     Returns:
         dictionary with updated microclimate variables
@@ -169,7 +172,7 @@ def run_microclimate(
     #  Ventilation rate above canopy, [s-1]
     ventilation_rate = wind.calculate_ventilation_rate(
         aerodynamic_resistance=aerodynamic_resistance_canopy[0],
-        characteristic_height=canopy_height,
+        characteristic_height=canopy_height + zero_plane_displacement,
     )
 
     # -------------------------------------------------------------------------
@@ -201,8 +204,10 @@ def run_microclimate(
     # If hourly input data is provided, iterate over day, else equilibrium assumption
     hourly_time_interval = max(int(time_interval / core_constants.seconds_to_hour), 1)
 
-    # TODO enable daily input in data object and select time index
-    if time_interval <= core_constants.seconds_to_day:
+    # TODO Run diurnal cycle
+    # enable hourly input in data object, select by time index, iterate over hours,
+    # and return averages/min/max over the VE time interval, similar to hydrology
+    if core_constants.seconds_to_hour <= time_interval <= core_constants.seconds_to_day:
         iteration = hourly_time_interval
 
     else:
@@ -340,9 +345,18 @@ def run_microclimate(
             layer_thickness=above_ground_layer_thickness,
             ventilation_rate=ventilation_rate,
             mixing_coefficient=mixing_coefficient,
-            time_interval=1.0,  # TODO core_constants.seconds_to_hour,
+            limits=abiotic_bounds.air_temperature[:2],
+            time_interval=core_constants.seconds_to_hour,
         )
 
+        # NOTE Advection not implemented as everything is removed with time interval>1h
+        # advection_rate = (
+        #   data["wind_speed_ref"].isel(time_index=time_index).to_numpy()
+        #   / np.sqrt(cell_area)
+        # )
+        # advected_fraction = np.clip(advection_rate * time_interval, 0, 1)
+        # all_air_temperature[0] -=all_air_temperature[0] *advected_fraction
+
         # Update atmospheric humidity/VPD
         # Saturated vapour pressure of air, [kPa]
         saturated_vapour_pressure_air = calc_vp_sat(
@@ -350,38 +364,45 @@ def run_microclimate(
             core_const=PyrealmConst(),
         )
 
-        #  Actual vapour pressure of air, [kPa]
-        actual_vapour_pressure_air = abiotic_tools.calculate_actual_vapour_pressure(
-            air_temperature=DataArray(all_air_temperature),
-            relative_humidity=DataArray(relative_humidity),
-            pyrealm_const=PyrealmConst,
+        # Specific humidity of air, [kg kg-1]
+        specific_humidity_air = abiotic_tools.calculate_specific_humidity(
+            air_temperature=all_air_temperature,
+            relative_humidity=relative_humidity,
+            atmospheric_pressure=atmospheric_pressure,
+            molecular_weight_ratio_water_to_dry_air=(
+                core_constants.molecular_weight_ratio_water_to_dry_air
+            ),
+            pyrealm_const=PyrealmConst(),
         )
 
-        # Specific humidity of air, [kg kg-1] TODO external function
-        specific_humidity_air = (
+        # Calculate specific humidity at saturation
+        mixing_ratio_saturation = (
             core_constants.molecular_weight_ratio_water_to_dry_air
-            * actual_vapour_pressure_air
-        ) / (atmospheric_pressure - actual_vapour_pressure_air)
+            * saturated_vapour_pressure_air
+            / (atmospheric_pressure - saturated_vapour_pressure_air)
+        )
+        max_specific_humidity = mixing_ratio_saturation / (1 + mixing_ratio_saturation)
 
+        # Update atmospheric humidity variables
         new_atmospheric_humidity_vars = energy_balance.update_humidity_vpd(
             evapotranspiration=evapotranspiration[
                 layer_structure.index_filled_canopy
             ].to_numpy(),
             soil_evaporation=data["soil_evaporation"].to_numpy(),
             saturated_vapour_pressure=saturated_vapour_pressure_air,
-            specific_humidity=specific_humidity_air.to_numpy(),
+            specific_humidity=specific_humidity_air,
             layer_thickness=above_ground_layer_thickness,
             atmospheric_pressure=atmospheric_pressure,
             density_air=density_air,
             mixing_coefficient=mixing_coefficient,
             ventilation_rate=ventilation_rate,
-            wind_speed=data["wind_speed_ref"].isel(time_index=time_index).to_numpy(),
             molecular_weight_ratio_water_to_dry_air=(
                 core_constants.molecular_weight_ratio_water_to_dry_air
             ),
             dry_air_factor=abiotic_constants.dry_air_factor,
             cell_area=cell_area,
-            time_interval=1.0,  # TODO core_constants.seconds_to_hour,
+            limits=(0, max_specific_humidity[0]),  # TODO make layer specific
+            time_interval=core_constants.seconds_to_hour,
         )
 
         relative_humidity = new_atmospheric_humidity_vars["relative_humidity"]
@@ -450,7 +471,6 @@ def run_microclimate(
     output["density_air"] = density_air_out
 
     # Combine longwave emission in one variable
-    # Assumption: accumulated emission in time interval based on accumulated input
     longwave_emission = layer_structure.from_template()
     longwave_emission[layer_structure.index_filled_canopy] = longwave_emission_canopy
     longwave_emission[layer_structure.index_topsoil_scalar] = longwave_emission_soil
@@ -472,18 +492,18 @@ def run_microclimate(
         latent_heat_vapourisation
     )
     output["latent_heat_vapourisation"] = latent_heat_vapourisation_out
-    # Combine sensible heat flux in one variable, TODO consider time interval
+
+    # Combine sensible heat flux in one variable
     sensible_heat_flux = layer_structure.from_template()
     sensible_heat_flux[layer_structure.index_filled_canopy] = sensible_heat_flux_canopy
     sensible_heat_flux[layer_structure.index_topsoil_scalar] = sensible_heat_flux_soil
-    output["sensible_heat_flux"] = sensible_heat_flux  # * time_interval
+    output["sensible_heat_flux"] = sensible_heat_flux
 
-    # Combine latent heat flux in one variable, TODO consider time interval
-    # TODO adjust to model timestep, currently per second
+    # Combine latent heat flux in one variable
     latent_heat_flux = layer_structure.from_template()
     latent_heat_flux[layer_structure.index_filled_canopy] = latent_heat_flux_canopy
     latent_heat_flux[layer_structure.index_topsoil_scalar] = latent_heat_flux_soil
-    output["latent_heat_flux"] = latent_heat_flux  # * time_interval
+    output["latent_heat_flux"] = latent_heat_flux
 
     soil_temperature_out = layer_structure.from_template()
     soil_temperature_out[layer_structure.index_all_soil] = soil_temperature
@@ -499,7 +519,7 @@ def run_microclimate(
     canopy_temperature_out[layer_structure.index_filled_canopy] = canopy_temperature
     output["canopy_temperature"] = canopy_temperature_out
 
-    # TODO check dimensions write humidity/VPD
+    # Write humidity/VPD
     for var in ["relative_humidity", "vapour_pressure", "vapour_pressure_deficit"]:
         var_out = layer_structure.from_template()
         var_out[layer_structure.index_filled_atmosphere] = (
diff --git a/virtual_ecosystem/models/abiotic/wind.py b/virtual_ecosystem/models/abiotic/wind.py
index a77e75b91..aeb51498c 100644
--- a/virtual_ecosystem/models/abiotic/wind.py
+++ b/virtual_ecosystem/models/abiotic/wind.py
@@ -211,7 +211,7 @@ def calculate_friction_velocity(
 def calculate_ventilation_rate(
     aerodynamic_resistance: float | NDArray[np.floating],
     characteristic_height: float | NDArray[np.floating],
-) -> float | NDArray[np.floating]:
+) -> NDArray[np.floating]:
     """Calculate ventilation rate from the top of the canopy to atmosphere above.
 
     This function calculates the rate of water and heat exchange between the top of the
@@ -219,7 +219,8 @@ def calculate_ventilation_rate(
 
     Args:
         aerodynamic_resistance: Aerodynamic resistance, [s m-1]
-        characteristic_height: Vertical scale of exchange, [m]
+        characteristic_height: Vertical scale of exchange, typically canopy height +
+            zero plane displacement height [m]
 
     Returns:
         Ventilation rate [s-1]
@@ -248,7 +249,7 @@ def calculate_mixing_coefficients_canopy(
 
     .. math::
 
-        k_{H,M}(z)=\kappa u_{*}z(1-z h_c)^{2}
+        k_{H,M}(z)=\kappa u_{*}z(\frac{1-z}{h_c})^{2}
 
     where :math:`\kappa` is the von Karman constant (dimensionless), :math:`u_{*}` is
     the friction velocity (m s-1), :math:`z` is the height (m) for which coefficients
@@ -278,49 +279,143 @@ def calculate_mixing_coefficients_canopy(
     return mixing_coefficients
 
 
+def clamp_variable_within_limits(
+    variable: NDArray[np.floating], limits: tuple[float, float]
+) -> NDArray[np.floating]:
+    """Clamp an array of canopy data within limits.
+
+    This function iterates from the bottom of the canopy, clamping the values of the
+    input array within the limits. When a value is altered by clamping, the residual is
+    added to the layer above to maintain the variable total within cells. Residual
+    values may be redistributed across multiple layers and empty values (representing
+    unoccupied canopy layers) are skipped.
+
+    Note:
+        If the vertical layers cannot absorb all of the accumulated residuals without
+        themselves being clamped, then the values in the top layer can still fall
+        outside the clamping limits.
+
+    Args:
+        variable: A numpy array containing canopy data.
+        limits: A tuple giving the upper and lower bounds within which to clamp the data
+    """
+
+    # Get a map of nan values and initialise the out_of_limits array
+    out_of_limits = np.zeros_like(variable[0])
+    nan_map = np.isnan(variable)
+    n_layers = variable.shape[0]
+
+    # Loop up from the row index of lowest layer, stopping before the top layer
+    for layer in np.arange(n_layers - 1, 0, -1):
+        # Calculate the clamped values for the current layer
+        in_limits = np.clip(variable[layer], *limits)
+
+        # Add under and overshoots to the out_of_limits array, trapping cells that
+        # contain no vegetation in the layer (np.nan)
+        out_of_limits += np.where(nan_map[layer], 0, variable[layer] - in_limits)
+
+        # Set the clamped data in the current layer
+        variable[layer] = in_limits
+
+        # Add out of limits to the layer above
+        variable[layer - 1] += out_of_limits
+        # Update out_of_limits
+        # - np.nan cells carry over the current out_of_limits total
+        # - otherwise the out_of_limits has been set into the layer above, so is zeroed
+        out_of_limits = np.where(nan_map[layer - 1], out_of_limits, 0)
+
+    return variable
+
+
 def mix_and_ventilate(
     input_variable: NDArray[np.floating],
     layer_thickness: NDArray[np.floating],
     mixing_coefficient: NDArray[np.floating],
-    ventilation_rate: float | NDArray[np.floating],
+    ventilation_rate: NDArray[np.floating],
+    limits: tuple[float, float],
     time_interval: float,
 ) -> NDArray[np.floating]:
-    """Mix and ventilate vertically.
+    """Apply vertical mixing and top-layer ventilation across multiple vertical layers.
+
+    This function simulates diffusion-like mixing between vertical layers based on local
+    gradients of atmospheric variables (e.g. temperature, relative humidity) and
+    layer-specific mixing coefficients. For each internal layer (excluding the top and
+    bottom), it computes upward and downward fluxes using the nearest valid
+    (finite) values above and below, respectively. The fluxes are scaled by the layer
+    thickness and applied to update the variable.
+
+    Additionally, the function applies a ventilation adjustment to the top layer of each
+    column, representing heat or water exchange with the  above the canopy. This is
+    based on the difference between the top and next valid layer, scaled by a
+    user-provided ventilation rate, with optional limits to prevent overcorrection or
+    negative concentrations.
 
-    This function takes an atmospheric variable such as temperature or specific
-    humidity, mixed vertically between layers and ventilates at the top of the canopy.
+    Advection is currently not implemented as everything is removed with time interval
+    > 1h.
 
     Args:
         input_variable: Input variable for all true atmospheric layers
         layer_thickness: Layer thickness, [m]
         mixing_coefficient: Turbulent mixing coefficients for canopy, [m2 s-1]
         ventilation_rate: Ventilation rate, [s-1]
+        limits: Upper and lower limit for input variable, avoid overshoot when mixing
         time_interval: Time interval, [s]
 
     Returns:
         Vertically mixed input variable
     """
 
+    # Copy the input to update with mixing
     input_variable_mixed = input_variable.copy()
-    for i in range(2, len(input_variable) - 1):
-        flux_up = (
-            mixing_coefficient[i - 1]
-            * (input_variable[i - 1] - input_variable[i])
-            / layer_thickness[i]
-        )
-        flux_down = (
-            mixing_coefficient[i + 1]
-            * (input_variable[i + 1] - input_variable[i])
-            / layer_thickness[i]
-        )
-        change_input_variable = (
-            (flux_up + flux_down) * time_interval / layer_thickness[i]
-        )
-        input_variable_mixed[i] += change_input_variable
-
-    # Ventilation at top
-    input_variable_mixed[0] += (
-        ventilation_rate * (input_variable[0] - input_variable[1]) * time_interval
+    n_layers, n_cells = input_variable_mixed.shape
+
+    # Extract the layer thickness and current value for the canopy layers, excluding the
+    # surface layer and above canopy values
+    value_canopy = input_variable[1:-1]
+    canopy_layer_thickness = layer_thickness[1:-1]
+
+    # Get the value and mixing coefficients from the layer above.
+    value_above = input_variable[0:-2]
+    mix_above = mixing_coefficient[0:-2]
+
+    # Get the value and mixing coefficients from the layer above.
+    value_below = input_variable[2:]
+    mix_below = mixing_coefficient[2:]
+
+    # In-fill missing below values from the surface layer
+    value_below = np.where(np.isnan(value_below), input_variable[-1], value_below)
+    mix_below = np.where(np.isnan(mix_below), mix_below[-1], mix_below)
+
+    # Calculate fluxes (positive upward, negative downward) and update variable
+    flux_up = mix_above * (value_above - value_canopy) / canopy_layer_thickness
+    flux_down = mix_below * (value_below - value_canopy) / canopy_layer_thickness
+
+    value_change = (flux_up + flux_down) * time_interval / canopy_layer_thickness
+    input_variable_mixed[1:-1] += value_change
+
+    # Calculate ventilation using the value from the highest layer, which is either the
+    # top canopy layer or the surface layer if no canopy is present
+    vent_below = np.where(
+        np.isnan(input_variable_mixed[1]),
+        input_variable_mixed[-1],
+        input_variable_mixed[1],
+    )
+
+    # Get the ventilation delta and change over time
+    delta_vent = input_variable_mixed[0] - vent_below
+    vent_change = ventilation_rate * delta_vent * time_interval
+
+    # Clip the maximum ventilation rate to 10%
+    vent_max_change = 0.1 * abs(input_variable[0])
+    vent_change = np.clip(vent_change, -vent_max_change, vent_max_change)
+
+    # Update the current values
+    input_variable_mixed[0] += vent_change
+    input_variable_mixed[1] -= vent_change
+
+    # Redistribute overshoot/undershoot
+    input_variable_mixed = clamp_variable_within_limits(
+        variable=input_variable_mixed, limits=limits
     )
 
     return input_variable_mixed
@@ -337,9 +432,9 @@ def advect_water_from_toplayer(
     """Remove water by advection from above canopy layer.
 
     Args:
-        specific_humidity: Specific humidity in each layer, [kg kg-1]
-        layer_thickness: Thickness of each layer, [m]
-        density_air: Air density in each layer, [kg m-3]
+        specific_humidity: Specific humidity in top layer, [kg kg-1]
+        layer_thickness: Thickness of top layer, [m]
+        density_air: Air density in top layer, [kg m-3]
         wind_speed: Horizontal wind speed above canopy, [m s-1]
         characteristic_length: Horizontal length scale of the grid cell, [m]
         time_interval: Time step, [s]