Atmosphere Thermodynamics

Breeze implements thermodynamic relations for moist atmospheres. By "moist", we mean that the atmosphere is a mixture of four components: two gas phases (i) "dry" air and (ii) "vapor", and two "condensed" phases (iii) "liquid", and (iv) "ice". Moisture makes life interesting, because vapor can condense or solidify (and liquid can freeze) into liquid droplets and ice particles - with major consequences.

On Earth, dry air is itself a mixture of gases with fixed composition constant molar mass. Dry air on Earth's is mostly nitrogen, oxygen, and argon, whose combination produces the typical (and Breeze's default) dry air molar mass

using Breezethermo = ThermodynamicConstants()thermo.dry_air.molar_mass

0.02897

The vapor, liquid, and ice components are $\mathrm{H_2 O}$, also known as "water". Water vapor, which in Breeze has the default molar mass

thermo.vapor.molar_mass

0.018015

is lighter than dry air. As a result, moist, humid air is lighter than dry air.

Liquid in Earth's atmosphere consists of falling droplets that range from tiny, nearly-suspended mist particles to careening fat rain drops. Ice in Earth's atmosphere consists of crystals, graupel, sleet, hail and, snow.

"Moist" thermodynamic relations for a four-component mixture

What does it mean that moist air is a mixture of four components? It means that the total mass $\mathcal{M}$ of air per volume, or density, $ρ$, can be expressed as the sum of the masses of the individual components over the total volume $V$,

\[ρ = \frac{\mathcal{M}}{V} = \frac{\mathcal{M}ᵈ + \mathcal{M}ᵛ + \mathcal{M}ˡ + \mathcal{M}ⁱ}{V} = ρᵈ + ρᵛ + ρˡ + ρⁱ\]

where $\mathcal{M}ᵈ$, $\mathcal{M}ᵛ$, $\mathcal{M}ˡ$, and $\mathcal{M}ⁱ$ are the masses of dry air, vapor, liquid, and ice, respectively, while $ρᵈ$, $ρᵛ$, $ρˡ$, and $ρⁱ$ denote their fractional densities. We likewise define the mass fractions of each component,

\[qᵈ ≡ \frac{\mathcal{M}ᵈ}{\mathcal{M}} = \frac{ρᵈ}{ρ} , \qquad qᵛ ≡ \frac{\mathcal{M}ᵛ}{\mathcal{M}} = \frac{ρᵛ}{ρ} , \qquad qˡ ≡ \frac{\mathcal{M}ˡ}{\mathcal{M}} = \frac{ρˡ}{ρ}, \qquad \text{and} \qquad qⁱ ≡ \frac{ρⁱ}{ρ} = \frac{\mathcal{M}ⁱ}{\mathcal{M}} .\]

The liquid and ice components are not always present. For example, a model with warm-phase microphysics does not have ice. With no microphysics at all, there is no liquid or ice.

By definition, all of the mass fractions sum up to unity,

\[1 = qᵈ + qᵛ + qˡ + qⁱ ,\]

so that, using $qᵗ = qᵛ + qˡ + qⁱ$, the dry air mass fraction can be diagnosed with $qᵈ = 1 - qᵗ$. The sometimes tedious bookkeeping required to correctly diagnose the effective mixture properties of moist air are facilitated by Breeze's handy MoistureMassFractions abstraction. For example,

q = Breeze.Thermodynamics.MoistureMassFractions(0.01, 0.002, 1e-5)

MoistureMassFractions{Float64}: 
├── vapor:  0.01
├── liquid: 0.002
└── ice:    1.0e-5

from which we can compute the total moisture mass fraction,

qᵗ = Breeze.Thermodynamics.total_specific_moisture(q)

0.01201

And the dry as well,

qᵈ = Breeze.Thermodynamics.dry_air_mass_fraction(q)

0.98799

To be sure,

qᵈ + qᵗ

1.0

Two laws for ideal gases

Both dry air and vapor are modeled as ideal gases, which means that the ideal gas law relates pressure $p$, temperature $T$, and density $ρ$,

\[p = ρ R T .\]

Above, $R ≡ ℛ / m$ is the specific gas constant given the molar or "universal" gas constant $ℛ ≈ 8.31 \; \mathrm{J} \, \mathrm{K}^{-1} \, \mathrm{mol}^{-1}$ and molar mass $m$ of the gas species under consideration.

The first law of thermodynamics, aka "conservation of energy", states that infinitesimal changes in "heat content"^[1] $\mathrm{d} \mathcal{H}$ are related to infinitesimal changes in temperature $\mathrm{d} T$ and pressure $\mathrm{d} p$ according to:^[2]

\[\mathrm{d} \mathcal{H} = cᵖ \mathrm{d} T - \frac{\mathrm{d} p}{ρ} ,\]

where $cᵖ$ is the specific heat capacity at constant pressure of the gas in question.

For example, to represent dry air typical for Earth, with molar mass $m = 0.029 \; \mathrm{kg} \, \mathrm{mol}^{-1}$ and constant-pressure heat capacity $c^p = 1005 \; \mathrm{J} \, \mathrm{kg}^{-1} \, \mathrm{K}^{-1}$, we write

using Breeze.Thermodynamics: IdealGasdry_air = IdealGas(molar_mass=0.029, heat_capacity=1005)

IdealGas{Float64}(molar_mass=0.029, heat_capacity=1005.0)

We can also change the properties of dry air by specifying new values when constructing ThermodynamicConstants,

weird_thermo = ThermodynamicConstants(dry_air_molar_mass=0.042, dry_air_heat_capacity=420)weird_thermo.dry_air

IdealGas{Float64}(molar_mass=0.042, heat_capacity=420.0)

Potential temperature and "adiabatic" transformations

Within adiabatic transformations, $\mathrm{d} \mathcal{H} = 0$. Then, combining the ideal gas law with conservation of energy yields

\[\frac{\mathrm{d} T}{T} = \frac{R}{cᵖ} \frac{\mathrm{d} p}{p} ,\]

which implies that $T ∼ ( p / p₀ )^{R / cᵖ}$, where $p₀$ is some reference pressure value.

As a result, the potential temperature, $θ$, defined as

\[θ ≡ T \big / \left ( \frac{p}{p₀} \right )^{R / cᵖ} = \frac{T}{Π} ,\]

remains constant under adiabatic transformations. Notice that above, we also defined the Exner function, $Π ≡ ( p / p₀ )^{R / cᵖ}$.

Hydrostatic balance

Next we consider a reference state that does not exchange energy with its environment (i.e., $\mathrm{d} \mathcal{H} = 0$) and thus has constant potential temperature

\[θ₀ = Tᵣ \left ( \frac{p₀}{pᵣ} \right )^{R / cᵖ} .\]

Hydrostatic balance requires

\[∂_z pᵣ = - ρᵣ g ,\]

where $g$ is gravitational acceleration, naturally by default

thermo.gravitational_acceleration

9.81

By combining the hydrostatic balance with the ideal gas law and the definition of potential temperature we get

\[\frac{pᵣ}{p₀} = \left (1 - \frac{g z}{cᵖ θ₀} \right )^{cᵖ / R} .\]

Thus

\[\begin{align*} Tᵣ(z) & = θ₀ \left ( \frac{pᵣ}{p₀} \right )^{R / cᵖ} \\ & = θ₀ - \frac{g}{cᵖ} z, \end{align*}\]

and

\[ρᵣ(z) = \frac{p₀}{Rᵈ θ₀} \left ( 1 - \frac{g z}{cᵖ θ₀} \right )^{cᵖ / R - 1} .\]

The quantity $g / cᵖ ≈ 9.76 \;\mathrm{K}\,\mathrm{km}^{-1}$ that appears above is also referred to as the "dry adiabatic lapse rate".

An example of a dry reference state in Breeze

We can visualise a hydrostatic reference profile evaluating Breeze's reference-state utilities (which assume a dry reference state) on a one-dimensional RectilinearGrid. In the following code, the superscript $d$ denotes dry air, e.g., an ideal gas with $Rᵈ = 286.71 \; \mathrm{J} \, \mathrm{K}^{-1}$:

using Breezeusing CairoMakiegrid = RectilinearGrid(size=160, z=(0, 12_000), topology=(Flat, Flat, Bounded))thermo = ThermodynamicConstants()reference_state = ReferenceState(grid, thermo, base_pressure=101325, potential_temperature=288)pᵣ = reference_state.pressureρᵣ = reference_state.densityRᵈ = Breeze.Thermodynamics.dry_air_gas_constant(thermo)cᵖᵈ = thermo.dry_air.heat_capacityp₀ = reference_state.base_pressureθ₀ = reference_state.potential_temperatureg = thermo.gravitational_acceleration# Verify that Tᵣ = θ₀ - (g / cᵖᵈ) * zz = KernelFunctionOperation{Center, Center, Center}(znode, grid, Center(), Center(), Center())Tᵣ₁ = Field(θ₀ * (pᵣ / p₀)^(Rᵈ / cᵖᵈ))Tᵣ₂ = Field(θ₀ - (g / cᵖᵈ) * z)fig = Figure()axT = Axis(fig[1, 1]; xlabel = "Temperature (ᵒK)", ylabel = "Height (m)")lines!(axT, Tᵣ₁)lines!(axT, Tᵣ₂, linestyle = :dash, color = :orange, linewidth = 2)axp = Axis(fig[1, 2]; xlabel = "Pressure (10⁵ Pa)", yticklabelsvisible = false)lines!(axp, pᵣ / 1e5)axρ = Axis(fig[1, 3]; xlabel = "Density (kg m⁻³)", yticklabelsvisible = false)lines!(axρ, ρᵣ)fig

The gaseous nature of moist air

To define the gaseous nature of moist air - that is, the equation of state relating density and pressure, we assume that the volume of liquid and ice components are negligible. As a result, moist air pressure is the sum of partial pressures of vapor and dry air with no contribution from liquid or ice phases,

\[p = pᵈ + pᵛ .\]

Because the dry air and vapor components are ideal gases, their densities are related to pressure through the ideal gas law,

\[pᵈ = ρᵈ Rᵈ T \qquad \text{and} \qquad pᵛ = ρᵛ Rᵛ T ,\]

where $T$ is temperature, $Rⁱ = ℛ / m^β$ is the specific gas constant for component $β$, $m^β$ is the molar mass of component $β$, and $ℛ$ is the molar or "universal" gas constant,

thermo = ThermodynamicConstants()thermo.molar_gas_constant

8.314462618

ThermodynamicConstants, which is central to Breeze's implementation of moist thermodynamics. holds constants like the molar gas constant and molar masses, latent heats, gravitational acceleration, and more,

thermo

ThermodynamicConstants{Float64}:
├── molar_gas_constant: 8.314462618
├── gravitational_acceleration: 9.81
├── energy_reference_temperature: 273.15
├── triple_point_temperature: 273.16
├── triple_point_pressure: 611.657
├── dry_air: IdealGas{Float64}(molar_mass=0.02897, heat_capacity=1005.0)
├── vapor: IdealGas{Float64}(molar_mass=0.018015, heat_capacity=1850.0)
├── liquid: CondensedPhase{Float64}(reference_latent_heat=2.5008e6, heat_capacity=4181.0, density=1000.0)
├── ice: CondensedPhase{Float64}(reference_latent_heat=2.834e6, heat_capacity=2108.0, density=917.0)
└── saturation_vapor_pressure: ClausiusClapeyron()

These default values evince basic facts about water vapor air typical to Earth's atmosphere: for example, the molar masses of dry air (itself a mixture of mostly nitrogen, oxygen, and argon), and water vapor are $mᵈ = 0.029 \; \mathrm{kg} \, \mathrm{mol}^{-1}$ and $mᵛ = 0.018 \; \mathrm{kg} \, \mathrm{mol}^{-1}$. And even more interesting, the triple point temperature and pressure of water vapor are

thermo.triple_point_temperature, thermo.triple_point_pressure

(273.16, 611.657)

not so far from the typical conditions we experience on Earth's surface - one of the reasons that things are so interesting down here. Also, that temperature is not a typo: the triple point temperature really is just $+0.01^\circ$C.

It's then convenient to introduce the "mixture" gas constant $Rᵐ(qᵛ)$ such that

\[p = ρ Rᵐ T, \qquad \text{where} \qquad Rᵐ ≡ qᵈ Rᵈ + qᵛ Rᵛ .\]

To illustrate, let's compute the mixture gas constant $Rᵐ$ for air with a small amount of water vapor. The contribution of vapor increases $Rᵐ$ above the dry air value:

q = Breeze.Thermodynamics.MoistureMassFractions(0.01, 0.0, 0.0) # 1% vapor by massRᵈ = Breeze.Thermodynamics.dry_air_gas_constant(thermo)Rᵐ = Breeze.Thermodynamics.mixture_gas_constant(q, thermo)Rᵐ - Rᵈ # shows the uplift from the vapor component

1.7452747490884803

A small increase in specific humidity increases the effective gas constant of air.

The thermal properties of moist air

Though we neglect the volume of liquid and ice, we do not neglect their mass or energy. The heat capacity of moist air thus includes contributions from all four components,

\[cᵖᵐ = qᵈ cᵖᵈ + qᵛ cᵖᵛ + qˡ cˡ + qⁱ cⁱ,\]

where the $cᵖᵝ$ denote the specific heat capacity at constant pressure of constituent $β$, and we have neglected the superscript $p$ for liquid and ice because they are assumed incompressible (their specific heats and constant pressure or volume are the same). We call $cᵖᵐ$ the "mixture heat capacity", and because with default parameters the heat capacity of dry air is the smallest of either vapor, liquid, or ice, any moisture at all tends to increase the mixture heat capacity,

q = Breeze.Thermodynamics.MoistureMassFractions(0.01, 0.0, 0.0)cᵖᵈ = thermo.dry_air.heat_capacitycᵖᵐ = Breeze.Thermodynamics.mixture_heat_capacity(q, thermo)cᵖᵐ - cᵖᵈ

8.450000000000045

The Clausius–Clapeyron relation and saturation vapor pressure

The Clausius–Clapeyron relation for an ideal gas describes how saturation vapor pressure changes with temperature:

\[\frac{\mathrm{d} pᵛ⁺}{\mathrm{d} T} = \frac{pᵛ⁺ ℒ^β(T)}{Rᵛ T^2} ,\]

where $pᵛ⁺$ is saturation vapor pressure over a surface of condensed phase $β$, $T$ is temperature, $Rᵛ$ is the specific gas constant for vapor, and $ℒ^β(T)$ is the latent heat of the phase transition from vapor to phase $β$. For atmospheric moist air, the relevant condensed phases are liquid water ($β = l$) and ice ($β = i$).

Temperature-dependent latent heat

For a thermodynamic formulation that uses constant (i.e. temperature-independent) specific heats, the latent heat of a phase transition is linear in temperature:

\[ℒ^β(T) = ℒ^β_0 + \Delta c^β \, T ,\]

where $ℒ^β_0 ≡ ℒ^β(T=0)$ is the latent heat at absolute zero and $\Delta c^β ≡ c_p^v - c^β$ is the constant difference between the vapor specific heat capacity at constant pressure and the specific heat capacity of the condensed phase $β$.

Note that we typically parameterize the latent heat in terms of a reference temperature $T_r$ that is well above absolute zero. In that case, the latent heat is written

\[ℒ^β(T) = ℒ^β_r + \Delta c^β (T - T_r) \qquad \text{and} \qquad ℒ^β_0 = ℒ^β_r - \Delta c^β T_r ,\]

where $ℒ^β_r$ is the latent heat at the reference temperature $T_r$.

Integration of the Clausius-Clapeyron relation

To find the saturation vapor pressure as a function of temperature, we integrate the Clausius-Clapeyron relation with the temperature-linear latent heat model from the triple point pressure and temperature $(p^{tr}, T^{tr})$ to a generic pressure $pᵛ⁺$ and temperature $T$:

\[\int_{p^{tr}}^{pᵛ⁺} \frac{\mathrm{d} p}{p} = \int_{T^{tr}}^{T} \frac{ℒ^β_0 + \Delta c^β T'}{Rᵛ T'^2} \, \mathrm{d} T' .\]

Evaluating the integrals yields

\[\log\left(\frac{pᵛ⁺}{p^{tr}}\right) = -\frac{ℒ^β_0}{Rᵛ T} + \frac{ℒ^β_0}{Rᵛ T^{tr}} + \frac{\Delta c^β}{Rᵛ} \log\left(\frac{T}{T^{tr}}\right) .\]

Exponentiating both sides gives the closed-form solution:

\[pᵛ⁺(T) = p^{tr} \left ( \frac{T}{T^{tr}} \right )^{\Delta c^β / Rᵛ} \exp \left [ \frac{ℒ^β_0}{Rᵛ} \left (\frac{1}{T^{tr}} - \frac{1}{T} \right ) \right ] .\]

Example: liquid water and ice parameters

Consider parameters for liquid water,

using Breeze.Thermodynamics: CondensedPhaseliquid_water = CondensedPhase(reference_latent_heat=2500800, heat_capacity=4181, density=1000)

CondensedPhase{Float64}(reference_latent_heat=2.5008e6, heat_capacity=4181.0, density=1000.0)

and water ice,

water_ice = CondensedPhase(reference_latent_heat=2834000, heat_capacity=2108, density=917)

CondensedPhase{Float64}(reference_latent_heat=2.834e6, heat_capacity=2108.0, density=917.0)

These represent the latent heat of vaporization at the reference temperature, specific heat capacity, and reference density of each condensed phase. We can compute the specific heat difference $\Delta c^β$ for liquid water:

using Breeze.Thermodynamics: vapor_gas_constantcᵖᵛ = thermo.vapor.heat_capacitycˡ = thermo.liquid.heat_capacityΔcˡ = cᵖᵛ - cˡ

-2331.0

This difference $\Delta c^l$ above is negative because water vapor has a lower heat capacity than liquid water.

Mixed-phase saturation vapor pressure

In atmospheric conditions near the freezing point, condensate may exist as a mixture of liquid and ice. Following Pressel et al. (2015), we model the saturation vapor pressure over a mixed-phase surface using a liquid fraction $λ$ that varies smoothly between 0 (pure ice) and 1 (pure liquid). The effective latent heat and specific heat difference for the mixture are computed as weighted averages:

\[ℒ^{li}_0 = λ \, ℒ^l_0 + (1 - λ) \, ℒ^i_0 ,\]

\[\Delta c^{li} = λ \, \Delta c^l + (1 - λ) \, \Delta c^i .\]

These effective properties are then used in the Clausius-Clapeyron formula to compute the saturation vapor pressure over the mixed-phase surface. This approach ensures thermodynamic consistency and smooth transitions between pure liquid and pure ice states.

We can illustrate this by computing the mixed-phase specific heat difference for a 50/50 mixture:

Δcⁱ = thermo.vapor.heat_capacity - thermo.ice.heat_capacityλ = 0.5Δcˡⁱ = λ * Δcˡ + (1 - λ) * Δcⁱ

-1294.5

Visualizing saturation vapor pressure

The saturation vapor pressure over liquid, ice, and mixed-phase surfaces can be computed and visualized:

using Breezeusing Breeze.Thermodynamics: saturation_vapor_pressure, PlanarMixedPhaseSurfacethermo = ThermodynamicConstants()T = collect(200:0.1:320)pᵛˡ⁺ = [saturation_vapor_pressure(Tⁱ, thermo, thermo.liquid) for Tⁱ in T]pᵛⁱ⁺ = [saturation_vapor_pressure(Tⁱ, thermo, thermo.ice) for Tⁱ in T]# Mixed-phase surface with 50% liquid, 50% icemixed_surface = PlanarMixedPhaseSurface(0.5)pᵛᵐ⁺ = [saturation_vapor_pressure(Tⁱ, thermo, mixed_surface) for Tⁱ in T]# Mask ice and mixed-phase pressures above the freezing pointfreezing_temperature = 273.15pᵛⁱ⁺[T .> freezing_temperature] .= NaNpᵛᵐ⁺[T .> 273.15] .= NaNusing CairoMakiefig = Figure()ax = Axis(fig[1, 1], xlabel="Temperature (ᵒK)", ylabel="Saturation vapor pressure pᵛ⁺ (Pa)",          yscale = log10, xticks=200:20:320)lines!(ax, T, pᵛˡ⁺, label="liquid", linewidth=2)lines!(ax, T, pᵛⁱ⁺, label="ice", linestyle=:dash, linewidth=2)lines!(ax, T, pᵛᵐ⁺, label="mixed (λ=0.5)", linestyle=:dot, linewidth=2, color=:purple)axislegend(ax, position=:rb)fig

The mixed-phase saturation vapor pressure lies between the liquid and ice curves, providing a smooth interpolation between the two pure phases.

The Flatau polynomial fits

For performance-sensitive configurations Breeze provides FlatauPolynomial: the eighth-order polynomial fits of Flatau et al. (1992) to the liquid and ice saturation curves (the relative-error-norm coefficient sets in operational use in WRF-family microphysics). The polynomial agrees with the integrated Clausius–Clapeyron formulation to within 0.2 % over 233–313 K (liquid) while replacing a ^ and an exp with a branch-free Horner evaluation — roughly 70× cheaper per call on CPU Float64 and free of the FP64 transcendental penalty on GPUs, which matters because SaturationAdjustment evaluates the saturation curve inside a secant iteration in every cell at every stage:

constants = ThermodynamicConstants(saturation_vapor_pressure = FlatauPolynomial())

The Tetens formula for saturation vapor pressure

In addition to the first-principles ClausiusClapeyron, Breeze also supports the empirical TetensFormula, which was used in early atmosphere models due to its simplicity. We include TetensFormula solely for the purpose of model intercomparisons.

The Tetens formula approximates saturation vapor pressure as:

\[pᵛ⁺(T) = pᵛ⁺_r \exp \left( a \frac{T - T_r}{T - δT} \right) ,\]

where $Tᵣ$ is a reference temperature, $δT$ is a temperature offset, and the coefficients $a$ and $δT$ differ for liquid and ice surfaces. Default values for liquid are from Monteith and Unsworth (2014), and for ice from Murray (1967).

Let's compare the two formulations over atmospheric temperatures. We use the Clausius-Clapeyron formulation for liquid, ice, and mixed-phase surfaces, and compare with the Tetens formula for liquid and ice:

using Breezeusing Breeze.Thermodynamics: saturation_vapor_pressure,                             PlanarLiquidSurface, PlanarIceSurface, PlanarMixedPhaseSurface,                             TetensFormulausing CairoMakieclausius_clapeyron = ThermodynamicConstants()tetens = ThermodynamicConstants(saturation_vapor_pressure=TetensFormula())liquid, ice = PlanarLiquidSurface(), PlanarIceSurface()T = collect(220:0.5:320)# Clausius-Clapeyron: liquid, ice, and mixed-phase (λ=0.5)pᵛˡ⁺_cc = [saturation_vapor_pressure(Tⁱ, clausius_clapeyron, liquid) for Tⁱ in T]pᵛⁱ⁺_cc = [saturation_vapor_pressure(Tⁱ, clausius_clapeyron, ice) for Tⁱ in T]# Tetens formula: liquid and icepᵛˡ⁺_tf = [saturation_vapor_pressure(Tⁱ, tetens, liquid) for Tⁱ in T]pᵛⁱ⁺_tf = [saturation_vapor_pressure(Tⁱ, tetens, ice) for Tⁱ in T]# Mask ice above triple point for clarityTᵗʳ = clausius_clapeyron.triple_point_temperaturepᵛⁱ⁺_cc[T .> Tᵗʳ] .= NaNpᵛⁱ⁺_tf[T .> Tᵗʳ] .= NaN# Phase colors: dark blue for liquid, orange for ice, green for mixedc_liquid = :darkbluec_ice = :darkorangec_mixed = :greenfig = Figure(size=(1000, 400))# Saturation vapor pressure comparisonax1 = Axis(fig[1, 1], xlabel="Temperature (K)", ylabel="Saturation vapor pressure (Pa)",           yscale=log10, title="Saturation vapor pressure comparison")# Clausius-Clapeyronlines!(ax1, T, pᵛˡ⁺_cc, linewidth=4, color=(c_liquid, 0.6), label="C-C liquid")lines!(ax1, T, pᵛⁱ⁺_cc, linewidth=4, color=(c_ice, 0.6), label="C-C ice")# Tetens formulalines!(ax1, T, pᵛˡ⁺_tf, linewidth=2, color=c_liquid, linestyle=:dash, label="Tetens liquid")lines!(ax1, T, pᵛⁱ⁺_tf, linewidth=2, color=c_ice, linestyle=:dash, label="Tetens ice")axislegend(ax1, position=:rb)# Relative difference (Tetens - C-C) / C-Cax2 = Axis(fig[1, 2], xlabel="Temperature (K)", ylabel="Relative difference (%)",           title="(Tetens - C-C) / C-C × 100")rel_diff_liquid = @. 100 * (pᵛˡ⁺_tf - pᵛˡ⁺_cc) / pᵛˡ⁺_ccrel_diff_ice = @. 100 * (pᵛⁱ⁺_tf - pᵛⁱ⁺_cc) / pᵛⁱ⁺_cclines!(ax2, T, rel_diff_liquid, linewidth=2, color=c_liquid, label="liquid")lines!(ax2, T, rel_diff_ice, linewidth=2, color=c_ice, label="ice")hlines!(ax2, [0], color=:gray, linestyle=:dot)axislegend(ax2, position=:rb)fig

The Tetens formula agrees well with the Clausius-Clapeyron relation over typical atmospheric temperatures (roughly 230–320 K), with relative differences typically less than 1%. The Tetens formula is calibrated for this range and may diverge at extreme temperatures. For most atmospheric applications, either formulation is suitable.

Saturation specific humidity

The saturation specific humidity $qᵛ⁺$ is the maximum amount of water vapor that can exist in equilibrium with a condensed phase at a given temperature and density. It is related to the saturation vapor pressure by:

\[qᵛ⁺ ≡ \frac{ρᵛ⁺}{ρ} = \frac{pᵛ⁺}{ρ Rᵛ T} ,\]

where $ρᵛ⁺$ is the vapor density at saturation, $ρ$ is the total air density, and $Rᵛ$ is the specific gas constant for water vapor.

Visualizing saturation vapor pressure and specific humidity

We can visualize how both saturation vapor pressure and saturation specific humidity vary with temperature for different liquid fractions, demonstrating the smooth interpolation provided by the mixed-phase model:

using Breezeusing Breeze.Thermodynamics: saturation_vapor_pressure, saturation_specific_humidity, PlanarMixedPhaseSurfacethermo = ThermodynamicConstants()# Temperature range covering typical atmospheric conditionsT = collect(250:0.1:320)p₀ = 101325  # Surface pressure (Pa)Rᵈ = Breeze.Thermodynamics.dry_air_gas_constant(thermo)# Liquid fractions to visualizeλ_values = [0.0, 0.25, 0.5, 0.75, 1.0]labels = ["ice (λ=0)", "λ=0.25", "λ=0.5", "λ=0.75", "liquid (λ=1)"]colors = [:blue, :cyan, :purple, :orange, :red]linestyles = [:solid, :dash, :dot, :dashdot, :solid]using CairoMakiefig = Figure(size=(1000, 400))# Panel 1: Saturation vapor pressureax1 = Axis(fig[1, 1], xlabel="Temperature (K)", ylabel="Saturation vapor pressure (Pa)",           yscale=log10, title="Saturation vapor pressure")for (i, λ) in enumerate(λ_values)    surface = PlanarMixedPhaseSurface(λ)    pᵛ⁺ = [saturation_vapor_pressure(Tⁱ, thermo, surface) for Tⁱ in T]    lines!(ax1, T, pᵛ⁺, label=labels[i], color=colors[i], linestyle=linestyles[i], linewidth=2)endaxislegend(ax1, position=:lt)# Panel 2: Saturation specific humidityax2 = Axis(fig[1, 2], xlabel="Temperature (K)", ylabel="Saturation specific humidity (kg/kg)",           title="Saturation specific humidity")for (i, λ) in enumerate(λ_values)    surface = PlanarMixedPhaseSurface(λ)    qᵛ⁺ = zeros(length(T))    for (j, Tⁱ) in enumerate(T)        ρ = p₀ / (Rᵈ * Tⁱ)  # Approximate density using dry air        qᵛ⁺[j] = saturation_specific_humidity(Tⁱ, ρ, thermo, surface)    end    lines!(ax2, T, qᵛ⁺, label=labels[i], color=colors[i], linestyle=linestyles[i], linewidth=2)endfig

This figure shows how the liquid fraction $λ$ smoothly interpolates between pure ice ($λ = 0$) and pure liquid ($λ = 1$). At lower temperatures, the differences between phases are more pronounced. The mixed-phase model allows for realistic representation of conditions near the freezing point where both liquid and ice may coexist.

Saturation specific humidity from temperature and pressure alone

The expression $qᵛ⁺ = pᵛ⁺ / (ρ Rᵛ T)$ is circular: the density $ρ$ depends on the moisture content of the air, which is itself what we are trying to compute. Recall the mixture gas constant $Rᵐ$ introduced earlier,

\[Rᵐ = qᵈ Rᵈ + qᵛ Rᵛ , \qquad \text{so} \qquad ρ = \frac{p}{Rᵐ T} .\]

Computing $qᵛ⁺$ therefore requires knowing $qᵛ$ in advance. We can break the circularity by making an assumption about the moisture state of the air. There are two useful assumptions, each yielding a closed-form expression in $T$ and $p$ alone.

"Total moisture" assumption: all moisture is vapor at saturation

If we assume the air is saturated and contains no condensate ($qˡ = qⁱ = 0$, hence $qᵗ = qᵛ = qᵛ⁺$), the mixture gas constant becomes $Rᵐ = (1 - qᵛ⁺) Rᵈ + qᵛ⁺ Rᵛ$. Substituting into the ideal-gas relation and into $qᵛ⁺ = pᵛ⁺ / (ρ Rᵛ T)$ and solving for $qᵛ⁺$,

\[qᵛ⁺ = \frac{ϵᵈᵛ \, pᵛ⁺(T)}{p + δᵈᵛ \, pᵛ⁺(T)} ,\]

where $ϵᵈᵛ ≡ Rᵈ / Rᵛ ≈ 0.622$ is the dry-air-to-vapor molar mass ratio and $δᵈᵛ ≡ ϵᵈᵛ - 1 ≈ -0.378$. This is what saturation_total_specific_moisture computes.

This expression uses density via the ideal gas law, but resolves the circular dependency between density and specific humidity by assuming saturation. The resulting formula coincides with the saturation specific humidity used in the COARE 3.6 Edson et al. (2013) air-sea bulk-flux algorithm.

"Equilibrium" assumption: condensate may be present

If the air is in thermodynamic equilibrium with possibly some condensate ($qˡ$ and $qⁱ$ not necessarily zero, $qᵗ ≥ qᵛ⁺$, $qᵛ = qᵛ⁺$), the mixture gas constant is $Rᵐ = (1 - qᵗ) Rᵈ + qᵛ⁺ Rᵛ$ and $qᵗ$ is taken as a known input. Solving as above,

\[qᵛ⁺ = \frac{ϵᵈᵛ \, (1 - qᵗ) \, pᵛ⁺(T)}{p - pᵛ⁺(T)} ,\]

which is equation (37) of Pressel et al. (2015) and the saturated branch of equilibrium_saturation_specific_humidity. When the air is unsaturated ($qᵗ < qᵛ⁺$), equilibrium_saturation_specific_humidity falls back to the original density-based formula with $ρ$ computed using $Rᵐ = (1 - qᵗ) Rᵈ + qᵗ Rᵛ$.

Moist static energy

For moist air, a convenient thermodynamic invariant that couples temperature, composition, and height is the moist static energy (MSE),

\[s ≡ cᵖᵐ \, T + g z - ℒˡ₀ \, qˡ - ℒⁱ₀ qⁱ ,\]

where $ℒˡ₀$ and $ℒⁱ₀$ are the latent heats of condensation and deposition extrapolated to absolute zero, $ℒ^β_0 = ℒ^β_r - \Delta c^β T_r$, as derived above. Because $cᵖᵐ T$ measures the sensible enthalpy of every constituent from absolute zero, the latent terms are referenced there too. The enthalpy of vapor then exceeds that of condensate $β$ by exactly $ℒ^β(T)$, so a phase change at fixed $s$ releases the temperature-dependent latent heat: condensing a mass fraction $δ$ of vapor warms the air by $δ ℒˡ(T) / cᵖᵐ$, with $cᵖᵐ$ evaluated after the change. Pairing $cᵖᵐ T$ with the reference latent heats $ℒ^β_r$ instead would misstate that heating by $\Delta c^β T_r$ per unit mass changing phase.

Liquid-ice potential temperature

Breeze's default thermodynamic formulation, LiquidIcePotentialTemperatureFormulation, advances the liquid-ice potential temperature density $ρθ$, where $θ$ is shorthand for $θˡⁱ$. The liquid-ice potential temperature is defined from the moist enthalpy so that it is exactly invariant under any phase change at fixed pressure, and so that it reduces to the ordinary potential temperature $T / Π$ wherever there is no condensate. The price of the first property is that $θˡⁱ$ is not exactly conserved when cloudy air expands: its prognostic equation carries an explicit expansion source, proportional to the condensate. This section derives the definition, its closed-form inversion for temperature, and every term of its prognostic equation from the first law of thermodynamics, starting from the moist enthalpy.

Moist enthalpy and Kirchhoff's law

With constant heat capacities and liquid water at absolute zero as the origin of enthalpy, the specific enthalpies of the four constituents of moist air are

\[hᵈ = cᵖᵈ T , \qquad hᵛ = cᵖᵛ T + ℒˡ₀ , \qquad hˡ = cˡ T , \qquad hⁱ = cⁱ T + ℒˡ₀ - ℒⁱ₀ ,\]

where $ℒˡ₀$ and $ℒⁱ₀$ are the latent heats of condensation and deposition extrapolated to absolute zero, $ℒ^β_0 = ℒ^β_r - Δc^β T_r$, as in the section on temperature-dependent latent heat above. The enthalpy differences between vapor and the condensed phases,

\[hᵛ - hˡ = ℒˡ₀ + (cᵖᵛ - cˡ) T = ℒˡ(T) , \qquad hᵛ - hⁱ = ℒⁱ₀ + (cᵖᵛ - cⁱ) T = ℒⁱ(T) ,\]

are the temperature-dependent latent heats of Kirchhoff's law — the same $ℒ^β(T)$ that Breeze integrates in the Clausius–Clapeyron relation. The specific enthalpy of moist air is the mass-weighted sum over the constituents $β ∈ \{d, v, l, i\}$,

\[h = \sum_β q^β h^β = cᵖᵐ T + ℒˡ₀ qᵛ + (ℒˡ₀ - ℒⁱ₀) qⁱ ,\]

so that $∂h / ∂T = cᵖᵐ$ at fixed composition, and $h$ does not depend on pressure because the gases are ideal and the condensate is incompressible.

The first law in enthalpy form, per unit volume, reads

\[ρ \frac{\mathrm{D} h}{\mathrm{D} t} = \frac{\mathrm{D} p}{\mathrm{D} t} + ρ (Q + ε) - \boldsymbol{\nabla \cdot} \Big( \boldsymbol{J}_q + \sum_β h^β \boldsymbol{J}^β \Big) ,\]

where $Q$ is the heating per unit mass (for radiation, $ρ Q = Q_{\mathscr{I}}$ of the radiative transfer section), $ε$ is the dissipation of kinetic energy into heat, $\boldsymbol{J}_q$ is the sensible heat flux (conduction and its subgrid counterpart), and $\boldsymbol{J}^β$ is the flux of constituent $β$ relative to the barycentric motion — subgrid transport and the fall of condensate — with $\sum_β \boldsymbol{J}^β = 0$. Two consequences of this equation are the yardsticks for everything that follows: in adiabatic expansion at fixed composition, $\mathrm{d} h = \mathrm{d} p / ρ = Rᵐ T \, \mathrm{d} \ln p$, and in a phase change at fixed pressure, $\mathrm{d} h = 0$, whether or not the change is reversible.

Because total water is unchanged by a phase change, it is convenient to subtract $ℒˡ₀ qᵗ$ from the enthalpy and work with

\[H ≡ h - ℒˡ₀ qᵗ = cᵖᵐ T - ℒˡ₀ qˡ - ℒⁱ₀ qⁱ ,\]

which is the moist static energy without its geopotential, $s = H + g z$. $H$ changes exactly when $h$ does in every process that conserves $qᵗ$, and it is invariant under every phase change at fixed pressure.

Vapor-equivalent properties

A phase change moves mass between $qᵛ$, $qˡ$ and $qⁱ$, so it changes $cᵖᵐ$ and $Rᵐ$. The vapor-equivalent heat capacity and gas constant — the values the parcel would have if all of its water were vapor — do not:

\[\hat{c}ᵖ ≡ qᵈ cᵖᵈ + qᵗ cᵖᵛ , \qquad \hat{R} ≡ qᵈ Rᵈ + qᵗ Rᵛ , \qquad \hat{κ} ≡ \frac{\hat{R}}{\hat{c}ᵖ} , \qquad \hat{Π} ≡ \left( \frac{p}{pˢᵗ} \right)^{\hat{κ}} .\]

Throughout the documentation a hat marks such a vapor-equivalent quantity; it depends on the composition only through $qᵗ$, and equals the actual mixture value wherever there is no condensate. The actual mixture properties are the vapor-equivalent ones plus condensate corrections,

\[cᵖᵐ = \hat{c}ᵖ + (cˡ - cᵖᵛ) qˡ + (cⁱ - cᵖᵛ) qⁱ , \qquad Rᵐ = \hat{R} - Rᵛ (qˡ + qⁱ) ,\]

and substituting the first of these into $H$ together with Kirchhoff's law gives

\[H = \hat{c}ᵖ T - ℒˡ(T) \, qˡ - ℒⁱ(T) \, qⁱ .\]

Definition and inversion

The liquid-ice potential temperature is the enthalpy $H$, measured in temperature units with the vapor-equivalent heat capacity, and referred to the standard pressure with the vapor-equivalent Exner function:

\[θˡⁱ ≡ \frac{H}{\hat{c}ᵖ \hat{Π}} = \frac{1}{\hat{Π}} \left[ T - \frac{ℒˡ(T) \, qˡ + ℒⁱ(T) \, qⁱ}{\hat{c}ᵖ} \right] .\]

The second form is the familiar linear liquid-ice potential temperature with the latent heats evaluated at the actual temperature and divided by $\hat{c}ᵖ$. The two ingredients only work together: the numerator $H$ is invariant under every phase change at fixed pressure, and the denominator $\hat{c}ᵖ \hat{Π}$ depends only on $qᵗ$ and $p$, so the ratio is invariant too.

Inverting the definition for temperature is a closed-form operation, since $\hat{c}ᵖ \hat{Π} θˡⁱ = cᵖᵐ T - ℒˡ₀ qˡ - ℒⁱ₀ qⁱ$ is linear in $T$:

\[T = \frac{\hat{c}ᵖ \hat{Π} \, θˡⁱ + ℒˡ₀ qˡ + ℒⁱ₀ qⁱ}{cᵖᵐ} .\]

On the anelastic core the pressure in $\hat{Π}$ is the reference pressure $pᵣ(z)$, and the inversion costs one power. On the compressible core the pressure is itself diagnosed from density and temperature, $p = ρ Rᵐ T$, so the same relation is implicit in $T$ and is solved iteratively; see the temperature_solver keyword of LiquidIcePotentialTemperatureFormulation.

In clear air $H = cᵖᵐ T$ and $\hat{c}ᵖ = cᵖᵐ$, so $θˡⁱ = T / \hat{Π}$ is the potential temperature of unsaturated moist air with the exact exponent $Rᵐ / cᵖᵐ$ for its composition. In cloudy air $θˡⁱ$ has an equally concrete meaning. Evaporate all of a parcel's condensate at fixed pressure and fixed enthalpy: its temperature becomes $T' = H / \hat{c}ᵖ$, because clear air with the same $qᵗ$ has $h = \hat{c}ᵖ T' + ℒˡ₀ qᵗ$. Then bring it to the standard pressure along its (now condensate-free) adiabat: its potential temperature is $T' / \hat{Π} = θˡⁱ$. The liquid-ice potential temperature is therefore the potential temperature the parcel would have if its condensate were evaporated isobarically, computed exactly rather than to first order in the condensate.

The prognostic equation of a thermodynamic state function

Every thermodynamic formulation in Breeze advances some function $χ(T, p, q)$ of the temperature, the pressure and the water mass fractions $q = (qᵛ, qˡ, qⁱ)$, in place of the temperature itself. Its prognostic equation follows from the first law and the water budgets, and every source it needs is a partial derivative of its definition. The derivation below is written once for a general $χ$; the liquid-ice potential temperature is then a special case.

Each water constituent $β ∈ \{v, l, i\}$ obeys

\[ρ \frac{\mathrm{D} q^β}{\mathrm{D} t} = - \boldsymbol{\nabla \cdot} \boldsymbol{J}^β + ρ \dot{S}^β , \qquad \dot{S}^v + \dot{S}^l + \dot{S}^i = 0 ,\]

where $\dot{S}^β$ is the net rate at which phase changes produce constituent $β$, and dry air carries the balancing relative flux $\boldsymbol{J}^d = -(\boldsymbol{J}^v + \boldsymbol{J}^l + \boldsymbol{J}^i)$. At fixed pressure, $\mathrm{d} h = cᵖᵐ \, \mathrm{d} T + \sum_β (h^β - hᵈ) \, \mathrm{d} q^β$ with the sum over the water constituents, since the dry-air mass fraction takes up whatever the water constituents give away. Substituting the enthalpy and water equations and using $\sum_β \boldsymbol{J}^β = 0$ and $\dot{S}^v = -\dot{S}^l - \dot{S}^i$ gives the temperature equation

\[ρ cᵖᵐ \frac{\mathrm{D} T}{\mathrm{D} t} = \underbrace{\frac{\mathrm{D} p}{\mathrm{D} t}}_{\text{expansion}} + \underbrace{ρ \big( ℒˡ \dot{S}^l + ℒⁱ \dot{S}^i \big)}_{\text{phase change}} + \underbrace{ρ (Q + ε)}_{\text{heating}} - \underbrace{\boldsymbol{\nabla \cdot} \boldsymbol{J}_q - \sum_β (c^β - cᵖᵈ) \, \boldsymbol{J}^β \boldsymbol{\cdot \nabla} T}_{\text{conduction and heat carried by moving water}} ,\]

where $c^β$ is $cᵖᵛ$, $cˡ$ or $cⁱ$. The last term arises because the enthalpy that a relative flux carries, $(h^β - hᵈ) \boldsymbol{J}^β$, is not simply $h^β - hᵈ$ times the mass it moves when the temperature varies along the flux.

Now let $χ(T, p, q)$ be any thermodynamic state function. By the chain rule, with each partial derivative holding the other two of $(T, p, q)$ fixed and $∂ / ∂q^β$ exchanging constituent $β$ with dry air,

\[ρ \frac{\mathrm{D} χ}{\mathrm{D} t} = \frac{∂χ}{∂T} ρ \frac{\mathrm{D} T}{\mathrm{D} t} + ρ \frac{∂χ}{∂ p} \frac{\mathrm{D} p}{\mathrm{D} t} + \sum_β \frac{∂χ}{∂q^β} ρ \frac{\mathrm{D} q^β}{\mathrm{D} t} .\]

Substituting the temperature and water equations, writing $\mathrm{D} p / \mathrm{D} t = p \, \mathrm{D} \ln p / \mathrm{D} t$ and using $p / ρ = Rᵐ T$, each group of terms collects into the derivative of $χ$ along the process that drives it:

\[ρ \frac{\mathrm{D} χ}{\mathrm{D} t} = ρ \, \mathcal{E}_χ \frac{\mathrm{D} \ln p}{\mathrm{D} t} + ρ \sum_β \mathcal{P}^β_χ \dot{S}^β + χ_h \Big[ ρ (Q + ε) - \boldsymbol{\nabla \cdot} \boldsymbol{J}_q - \sum_β (c^β - cᵖᵈ) \, \boldsymbol{J}^β \boldsymbol{\cdot \nabla} T \Big] - \sum_β χ_{q^β} \boldsymbol{\nabla \cdot} \boldsymbol{J}^β ,\]

with four coefficients, each a derivative of the definition of $χ$:

The expansion coefficient is a derivative at fixed entropy $η$ because, at fixed composition, $T \, \mathrm{d} η = cᵖᵐ \, \mathrm{d} T - Rᵐ T \, \mathrm{d} \ln p$, so that $(∂T / ∂\ln p)_{η, q} = Rᵐ T / cᵖᵐ$; the phase-change coefficient is a derivative at fixed enthalpy because $0 = cᵖᵐ \, \mathrm{d} T + (h^β - hᵈ) \, \mathrm{d} q^β$ there. A conversion of constituent $β$ into $γ$ at rate $C$ contributes $(\mathcal{P}^γ_χ - \mathcal{P}^β_χ) C$, so $χ$ needs no phase-change source exactly when all of its $\mathcal{P}^β_χ$ are equal, and no expansion source exactly when $\mathcal{E}_χ = 0$.

No choice of $χ$ achieves both. Invariance under adiabatic expansion at fixed composition requires $χ$ to be a function of the mixture potential temperature $T (pˢᵗ / p)^{Rᵐ / cᵖᵐ}$ and $q$; invariance under isobaric phase change requires it to be a function of $H$, $p$ and $qᵗ$; and no function is both. Every formulation therefore computes at least one explicit source, and the choice of $χ$ decides which. The cost of omitting a term $X$ from the right side is an energy error: $χ$ is then wrong by $X \, \mathrm{d} t$, which at fixed pressure and composition is an enthalpy error of $X \, \mathrm{d} t / χ_h$, or $ρ X / χ_h$ per unit volume and time.

The equation for the liquid-ice potential temperature

For $χ = θˡⁱ = H / (\hat{c}ᵖ \hat{Π})$, the derivatives are $∂χ / ∂T = cᵖᵐ / (\hat{c}ᵖ \hat{Π})$ and $∂χ / ∂\ln p = - \hat{κ} χ$, so the four coefficients are

\[χ_h = \frac{1}{\hat{c}ᵖ \hat{Π}} , \qquad \mathcal{E}_θ = \frac{Rᵐ T - \hat{R} \hat{Π} θˡⁱ}{\hat{c}ᵖ \hat{Π}} = \frac{\hat{κ} \big( ℒˡ qˡ + ℒⁱ qⁱ \big) - Rᵛ (qˡ + qⁱ) \, T}{\hat{c}ᵖ \hat{Π}} , \qquad \mathcal{P}^v_θ = \mathcal{P}^l_θ = \mathcal{P}^i_θ ,\]

together with the composition derivatives $χ_{q^β}$ evaluated below for sedimentation. The second form of $\mathcal{E}_θ$ follows from $\hat{R} \hat{Π} θˡⁱ = \hat{κ} H$ and the vapor-equivalent decompositions of $H$ and $Rᵐ$. Both of its terms are proportional to the condensate: the expansion source vanishes in clear air, where $θˡⁱ$ is exactly conserved, and in cloud it is about $6 × 10^5 \, qˡ$ J kg⁻¹ per unit $\ln p$ — between one and two percent of the adiabatic cooling $Rᵐ T \, \mathrm{d} \ln p$ for 2 g kg⁻¹ of liquid — which is exactly the expansion work that the linearized definitions get wrong.

The equality of the phase-change coefficients is immediate from the definition rather than from the formula in the table: a phase change at fixed pressure holds $p$ and $qᵗ$ fixed, hence $\hat{c}ᵖ \hat{Π}$, and it holds $H$ fixed by the first law, hence $θˡⁱ$. This is exact for condensation, evaporation, deposition, sublimation, freezing and melting, at every temperature, in equilibrium or not. Microphysics schemes therefore never add a tendency to $ρθ$: they act on the temperature at fixed $θˡⁱ$, and the temperature after any microphysical step satisfies $cᵖᵐ \, \mathrm{d} T = ℒ^β(T) \, \mathrm{d} q^β$ automatically.

Collecting the terms, the prognostic equation for the density-weighted variable is

\[∂_t (ρ θˡⁱ) + \boldsymbol{\nabla \cdot} (ρ θˡⁱ \boldsymbol{u}) = \underbrace{ρ \, \mathcal{E}_θ \frac{\mathrm{D} \ln p}{\mathrm{D} t}}_{\text{expansion}} + \underbrace{\frac{ρ Q}{\hat{c}ᵖ \hat{Π}}}_{\text{heating}} - \boldsymbol{\nabla \cdot} \boldsymbol{J}_θ + \underbrace{\mathcal{C}_{sed}}_{\text{sedimentation}} + \underbrace{0}_{\text{phase change}} ,\]

where $\boldsymbol{J}_θ$ is the subgrid flux of $θˡⁱ$ that stands in for the conduction and water-flux terms (see What remains unconserved), and the sedimentation term is given below.

The expansion source on the two dynamical cores

The expansion source needs the material derivative of the pressure, which the two dynamical cores supply differently.

On the anelastic core the thermodynamic pressure is the hydrostatic reference pressure $pᵣ(z)$, so its material derivative is exact and local:

\[\frac{\mathrm{D} \ln pᵣ}{\mathrm{D} t} = w \, ∂_z \ln pᵣ = - \frac{g ρᵣ w}{pᵣ} .\]

The source is then $- ρᵣ \mathcal{E}_θ \, g ρᵣ w / pᵣ$: proportional to the vertical velocity, with a coefficient that depends only on the cell's own state. The $Rᵐ T / pᵣ$ inside $\mathcal{E}_θ$ is the actual specific volume, the same one that the anelastic buoyancy uses, which is what makes the sum of kinetic, potential and internal energy conserved in the sense of Pauluis (2008).

On the compressible core the pressure is prognostic in effect, and its exact material derivative follows from the equation of state $p = ρ Rᵐ T$, the continuity equation and the temperature equation:

\[\frac{\mathrm{D} \ln p}{\mathrm{D} t} = \frac{cᵖᵐ}{cᵛᵐ} \left[ - \boldsymbol{\nabla \cdot u} + \frac{\mathrm{D} \ln Rᵐ}{\mathrm{D} t} + \frac{Q^{\star}}{cᵖᵐ T} \right] , \qquad cᵛᵐ ≡ cᵖᵐ - Rᵐ ,\]

where $Q^{\star}$ collects every diabatic term of the temperature equation per unit mass: heating, latent heating, and the flux terms. This includes acoustic compression. Breeze does not evaluate it; the expansion source on the compressible core uses the hydrostatic estimate

\[\frac{\mathrm{D} \ln p}{\mathrm{D} t} ≈ - \frac{g ρ w}{p} ,\]

which drops the local pressure tendency and its horizontal advection. The alternative of keeping only the adiabatic part of the exact tendency, $-(cᵖᵐ / cᵛᵐ) \boldsymbol{\nabla \cdot u}$, is not better: heating also makes air expand, and that expansion is already accounted for, because heating is applied through the isobaric coefficient $χ_h = (∂χ / ∂h)_p$, which includes the work done during isobaric expansion. Charging the heating-driven divergence to the adiabatic term as well would count it twice, which the $Q^{\star}$ and $\mathrm{D} \ln Rᵐ / \mathrm{D} t$ terms of the exact tendency are there to subtract back out. A variable with $\mathcal{E}_χ = 0$ — the mixture potential temperature, the moist entropy — or a total energy formulation avoids the pressure tendency on this core altogether.

Heating, boundary fluxes and sedimentation

Heating at fixed composition changes $H$ by $Q$ per unit mass, so it changes $θˡⁱ$ by $Q / (\hat{c}ᵖ \hat{Π})$. The radiative heating $Q_{\mathscr{I}}$ in W m⁻³, an interior forcing supplied under the ρE key, and an energy flux in W m⁻² supplied as a boundary condition under ρE are all converted to the $ρθ$ tendency or flux by the same factor $1 / (\hat{c}ᵖ \hat{Π})$; see Boundary conditions and forcing. Water fluxes need no conversion. Adding vapor at fixed $θˡⁱ$ and fixed pressure changes the enthalpy by $θˡⁱ ∂(\hat{c}ᵖ \hat{Π}) / ∂qᵗ + ℒˡ₀$ per unit mass, which at $p = pˢᵗ$, where $\hat{Π} = 1$ and $θˡⁱ = T$ in clear air, is exactly $(cᵖᵛ - cᵖᵈ) T + ℒˡ₀ = hᵛ - hᵈ$: an evaporative surface flux delivers the enthalpy of the vapor relative to the dry air it displaces.

When condensate falls, the thermodynamic variable must account for the enthalpy that leaves with it and arrives from the cell above; the mechanism is described in the microphysics developer documentation. The content it requires is the change of $θˡⁱ$ per unit departing mass of phase $β$ at fixed temperature and pressure, along the composition increment $Δq$ that the dynamics prescribes (condensate out and dry air in on the anelastic core). Differentiating $\hat{c}ᵖ \hat{Π} θˡⁱ = H$ along that increment,

\[χ^β = \frac{ΔH|_T - θˡⁱ \, Δ(\hat{c}ᵖ \hat{Π})}{\hat{c}ᵖ \hat{Π}} , \qquad ΔH|_T = T \, Δcᵖᵐ - ℒˡ₀ Δqˡ - ℒⁱ₀ Δqⁱ , \qquad Δ(\hat{c}ᵖ \hat{Π}) = \hat{Π} \left[ Δ\hat{c}ᵖ + \hat{c}ᵖ \ln \frac{p}{pˢᵗ} \, Δ\hat{κ} \right] ,\]

with $Δ\hat{c}ᵖ = (cᵖᵛ - cᵖᵈ) Δqᵗ$ and $Δ\hat{κ} = \hat{κ} \big[ (Rᵛ - Rᵈ) / \hat{R} - (cᵖᵛ - cᵖᵈ) / \hat{c}ᵖ \big] Δqᵗ$. For liquid replaced by dry air, $ΔH|_T = (cˡ - cᵖᵈ) T - ℒˡ₀$ is the same content that the static energy transports, and $χ^l ≈ - ℒˡ(T) / (\hat{c}ᵖ \hat{Π})$ to leading order. The accompanying thermal response is $∂θˡⁱ / ∂h = 1 / (\hat{c}ᵖ \hat{Π})$. Because this Jacobian varies with the Exner function, the sedimentation tendency of $ρθ$ is a cell-local response rather than the divergence of a conservative flux: energy is consistent, but $∫ρθ$ is not conserved by precipitation, nor should it be, since heat released at one pressure is absorbed at another.

Choosing the Exner exponent

The derivation above holds for any exponent $κ$ that depends on the composition at most through $qᵗ$, with $\hat{Π}$ replaced by $Π = (p / pˢᵗ)^κ$ and $\mathcal{E}_θ = (Rᵐ T - κ \hat{c}ᵖ Π θˡⁱ) / (\hat{c}ᵖ Π)$. Phase change remains exact, since neither $\hat{c}ᵖ$ nor $Π$ changes, and expansion remains exact as long as the source is applied. LiquidIcePotentialTemperatureFormulation offers two choices through its exner_exponent keyword:

• VaporEquivalentExponent() (the default), $κ = \hat{κ}$. The expansion source vanishes in clear air, so unsaturated moist ascent conserves $θˡⁱ$ exactly, and only cloudy air carries a source. The inversion for temperature costs one power per cell, since $\hat{κ}$ varies with $qᵗ$.

• DryExponent(), $κ = κᵈ = Rᵈ / cᵖᵈ$. On the anelastic core $Πᵈ(z) = (pᵣ / pˢᵗ)^{κᵈ}$ is then a function of height alone, which can be computed once per level, and the inversion needs no transcendental function. The expansion source becomes

  \[\mathcal{E}_θ = \frac{(Rᵐ - κᵈ cᵖᵐ) \, T + κᵈ (ℒˡ₀ qˡ + ℒⁱ₀ qⁱ)}{\hat{c}ᵖ Πᵈ} ,\]

  whose first term no longer vanishes in clear air, where it reduces to $(Rᵐ - κᵈ cᵖᵐ) T = qᵛ (Rᵛ - κᵈ cᵖᵛ) T ≈ -0.23 \, qᵛ Rᵈ T$: about 0.4% of the adiabatic cooling at 15 g kg⁻¹ of vapor. This is the exponent of Betts, Tripoli and Cotton, and most cloud models, which omit the corresponding source; with the source applied it is as exact as the vapor-equivalent exponent, at the price of a source that acts everywhere.

What remains unconserved

With the exact definition and the expansion source, resolved adiabatic motion, phase change, heating and sedimentation are all energetically consistent in the continuum. Three things are not.

• Isobaric mixing. Physical mixing at one pressure conserves the mass-weighted enthalpy and total water. Since $h = \hat{c}ᵖ \hat{Π} θˡⁱ + ℒˡ₀ qᵗ$ is bilinear in $(θˡⁱ, qᵗ)$ through $\hat{c}ᵖ \hat{Π}$, mixing $θˡⁱ$ as a conserved scalar errs by $-\tfrac{1}{4} Δ(\hat{c}ᵖ \hat{Π}) Δθˡⁱ$ per equal-mass mixing event: of order 10 J kg⁻¹ for a cloudy parcel mixing with clear air, and zero for parcels of equal composition. This is negligible.

• Mixing between pressure levels. Numerical mixing by the advection scheme and subgrid mixing by the closure both move $θˡⁱ$ down its gradient between levels. Physically, such an exchange is an isentropic displacement followed by isobaric mixing, and in a stable layer the displacement costs kinetic energy. Mixing $θˡⁱ$ reaches approximately the right final thermodynamic state but spends no kinetic energy on it, so the potential energy it creates appears from nothing. Since $ρ \, \mathrm{D} H / \mathrm{D} t = - \hat{c}ᵖ \hat{Π} \, \boldsymbol{\nabla \cdot} \boldsymbol{J}_θ$ for a flux of $θˡⁱ$, the subgrid flux conserves $∫ ρ H$ only with the additional exchange term

  \[- \frac{\boldsymbol{J}_θ \boldsymbol{\cdot \nabla} (\hat{c}ᵖ \hat{Π})}{\hat{c}ᵖ \hat{Π}}\]

  in the $ρθ$ tendency, which an explicit closure could supply (a TKE closure already represents the kinetic energy that the exchange consumes through its buoyancy flux). Breeze does not apply it, and the implicit mixing of the advection scheme cannot be corrected this way at all. In the radiative-convective equilibrium simulation mentioned above, numerical mixing of $θˡⁱ$ creates about 3 W m⁻² and subgrid mixing about 0.5 W m⁻².

• Dissipative heating. Breeze omits $ε$ in every formulation: the kinetic energy that the advection scheme and the closure dissipate is not returned as heat. In a statistically steady state this loses the buoyancy work, 1–2 W m⁻² in the same simulation, from every thermodynamic variable alike.

Discretely, the expansion source uses a cell-centered vertical velocity while advection is in flux form, so the two do not cancel to machine precision in a closed domain; in a closed convecting box the remaining imbalance is about one percent of what the previous definition created.