From 182d573022704d51599d486f3bd3faa56cba0a43 Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Tue, 28 Apr 2026 11:36:49 +0200
Subject: [PATCH 01/17] Add Thermodynamics dep

---
 Project.toml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/Project.toml b/Project.toml
index eb46b53ba..cfdb548db 100644
--- a/Project.toml
+++ b/Project.toml
@@ -18,6 +18,7 @@ KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 RingGrids = "d1845624-ad4f-453b-8ff4-a8db365bf3a7"
 SpeedyWeatherInternals = "34489162-d270-4603-8b96-37b04f830d73"
+Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [weakdeps]
@@ -44,5 +45,6 @@ KernelAbstractions = "0.9"
 Oceananigans = "0.100, 0.101, 0.102, 0.103, 0.104, 0.105, 0.106"
 RingGrids = "0.1"
 SpeedyWeatherInternals = "0.1"
+Thermodynamics = "1"
 Unitful = "1"
 julia = "1.10"

From fc95052ed79b9e2e2bc1d99c97c670d430dc0453 Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Tue, 28 Apr 2026 11:37:13 +0200
Subject: [PATCH 02/17] Add Thermodynamics dep examples

---
 examples/Project.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/examples/Project.toml b/examples/Project.toml
index ee395dac8..03b4fffd9 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -19,6 +19,7 @@ SpeedyWeatherInternals = "34489162-d270-4603-8b96-37b04f830d73"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Terrarium = "80418d68-07fa-499d-ae2b-c12e531f5cd8"
+Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 
 [sources.Terrarium]
 path = ".."

From 9a79a1a5ea3d0a56ac3e56024f940d1c614eb538 Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Tue, 28 Apr 2026 11:38:32 +0200
Subject: [PATCH 03/17] Add WIP extension Thermodynamics

Co-authored-by: Copilot <copilot@github.com>
---
 src/processes/physical_constants.jl | 63 +++++++++++++++++++++++++----
 1 file changed, 56 insertions(+), 7 deletions(-)

diff --git a/src/processes/physical_constants.jl b/src/processes/physical_constants.jl
index 630b5bf8e..5bfb85607 100644
--- a/src/processes/physical_constants.jl
+++ b/src/processes/physical_constants.jl
@@ -1,3 +1,18 @@
+import Thermodynamics.Parameters:
+    AbstractThermodynamicsParameters,
+    R_d,
+    R_v,
+    cp_d,
+    cp_i,
+    cp_l,
+    cp_v,
+    LH_v0,
+    LH_s0,
+    T_0,
+    T_freeze,
+    T_triple,
+    press_triple
+
 """
     $TYPEDEF
 
@@ -6,18 +21,27 @@ A collection of general physical constants that do not (usually) need to be vari
 Properties:
 $FIELDS
 """
-@kwdef struct PhysicalConstants{NF}
+@kwdef struct PhysicalConstants{NF} <: AbstractThermodynamicsParameters{NF}
     "Density of water in kg/m^3"
     ρw::NF = 1000.0
 
     "Density of ice in kg/m^3"
-    ρi::NF = 916.2
+    ρi::NF = 916.7
 
     "Density of air at standard pressure and 0°C in kg/m^3"
     ρₐ::NF = 1.293
 
-    "Specific heat capacity of dry air at standard pressure and 0°C in J/(m^3*K)"
-    cₐ::NF = 1005.7
+    "Isobaric specific heat capacity of dry air at standard pressure and 0°C in J/(m^3*K)"
+    cp_d::NF = 1004.5
+
+    "Isobaric specific heat capacity of ice at standard pressure and 0°C in J/(m^3*K)"
+    cp_i::NF = 2070.0
+
+    "Isobaric specific heat capacity of liquid water at standard pressure and 0°C in J/(m^3*K)"
+    cp_l::NF = 4186.0
+
+    "Isobaric specific heat capacity of water vapor at standard pressure and 0°C in J/(m^3*K)"
+    cp_v::NF = 1846.0
 
     "Sepcific latent heat of fusion of water in J/kg"
     Lsl::NF = 3.34e5
@@ -32,7 +56,16 @@ $FIELDS
     g::NF = 9.80665
 
     "Reference temperature (0°C in Kelvin)"
-    Tref::NF = 273.15
+    T_ref::NF = 273.16
+
+    "Freezing temperature of water in Kelvin"
+    T_freeze::NF = 273.16
+
+    "Triple point temperature of water in Kelvin"
+    T_triple::NF = 273.16
+
+    "Triple point pressure of water in Pa"
+    press_triple::NF = 611.657
 
     "Stefan-Boltzmann constant in J/(s*m^2*K^4)"
     σ::NF = 5.6704e-8
@@ -43,8 +76,11 @@ $FIELDS
     "Ratio of molecular weight of water vapor to dry air"
     ε::NF = 0.622
 
-    "Specific gas constant of air in J/(kg*K)"
-    Rₐ::NF = 287.058
+    "Specific gas constant of dry air in J/(kg*K)"
+    R_d::NF = 287.058
+
+    "Specific gas constant of water vapor in J/(kg*K)"
+    R_v::NF = 461.5
 
     "Atomic mass of carbon [gC/mol]"
     C_mass::NF = 12.0
@@ -52,6 +88,19 @@ end
 
 PhysicalConstants(::Type{NF}; kwargs...) where {NF} = PhysicalConstants{NF}(; kwargs...)
 
+@inline R_d(c::PhysicalConstants) = c.R_d
+@inline R_v(c::PhysicalConstants) = c.R_v
+@inline cp_d(c::PhysicalConstants) = c.cp_d
+@inline cp_i(c::PhysicalConstants) = c.cp_i
+@inline cp_l(c::PhysicalConstants) = c.cp_l
+@inline cp_v(c::PhysicalConstants) = c.cp_v
+@inline LH_v0(c::PhysicalConstants) = c.Llg
+@inline LH_s0(c::PhysicalConstants) = c.Lsg
+@inline T_0(c::PhysicalConstants) = c.T_ref
+@inline T_freeze(c::PhysicalConstants) = c.T_freeze
+@inline T_triple(c::PhysicalConstants) = c.T_triple
+@inline press_triple(c::PhysicalConstants) = c.press_triple
+
 """
     celsius_to_kelvin(c::PhysicalConstants, T)
 

From ab3b7e5e152cfc7d3a260946b3c8f54302bcc59f Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Tue, 28 Apr 2026 15:10:07 +0200
Subject: [PATCH 04/17] Add Rv_over_Rd method

---
 src/processes/physical_constants.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/processes/physical_constants.jl b/src/processes/physical_constants.jl
index 5bfb85607..8af6763e3 100644
--- a/src/processes/physical_constants.jl
+++ b/src/processes/physical_constants.jl
@@ -2,6 +2,7 @@ import Thermodynamics.Parameters:
     AbstractThermodynamicsParameters,
     R_d,
     R_v,
+    Rv_over_Rd,
     cp_d,
     cp_i,
     cp_l,
@@ -101,6 +102,9 @@ PhysicalConstants(::Type{NF}; kwargs...) where {NF} = PhysicalConstants{NF}(; kw
 @inline T_triple(c::PhysicalConstants) = c.T_triple
 @inline press_triple(c::PhysicalConstants) = c.press_triple
 
+# Derived parameters
+@inline Rv_over_Rd(c::PhysicalConstants) = R_v(c) / R_d(c)
+
 """
     celsius_to_kelvin(c::PhysicalConstants, T)
 

From 38a2a2e0bef41e825274d2150cba51be85b2b0a4 Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Tue, 28 Apr 2026 15:16:53 +0200
Subject: [PATCH 05/17] Remove constant air density

---
 src/processes/physical_constants.jl | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/processes/physical_constants.jl b/src/processes/physical_constants.jl
index 8af6763e3..715362f9d 100644
--- a/src/processes/physical_constants.jl
+++ b/src/processes/physical_constants.jl
@@ -13,6 +13,7 @@ import Thermodynamics.Parameters:
     T_freeze,
     T_triple,
     press_triple
+import Thermodynamics: air_density
 
 """
     $TYPEDEF
@@ -29,9 +30,6 @@ $FIELDS
     "Density of ice in kg/m^3"
     ρi::NF = 916.7
 
-    "Density of air at standard pressure and 0°C in kg/m^3"
-    ρₐ::NF = 1.293
-
     "Isobaric specific heat capacity of dry air at standard pressure and 0°C in J/(m^3*K)"
     cp_d::NF = 1004.5
 

From 9be5386ef7ef8009c7b1917d7f8f75150f26063d Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Tue, 28 Apr 2026 15:17:19 +0200
Subject: [PATCH 06/17] Add utils from Thermodynamics

Co-authored-by: Copilot <copilot@github.com>
---
 src/processes/physics_utils.jl | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/processes/physics_utils.jl b/src/processes/physics_utils.jl
index dfc70e920..5e953aebb 100644
--- a/src/processes/physics_utils.jl
+++ b/src/processes/physics_utils.jl
@@ -1,3 +1,9 @@
+import Thermodynamics:
+    partial_pressure_vapor,
+    saturation_vapor_pressure,
+    q_vap_from_RH,
+    q_vap_saturation
+
 """
 Return the number of seconds per day in the given number format.
 """
@@ -36,6 +42,7 @@ Convert the vapor pressure `e` to specific humidity at the given pressure `p` ba
 molecular weight ratio ε.
 """
 @inline vapor_pressure_to_specific_humidity(e, p, ε) = ε * e / p
+# TODO: Replace by q_vap_from_p_vap
 
 """
     specific_humidity_to_vapor_pressure(q, p, ε)
@@ -47,7 +54,7 @@ molecular weight ratio ε.
     e = q * p / (ε + (1 - ε) * q)
     return e
 end
-
+# TODO: replace by partial_pressure_vapor
 """
     relative_to_specific_humidity(r_h, pr, T, ε)
 
@@ -74,10 +81,10 @@ Coefficients of August-Roche-Magnus equation taken from [alduchovImprovedMagnusF
 # References
 * [alduchovImprovedMagnusForm1996](@cite) Alduchov and Eskridge, Journal of Applied Meteorology and Climatology (1996)
 """
-@inline function saturation_vapor_pressure(T::NF) where {NF}
+@inline function saturation_vapor_pressure(c::PhysicalConstants, T::NF) where {NF}
     return if T <= zero(T)
-        saturation_vapor_pressure(T, NF(611.0), NF(22.46), NF(272.62))
+        saturation_vapor_pressure(c, T, Ice())
     else
-        saturation_vapor_pressure(T, NF(611.0), NF(17.62), NF(243.12))
+        saturation_vapor_pressure(c, T, Liquid())
     end
 end

From a7f81d7d31892bd37b35b510718d4171cfa33b9c Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Tue, 28 Apr 2026 17:48:17 +0200
Subject: [PATCH 07/17] Add WIP turbulent fluxes refactor

Co-authored-by: Copilot <copilot@github.com>
---
 .../atmosphere/prescribed_atmosphere.jl       | 22 +-------
 src/processes/physical_constants.jl           | 14 +++--
 src/processes/physics_utils.jl                | 51 ++++++-------------
 .../surface_energy/turbulent_fluxes.jl        | 36 +++++++++----
 4 files changed, 48 insertions(+), 75 deletions(-)

diff --git a/src/processes/atmosphere/prescribed_atmosphere.jl b/src/processes/atmosphere/prescribed_atmosphere.jl
index 42ac0cef6..787ac173c 100644
--- a/src/processes/atmosphere/prescribed_atmosphere.jl
+++ b/src/processes/atmosphere/prescribed_atmosphere.jl
@@ -102,26 +102,6 @@ variables(atmos::PrescribedAtmosphere{NF}) where {NF} = (
 
 @inline compute_tendencies!(state, grid, atmos::PrescribedAtmosphere) = nothing
 
-"""
-    $SIGNATURES
-
-Computes the vapor pressure deficit for an air parcel at temperature `T` with surface pressure `pres`
-and specific humidity of air `q_air`. Assumes that air parcel is over water when `T > 0°C` and over 
-ice when `T < 0°C`.
-"""
-@inline function vapor_pressure_deficit(c::PhysicalConstants{NF}, pres, q_air, T) where {NF}
-    # Compute saturation vapor pressure of air parcel at temperature T
-    e_sat = saturation_vapor_pressure(T)
-
-    # Convert air specific humidity to vapor pressure [Pa]
-    e_air = specific_humidity_to_vapor_pressure(q_air, pres, c.ε)
-
-    # Compute vapor pressure deficit [Pa]
-    vpd = max(e_sat - e_air, NF(0.1))
-
-    return vpd
-end
-
 """
     aerodynamic_resistance(i, j, grid, fields, atmos::PrescribedAtmosphere)
 
@@ -184,7 +164,7 @@ Computes the vapor pressure deficit (VPD) at atmospheric reference level given t
     T_air = air_temperature(i, j, grid, fields, atmos)
     q_air = specific_humidity(i, j, grid, fields, atmos)
     p = air_pressure(i, j, grid, fields, atmos)
-    vpd = vapor_pressure_deficit(c, p, q_air, T_air)
+    vpd = vapor_pressure_deficit(c, celsius_to_kelvin(c, T_air), p, q_air)
     return vpd
 end
 
diff --git a/src/processes/physical_constants.jl b/src/processes/physical_constants.jl
index 715362f9d..5ecfafd91 100644
--- a/src/processes/physical_constants.jl
+++ b/src/processes/physical_constants.jl
@@ -42,13 +42,13 @@ $FIELDS
     "Isobaric specific heat capacity of water vapor at standard pressure and 0°C in J/(m^3*K)"
     cp_v::NF = 1846.0
 
-    "Sepcific latent heat of fusion of water in J/kg"
+    "Specific latent heat of fusion of water in J/kg at 0°C"
     Lsl::NF = 3.34e5
 
-    "Specific latent heat of vaporization of water in J/kg"
+    "Specific latent heat of vaporization of water in J/kg at 0°C"
     Llg::NF = 2.257e6
 
-    "Specific latent heat of sublimation of water in J/kg"
+    "Specific latent heat of sublimation of water in J/kg at 0°C"
     Lsg::NF = 2.834e6
 
     "Gravitational constant in m/s^2"
@@ -72,9 +72,6 @@ $FIELDS
     "von Kármán constant"
     κ::NF = 0.4
 
-    "Ratio of molecular weight of water vapor to dry air"
-    ε::NF = 0.622
-
     "Specific gas constant of dry air in J/(kg*K)"
     R_d::NF = 287.058
 
@@ -102,13 +99,14 @@ PhysicalConstants(::Type{NF}; kwargs...) where {NF} = PhysicalConstants{NF}(; kw
 
 # Derived parameters
 @inline Rv_over_Rd(c::PhysicalConstants) = R_v(c) / R_d(c)
+@inline ε(c::PhysicalConstants) = R_d(c) / R_v(c)
 
 """
     celsius_to_kelvin(c::PhysicalConstants, T)
 
 Convert the given temperature in °C to Kelvin based on the constant `Tref`.
 """
-@inline celsius_to_kelvin(c::PhysicalConstants, T) = T + c.Tref
+@inline celsius_to_kelvin(c::PhysicalConstants, T) = T + c.T_ref
 
 """
     stefan_boltzmann(c::PhysicalConstants, T, ϵ)
@@ -123,4 +121,4 @@ and ϵ is the emissivity.
 
 Calcualte the psychrometric constant at the given atmospheric pressure `p`.
 """
-@inline psychrometric_constant(c::PhysicalConstants, p) = c.cₐ * p / (c.Llg * c.ε)
+@inline psychrometric_constant(c::PhysicalConstants, p) = cp_d(c) * p / (LH_v0(c) * ε(c))
diff --git a/src/processes/physics_utils.jl b/src/processes/physics_utils.jl
index 5e953aebb..9fc14b7a2 100644
--- a/src/processes/physics_utils.jl
+++ b/src/processes/physics_utils.jl
@@ -2,7 +2,10 @@ import Thermodynamics:
     partial_pressure_vapor,
     saturation_vapor_pressure,
     q_vap_from_RH,
-    q_vap_saturation
+    q_vap_saturation,
+    vapor_pressure_deficit,
+    Ice,
+    Liquid
 
 """
 Return the number of seconds per day in the given number format.
@@ -36,45 +39,22 @@ Compute partial pressure of CO2 from surface pressure and CO2 concentration in P
 end
 
 """
-    vapor_pressure_to_specific_humidity(e, p, ε)
+    relative_to_specific_humidity(r_h, pr, T, c::PhysicalConstants)
 
-Convert the vapor pressure `e` to specific humidity at the given pressure `p` based on the
-molecular weight ratio ε.
+Derives specific humidity from measured relative humidity `r_h` [%], air pressure `pr` [Pa],
+air temperature `T` [°C], and physical constants `c`. Assumes saturation over ice for
+`T <= 0°C` and over liquid water otherwise.
 """
-@inline vapor_pressure_to_specific_humidity(e, p, ε) = ε * e / p
-# TODO: Replace by q_vap_from_p_vap
-
-"""
-    specific_humidity_to_vapor_pressure(q, p, ε)
-
-Convert the specific humidity `q` to vapor pressure at the given pressure `p` based on the
-molecular weight ratio ε.
-"""
-@inline function specific_humidity_to_vapor_pressure(q, p, ε)
-    e = q * p / (ε + (1 - ε) * q)
-    return e
+@inline function relative_to_specific_humidity(r_h, pr, T, c::PhysicalConstants)
+    T_K = celsius_to_kelvin(c, T)
+    phase = T <= zero(T) ? Ice() : Liquid()
+    return q_vap_from_RH(c, pr, T_K, r_h / 100, phase)
 end
-# TODO: replace by partial_pressure_vapor
-"""
-    relative_to_specific_humidity(r_h, pr, T, ε)
-
-Derives specific humidity from measured relative humidity, air pressure, air temperature, and molecular weight ratio.
-"""
-@inline relative_to_specific_humidity(r_h, pr, T, ε) = ε * (r_h / 100) * saturation_vapor_pressure(T) / pr
-
-# saturation vapor pressure
-"""
-    saturation_vapor_pressure(T, a₁, a₂, a₃)
-
-August-Roche-Magnus equation for saturation vapor pressure at temperature `T` with empirical
-coefficients a₁, a₂, and a₃.
-"""
-@inline saturation_vapor_pressure(T, a₁, a₂, a₃) = a₁ * exp(a₂ * T / (T + a₃))
 
 """
     saturation_vapor_pressure(T)
 
-Saturation vapor pressure of an air parcel at the given temperature `T`. By default, the saturation vapor
+Saturation vapor pressure of an air parcel at the given temperature `T` in °C. By default, the saturation vapor
 pressure is computed over ice for `T <= 0°C` and over water for `T > 0°C`
 Coefficients of August-Roche-Magnus equation taken from [alduchovImprovedMagnusForm1996](@cite).
 
@@ -82,9 +62,10 @@ Coefficients of August-Roche-Magnus equation taken from [alduchovImprovedMagnusF
 * [alduchovImprovedMagnusForm1996](@cite) Alduchov and Eskridge, Journal of Applied Meteorology and Climatology (1996)
 """
 @inline function saturation_vapor_pressure(c::PhysicalConstants, T::NF) where {NF}
+    T_K = celsius_to_kelvin(c, T)
     return if T <= zero(T)
-        saturation_vapor_pressure(c, T, Ice())
+        saturation_vapor_pressure(c, T_K, Ice())
     else
-        saturation_vapor_pressure(c, T, Liquid())
+        saturation_vapor_pressure(c, T_K, Liquid())
     end
 end
diff --git a/src/processes/surface_energy/turbulent_fluxes.jl b/src/processes/surface_energy/turbulent_fluxes.jl
index 0c8f44ecb..a3dd6fd72 100644
--- a/src/processes/surface_energy/turbulent_fluxes.jl
+++ b/src/processes/surface_energy/turbulent_fluxes.jl
@@ -56,8 +56,8 @@ Computes the difference in vapor pressure between a saturated surface at tempera
 and the atmosphere, defined by its specific humidity `q_air`.
 """
 function vapor_pressure_difference(c::PhysicalConstants, p, q_air, T)
-    e_air = specific_humidity_to_vapor_pressure(q_air, p, c.ε)
-    e_sat_s = saturation_vapor_pressure(T)
+    e_air = partial_pressure_vapor(c, p, q_air)
+    e_sat_s = saturation_vapor_pressure(c, T)
     Δe = e_sat_s - e_air
     return Δe
 end
@@ -69,9 +69,11 @@ Computes the difference in specific humidity between a saturated surface at temp
 and the atmosphere, defined by its specific humidity `q_air`.
 """
 function specific_humidity_difference(c::PhysicalConstants, p, q_air, T)
-    e_sat_s = saturation_vapor_pressure(T)
-    q_sat_s = vapor_pressure_to_specific_humidity(e_sat_s, p, c.ε)
-    Δq = q_sat_s - q_air
+    T_K = celsius_to_kelvin(c, T)
+    # TODO: should use surface air density for better accuracy
+    ρₐ = air_density(c, T_K, p, q_air)
+    q_sat = q_vap_saturation(c, T_K, ρₐ)
+    Δq = q_sat - q_air
     return Δq
 end
 
@@ -139,11 +141,14 @@ Compute the sensible heat flux at `i, j` based on the current skin temperature a
         constants::PhysicalConstants,
         atmos::AbstractAtmosphere
     )
-    let ρₐ = constants.ρₐ, # density of air
-            cₐ = constants.cₐ, # specific heat capacity of air
+    let cₐ = cp_d(constants), # specific heat capacity of air
             rₐ = aerodynamic_resistance(i, j, grid, fields, atmos), # aerodynamic resistance
             Tₛ = skin_temperature(i, j, grid, fields, skinT), # skin temperature
             Tₐ = air_temperature(i, j, grid, fields, atmos), # air temperature
+            pres = air_pressure(i, j, grid, fields, atmos),
+            q_air = specific_humidity(i, j, grid, fields, atmos),
+            # TODO: should be evaluated at surface temperature for better accuracy
+            ρₐ = air_density(constants, celsius_to_kelvin(constants, Tₐ), pres, q_air),
             Q_T = (Tₛ - Tₐ) / rₐ  # bulk aerodynamic temperature-gradient
         # Calculate sensible heat flux (positive upwards)
         Hₛ = compute_sensible_heat_flux(tur, Q_T, ρₐ, cₐ)
@@ -166,8 +171,12 @@ to the latent heat flux.
         atmos::AbstractAtmosphere
     )
     let L = constants.Llg, # specific latent heat of vaporization of water
-            ρₐ = constants.ρₐ, # density of air
             Tₛ = skin_temperature(i, j, grid, fields, skinT),
+            Tₐ = air_temperature(i, j, grid, fields, atmos), # air temperature
+            pres = air_pressure(i, j, grid, fields, atmos),
+            q_air = specific_humidity(i, j, grid, fields, atmos),
+            # TODO: should be evaluated at surface temperature for better accuracy
+            ρₐ = air_density(constants, celsius_to_kelvin(constants, Tₐ), pres, q_air),
             rₐ = aerodynamic_resistance(i, j, grid, fields, atmos), # aerodynamic resistance
             Δq = compute_specific_humidity_difference(i, j, grid, fields, atmos, constants, Tₛ),
             Q_h = Δq / rₐ  # humidity flux
@@ -188,10 +197,15 @@ defined by `evtr` which is assumed to be already computed.
         i, j, grid, fields,
         tur::DiagnosedTurbulentFluxes,
         evtr::AbstractEvapotranspiration,
-        constants::PhysicalConstants
+        constants::PhysicalConstants,
+        atmos::AbstractAtmosphere
     )
     let L = constants.Llg, # specific latent heat of vaporization of water
-            ρₐ = constants.ρₐ, # density of air
+            Tₐ = air_temperature(i, j, grid, fields, atmos), # air temperature
+            pres = air_pressure(i, j, grid, fields, atmos),
+            q_air = specific_humidity(i, j, grid, fields, atmos),
+            # TODO: should be evaluated at surface temperature for better accuracy
+            ρₐ = air_density(constants, celsius_to_kelvin(constants, Tₐ), pres, q_air),
             Q_h = surface_humidity_flux(i, j, grid, fields, evtr)   # humidity flux
         # Calculate latent heat flux (positive upwards)
         Hₗ = compute_latent_heat_flux(tur, Q_h, ρₐ, L)
@@ -228,5 +242,5 @@ end
     # compute sensible heat flux
     out.sensible_heat_flux[i, j, 1] = compute_sensible_heat_flux(i, j, grid, fields, tur, skinT, constants, atmos)
     # compute latent heat flux
-    out.latent_heat_flux[i, j, 1] = compute_latent_heat_flux(i, j, grid, fields, tur, evtr, constants)
+    out.latent_heat_flux[i, j, 1] = compute_latent_heat_flux(i, j, grid, fields, tur, evtr, constants, atmos)
 end

From 89cbc7c993e654826cfce9a5b6c462d17cc72018 Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Tue, 28 Apr 2026 18:41:19 +0200
Subject: [PATCH 08/17] Fix tests

Co-authored-by: Copilot <copilot@github.com>
---
 src/processes/physics_utils.jl                    | 15 +++++++++++++++
 .../surface_energy/surface_energy_balance.jl      |  2 +-
 src/processes/surface_energy/turbulent_fluxes.jl  |  2 +-
 test/surface_energy/radiative_fluxes.jl           |  3 ++-
 4 files changed, 19 insertions(+), 3 deletions(-)

diff --git a/src/processes/physics_utils.jl b/src/processes/physics_utils.jl
index 9fc14b7a2..890abcf1e 100644
--- a/src/processes/physics_utils.jl
+++ b/src/processes/physics_utils.jl
@@ -69,3 +69,18 @@ Coefficients of August-Roche-Magnus equation taken from [alduchovImprovedMagnusF
         saturation_vapor_pressure(c, T_K, Liquid())
     end
 end
+
+"""
+    q_vap_saturation(c::PhysicalConstants, T, ρ)
+
+Saturation specific humidity at temperature `T` [°C] and density `ρ` [kg/m³]. Dispatches
+over ice for `T <= 0°C` and over liquid water otherwise.
+"""
+@inline function q_vap_saturation(c::PhysicalConstants, T, ρ)
+    T_K = celsius_to_kelvin(c, T)
+    return if T <= zero(T)
+        q_vap_saturation(c, T_K, ρ, Ice())
+    else
+        q_vap_saturation(c, T_K, ρ, Liquid())
+    end
+end
diff --git a/src/processes/surface_energy/surface_energy_balance.jl b/src/processes/surface_energy/surface_energy_balance.jl
index b409ebdee..a1db2b20e 100644
--- a/src/processes/surface_energy/surface_energy_balance.jl
+++ b/src/processes/surface_energy/surface_energy_balance.jl
@@ -136,7 +136,7 @@ Fused kernel function that computes the radiative and turbulent fluxes, as well
         out.latent_heat_flux[i, j, 1] = compute_latent_heat_flux(i, j, grid, fields, seb.turbulent_fluxes, seb.skin_temperature, constants, atmos)
     else
         # Coupling with surface hydrology ET scheme
-        out.latent_heat_flux[i, j, 1] = compute_latent_heat_flux(i, j, grid, fields, seb.turbulent_fluxes, evtr, constants)
+        out.latent_heat_flux[i, j, 1] = compute_latent_heat_flux(i, j, grid, fields, seb.turbulent_fluxes, evtr, constants, atmos)
     end
     # Compute ground heat flux
     out.ground_heat_flux[i, j, 1] = compute_ground_heat_flux(i, j, grid, fields, seb.skin_temperature, seb.radiative_fluxes, seb.turbulent_fluxes)
diff --git a/src/processes/surface_energy/turbulent_fluxes.jl b/src/processes/surface_energy/turbulent_fluxes.jl
index a3dd6fd72..149bb9061 100644
--- a/src/processes/surface_energy/turbulent_fluxes.jl
+++ b/src/processes/surface_energy/turbulent_fluxes.jl
@@ -72,7 +72,7 @@ function specific_humidity_difference(c::PhysicalConstants, p, q_air, T)
     T_K = celsius_to_kelvin(c, T)
     # TODO: should use surface air density for better accuracy
     ρₐ = air_density(c, T_K, p, q_air)
-    q_sat = q_vap_saturation(c, T_K, ρₐ)
+    q_sat = q_vap_saturation(c, T, ρₐ)
     Δq = q_sat - q_air
     return Δq
 end
diff --git a/test/surface_energy/radiative_fluxes.jl b/test/surface_energy/radiative_fluxes.jl
index ec9e0ada4..f6201b035 100644
--- a/test/surface_energy/radiative_fluxes.jl
+++ b/test/surface_energy/radiative_fluxes.jl
@@ -34,6 +34,7 @@ end
     set!(state.surface_longwave_down, surface_longwave_down)
     compute_auxiliary!(state, grid, radiative_fluxes, seb, model.constants, model.atmosphere)
     @test all(state.surface_shortwave_up .≈ 0.5 * surface_shortwave_down)
-    @test all(state.surface_longwave_up .≈ (1 - 0.9) * surface_longwave_down + Terrarium.stefan_boltzmann(model.constants, 273.15, 0.9))
+    T_skin_K = Terrarium.celsius_to_kelvin(model.constants, 0.0)
+    @test all(state.surface_longwave_up .≈ (1 - 0.9) * surface_longwave_down + Terrarium.stefan_boltzmann(model.constants, T_skin_K, 0.9))
     @test all(state.surface_net_radiation .≈ state.surface_shortwave_up - surface_shortwave_down + state.surface_longwave_up - surface_longwave_down)
 end

From fa1e0ccb7bc53324495120554c77c0f0f4029376 Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Wed, 29 Apr 2026 15:50:05 +0200
Subject: [PATCH 09/17] Add explicit vpd wrapper

Co-authored-by: Copilot <copilot@github.com>
---
 .../atmosphere/prescribed_atmosphere.jl          | 13 +++++++++++++
 src/processes/physics_utils.jl                   |  1 -
 test/Project.toml                                |  3 ++-
 test/atmosphere/prescribed_atmosphere.jl         | 16 ++++++++++++++++
 4 files changed, 31 insertions(+), 2 deletions(-)
 create mode 100644 test/atmosphere/prescribed_atmosphere.jl

diff --git a/src/processes/atmosphere/prescribed_atmosphere.jl b/src/processes/atmosphere/prescribed_atmosphere.jl
index 787ac173c..090fadbbc 100644
--- a/src/processes/atmosphere/prescribed_atmosphere.jl
+++ b/src/processes/atmosphere/prescribed_atmosphere.jl
@@ -1,3 +1,5 @@
+import Thermodynamics
+
 """
 Generic type representing the concentration of a particular tracer gas in the atmosphere.
 """
@@ -102,6 +104,17 @@ variables(atmos::PrescribedAtmosphere{NF}) where {NF} = (
 
 @inline compute_tendencies!(state, grid, atmos::PrescribedAtmosphere) = nothing
 
+"""
+    $SIGNATURES
+Computes the vapor pressure deficit for an air parcel at temperature `T` [°C] with 
+surface pressure `pres` [Pa] and specific humidity of air `q_air` [kg/kg].
+Assumes that air parcel is over water when `T > 0°C` and over ice when `T < 0°C`.
+"""
+@inline function vapor_pressure_deficit(c::PhysicalConstants, T, pres, q_air)
+    T_K = celsius_to_kelvin(c, T)
+    vpd = Thermodynamics.vapor_pressure_deficit(c, T_K, pres, q_air)
+    return vpd
+end
 """
     aerodynamic_resistance(i, j, grid, fields, atmos::PrescribedAtmosphere)
 
diff --git a/src/processes/physics_utils.jl b/src/processes/physics_utils.jl
index 890abcf1e..a5c8befaf 100644
--- a/src/processes/physics_utils.jl
+++ b/src/processes/physics_utils.jl
@@ -3,7 +3,6 @@ import Thermodynamics:
     saturation_vapor_pressure,
     q_vap_from_RH,
     q_vap_saturation,
-    vapor_pressure_deficit,
     Ice,
     Liquid
 
diff --git a/test/Project.toml b/test/Project.toml
index 321141a56..3dd8dac9f 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -3,15 +3,16 @@ Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
-NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 FreezeCurves = "71e4ad71-e4f2-45a3-aa0a-91ffaa9676be"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
+NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 RingGrids = "d1845624-ad4f-453b-8ff4-a8db365bf3a7"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
diff --git a/test/atmosphere/prescribed_atmosphere.jl b/test/atmosphere/prescribed_atmosphere.jl
new file mode 100644
index 000000000..ad61ce0e7
--- /dev/null
+++ b/test/atmosphere/prescribed_atmosphere.jl
@@ -0,0 +1,16 @@
+using Terrarium
+using Test
+import Thermodynamics as TD
+
+
+@testset "Vapor pressure deficit" begin
+    FT = Float64
+    physical_constants = PhysicalConstants(FT)
+    temperature_degC = 25.0
+    temperature_K = Terrarium.celsius_to_kelvin(physical_constants, temperature_degC)
+    pressure = 101325.0 # Pa
+    total_specific_humidity = 0.01 # kg/kg
+    vpd_TD = TD.vapor_pressure_deficit(physical_constants, temperature_K, pressure, total_specific_humidity)
+    vpd_terrarium = Terrarium.vapor_pressure_deficit(physical_constants, temperature_degC, pressure, total_specific_humidity)
+    @test vpd_TD ≈ vpd_terrarium
+end

From 19ceba357fc8c1979ad7dc71dcca376433cdfcef Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Thu, 30 Apr 2026 11:43:05 +0200
Subject: [PATCH 10/17] Fix units temperature in VPD kernel function

---
 src/processes/atmosphere/prescribed_atmosphere.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/processes/atmosphere/prescribed_atmosphere.jl b/src/processes/atmosphere/prescribed_atmosphere.jl
index 090fadbbc..928806f57 100644
--- a/src/processes/atmosphere/prescribed_atmosphere.jl
+++ b/src/processes/atmosphere/prescribed_atmosphere.jl
@@ -177,7 +177,7 @@ Computes the vapor pressure deficit (VPD) at atmospheric reference level given t
     T_air = air_temperature(i, j, grid, fields, atmos)
     q_air = specific_humidity(i, j, grid, fields, atmos)
     p = air_pressure(i, j, grid, fields, atmos)
-    vpd = vapor_pressure_deficit(c, celsius_to_kelvin(c, T_air), p, q_air)
+    vpd = vapor_pressure_deficit(c, T_air, p, q_air)
     return vpd
 end
 

From 4f460b9b7a4cd75a33b5b2d789758bef3c455de7 Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Thu, 30 Apr 2026 11:58:33 +0200
Subject: [PATCH 11/17] clarify to do

Co-authored-by: Copilot <copilot@github.com>
---
 src/processes/surface_energy/turbulent_fluxes.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/processes/surface_energy/turbulent_fluxes.jl b/src/processes/surface_energy/turbulent_fluxes.jl
index 149bb9061..2465c9f7a 100644
--- a/src/processes/surface_energy/turbulent_fluxes.jl
+++ b/src/processes/surface_energy/turbulent_fluxes.jl
@@ -147,7 +147,7 @@ Compute the sensible heat flux at `i, j` based on the current skin temperature a
             Tₐ = air_temperature(i, j, grid, fields, atmos), # air temperature
             pres = air_pressure(i, j, grid, fields, atmos),
             q_air = specific_humidity(i, j, grid, fields, atmos),
-            # TODO: should be evaluated at surface temperature for better accuracy
+            # TODO: density should be evaluated at surface temperature for better accuracy
             ρₐ = air_density(constants, celsius_to_kelvin(constants, Tₐ), pres, q_air),
             Q_T = (Tₛ - Tₐ) / rₐ  # bulk aerodynamic temperature-gradient
         # Calculate sensible heat flux (positive upwards)
@@ -175,7 +175,7 @@ to the latent heat flux.
             Tₐ = air_temperature(i, j, grid, fields, atmos), # air temperature
             pres = air_pressure(i, j, grid, fields, atmos),
             q_air = specific_humidity(i, j, grid, fields, atmos),
-            # TODO: should be evaluated at surface temperature for better accuracy
+            # TODO: density should be evaluated at surface temperature for better accuracy
             ρₐ = air_density(constants, celsius_to_kelvin(constants, Tₐ), pres, q_air),
             rₐ = aerodynamic_resistance(i, j, grid, fields, atmos), # aerodynamic resistance
             Δq = compute_specific_humidity_difference(i, j, grid, fields, atmos, constants, Tₛ),
@@ -204,7 +204,7 @@ defined by `evtr` which is assumed to be already computed.
             Tₐ = air_temperature(i, j, grid, fields, atmos), # air temperature
             pres = air_pressure(i, j, grid, fields, atmos),
             q_air = specific_humidity(i, j, grid, fields, atmos),
-            # TODO: should be evaluated at surface temperature for better accuracy
+            # TODO: density should be evaluated at surface temperature for better accuracy
             ρₐ = air_density(constants, celsius_to_kelvin(constants, Tₐ), pres, q_air),
             Q_h = surface_humidity_flux(i, j, grid, fields, evtr)   # humidity flux
         # Calculate latent heat flux (positive upwards)

From 6ac12641f64a39fc137c702e136bfca040c5ffbe Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Thu, 30 Apr 2026 14:15:36 +0200
Subject: [PATCH 12/17] Add wrapper specific heat capacity moist air

Co-authored-by: Copilot <copilot@github.com>
---
 src/processes/physical_constants.jl | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/src/processes/physical_constants.jl b/src/processes/physical_constants.jl
index 5ecfafd91..218e681a4 100644
--- a/src/processes/physical_constants.jl
+++ b/src/processes/physical_constants.jl
@@ -13,8 +13,7 @@ import Thermodynamics.Parameters:
     T_freeze,
     T_triple,
     press_triple
-import Thermodynamics: air_density
-
+import Thermodynamics: air_density, cp_m
 """
     $TYPEDEF
 
@@ -101,6 +100,13 @@ PhysicalConstants(::Type{NF}; kwargs...) where {NF} = PhysicalConstants{NF}(; kw
 @inline Rv_over_Rd(c::PhysicalConstants) = R_v(c) / R_d(c)
 @inline ε(c::PhysicalConstants) = R_d(c) / R_v(c)
 
+"""
+    specific_heat_capacity_moist_air(c::PhysicalConstants, q)
+
+Compute the isobaric specific heat capacity [J/(kg*K)] of moist air as a function 
+of the total specific humidity `q` [kg/kg]
+"""
+@inline specific_heat_capacity_moist_air(c::PhysicalConstants, q) = cp_m(c, q)
 """
     celsius_to_kelvin(c::PhysicalConstants, T)
 

From e7e5742b27740ee72014012e407600d93a093c62 Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Thu, 30 Apr 2026 14:16:48 +0200
Subject: [PATCH 13/17] Add moist air specific heat capacity for sensible heat
 flux

Co-authored-by: Copilot <copilot@github.com>
---
 src/processes/surface_energy/turbulent_fluxes.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/processes/surface_energy/turbulent_fluxes.jl b/src/processes/surface_energy/turbulent_fluxes.jl
index 2465c9f7a..a500f402b 100644
--- a/src/processes/surface_energy/turbulent_fluxes.jl
+++ b/src/processes/surface_energy/turbulent_fluxes.jl
@@ -141,12 +141,12 @@ Compute the sensible heat flux at `i, j` based on the current skin temperature a
         constants::PhysicalConstants,
         atmos::AbstractAtmosphere
     )
-    let cₐ = cp_d(constants), # specific heat capacity of air
-            rₐ = aerodynamic_resistance(i, j, grid, fields, atmos), # aerodynamic resistance
+    let rₐ = aerodynamic_resistance(i, j, grid, fields, atmos), # aerodynamic resistance
             Tₛ = skin_temperature(i, j, grid, fields, skinT), # skin temperature
             Tₐ = air_temperature(i, j, grid, fields, atmos), # air temperature
             pres = air_pressure(i, j, grid, fields, atmos),
             q_air = specific_humidity(i, j, grid, fields, atmos),
+            cₐ = specific_heat_capacity_moist_air(constants, q_air), # specific heat capacity of moist air
             # TODO: density should be evaluated at surface temperature for better accuracy
             ρₐ = air_density(constants, celsius_to_kelvin(constants, Tₐ), pres, q_air),
             Q_T = (Tₛ - Tₐ) / rₐ  # bulk aerodynamic temperature-gradient

From 955e575f29505ec14daa76c8050bf0413681a7d5 Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Thu, 30 Apr 2026 14:55:44 +0200
Subject: [PATCH 14/17] Update docstrings

Co-authored-by: Copilot <copilot@github.com>
---
 docs/make.jl                                  |  1 +
 .../atmosphere/prescribed_atmosphere.jl       |  5 +++--
 src/processes/physical_constants.jl           |  3 ++-
 src/processes/physics_utils.jl                | 13 ++++++------
 .../surface_energy/turbulent_fluxes.jl        | 21 +++++++++++--------
 5 files changed, 24 insertions(+), 19 deletions(-)

diff --git a/docs/make.jl b/docs/make.jl
index bcc9c6879..09c7e9bcb 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -146,6 +146,7 @@ links = InterLinks(
     "KernelAbstractions" => "https://juliagpu.github.io/KernelAbstractions.jl/stable/",
     "SpeedyWeather" => "https://speedyweather.github.io/SpeedyWeatherDocumentation/stable/",
     "FreezeCurves" => "https://cryogrid.github.io/FreezeCurves.jl/stable/",
+    "Thermodynamics" => "https://clima.github.io/Thermodynamics.jl/stable/",
 )
 
 
diff --git a/src/processes/atmosphere/prescribed_atmosphere.jl b/src/processes/atmosphere/prescribed_atmosphere.jl
index 928806f57..88010983f 100644
--- a/src/processes/atmosphere/prescribed_atmosphere.jl
+++ b/src/processes/atmosphere/prescribed_atmosphere.jl
@@ -107,8 +107,9 @@ variables(atmos::PrescribedAtmosphere{NF}) where {NF} = (
 """
     $SIGNATURES
 Computes the vapor pressure deficit for an air parcel at temperature `T` [°C] with 
-surface pressure `pres` [Pa] and specific humidity of air `q_air` [kg/kg].
+pressure `pres` [Pa] and specific humidity `q_air` [kg/kg].
 Assumes that air parcel is over water when `T > 0°C` and over ice when `T < 0°C`.
+Wrapper around [`vapor_pressure_deficit`](@extref Thermodynamics.vapor_pressure_deficit).
 """
 @inline function vapor_pressure_deficit(c::PhysicalConstants, T, pres, q_air)
     T_K = celsius_to_kelvin(c, T)
@@ -171,7 +172,7 @@ Retrieve or compute the specific_humidity at the current time step.
 """
     $TYPEDSIGNATURES
 
-Computes the vapor pressure deficit (VPD) at atmospheric reference level given the current atmospheric fields
+Computes the vapor pressure deficit (VPD) [Pa] at atmospheric reference level given the current atmospheric fields
 """
 @propagate_inbounds function compute_vapor_pressure_deficit(i, j, grid, fields, atmos::AbstractAtmosphere, c::PhysicalConstants)
     T_air = air_temperature(i, j, grid, fields, atmos)
diff --git a/src/processes/physical_constants.jl b/src/processes/physical_constants.jl
index 218e681a4..77b9b71a2 100644
--- a/src/processes/physical_constants.jl
+++ b/src/processes/physical_constants.jl
@@ -104,7 +104,8 @@ PhysicalConstants(::Type{NF}; kwargs...) where {NF} = PhysicalConstants{NF}(; kw
     specific_heat_capacity_moist_air(c::PhysicalConstants, q)
 
 Compute the isobaric specific heat capacity [J/(kg*K)] of moist air as a function 
-of the total specific humidity `q` [kg/kg]
+of the total specific humidity `q` [kg/kg]. Wrapper around 
+[`cp_m`](@extref Thermodynamics.cp_m). 
 """
 @inline specific_heat_capacity_moist_air(c::PhysicalConstants, q) = cp_m(c, q)
 """
diff --git a/src/processes/physics_utils.jl b/src/processes/physics_utils.jl
index a5c8befaf..b7c8439bf 100644
--- a/src/processes/physics_utils.jl
+++ b/src/processes/physics_utils.jl
@@ -42,7 +42,8 @@ end
 
 Derives specific humidity from measured relative humidity `r_h` [%], air pressure `pr` [Pa],
 air temperature `T` [°C], and physical constants `c`. Assumes saturation over ice for
-`T <= 0°C` and over liquid water otherwise.
+`T <= 0°C` and over liquid water otherwise. Wrapper around
+[`q_vap_from_RH`](@extref Thermodynamics.q_vap_from_RH).
 """
 @inline function relative_to_specific_humidity(r_h, pr, T, c::PhysicalConstants)
     T_K = celsius_to_kelvin(c, T)
@@ -54,11 +55,8 @@ end
     saturation_vapor_pressure(T)
 
 Saturation vapor pressure of an air parcel at the given temperature `T` in °C. By default, the saturation vapor
-pressure is computed over ice for `T <= 0°C` and over water for `T > 0°C`
-Coefficients of August-Roche-Magnus equation taken from [alduchovImprovedMagnusForm1996](@cite).
-
-# References
-* [alduchovImprovedMagnusForm1996](@cite) Alduchov and Eskridge, Journal of Applied Meteorology and Climatology (1996)
+pressure is computed over ice for `T <= 0°C` and over water for `T > 0°C`. Wrapper around
+[`saturation_vapor_pressure`](@extref Thermodynamics.saturation_vapor_pressure).
 """
 @inline function saturation_vapor_pressure(c::PhysicalConstants, T::NF) where {NF}
     T_K = celsius_to_kelvin(c, T)
@@ -73,7 +71,8 @@ end
     q_vap_saturation(c::PhysicalConstants, T, ρ)
 
 Saturation specific humidity at temperature `T` [°C] and density `ρ` [kg/m³]. Dispatches
-over ice for `T <= 0°C` and over liquid water otherwise.
+over ice for `T <= 0°C` and over liquid water otherwise. Wrapper around
+[`q_vap_saturation`](@extref Thermodynamics.q_vap_saturation).
 """
 @inline function q_vap_saturation(c::PhysicalConstants, T, ρ)
     T_K = celsius_to_kelvin(c, T)
diff --git a/src/processes/surface_energy/turbulent_fluxes.jl b/src/processes/surface_energy/turbulent_fluxes.jl
index a500f402b..e9f241a51 100644
--- a/src/processes/surface_energy/turbulent_fluxes.jl
+++ b/src/processes/surface_energy/turbulent_fluxes.jl
@@ -30,7 +30,7 @@ DiagnosedTurbulentFluxes(::Type{NF}) where {NF} = DiagnosedTurbulentFluxes{NF}()
 """
     $TYPEDSIGNATURES
 
-Compute the sensible heat flux as a function of the bulk aerodynamic temperature gradient `Q_T` [K m/s]
+Compute the sensible heat flux [W/m²] as a function of the bulk aerodynamic temperature gradient `Q_T` [K m/s]
 and the density `ρₐ` [kg/m³] and specific heat capacity `cₐ` [J/kg K] of air.
 """
 function compute_sensible_heat_flux(::DiagnosedTurbulentFluxes, Q_T, ρₐ, cₐ)
@@ -41,7 +41,7 @@ end
 """
     $TYPEDSIGNATURES
 
-Compute the latent heat flux as a function of the humidity flux `Q_h` [m/s], the density `ρₐ` [kg/m³] of air,
+Compute the latent heat flux [W/m²] as a function of the humidity flux `Q_h` [m/s], the density `ρₐ` [kg/m³] of air,
 and the specific latent heat of fusion `Lsl` [J/kg].
 """
 function compute_latent_heat_flux(::DiagnosedTurbulentFluxes, Q_h, ρₐ, Lsl)
@@ -52,8 +52,11 @@ end
 """
     $TYPEDSIGNATURES
     
-Computes the difference in vapor pressure between a saturated surface at temperature `T`
-and the atmosphere, defined by its specific humidity `q_air`.
+Computes the difference in vapor pressure between a saturated surface at temperature `T` [°C]
+and the atmosphere, defined by its specific humidity `q_air` [kg/kg] and pressure `p` [Pa]. 
+Relies on [Thermodynamics.jl](https://github.com/CliMA/Thermodynamics.jl) via
+[`partial_pressure_vapor`](@extref Thermodynamics.partial_pressure_vapor) and
+[`saturation_vapor_pressure`](@ref saturation_vapor_pressure).
 """
 function vapor_pressure_difference(c::PhysicalConstants, p, q_air, T)
     e_air = partial_pressure_vapor(c, p, q_air)
@@ -65,8 +68,8 @@ end
 """
     $TYPEDSIGNATURES
 
-Computes the difference in specific humidity between a saturated surface at temperature `T`
-and the atmosphere, defined by its specific humidity `q_air`.
+Computes the difference in specific humidity between a saturated surface at temperature `T` [°C]
+and the atmosphere, defined by its specific humidity `q_air` [kg/kg] and pressure `p` [Pa].
 """
 function specific_humidity_difference(c::PhysicalConstants, p, q_air, T)
     T_K = celsius_to_kelvin(c, T)
@@ -108,7 +111,7 @@ end
 """
     $TYPEDSIGNATURES
 
-Computes the specific humidity difference between a saturated surface at temperature `T` and the current atmospheric fields.
+Computes the specific humidity difference [kg/kg] between a saturated surface at temperature `T` [°C] and the current atmospheric fields.
 """
 @propagate_inbounds function compute_specific_humidity_difference(i, j, grid, fields, atmos::AbstractAtmosphere, c::PhysicalConstants, T)
     q_air = specific_humidity(i, j, grid, fields, atmos)
@@ -120,7 +123,7 @@ end
 """
     $TYPEDSIGNATURES
 
-Computes the vapor pressure difference between a saturated surface at temperature `T` and the current atmospheric fields.
+Computes the vapor pressure difference [Pa] between a saturated surface at temperature `T` [°C] and the current atmospheric fields.
 """
 @propagate_inbounds function compute_vapor_pressure_difference(i, j, grid, fields, atmos::AbstractAtmosphere, c::PhysicalConstants, T)
     q_air = specific_humidity(i, j, grid, fields, atmos)
@@ -160,7 +163,7 @@ end
     $TYPEDSIGNATURES
 
 Compute the bare ground latent heat flux at `i, j` based on the current skin temperature
-and atmospheric conditions. This imlementation assumes that evaporation is the only contributor
+and atmospheric conditions. This implementation assumes that evaporation is the only contributor
 to the latent heat flux.
 """
 @inline function compute_latent_heat_flux(

From 1162689933ddd97739c1b1bbedb5a859150ba81e Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Thu, 30 Apr 2026 15:07:47 +0200
Subject: [PATCH 15/17] Rename q_vap_saturation wrapper

---
 src/processes/physics_utils.jl                   | 4 ++--
 src/processes/surface_energy/turbulent_fluxes.jl | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/processes/physics_utils.jl b/src/processes/physics_utils.jl
index b7c8439bf..2ec0e4442 100644
--- a/src/processes/physics_utils.jl
+++ b/src/processes/physics_utils.jl
@@ -68,13 +68,13 @@ pressure is computed over ice for `T <= 0°C` and over water for `T > 0°C`. Wra
 end
 
 """
-    q_vap_saturation(c::PhysicalConstants, T, ρ)
+    saturation_specific_humidity_vapor(c::PhysicalConstants, T, ρ)
 
 Saturation specific humidity at temperature `T` [°C] and density `ρ` [kg/m³]. Dispatches
 over ice for `T <= 0°C` and over liquid water otherwise. Wrapper around
 [`q_vap_saturation`](@extref Thermodynamics.q_vap_saturation).
 """
-@inline function q_vap_saturation(c::PhysicalConstants, T, ρ)
+@inline function saturation_specific_humidity_vapor(c::PhysicalConstants, T, ρ)
     T_K = celsius_to_kelvin(c, T)
     return if T <= zero(T)
         q_vap_saturation(c, T_K, ρ, Ice())
diff --git a/src/processes/surface_energy/turbulent_fluxes.jl b/src/processes/surface_energy/turbulent_fluxes.jl
index e9f241a51..d50695379 100644
--- a/src/processes/surface_energy/turbulent_fluxes.jl
+++ b/src/processes/surface_energy/turbulent_fluxes.jl
@@ -75,7 +75,7 @@ function specific_humidity_difference(c::PhysicalConstants, p, q_air, T)
     T_K = celsius_to_kelvin(c, T)
     # TODO: should use surface air density for better accuracy
     ρₐ = air_density(c, T_K, p, q_air)
-    q_sat = q_vap_saturation(c, T, ρₐ)
+    q_sat = saturation_specific_humidity_vapor(c, T, ρₐ)
     Δq = q_sat - q_air
     return Δq
 end

From 6bed9adb33a2c76946e2c97f6816e4ad25ce1d78 Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Thu, 30 Apr 2026 16:33:36 +0200
Subject: [PATCH 16/17] Add vapor pressure to specific humidity util

Co-authored-by: Copilot <copilot@github.com>
---
 src/processes/physics_utils.jl | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/src/processes/physics_utils.jl b/src/processes/physics_utils.jl
index 2ec0e4442..e3fef5032 100644
--- a/src/processes/physics_utils.jl
+++ b/src/processes/physics_utils.jl
@@ -51,6 +51,15 @@ air temperature `T` [°C], and physical constants `c`. Assumes saturation over i
     return q_vap_from_RH(c, pr, T_K, r_h / 100, phase)
 end
 
+"""
+    vapor_pressure_to_specific_humidity(c::PhysicalConstants, e, pr)
+
+Derives specific humidity from measured vapor pressure `e` [Pa] and air pressure `pr` [Pa]. 
+"""
+@inline function vapor_pressure_to_specific_humidity(c::PhysicalConstants, e, pr)
+    return ε(c) * e / (pr - e * (1 - ε(c)))
+end
+
 """
     saturation_vapor_pressure(T)
 

From c463526d5bf8e695ad91f09cc7f2b2837184a65d Mon Sep 17 00:00:00 2001
From: Olivier Bonte <olivier.bonte@hotmail.com>
Date: Thu, 30 Apr 2026 16:55:01 +0200
Subject: [PATCH 17/17] Adapt docs

Co-authored-by: Copilot <copilot@github.com>
---
 .../surface_energy/turbulent_fluxes.md        |  7 ++--
 .../src/processes/utils/physical_constants.md | 28 +++++++++------
 docs/src/processes/utils/physics_utils.md     | 35 +++++--------------
 docs/src/references.bib                       | 14 ++++++++
 src/processes/physics_utils.jl                |  2 +-
 5 files changed, 46 insertions(+), 40 deletions(-)

diff --git a/docs/src/processes/surface_energy/turbulent_fluxes.md b/docs/src/processes/surface_energy/turbulent_fluxes.md
index d7b720b29..e65bdd93e 100644
--- a/docs/src/processes/surface_energy/turbulent_fluxes.md
+++ b/docs/src/processes/surface_energy/turbulent_fluxes.md
@@ -37,7 +37,7 @@ H_s = c_a \rho_a \frac{\Delta T}{r_a}
 \end{equation}
 ```
 
-where $c_a$ is the specific heat capacity of air (J/kg/K), $\rho_a$ is the air density (kg/m³), $\Delta T = T_s - T_a$ is the temperature difference (K) (positive if surface warmer than air), and $r_a$ is the aerodynamic resistance (s/m). $H_s$ is **positive when surface is warmer than air** (heat flows upward), and **negative when surface is cooler** (heat flows downward).
+where $c_a$ is the specific heat capacity of moist air (J/kg/K) and $\rho_a$ is the moist air density (kg/m³), both computed dynamically based on the local atmosphere state using `Thermodynamics.jl`. $\Delta T = T_s - T_a$ is the temperature difference (K) (positive if surface warmer than air), and $r_a$ is the aerodynamic resistance (s/m). $H_s$ is **positive when surface is warmer than air** (heat flows upward), and **negative when surface is cooler** (heat flows downward).
 
 
 ### Latent heat flux
@@ -50,13 +50,14 @@ Evaporation and sublimation remove heat from the surface through the latent heat
 \end{equation}
 ```
 
-with $T_s$ the skin temperature. The corresponding specific humidity difference is
+with $T_s$ the skin temperature. The specific humidity difference $\Delta q$ is 
 
 ```math
 \begin{equation}
-\Delta q = \frac{\varepsilon \Delta e}{p}.
+\Delta q = q_s - q_a = q_{\text{sat}}(T_s) - q_a
 \end{equation}
 ```
+with $q_{\text{sat}}$ the saturation specific humidity at the surface temperature and $q_a$ the specific humidity of the atmosphere. 
 
 The latent heat flux is then computed as:
 
diff --git a/docs/src/processes/utils/physical_constants.md b/docs/src/processes/utils/physical_constants.md
index dcdf35135..8d2e3cc60 100644
--- a/docs/src/processes/utils/physical_constants.md
+++ b/docs/src/processes/utils/physical_constants.md
@@ -14,7 +14,7 @@ using InteractiveUtils
 
 ## Overview
 
-[`PhysicalConstants`](@ref) collects fundamental physical constants used throughout Terrarium's process implementations. All constants are stored as fields of a single struct so that they are passed explicitly through the call graph — avoiding global state and keeping the code fully differentiable with Enzyme.jl. The struct is parametrically typed so that constants are automatically promoted to the model's numeric precision `NF`.
+[`PhysicalConstants`](@ref) collects fundamental physical constants used throughout Terrarium's process implementations. All constants are stored as fields of a single struct so that they are passed explicitly through the call graph — avoiding global state and keeping the code fully differentiable with Enzyme.jl. The struct is parametrically typed so that constants are automatically promoted to the model's numeric precision `NF`. It also subtypes `AbstractThermodynamicsParameters` to integrate directly with `Thermodynamics.jl`.
 
 ```@docs; canonical = false
 PhysicalConstants
@@ -26,18 +26,23 @@ construction to support unit-testing or sensitivity studies.
 | Field | Symbol | Default | Units | Description |
 |---|---|---|---|---|
 | `ρw` | $\rho_w$ | 1000.0 | kg/m³ | Density of liquid water |
-| `ρi` | $\rho_i$ | 916.2 | kg/m³ | Density of ice |
-| `ρₐ` | $\rho_a$ | 1.293 | kg/m³ | Density of dry air at 0°C, 1 atm |
-| `cₐ` | $c_a$ | 1005.7 | J/(kg·K) | Specific heat of dry air |
-| `Lsl` | $L_{sl}$ | 3.34×10⁵ | J/kg | Latent heat of fusion |
-| `Llg` | $L_{lv}$ | 2.257×10⁶ | J/kg | Latent heat of vaporization |
-| `Lsg` | $L_{sg}$ | 2.834×10⁶ | J/kg | Latent heat of sublimation |
+| `ρi` | $\rho_i$ | 916.7 | kg/m³ | Density of ice |
+| `cp_d` | $c_{p,d}$ | 1004.5 | J/(kg·K) | Isobaric specific heat capacity of dry air at 0°C |
+| `cp_i` | $c_{p,i}$ | 2070.0 | J/(kg·K) | Isobaric specific heat capacity of ice at 0°C |
+| `cp_l` | $c_{p,l}$ | 4186.0 | J/(kg·K) | Isobaric specific heat capacity of liquid water at 0°C |
+| `cp_v` | $c_{p,v}$ | 1846.0 | J/(kg·K) | Isobaric specific heat capacity of water vapor at 0°C|
+| `Lsl` | $L_{sl}$ | 3.34×10⁵ | J/kg | Latent heat of fusion at 0°C |
+| `Llg` | $L_{lv}$ | 2.257×10⁶ | J/kg | Latent heat of vaporization at 0°C |
+| `Lsg` | $L_{sg}$ | 2.834×10⁶ | J/kg | Latent heat of sublimation at 0°C |
 | `g` | $g$ | 9.80665 | m/s² | Gravitational acceleration |
-| `Tref` | $T_{\text{ref}}$ | 273.15 | K | Reference temperature (0°C) |
+| `T_ref` | $T_{\text{ref}}$ | 273.16 | K | Reference temperature |
+| `T_freeze` | $T_{\text{freeze}}$ | 273.16 | K | Freezing temperature of water |
+| `T_triple` | $T_{\text{triple}}$ | 273.16 | K | Triple point temperature of water |
+| `press_triple` | $p_{\text{triple}}$ | 611.657 | Pa | Triple point pressure of water |
 | `σ` | $\sigma$ | 5.6704×10⁻⁸ | W/(m²·K⁴) | Stefan-Boltzmann constant |
 | `κ` | $\kappa$ | 0.4 | — | von Kármán constant |
-| `ε` | $\varepsilon$ | 0.622 | — | Ratio of molecular weights $M_v / M_d$ |
-| `Rₐ` | $R_a$ | 287.058 | J/(kg·K) | Specific gas constant of dry air |
+| `R_d` | $R_d$ | 287.058 | J/(kg·K) | Specific gas constant of dry air |
+| `R_v` | $R_v$ | 461.5 | J/(kg·K) | Specific gas constant of water vapor |
 | `C_mass` | — | 12.0 | gC/mol | Molar mass of carbon |
 
 ## Methods
@@ -46,6 +51,9 @@ construction to support unit-testing or sensitivity studies.
 celsius_to_kelvin
 ```
 
+```@docs; canonical = false
+specific_heat_capacity_moist_air
+```
 ```@docs; canonical = false
 stefan_boltzmann
 ```
diff --git a/docs/src/processes/utils/physics_utils.md b/docs/src/processes/utils/physics_utils.md
index f7786b36d..161122a58 100644
--- a/docs/src/processes/utils/physics_utils.md
+++ b/docs/src/processes/utils/physics_utils.md
@@ -14,42 +14,25 @@ using InteractiveUtils
 
 ## Overview
 
-This module provides small, self-contained thermodynamic and atmospheric utility functions that are shared across multiple process implementations. All functions are `@inline`d and scalar-valued; they are intended to be called from within kernel functions.
+This module provides small, self-contained thermodynamic and atmospheric utility functions that are shared across multiple process implementations. All functions are `@inline`d and scalar-valued; they are intended to be called from within kernel functions. Note that these functions are mostly intended for internal use within the `Terrarium.jl` codebase, and are therefore not exported as part of the public API. However, once can still access them via `Terrarium.function_name` if needed.
 
-### Saturation vapor pressure
-
-The saturation vapor pressure $e_{\text{sat}}$ is computed using the August-Roche-Magnus empirical formula [alduchovImprovedMagnusForm1996](@cite):
-
-```math
-\begin{equation}
-e_{\text{sat}}(T) = a_1 \exp\!\left(\frac{a_2 T}{T + a_3}\right)
-\end{equation}
-```
+Some of the functionality here builds upon [Thermodynamics.jl](@extref `Thermodynamics.jl`), which provides the basic thermodynamic building blocks, based on Rankine-Kirchhoff approximations. For a mathematical overview of these approximations, see [here](@extref Thermodynamics `Formulation`).
 
-where $T$ is temperature in °C and the coefficients differ for liquid water ($T \geq 0$°C) and ice ($T < 0$°C):
+### Saturation vapor pressure
 
-| Phase | $a_1$ (Pa) | $a_2$ | $a_3$ (°C) |
-|---|---|---|---|
-| Liquid water ($T \geq 0$°C) | 611.0 | 17.62 | 243.12 |
-| Ice ($T < 0$°C) | 611.0 | 22.46 | 272.62 |
+The saturation vapor pressure $e_{\text{sat}}$ is computed using a wrapper around [`saturation_vapor_pressure`](@extref Thermodynamics.saturation_vapor_pressure) from `Thermodynamics.jl` and is based on the integration of the Clausius-Clapeyron relation (see [here](@extref Thermodynamics 9.-Saturation-Vapor-Pressure) for details). 
 
 ### Vapor pressure and humidity conversions
 
-Specific humidity $q$ and vapor pressure $e$ are related through the molecular weight ratio $\varepsilon = M_v / M_d \approx 0.622$ and the total atmospheric pressure $p$:
-
+Specific humidity $q$ and vapor pressure $e$ conversions (including related variables like relative humidity) are handled by `Thermodynamics.jl` functions (or wrappers around these functions). Specifically:
+- Transform $q$ to $e$ via [`partial_pressure_vapor`](@extref Thermodynamics.partial_pressure_vapor).
+- Transform $e$ to $q$ via [`vapor_pressure_to_specific_humidity`](@ref). With $\varepsilon = R_d / R_v$, the conversion is given by [shuttleworthTerrestrialHydrometWaterVapor2012; Eq. (2.8)](@cite) using the total air pressure $p$:
 ```math
 \begin{equation}
-q = \frac{\varepsilon \, e}{p}
+q = \frac{\varepsilon e}{p - e (1 - \varepsilon)} .
 \end{equation}
 ```
 
-The inverse conversion (vapor pressure from specific humidity) accounts for the partial pressure of dry air:
-
-```math
-\begin{equation}
-e = \frac{q \, p}{\varepsilon + (1 - \varepsilon) q}
-\end{equation}
-```
 
 ### Partial pressures of trace gases
 
@@ -77,7 +60,7 @@ saturation_vapor_pressure
 ```
 
 ```@docs; canonical = false
-vapor_pressure_deficit
+saturation_specific_humidity_vapor
 ```
 
 ```@docs; canonical = false
diff --git a/docs/src/references.bib b/docs/src/references.bib
index 07084f036..99b34efb7 100644
--- a/docs/src/references.bib
+++ b/docs/src/references.bib
@@ -302,6 +302,20 @@ @article{schaphoffLPJmL4DynamicGlobal2018
   langid = {english}
 }
 
+
+@inbook{shuttleworthTerrestrialHydrometWaterVapor2012,
+  publisher = {John Wiley \& Sons, Ltd},
+  isbn = {9781119951933},
+  author = {Shuttleworth, W. J.},
+  title = {Water Vapor in the Atmosphere},
+  booktitle = {Terrestrial Hydrometeorology},
+  chapter = {2},
+  pages = {14-24},
+  doi = {10.1002/9781119951933.ch2},
+  year = {2012},
+}
+
+
 @article{vandijkModellingInterception2001,
   title = {Modelling rainfall interception by vegetation of variable density using an adapted analytical model. Part 1. Model description},
   journal = {Journal of Hydrology},
diff --git a/src/processes/physics_utils.jl b/src/processes/physics_utils.jl
index e3fef5032..ef6610f63 100644
--- a/src/processes/physics_utils.jl
+++ b/src/processes/physics_utils.jl
@@ -1,5 +1,5 @@
 import Thermodynamics:
-    partial_pressure_vapor,
+    partial_pressure_vapor, # used as if for conversion from specific humidity to vapor pressure
     saturation_vapor_pressure,
     q_vap_from_RH,
     q_vap_saturation,