Skip to content

Thermodynamics

JAX implementations of the aerosol/atmospheric thermodynamics functions. All functions are elementwise, broadcast naturally over arrays, and are differentiable via jax.grad. They are faithful translations of the NumPy reference implementations in pyrcel.legacy.thermo.

Legacy reference

The original NumPy implementations remain available in pyrcel.legacy.thermo for cross-checking. They are not part of the v2 numerical path and are not differentiable.

sigma_w 

sigma_w(T: ArrayLike) -> Array

Surface tension of water for a given temperature.

Parameters:

Name Type Description Default
T float

Ambient temperature, K.

 required

Returns:

Type Description
float

Surface tension, J/m².
Notes

Linear fit to experimental data:

\[\sigma_w(T) = 0.0761 - 1.55 \times 10^{-4}(T - 273.15)\]
See Also

pyrcel.legacy.thermo.sigma_w

Source code in pyrcel/thermo.py 
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
def sigma_w(T: ArrayLike) -> Array:
    r"""Surface tension of water for a given temperature.

    Parameters
    ----------
    T : float
        Ambient temperature, K.

    Returns
    -------
    float
        Surface tension, J/m².

    Notes
    -----
    Linear fit to experimental data:

    $$\sigma_w(T) = 0.0761 - 1.55 \times 10^{-4}(T - 273.15)$$

    See Also
    --------
    pyrcel.legacy.thermo.sigma_w
    """
    return 0.0761 - 1.55e-4 * (T - 273.15)

es 

es(T_c: ArrayLike) -> Array

Saturation vapor pressure over water for a given temperature.

Parameters:

Name Type Description Default
T_c float

Ambient temperature, °C.

 required

Returns:

Type Description
float

Saturation vapor pressure, Pa.
Notes

Magnus formula:

\[e_s(T_c) = 611.2 \exp\!\left(\frac{17.67\, T_c}{T_c + 243.5}\right)\]
See Also

pyrcel.legacy.thermo.es

Source code in pyrcel/thermo.py 
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def es(T_c: ArrayLike) -> Array:
    r"""Saturation vapor pressure over water for a given temperature.

    Parameters
    ----------
    T_c : float
        Ambient temperature, °C.

    Returns
    -------
    float
        Saturation vapor pressure, Pa.

    Notes
    -----
    Magnus formula:

    $$e_s(T_c) = 611.2 \exp\!\left(\frac{17.67\, T_c}{T_c + 243.5}\right)$$

    See Also
    --------
    pyrcel.legacy.thermo.es
    """
    return 611.2 * jnp.exp(17.67 * T_c / (T_c + 243.5))

Seq 

Seq(r: ArrayLike, r_dry: ArrayLike, T: ArrayLike, kappa: ArrayLike) -> Array

κ-Köhler equilibrium supersaturation over an aerosol particle.

Two numerical stability improvements over the naïve formulation:

  1. Difference-of-cubes: r³ - r_dry³ is rewritten as (r - r_dry)(r² + r·r_dry + r_dry²) to avoid catastrophic cancellation when r ≈ r_dry (near the dry particle limit).
  2. Compensated final step: exp(A)·B - 1 is replaced by B·expm1(A) + (B - 1) where B - 1 = -κ·r_dry³ / denom is computed without cancellation, so the result is accurate even when A ≪ 1 (large droplets) and B ≈ 1 (dilute solution).

For kappa <= 0 the solute term is suppressed (B = 1, B - 1 = 0), so the result reduces to the pure-curvature Kelvin term expm1(A).

Parameters:

Name Type Description Default
r array or float

Droplet radius, m.

 required
r_dry array or float

Dry particle radius, m.

 required
T float

Ambient temperature, K.

 required
kappa array or float

Particle hygroscopicity parameter. For kappa <= 0 the solute term is zeroed out and only the Kelvin curvature term remains.

 required

Returns:

Type Description
array or float

Equilibrium supersaturation \(S_\mathrm{eq}\).
Notes

The full κ-Köhler equation [Petters2007]:

\[S_\mathrm{eq}(r) = \exp\!\left(\frac{A}{r}\right) \frac{r^3 - r_d^3}{r^3 - r_d^3(1 - \kappa)} - 1\]

where the Kelvin parameter is

\[A = \frac{2 M_w \sigma_w(T)}{R T \rho_w}\]
References

[Petters2007] Petters, M. D., & Kreidenweis, S. M. (2007). A single parameter representation of hygroscopic growth and cloud condensation nucleus activity. Atmos. Chem. Phys., 7, 1961–1971. doi:10.5194/acp-7-1961-2007

See Also

pyrcel.legacy.thermo.Seq

Source code in pyrcel/thermo.py 
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
def Seq(r: ArrayLike, r_dry: ArrayLike, T: ArrayLike, kappa: ArrayLike) -> Array:
    r"""κ-Köhler equilibrium supersaturation over an aerosol particle.

    Two numerical stability improvements over the naïve formulation:

    1. **Difference-of-cubes**: ``r³ - r_dry³`` is rewritten as
       ``(r - r_dry)(r² + r·r_dry + r_dry²)`` to avoid catastrophic
       cancellation when ``r ≈ r_dry`` (near the dry particle limit).
    2. **Compensated final step**: ``exp(A)·B - 1`` is replaced by
       ``B·expm1(A) + (B - 1)`` where ``B - 1 = -κ·r_dry³ / denom`` is
       computed without cancellation, so the result is accurate even when
       ``A ≪ 1`` (large droplets) and ``B ≈ 1`` (dilute solution).

    For ``kappa <= 0`` the solute term is suppressed (``B = 1``,
    ``B - 1 = 0``), so the result reduces to the pure-curvature Kelvin
    term ``expm1(A)``.

    Parameters
    ----------
    r : array or float
        Droplet radius, m.
    r_dry : array or float
        Dry particle radius, m.
    T : float
        Ambient temperature, K.
    kappa : array or float
        Particle hygroscopicity parameter.  For ``kappa <= 0`` the solute
        term is zeroed out and only the Kelvin curvature term remains.

    Returns
    -------
    array or float
        Equilibrium supersaturation $S_\mathrm{eq}$.

    Notes
    -----
    The full κ-Köhler equation [Petters2007]:

    $$S_\mathrm{eq}(r) = \exp\!\left(\frac{A}{r}\right)
      \frac{r^3 - r_d^3}{r^3 - r_d^3(1 - \kappa)} - 1$$

    where the Kelvin parameter is

    $$A = \frac{2 M_w \sigma_w(T)}{R T \rho_w}$$

    References
    ----------
    [Petters2007] Petters, M. D., & Kreidenweis, S. M. (2007). A single
    parameter representation of hygroscopic growth and cloud condensation
    nucleus activity. *Atmos. Chem. Phys.*, **7**, 1961–1971.
    doi:[10.5194/acp-7-1961-2007](https://doi.org/10.5194/acp-7-1961-2007)

    See Also
    --------
    pyrcel.legacy.thermo.Seq
    """
    A = (2.0 * c.Mw * sigma_w(T)) / (c.R * T * c.rho_w * r)

    # r³ - r_dry³ = (r - r_dry)(r² + r·r_dry + r_dry²); accurate near r ≈ r_dry.
    delta = r - r_dry
    sum_sq = r**2 + r * r_dry + r_dry**2
    num = delta * sum_sq  # r³ - r_dry³
    krd3 = kappa * r_dry**3
    denom = num + krd3  # r³ - r_dry³·(1 - κ)
    B_full = num / denom
    Bm1_full = -krd3 / denom  # B - 1, avoids cancellation when B ≈ 1

    B = jnp.where(kappa > 0.0, B_full, 1.0)
    Bm1 = jnp.where(kappa > 0.0, Bm1_full, 0.0)

    # exp(A)·B - 1 = B·expm1(A) + (B - 1); accurate when A ≪ 1.
    return B * jnp.expm1(A) + Bm1

rho_air 

rho_air(T: ArrayLike, P: ArrayLike, RH: ArrayLike = 1.0) -> Array

Density of moist air for a given temperature, pressure, and relative humidity.

Parameters:

Name Type Description Default
T float

Ambient temperature, K.

 required
P float

Ambient pressure, Pa.

 required
RH float

Relative humidity, decimal (default 1.0).

 1.0

Returns:

Type Description
float

Air density, kg/m³.
Notes

Uses the virtual temperature \(T_v = T(1 + 0.61\, q_\mathrm{sat})\) with

\[q_\mathrm{sat} = \mathrm{RH} \cdot 0.622 \cdot \frac{e_s(T)}{P}\]

so that

\[\rho = \frac{P}{R_d\, T_v}\]
See Also

pyrcel.legacy.thermo.rho_air

Source code in pyrcel/thermo.py 
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
def rho_air(T: ArrayLike, P: ArrayLike, RH: ArrayLike = 1.0) -> Array:
    r"""Density of moist air for a given temperature, pressure, and relative humidity.

    Parameters
    ----------
    T : float
        Ambient temperature, K.
    P : float
        Ambient pressure, Pa.
    RH : float, optional
        Relative humidity, decimal (default 1.0).

    Returns
    -------
    float
        Air density, kg/m³.

    Notes
    -----
    Uses the virtual temperature $T_v = T(1 + 0.61\, q_\mathrm{sat})$ with

    $$q_\mathrm{sat} = \mathrm{RH} \cdot 0.622 \cdot \frac{e_s(T)}{P}$$

    so that

    $$\rho = \frac{P}{R_d\, T_v}$$

    See Also
    --------
    pyrcel.legacy.thermo.rho_air
    """
    qsat = RH * 0.622 * (es(T - 273.15) / P)
    Tv = T * (1.0 + 0.61 * qsat)
    return P / c.Rd / Tv

ka 

ka(T: ArrayLike, rho: ArrayLike, r: ArrayLike) -> Array

Thermal conductivity of air, modified for non-continuum effects.

Parameters:

Name Type Description Default
T float

Ambient temperature, K.

 required
rho float

Ambient air density, kg/m³.

 required
r array or float

Droplet/particle radius, m.

 required

Returns:

Type Description
array or float

Non-continuum-corrected thermal conductivity, J m⁻¹ s⁻¹ K⁻¹.
Notes
\[k_a'(T, \rho, r) = \frac{k_a(T)}{1 + \dfrac{k_a(T)}{\alpha_t\, r\, \rho\, C_p} \sqrt{\dfrac{2\pi M_a}{R T}}}\]

where \(\alpha_t\) is the thermal accommodation coefficient.

See Also

pyrcel.legacy.thermo.ka

Source code in pyrcel/thermo.py 
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
def ka(T: ArrayLike, rho: ArrayLike, r: ArrayLike) -> Array:
    r"""Thermal conductivity of air, modified for non-continuum effects.

    Parameters
    ----------
    T : float
        Ambient temperature, K.
    rho : float
        Ambient air density, kg/m³.
    r : array or float
        Droplet/particle radius, m.

    Returns
    -------
    array or float
        Non-continuum-corrected thermal conductivity, J m⁻¹ s⁻¹ K⁻¹.

    Notes
    -----
    $$k_a'(T, \rho, r) = \frac{k_a(T)}{1 + \dfrac{k_a(T)}{\alpha_t\, r\, \rho\, C_p}
      \sqrt{\dfrac{2\pi M_a}{R T}}}$$

    where $\alpha_t$ is the thermal accommodation coefficient.

    See Also
    --------
    pyrcel.legacy.thermo.ka
    """
    ka_t = ka_cont(T)
    denom = 1.0 + (ka_t / (c.at * r * rho * c.Cp)) * jnp.sqrt((2.0 * PI * c.Ma) / (c.R * T))
    return ka_t / denom

dv 

dv(T: ArrayLike, r: ArrayLike, P: ArrayLike, accom: ArrayLike = c.ac) -> Array

Diffusivity of water vapor in air, modified for non-continuum effects.

Parameters:

Name Type Description Default
T float

Ambient temperature, K.

 required
r array or float

Droplet/particle radius, m.

 required
P float

Ambient pressure, Pa.

 required
accom float

Condensation coefficient (default pyrcel.constants.ac).

 ac

Returns:

Type Description
array or float

Non-continuum-corrected water vapor diffusivity, m²/s.
Notes
\[D_v'(T, r, P) = \frac{D_v(T,P)}{1 + \dfrac{D_v(T,P)}{\alpha_c\, r} \sqrt{\dfrac{2\pi M_w}{R T}}}\]

where \(\alpha_c\) is the condensation (accommodation) coefficient.

See Also

pyrcel.legacy.thermo.dv

Source code in pyrcel/thermo.py 
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
def dv(T: ArrayLike, r: ArrayLike, P: ArrayLike, accom: ArrayLike = c.ac) -> Array:
    r"""Diffusivity of water vapor in air, modified for non-continuum effects.

    Parameters
    ----------
    T : float
        Ambient temperature, K.
    r : array or float
        Droplet/particle radius, m.
    P : float
        Ambient pressure, Pa.
    accom : float, optional
        Condensation coefficient (default `pyrcel.constants.ac`).

    Returns
    -------
    array or float
        Non-continuum-corrected water vapor diffusivity, m²/s.

    Notes
    -----
    $$D_v'(T, r, P) = \frac{D_v(T,P)}{1 + \dfrac{D_v(T,P)}{\alpha_c\, r}
      \sqrt{\dfrac{2\pi M_w}{R T}}}$$

    where $\alpha_c$ is the condensation (accommodation) coefficient.

    See Also
    --------
    pyrcel.legacy.thermo.dv
    """
    dv_t = dv_cont(T, P)
    denom = 1.0 + (dv_t / (accom * r)) * jnp.sqrt((2.0 * PI * c.Mw) / (c.R * T))
    return dv_t / denom