diff --git a/docs/source/refs.bib b/docs/source/refs.bib
index c915bc45c..40a07bb00 100644
--- a/docs/source/refs.bib
+++ b/docs/source/refs.bib
@@ -1,4 +1,30 @@
 
+@book{holmes_wind_2019,
+author = {Holmes, John D, and Seifu Bekele},
+year = {2020},
+title = {Wind Loading of Structures. Fourth edition.},
+publisher={Boca Raton: CRC Press}
+}
+@article{jansson_coupled_2004,
+author = {Jansson, Per-Erik and Karlberg, Lisa},
+year = {2004},
+month = {01},
+pages = {},
+title = {Coupled heat and mass transfer model for soil-plant-atmosphere systems. Royal Institute of Technology, Department of Civil and Environmental Engineering},
+journal = {Stockholm}
+}
+@Article{wolfe_forest_2011,
+AUTHOR = {Wolfe, G. M. and Thornton, J. A. and McKay, M. and Goldstein, A. H.},
+TITLE = {Forest-atmosphere exchange of ozone: sensitivity to very reactive biogenic VOC emissions and implications for in-canopy photochemistry},
+JOURNAL = {Atmospheric Chemistry and Physics},
+VOLUME = {11},
+YEAR = {2011},
+NUMBER = {15},
+PAGES = {7875--7891},
+URL = {https://acp.copernicus.org/articles/11/7875/2011/},
+DOI = {10.5194/acp-11-7875-2011}
+}
+
 @article{raupach_coherent_1996,
 author = {Raupach, Michael and Finnigan, John and Brunet, Y.},
 year = {1996},
diff --git a/docs/source/virtual_ecosystem/implementation/abiotic_implementation.md b/docs/source/virtual_ecosystem/implementation/abiotic_implementation.md
index fe786dc03..ef3018d4b 100644
--- a/docs/source/virtual_ecosystem/implementation/abiotic_implementation.md
+++ b/docs/source/virtual_ecosystem/implementation/abiotic_implementation.md
@@ -57,13 +57,13 @@ display_markdown(
 ## Model overview
 
 The exchange of energy between the Earth's surface or canopy and the surrounding
-atmosphere involves five important processes:
+atmosphere involves five important categories of processes:
 
 - *Absorption* and *emission* of electromagnetic radiation by the surface/canopy
 - *Thermal conduction* of heat energy within the ground
 - *Turbulent transfer* of heat energy towards or away from the surface within the
   atmosphere
-- *Evaporation* and *condensation* of water
+- *Evaporation*, *transpiration*, and *condensation* of water
 - *Primary productivity*
 
 Each of these processes can be associated with an energy flux density, which is the rate
@@ -72,14 +72,14 @@ of transfer of energy normal to a surface of unit area (in $\mathrm{W\,m^{-2}}$)
 The energy balance of a surface layer of finite depth and unit horizontal area can be
 written as:
 
-$$\frac{dQ}{dt} = R_N - G - H - \lambda E (- PP)$$
+$$\frac{dQ}{dt} = R_n - G - H - \lambda E (- PP)$$
 
 where:
 
 $Q$:
 Total heat energy stored in the surface layer.
 
-$R_N$:
+$R_n$:
 Net surface irradiance (commonly referred to as the net radiation). It
 represents the gain of energy by the surface from radiation. It is a positive number
 when it is towards the surface.
@@ -87,7 +87,7 @@ when it is towards the surface.
 $G$:
 Ground Heat Flux. It is the loss of energy by heat conduction through the
 lower boundary. It is a positive number when it is directed away from the surface into
-ground. The value at the surface is denoted G0.
+ground. The value at the surface is denoted $G_{0}$.
 
 $H$:
 Sensible Heat Flux. It represents the loss of energy by the
@@ -96,23 +96,36 @@ away from the surface into the atmosphere.
 
 $\lambda E$:
 Latent Heat Flux. It represents a loss of energy from the
-surface due to evaporation. ($\lambda$  is the specific latent heat of evaporation,
+surface due to evaporation and/or transpiration. ($\lambda$  is the specific latent heat
+of evaporation,
 units $\mathrm{J\,kg^{-1}}$ and E is the evaporation rate, with units
 $\mathrm{kg\,m^{-2}\,s^{-1}}$).
 
 $PP$:
 Primary productivity, represents the energy that plants use to photosynthesize.
 
+```{note}
+Calculating abiotic processes at coarse time scales can lead to inaccuracies, so the
+abiotic model uses an internal hourly time step. If the Virtual Ecosystem is run with a
+coarser update interval, the abiotic model first partitions the input data into hourly
+values where required, then simulates a single representative hour for the entire
+interval.
+
+The outputs returned by the abiotic model are therefore equilibrium values for that
+representative hour. A planned future improvement is to allow true hourly input, so the
+model can capture full diurnal cycles and return time-averaged values of key variables.
+```
+
 ### Net radiation
 
 The current representation of the radiation balance is limited to the reflection and
 absorption of direct downward shortwave radiation and the emission of longwave radiation
 as part of the surface energy balance.
 
-The net radiation $R_N$ ($\mathrm{W\,m^{-2}}$) at the leaf or soil surface is
+The net radiation $R_n$ ($\mathrm{W\,m^{-2}}$) at the leaf or soil surface is
 calculated as:
 
-$$R_N = S_0 \cdot (1 - \alpha) - \epsilon_{s} \sigma T^{4}$$
+$$R_n = S_0 \cdot (1 - \alpha) - \epsilon_{s} \sigma T^{4}$$
 
 where:
 
@@ -147,28 +160,29 @@ soil surface by partitioning net radiation $R_N$ into different fluxes.
 
 The **sensible heat flux** from the soil surface is given by:
 
-$$H_{S} = \frac {\rho_{air} C_{air} (T_{S} - T_{A})}{r_{A}}$$
+$$H_{s} = \frac {\rho_{a} c_{p} (T_{s} - T_{a})}{r_{a}}$$
 
 where:
 
-$T_S$:
+$T_s$:
 Soil surface temperature (°C)
 
-$T_A$:
+$T_a$:
 Air temperature in the bottom atmospheric layer (°C)
 
-$r_A$:
+$r_a$:
 Aerodynamic resistance of the soil surface ($\mathrm{s\,m^{-1}}$)
 
-$\rho_{air}$:
+$\rho_{a}$:
 Air density ($\mathrm{kg\,m^{-3}}$)
 
-$C_{air}$:
-Specific heat capacity of air ($\mathrm{J\,kg^{-1}\,K^{-1}}$)
+$c_{p}$:
+Specific heat capacity of air at constant pressure ($\mathrm{J\,kg^{-1}\,K^{-1}}$)
 
-The aerodynamic resistance of the soil surface is given by:
+The aerodynamic resistance of the soil surface is given by
+{cite:p}`barton_parameterization_1979`:
 
-$$r_{A} = \frac {C_{E}}{u}$$
+$$r_{a} = \frac{1}{C_{E} u}$$
 
 where:
 
@@ -181,9 +195,10 @@ Drag coefficient for evaporation (–)
 The **latent heat flux** is derived by conversion of surface evaporation as
 calculated by the hydrology model.
 
-The **ground heat flux** is calculated as the residual of the energy balance:
+The **ground heat flux** is calculated as the residual of the energy balance at the
+soil surface:
 
-$$G = R_N - H_S - \lambda E_S$$
+$$G = R_n - H_s - \lambda E_s$$
 
 ### Soil temperature update
 
@@ -197,7 +212,7 @@ each soil depth.
 The **soil thermal diffusivity** $\alpha$ ($\mathrm{m^{2}\,s^{-1}}$) determines the rate
 at which heat is conducted through the soil. It is defined as:
 
-$$\alpha = \frac{k}{\rho c}$$
+$$\alpha = \frac{k}{\rho_s c_s}$$
 
 where:
 
@@ -205,16 +220,17 @@ $k$:
 Soil thermal conductivity ($\mathrm{W\,m^{-1}\,K^{-1}}$), indicating how
   easily heat moves through soil
 
-$\rho$:
-Soil bulk density ($\mathrm{kg\,m^{-3}}$), including solids and pore spaces
+$\rho_s$:
+Soil bulk density ($\mathrm{kg\,m^{-3}}$), including solids and pore spaces, currently
+constant across all grid cells and layers
 
-$c$:
+$c_s$:
 Soil specific heat capacity ($\mathrm{J\,kg^{-1}\,K^{-1}}$), the energy required to
 raise the temperature of 1 kg of soil by 1 K.
 
 #### Temperature Update Scheme
 
-Let $T_i^t$ represent the temperature (°C or K) of the $i^{\text{th}}$ soil layer at time
+Let $T_i^t$ represent the temperature (°C) of the $i^{\text{th}}$ soil layer at time
 $t$. The soil column is discretized into $n$ layers, each of thickness $\Delta z$ (m),
 and time advances in steps of $\Delta t$ (s).
 
@@ -268,10 +284,10 @@ layer is as follows:
 
 ```{math}
     & \frac{dQ}{dt} \\
-    & = R_{N} - H - \lambda E (- PP)\\
-    & = R_{\text{abs}} - \epsilon_{s} \sigma T_{L}^{4} -
-    \frac{\rho_a c_p}{r_a}(T_{L} - T_{A})
-    - \lambda g_{v} \frac {e_{L} - e_{A}}{p_{A}} (- PP)\\
+    & = R_{n} - H_l - \lambda E_l (- PP)\\
+    & = R_{\text{abs}} - \epsilon_{l} \sigma T_{l}^{4} -
+    \frac{\rho_a c_p}{r_a}(T_{l} - T_{a})
+    - \lambda g_{v} \frac {e_{l} - e_{a}}{p_{a}} (- PP)\\
     & = 0
 ```
 
@@ -284,35 +300,35 @@ Shortwave radiation absorbed by the canopy, equivalent to $S_0 (1-\alpha)$
 $R_{\text{em}}$:
 Emitted longwave radiation from the canopy ($\mathrm{W\,m^{-2}}$)
 
-$H$:
+$H_{l}$:
 Sensible heat flux from the canopy to the air ($\mathrm{W\,m^{-2}}$)
 
-$\lambda E$:
+$\lambda E_{l}$:
 Latent heat flux associated with transpiration from the canopy to the air
 ($\mathrm{W\,m^{-2}}$)
 
-$\epsilon_{s}$:
+$\epsilon_{l}$:
 Emissivity of the leaf (-), typically close to 1
 
 $\sigma$:
 Stefan–Boltzmann constant ($5.67 \times 10^{-8}\,\mathrm{W\,m^{-2}\,K^{-4}}$)
 
-$T_{L}$:
+$T_{l}$:
 Temperature of the leaf (°C)
 
-$T_{A}$:
+$T_{a}$:
 Temperature of the air surrounding the leaf (°C)
 
 $\lambda$:
 Latent heat of vapourisation of water ($\mathrm{kJ\,kg^{-1}}$)
 
-$e_{L}$:
+$e_{l}$:
 Effective vapour pressure of the leaf (kPa)
 
-$e_{A}$:
+$e_{a}$:
 Vapour pressure of air (kPa)
 
-$p_{A}$:
+$p_{a}$:
 Atmospheric pressure (kPa)
 
 $g_{v}$:
@@ -338,22 +354,22 @@ To iteratively solve for the leaf temperature that satisfies the energy balance
 $\frac{dQ}{dt}$ = 0, we use the Newton method:
 
 ```{math}
-T_L^{\text{new}} = T_L^{\text{old}} + W \cdot
-\frac{\frac{dQ}{dt}}{\frac{\partial \frac{dQ}{dt}}{\partial T_L^{\text{old}}}}
+T_l^{\text{new}} = T_l^{\text{old}} + W \cdot
+\frac{\frac{dQ}{dt}}{\frac{\partial \frac{dQ}{dt}}{\partial T_l^{\text{old}}}}
 ```
 
 where:
 
-$T_L^{\text{old}}$:
+$T_l^{\text{old}}$:
 Current estimate of leaf temperature (°C)
 
-$T_L^{\text{new}}$:
+$T_l^{\text{new}}$:
 Updated estimate of leaf temperature (°C)
 
 $W$:
 Step-size weighting factor (–), typically between 0.1 and 1
 
-$\frac{\partial \frac{dQ}{dt}}{\partial T_L^{\text{old}}}$:
+$\frac{\partial \frac{dQ}{dt}}{\partial T_l^{\text{old}}}$:
 The first derivative of the energy balance with respect to temperature
 
 This update adjusts the leaf temperature proportionally to the energy imbalance, scaled
@@ -367,10 +383,10 @@ The temperature derivative of the energy balance as formulated above
 is calculated analytically as:
 
 ```{math}
-\frac{\partial \frac{dQ}{dt}}{\partial T_L^{\text{old}}} =
+\frac{\partial \frac{dQ}{dt}}{\partial T_l^{\text{old}}} =
 \frac{\rho_a c_p}{r_a} +
 \frac{\rho_a \Delta_v}{r_a + r_s} \lambda +
-4 \epsilon \sigma (T_L^{\text{old}} + 273.15)^3
+4 \epsilon_l \sigma (T_l^{\text{old}} + 273.15)^3
 ```
 
 where:
@@ -393,13 +409,13 @@ Slope of the saturation vapour pressure curve ($\mathrm{kPa\, K^{-1}}$)
 $\lambda$:
 Latent heat of vapourisation of water ($\mathrm{kJ\, kg^{-1}}$)
 
-$\epsilon$:
+$\epsilon_l$:
 Leaf emissivity (-)
 
 $\sigma$:
 Stefan–Boltzmann constant ($5.67 \times 10^{-8}\,\mathrm{W\,m^{-2}\,K^{-4}}$)
 
-$T_L^{\text{old}}$:
+$T_l^{\text{old}}$:
 Previous estimate of leaf temperature (°C, converted to K in the radiation term)
 
 This derivative represents the rate at which each energy loss term changes with leaf
@@ -412,15 +428,15 @@ After updating the canopy temperature, we update the air temperature in the
 adjacent canopy layer to reflect its coupling with the leaf temperature following
 {cite:t}`bonan_climate_2019`:
 
-$$H = \frac{\rho_a c_p}{r_a}(T_{L} - T_{A})$$
+$$H = \frac{\rho_a c_p}{r_a}(T_{l} - T_{a})$$
 
 and
 
-$$T_{A}^{\text{new}} = T_{A}^{\text{old}} + \frac{H \Delta t}{\rho_a c_p z}$$
+$$T_{a}^{\text{new}} = T_{a}^{\text{old}} + \frac{H \Delta t}{\rho_a c_p z}$$
 
 where:
 
-$T_A$:
+$T_a$:
 Air temperature, (°C)
 
 $z$:
@@ -429,41 +445,31 @@ Thickness of the air layer we are updating, (m)
 Finally, we consider vertical mixing between layers and heat is transferred to the
 air above the canopy.
 
-#### Update of atmospheric moisture
+```{note}
+Advection of heat above the canopy is currently not implemented as everything is
+removed with time interval >= 1h and horizontal transfer is not considered.
+```
 
-To account for moisture added to the atmosphere from canopy transpiration and soil
-evaporation, the model updates key atmospheric humidity variables in each vertical
-layer. This ensures consistency in the representation of atmospheric water content
-across the grid.
+#### Update of atmospheric moisture
 
 Evapotranspiration and soil evaporation are initially provided in millimetres of water
-depth. These values are converted to a mass of water per unit volume of air (kg m⁻³)
-using the grid cell area. The evaporated water is then added to the relevant atmospheric
+depth. These values are converted to a mass of water per unit volume of air
+($\mathrm{kg\, m^{-3}}$) then added to the relevant atmospheric
 layers: canopy evapotranspiration is distributed across the layers surrounding the
 vegetation, while soil evaporation is added to the lowest layer near the surface.
 
 Using the updated water mass, specific humidity is recalculated for each layer by
 dividing the total water mass by the volume of air in that layer. Then, the new specific
 humidity is vertically mixed between layers and ventilated at the top of the canopy to
-make sure that water does not accumulate unrealistcaly in the canopy but stays connected
-to the atmosphere above. The resulting change in
-specific humidity is then used to compute the new vapour pressure, taking into account
-the atmospheric pressure and the molecular weight difference between water vapour and
-dry air. To maintain physical realism, the vapour pressure is capped at the saturated
-vapour pressure, avoiding supersaturation.
-
-Finally, the model derives relative humidity as the ratio of vapour pressure to
-saturated vapour pressure, expressed as a percentage. The vapour pressure deficit (VPD)
-is then calculated as the difference between saturated and actual vapour pressure,
-indicating the remaining atmospheric demand for water.
-
-This update step ensures that changes in canopy and soil water fluxes are accurately
-reflected in the atmospheric humidity profile, which in turn affects subsequent energy
-and water balance calculations.
+make sure that water does not accumulate unrealistcally in the canopy but stays connected
+to the atmosphere above. To maintain physical realism, additional redistribution steps
+are taken where necessary until all layers in the canopy are within realistic bounds.
+The resulting change in specific humidity is then used to compute the new vapour pressure
+, relative humidity, and vapour pressure deficit.
 
 ```{note}
-At the moment we get 100% relative humidity in the canopy,and 0% relative humidity above
-the canopy. This will be addressed during the model calibration.
+Advection of water above the canopy is currently not implemented as everything is
+removed with time interval >= 1h and horizontal transfer is not considered.
 ```
 
 ### Wind
@@ -506,12 +512,13 @@ with
 $$R = \sqrt{C_s + \frac{C_r LAI}{2}}$$
 
 where $C_{s}$ is the substrate surface drag coefficient, $C_{r}$ is the roughness
-element (vegetation) drag coefficient, $C_d$ is the roughness sublayer depth parameter,
+element (vegetation) drag coefficient, $C_{d}$ is the roughness sublayer depth parameter,
 $\kappa$ is the von Karman constant, and $LAI$ is the leaf area index
 ($\mathrm{m\,m^{-1}}$).
 
 The **wind speed** ($\mathrm{m\,s^{-1}}$) at any height $z$ (m) is computed using the
-logarithmic wind profile under neutral conditions:
+logarithmic wind profile under neutral conditions (based on
+{cite:t}`holmes_wind_2019`):
 
 ```{math}
 u(z) = u_{\text{ref}} \cdot \frac{\ln\left( \frac{z - d}{z_0} \right)}
@@ -526,14 +533,16 @@ Minimum wind speed is enforced below the canopy to avoid unrealistically low tur
 transport.
 
 **Friction velocity** $u_{*}$ ($\mathrm{m\,s^{-1}}$) quantifies the shear stress
-imposed by wind near the surface and is calculated from the wind speed profile:
+imposed by wind near the surface and is calculated from the wind speed profile
+(based on {cite:t}`holmes_wind_2019`):
 
 $$u_* = \frac{\kappa \cdot u(z)}{\ln\left( \frac{z - d}{z_0} \right)}$$
 
 Friction velocity is used to estimate turbulence strength and mixing coefficients.
 
 The **aerodynamic resistance** $r_a$ ($\mathrm{s\,m^{-1}}$) quantifies the resistance to
-vertical transfer of scalars (heat, water vapour) between surface and air:
+vertical transfer of scalars (heat, water vapour) between surface and air
+(based on {cite:t}`jansson_coupled_2004`):
 
 ```{math}
 r_a = \frac{1}{g_a} = \frac{\left[ \ln\left( \frac{z - d}{z_0} \right) \right]^2}
@@ -563,7 +572,7 @@ This particular form goes to zero at both z=0 and z=h and peaks somewhere within
 canopy.
 
 The **ventilation rate** $v$ represents the rate of air exchange above the
-canopy and is defined as:
+canopy and is defined as (after {cite:t}`wolfe_forest_2011`):
 
 $$v = \frac{1}{r_a \cdot h}$$
 
diff --git a/virtual_ecosystem/models/abiotic/abiotic_model.py b/virtual_ecosystem/models/abiotic/abiotic_model.py
index ed916c150..bb09b5b6c 100644
--- a/virtual_ecosystem/models/abiotic/abiotic_model.py
+++ b/virtual_ecosystem/models/abiotic/abiotic_model.py
@@ -220,10 +220,6 @@ def spinup(self) -> None:
     def _update(self, time_index: int, **kwargs: Any) -> None:
         """Function to update the abiotic model.
 
-        TODO the units of fluxes are in W m-2 and we need to make sure that the input
-        energy over a time interval is coherent with the calculations of fluxes in that
-        time interval.
-
         Args:
             time_index: The index of the current time step in the data object.
             **kwargs: Further arguments to the update method.
diff --git a/virtual_ecosystem/models/abiotic/abiotic_tools.py b/virtual_ecosystem/models/abiotic/abiotic_tools.py
index 598b59c9d..cf02773d6 100644
--- a/virtual_ecosystem/models/abiotic/abiotic_tools.py
+++ b/virtual_ecosystem/models/abiotic/abiotic_tools.py
@@ -59,6 +59,7 @@ def calculate_air_density(
         specific_gas_constant_dry_air: Specific gas constant for dry air, [J kg-1 K-1]
         celsius_to_kelvin: Factor to convert temperature in Celsius to absolute
             temperature in Kelvin
+
     Returns:
         density of air, [kg m-3].
     """
@@ -99,8 +100,8 @@ def find_last_valid_row(array: NDArray[np.floating]) -> NDArray[np.floating]:
     """Find last valid value in array for each column.
 
     This function looks for the last valid value in each column of a 2-dimensional
-    array. If the previous value is nan, it moved up the array. If all values are nan,
-    the value is set to nan, too.
+    array. If the previous value is nan, it moved up the array. If all values are NaN,
+    the value is set to NaN, too.
 
     Args:
         array: Two-dimesional array for which last valid values should be found
@@ -198,7 +199,7 @@ def compute_layer_thickness_for_varying_canopy(
     .
 
     Args:
-        heights: 2D array (n_layers, n_columns) of layer heights, [m]
+        heights: 2D array of layer heights, [m]
 
     Returns:
         2D array of layer thickness, [m], same shape as input
@@ -246,7 +247,7 @@ def calculate_specific_humidity(
         pyrealm_const: Pyrealm constants
 
     Returns:
-        Specific humidity (kg/kg)
+        Specific humidity, [kg kg-1]
     """
     # Saturation vapor pressure
     saturation_vapour_pressure = calc_vp_sat(
diff --git a/virtual_ecosystem/models/abiotic/constants.py b/virtual_ecosystem/models/abiotic/constants.py
index 104de5688..9ca54307a 100644
--- a/virtual_ecosystem/models/abiotic/constants.py
+++ b/virtual_ecosystem/models/abiotic/constants.py
@@ -114,12 +114,6 @@ class AbioticConsts(ConstantsDataclass):
     roughness significantly affects the wind flow over a particular terrain or
     surface. Implementation and value from :cite:t:`maclean_microclimc_2021`."""
 
-    canopy_temperature_ini_factor: float = 0.01
-    """Factor used to initialise canopy temperature, dimensionless.
-    
-    This factor is used to initialise canopy temperature in the model setup as a
-    function of air temperature and absorbed shortwave radiation."""
-
     light_extinction_coefficient: float = 0.01
     """Light extinction coefficient for canopy, unitless.
     
diff --git a/virtual_ecosystem/models/abiotic/energy_balance.py b/virtual_ecosystem/models/abiotic/energy_balance.py
index 7cf1bed31..47c3d537a 100644
--- a/virtual_ecosystem/models/abiotic/energy_balance.py
+++ b/virtual_ecosystem/models/abiotic/energy_balance.py
@@ -7,33 +7,37 @@
 and the exchange of heat between different layers of the canopy, must be considered
 explicitly, see :cite:t:`maclean_microclimc_2021`. This is currently not implemented.)
 
-Under steady-state, the balance equation for the leaves in each canopy layer is as
+Under steady-state, the balance equation :math:`\frac{dQ}{dt}` for the leaves in each
+canopy layer is as
 follows (after :cite:t:`maclean_microclimc_2021`):
 
 .. math::
-    R_{abs} - R_{em} - H - \lambda E - PP
-    = R_{abs} - \epsilon_{s} \sigma T_{L}^{4} - \frac{\rho_a c_p}{r_a}(T_{L} - T_{A})
-    - \lambda g_{v} \frac {e_{L} - e_{A}}{p_{a}} - PP = 0
+    \frac{dQ}{dt}
+    = R_{abs} - R_{em} - H - \lambda E - PP
+    = R_{abs} - \epsilon_{l} \sigma T_{l}^{4} - \frac{\rho_{a} c_p}{r_a}(T_{l} - T_{a})
+    - \lambda g_{v} \frac {e_{l} - e_{a}}{p_{a}} - PP = 0
 
 where :math:`R_{abs}` is absorbed radiation, :math:`R_{em}` emitted radiation, :math:`H`
-the sensible heat flux, :math:`\lambda E` the latent heat flux, :math:`\epsilon_{s}` the
-emissivity of the leaf, :math:`\sigma` the Stefan-Boltzmann constant, :math:`T_{L}` the
-absolute temperature of the leaf, :math:`T_{A}` the absolute temperature of the air
+the sensible heat flux, :math:`\lambda E` the latent heat flux, :math:`\epsilon_{l}` the
+emissivity of the leaf, :math:`\sigma` the Stefan-Boltzmann constant, :math:`T_{l}` the
+absolute temperature of the leaf, :math:`T_{a}` the absolute temperature of the air
 surrounding the leaf, :math:`\lambda` the latent heat of vapourisation of water,
-:math:`e_{L}` the effective vapour pressure of the leaf, :math:`e_{A}` the vapour
+:math:`e_{l}` the effective vapour pressure of the leaf, :math:`e_{a}` the vapour
 pressure of air and :math:`p_{a}` atmospheric pressure. :math:`\rho_a` is the density of
 air, :math:`c_{p}` is the specific heat capacity of air
-at constant pressure, :math:`r_a` is the aerodynamic resistance of the surface (leaf or
-soil), :math:`g_{v}` represents the conductivity for vapour loss from the leaves as a
+at constant pressure, :math:`r_{a}` is the aerodynamic resistance of the surface (leaf
+or soil), :math:`g_{v}` represents the conductivity for vapour loss from the leaves as a
 function of the stomatal conductivity, :math:`PP` stands for primary productivity.
 
 A challenge in solving this equation is the dependency of latent heat and emitted
 radiation on leaf temperature. We use a Newton approximation to update
 leaf temperature and air temperature iteratively.
 
-After updating each layer, temperature and vapor are mixed vertically between layers.
-Ventilation and advection are considered at the top of the canopy to remove some of the
-water and heat from the system.
+After updating each layer, temperature and vapor are mixed vertically between
+atmospheric layers.
+Advection at the top of the canopy is currently not considered as we don't have
+have horizontal exchange between grid cells and air above canopy values would be
+unrealistic.
 
 TODO plants use a fraction of the absorbed radiation of photosynthesis, this needs to be
 subtracted from the energy balance
@@ -64,12 +68,12 @@ def initialise_canopy_and_soil_fluxes(
     Args:
         air_temperature: Air temperature, [C]
         layer_structure: Instance of LayerStructure
-        light_extinction_coefficient: Light extinction coefficient for canopy
+        light_extinction_coefficient: Light extinction coefficient for canopy, unitless
         initial_flux_value: Initial non-zero flux, [W m-2]
 
     Returns:
         Dictionary with canopy temperature, [C], sensible and latent heat flux (canopy
-        and soil), [W m-2], and ground heat flux [W m-2].
+        and soil), [W m-2], and ground heat flux, [W m-2].
     """
 
     output = {}
@@ -103,7 +107,7 @@ def calculate_longwave_emission(
     emissivity: float | NDArray[np.floating],
     stefan_boltzmann: float,
 ) -> NDArray[np.floating]:
-    """Calculate longwave emission using the Stefan Boltzmann law, [W m-2].
+    """Calculate longwave emission using the Stefan Boltzmann law.
 
     According to the Stefan Boltzmann law, the amount of radiation emitted per unit time
     from the area of a black body at absolute temperature is directly proportional to
@@ -128,16 +132,16 @@ def calculate_sensible_heat_flux(
     surface_temperature: NDArray[np.floating],
     aerodynamic_resistance: float | NDArray[np.floating],
 ) -> NDArray[np.floating]:
-    r"""Calculate sensible heat flux, [W m-2].
+    r"""Calculate sensible heat flux.
 
     The sensible heat flux :math:`H` is calculated using the following equation:
 
     .. math::
-        H = \frac{\rho_{a} c_{p}}{r_{a}} (T_{S} - T_{A})
+        H = \frac{\rho_{a} c_{p}}{r_{a}} (T_{s} - T_{a})
 
     where :math:`\rho_{a}` is the density of air, :math:`c_{p}` is the specific heat
     capacity of air at constant pressure, :math:`r_{a}` is the aerodynamic resistance of
-    the surface, :math:`T_{S}` is the surface temperature, and :math:`T_{A}` is the air
+    the surface, :math:`T_{s}` is the surface temperature, and :math:`T_{a}` is the air
     temperature.
 
     Args:
@@ -174,10 +178,10 @@ def update_soil_temperature(
     Soil thermal diffusivity:
 
     .. math::
-        \alpha = \frac{\lambda}{\rho c}
+        \alpha = \frac{\lambda}{\rho_s c_s}
 
     where :math:`\lambda` is the soil thermal conductivity [W m-1 K-1],
-    :math:`\rho` is the soil bulk density [kg m-3], :math:`c` is the specific heat
+    :math:`\rho_s` is the soil bulk density [kg m-3], :math:`c_s` is the specific heat
     capacity of soil [J kg-1 K-1].
 
     Internal layer update:
@@ -189,7 +193,7 @@ def update_soil_temperature(
     Top layer update with ground heat flux:
 
     .. math::
-        T_0^{t+\Delta t} = T_0^t + (\Delta t / (\rho c \Delta z)) * G
+        T_0^{t+\Delta t} = T_0^t + (\Delta t / (\rho_s c_s \Delta z)) * G
 
     No-heat-flux bottom boundary condition:
 
@@ -200,9 +204,9 @@ def update_soil_temperature(
     Args:
         ground_heat_flux: Ground heat flux at top soil, [W m-2]
         soil_temperature: Soil temperature for each soil layer, [C]
-        soil_thermal_conductivity: Thermal conductivity of soil [W m-2 K-1]
+        soil_thermal_conductivity: Thermal conductivity of soil, [W m-2 K-1]
         soil_bulk_density: Soil bulk density, [kg m-3]
-        specific_heat_capacity_soil: Specific heat capacity of soil [J kg-1 K-1]
+        specific_heat_capacity_soil: Specific heat capacity of soil, [J kg-1 K-1]
         soil_layer_thickness: Thickness of each soil layer, [m]
         time_interval: Time interval, [s]
 
@@ -266,11 +270,11 @@ def calculate_energy_balance_residual(
     The energy balance residual (:math:`\frac{dQ}{dt}`) for the canopy is given by:
 
     .. math::
-        \frac{dQ}{dt} = R_{abs} - \epsilon \sigma T_{L}^{4} - H - \lambda E - PP
+        \frac{dQ}{dt} = R_{abs} - \epsilon_{l} \sigma T_{l}^{4} - H - \lambda E - PP
 
     Where :math:`R_abs` is the absorbed shortwave radiation by the canopy,
-    :math:`\epsilon` is the leaf emissivity, :math:`\sigma` is the Stefan-Boltzmann
-    constant, :math:`T_{L}` is the leaf temperature, :math:`H` is the sensible heat
+    :math:`\epsilon_{l}` is the leaf emissivity, :math:`\sigma` is the Stefan-Boltzmann
+    constant, :math:`T_{l}` is the leaf temperature, :math:`H` is the sensible heat
     flux from the canopy, :math:`\lambda E` is the latent heat flux from the canopy,
     :math:`PP` is a fraction of the absorbed light is used in photosynthesis (PAR).
 
@@ -371,12 +375,12 @@ def solve_canopy_temperature(
     The energy balance for the canopy is given by:
 
     .. math::
-        R_{abs} - \epsilon \sigma T_{L}^{4} - H - Q_{LE} - PP = 0
+        R_{abs} - \epsilon_{l} \sigma T_{l}^{4} - H - \lambda E - PP = 0
 
     Where :math:`R_abs` is the absorbed shortwave radiation by the canopy,
-    :math:`\epsilon` is the leaf emissivity, :math:`\sigma` is the Stefan-Boltzmann
-    constant, :math:`T_{L}` is the leaf temperature, :math:`H` is the sensible heat
-    flux from the canopy, :math:`Q_{LE}` is the latent heat flux from the canopy, and
+    :math:`\epsilon_{l}` is the leaf emissivity, :math:`\sigma` is the Stefan-Boltzmann
+    constant, :math:`T_{l}` is the leaf temperature, :math:`H` is the sensible heat
+    flux from the canopy, :math:`\lambda E` is the latent heat flux from the canopy, and
     :math:`PP` is a fraction of the absorbed light is used in photosynthesis (PAR).
 
     Note that the latent heat flux is currently a constant given by the plant model.
@@ -385,18 +389,18 @@ def solve_canopy_temperature(
     The Newton linearization for canopy temperature update is:
 
     .. math::
-        T_{L}^{new} =
-        T_{L}^{old} + W \cdot \frac{EB} {\frac{\delta EB}{\delta T_{L}^{old}}}
+        T_{l}^{new} =
+        T_{l}^{old} + W \cdot \frac{EB} {\frac{\delta EB}{\delta T_{l}^{old}}}
 
-    where :math:`\frac{\delta EB}{\delta T_{L}^{old}}` is the first derivative of the
+    where :math:`\frac{\delta EB}{\delta T_{l}^{old}}` is the first derivative of the
     energy balance closure error to temperature, and :math:`W` is a weighting for the
     step size to ensure numerical stability. The derivative is estimated analytically:
 
     .. math::
-        \frac{\delta EB}{\delta T_{L}^{old}}
+        \frac{\delta EB}{\delta T_{l}^{old}}
         = \frac{\rho_{a} c_{p}} {r_{a}}
         + \frac{\rho_{a} \Delta_{v}}{(r_{a} + r_{s})} \lambda
-        + 4 \epsilon \sigma (T_{L}^{old} + 273.15)^{3}
+        + 4 \epsilon_{l} \sigma (T_{l}^{old} + 273.15)^{3}
 
     Where :math:`c_{p}` is the specific heat capacity of air, [J kg-1 K-1],
     :math:`\rho_{a}` is the density of air, [kg m-3], :math:`\Delta_{v}` is the slope of
@@ -521,20 +525,20 @@ def update_air_temperature(
 ) -> NDArray[np.floating]:
     r"""Update air temperature in steady state.
 
-    The new air temperature :math:`T_{A}^{new}` is updated following
+    The new air temperature :math:`T_{a}^{new}` is updated following
     :cite:t:`bonan_climate_2019`:
 
     .. math ::
-        H = \frac{\rho_a c_p}{r_a}(T_{L} - T_{A})
+        H = \frac{\rho_a c_p}{r_a}(T_{l} - T_{a})
 
     and
 
     .. math::
-        T_{A}^{new} = T_{A}^{old} + \frac{H}{\rho_a c_p z}
+        T_{a}^{new} = T_{a}^{old} + \frac{H}{\rho_a c_p z}
 
     where :math:`\rho_{a}` is the density of air, :math:`c_{p}` is the specific heat
     capacity of air at constant pressure, :math:`r_{a}` is the aerodynamic resistance of
-    the surface, :math:`T_{S}` is the surface temperature, :math:`T_{A}` is the air
+    the surface, :math:`T_{s}` is the surface temperature, :math:`T_{a}` is the air
     temperature, and :math:`z` is the thickness of the air layer we are updating.
 
     Args:
@@ -547,7 +551,7 @@ def update_air_temperature(
         mixing_layer_thickness: thickness of the air layer we are updating, [m]
 
     Returns:
-        Updated air temperatures, [C]
+        updated air temperatures, [C]
     """
 
     # Update temperatures
@@ -637,6 +641,7 @@ def update_humidity_vpd(
     )
 
     # NOTE Advection not implemented as everything is removed with time interval > 1h
+    # and horizontal transfer is not implemented
     # specific_humidity_advected = wind.advect_water_from_toplayer(
     #     specific_humidity=specific_humidity_updated[0],
     #     layer_thickness=layer_thickness[0],
diff --git a/virtual_ecosystem/models/abiotic/microclimate.py b/virtual_ecosystem/models/abiotic/microclimate.py
index 323c0c6dc..b68461c99 100644
--- a/virtual_ecosystem/models/abiotic/microclimate.py
+++ b/virtual_ecosystem/models/abiotic/microclimate.py
@@ -350,6 +350,7 @@ def run_microclimate(
         )
 
         # NOTE Advection not implemented as everything is removed with time interval>1h
+        # and horizontal transfer is not implemented
         # advection_rate = (
         #   data["wind_speed_ref"].isel(time_index=time_index).to_numpy()
         #   / np.sqrt(cell_area)
diff --git a/virtual_ecosystem/models/abiotic/wind.py b/virtual_ecosystem/models/abiotic/wind.py
index aeb51498c..b64425ea3 100644
--- a/virtual_ecosystem/models/abiotic/wind.py
+++ b/virtual_ecosystem/models/abiotic/wind.py
@@ -13,7 +13,7 @@ def calculate_zero_plane_displacement(
     leaf_area_index: NDArray[np.floating],
     zero_plane_scaling_parameter: float,
 ) -> NDArray[np.floating]:
-    """Calculate zero plane displacement height, [m].
+    """Calculate zero plane displacement height.
 
     The zero plane displacement height is a concept used in micrometeorology to describe
     the flow of air near the ground or over surfaces like a forest canopy or crops. It
@@ -55,7 +55,7 @@ def calculate_roughness_length_momentum(
     min_roughness_length: float,
     von_karman_constant: float,
 ) -> NDArray[np.floating]:
-    """Calculate roughness length governing momentum transfer, [m].
+    """Calculate roughness length governing momentum transfer.
 
     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
@@ -127,9 +127,10 @@ def calculate_wind_profile(
     zero_plane_displacement: NDArray[np.floating],
     min_wind_speed: float,
 ) -> NDArray[np.floating]:
-    r"""Calculate wind speed profile, [m s-1].
+    r"""Calculate wind speed profile.
 
-    The wind speed at different heights is calculated using the following equation:
+    The wind speed at different heights is calculated using the following equation
+    (based on :cite:t:`holmes_wind_2019`):
 
     .. math::
         u(z) = u_{ref} \times \frac{ \ln \left( \frac{z - d}{z_0} \right) }
@@ -174,13 +175,14 @@ def calculate_friction_velocity(
     zero_plane_displacement: NDArray[np.floating],
     von_karman_constant: float,
 ) -> NDArray[np.floating]:
-    r"""Calculate friction velocity, [m s-1].
+    r"""Calculate friction velocity.
 
     Friction velocity is a measure of the shear stress exerted by the wind on the
     Earth's surface, representing the velocity scale that relates to turbulent energy
     transfer near the surface.
 
-    The friction velocity (:math:`u_{*}`, [m s-1]) is calculated as
+    The friction velocity (:math:`u_{*}`, [m s-1]) is calculated as (based on
+    :cite:t:`holmes_wind_2019`):
 
     :math:`u_{*} = \frac{\kappa u}{\ln{(\frac{z - d}{z_0})}}`
 
@@ -215,7 +217,7 @@ def calculate_ventilation_rate(
     """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
-    canopy and the atmosphere above.
+    canopy and the atmosphere above after :cite:t:`wolfe_forest_2011`.
 
     Args:
         aerodynamic_resistance: Aerodynamic resistance, [s m-1]
@@ -351,7 +353,7 @@ def mix_and_ventilate(
     negative concentrations.
 
     Advection is currently not implemented as everything is removed with time interval
-    > 1h.
+    > 1h and horizontal transfer is not implemented.
 
     Args:
         input_variable: Input variable for all true atmospheric layers
@@ -470,9 +472,10 @@ def calculate_aerodynamic_resistance(
     wind_speed: NDArray[np.floating],
     von_karman_constant: float,
 ) -> NDArray[np.floating]:
-    r"""Calculate aerodynamic resistance in canopy, [s m-1].
+    r"""Calculate aerodynamic resistance in canopy.
 
-    The aerodynamic resistance :math:`r_{a}` is calculated as:
+    The aerodynamic resistance :math:`r_{a}` is calculated as (based on
+    :cite:t:`jansson_coupled_2004`):
 
     .. math::
         r_{a} = \frac{ln(\frac{z-d}{z_{m}})^{2}}{\kappa ^{2} u(z)}
@@ -483,7 +486,7 @@ def calculate_aerodynamic_resistance(
     :math:`u(z)` is the wind speed at height :math:`z`.
 
     Args:
-        wind_heights: Heights where wind speed is to be calculated [m].
+        wind_heights: Heights where wind speed is to be calculated, [m].
         roughness_length: Momentum roughness length, [m]
         zero_plane_displacement: Height above the actual ground where the wind speed is
             theoretically reduced to zero due to the obstruction caused by the roughness