Adding Clima Thermodynamics

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"

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/",
 )
 
 

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:
 

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
 ```

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

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},

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 = ".."

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.
 """
@@ -104,24 +106,16 @@ variables(atmos::PrescribedAtmosphere{NF}) where {NF} = (
 
 """
     $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))
-
+Computes the vapor pressure deficit for an air parcel at temperature `T` [°C] with 
+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)
+    vpd = Thermodynamics.vapor_pressure_deficit(c, T_K, pres, q_air)
     return vpd
 end
-
 """
     aerodynamic_resistance(i, j, grid, fields, atmos::PrescribedAtmosphere)
 
@@ -178,13 +172,13 @@ 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)
     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, T_air, p, q_air)
     return vpd
 end
 

src/processes/physical_constants.jl

@@ -1,3 +1,19 @@
+import Thermodynamics.Parameters:
+    AbstractThermodynamicsParameters,
+    R_d,
+    R_v,
+    Rv_over_Rd,
+    cp_d,
+    cp_i,
+    cp_l,
+    cp_v,
+    LH_v0,
+    LH_s0,
+    T_0,
+    T_freeze,
+    T_triple,
+    press_triple
+import Thermodynamics: air_density, cp_m
 """
     $TYPEDEF
 
@@ -6,58 +22,98 @@ 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
+    "Isobaric specific heat capacity of dry air at standard pressure and 0°C in J/(m^3*K)"
+    cp_d::NF = 1004.5
 
-    "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 ice at standard pressure and 0°C in J/(m^3*K)"
+    cp_i::NF = 2070.0
 
-    "Sepcific latent heat of fusion of water in J/kg"
+    "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
+
+    "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"
     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
 
     "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
 
-    "Specific gas constant of air in J/(kg*K)"
-    Rₐ::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
 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
+
+# Derived parameters
+@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]. Wrapper around 
+[`cp_m`](@extref Thermodynamics.cp_m). 
+"""
+@inline specific_heat_capacity_moist_air(c::PhysicalConstants, q) = cp_m(c, q)
 """
     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, ϵ)
@@ -72,4 +128,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))