# Effects of Turbulence on Upper-Tropospheric Ice Supersaturation

B. Kärcher, F. Hoffmann, A. Podglajen, A. Hertzog, R. Pluogonven, R. Atlas, M. Corcos, W. W. Grabowski, B. Gasparini

► To cite this version:

B. Kärcher, F. Hoffmann, A. Podglajen, A. Hertzog, R. Pluogonven, et al. Effects of Turbulence on Upper-Tropospheric Ice Supersaturation. Journal of the Atmospheric Sciences, 2024, 81, pp.1589-1604. 10.1175/JAS-D-23-0217.1 . insu-04729833

HAL Id: insu-04729833

https://insu.hal.science/insu-04729833v1

Submitted on 10 Oct 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Effects of Turbulence on Upper-Tropospheric Ice Supersaturation

B. KÄRCHER,ª F. HOFFMANN,ᵇ A. PODGLAJEN,ᶜ A. HERTZOG,ᶜ R. PLUOGONVEN,ᶜ R. ATLAS,ᶜ M. CORCOS,ᵈ
W. W. GRABOWSKI,ᶜ AND B. GASPARINIf

(Manuscript received 16 November 2023, in final form 23 May 2024, accepted 24 June 2024)

ABSTRACT: Effects of turbulence on ice supersaturation at cirrus heights (>8 km) remain unexplored. Small-scale mixing processes become important for high Reynolds number flows, which may develop below the buoyancy length scale (10–100 m). The current study couples a stochastic turbulent mixing model with reduced dimensionality to an entraining parcel model to investigate, in large-ensemble simulations, how supersaturation evolves due to homogeneous turbulence in the stably stratified, cloud-free upper troposphere. The rising parcel is forced by a mesoscale updraft. The perturbation of an initially homogeneous vertical distribution of supersaturation is studied after a 36-m ascent in a baseline case and several sensitivity scenarios. Turbulent mixing and associated temperature fluctuations alone lead to changes in ensemble-mean distributions with standard deviations in the range 0.001–0.006, while mean values are hardly affected. Large case-to-case variability in the supersaturation field is predicted with fluctuation amplitudes of up to 0.03, although such large values are rare. A vertical gradient of supersaturation (≈10⁻³ m⁻¹) is generated for high turbulence intensities due to the development of a dry-adiabatic lapse rate. Entrainment of slightly warmer (less than 0.1 K) environmental air into the parcel decreases the mean supersaturation by less than 0.01. Supersaturation fluctuations are substantially larger after entrainment events with an additional small offset in absolute humidity (by ±3.5%) between the parcel and environmental air. The predicted perturbations of ice supersaturation are significant enough to motivate studies of turbulence–ice nucleation interactions during cirrus formation that abandon the assumption of instantaneous mixing inherent to traditional parcel models.

SIGNIFICANCE STATEMENT: The purpose of our study is to investigate the effects of microscale turbulence on ice supersaturation in the upper troposphere. The associated variability in temperature and moisture fields is not resolved in cloud models and cannot easily be represented in terms of large-scale flow variables. We specify the conditions in which turbulent mixing and entrainment cause substantial variations in distributions of supersaturation. These include high turbulence intensity, strong atmospheric stability, and large moisture gradients. Our results suggest that turbulence may affect the strongly supersaturation-dependent ice formation processes in high-altitude clouds, pointing to the need to investigate cirrus formation in the presence of turbulence.

KEYWORDS: Cirrus clouds; Stochastic models; Subgrid-scale processes; Turbulence; Upper troposphere; Microscale processes/variability

**1. Introduction**

Ice supersaturation *s*, with respect to the cloud ice phase, is defined as the ratio of the partial pressure of atmospheric water vapor molecules (H₂O), *p*<sub>v</sub>, and the ice saturation vapor pressure, *p*<sub>s</sub>(*T*), where *T* is the air temperature, as defined by Murphy and Koop (2005), relative to saturation: *s* = (p<sub>v</sub>/p<sub>s</sub>) - 1.

Variability in ice supersaturation in cloud-free air is brought
about by vertical temperature variability and by horizontal vari-
ability in H₂O mixing ratios due to advective transport, convec-
tion, and gravity wave (GW) activity. Variability is enhanced
within cirrus due to ice nucleation, deposition, and sublimation

processes (Jensen et al. 2022) and remains for a time after cloud
dissipation. High-frequency aircraft measurements in the tropical
tropopause layer (TTL) point to variability in *s* down to 2-m
length scales in the vertical (Jensen et al. 2013).

GW activity in the upper troposphere (UT) is often associ-
ated with turbulence (Corcos et al. 2021; Atlas and Bretherton
2023). Turbulence sources include wind shear generation, lead-
ing to GW trapping or breaking (Fritts and Alexander 2003;
Vicari et al. 2024) and radiative-heating-driven convection
(Quante and Starr 2002; Hartmann et al. 2018), inducing an in-
ertial subrange turbulence kinetic energy cascade. While high-
resolution large-eddy simulations investigated flow topology
and turbulence statistics of stably stratified air motions in the
upper troposphere and lower stratosphere (Paoli et al. 2014),
obervational studies of ice supersaturation in the UT on even
smaller scales are not available.

A steady-state and isotropic turbulent flow field may de-
velop between the inner and outer length scales, bounding the

<footnote>ª Institut für Physik der Atmosphäre, DLR Oberpfaffenhofen, Wessling, Germany</footnote>

<footnote>ᵇMeteorologisches Institut, Ludwig-Maximilians-Universität München, Munich, Germany</footnote>

<footnote>ᶜLaboratoire de Météorologie Dynamique/IPSL, CNRS, École Polytechnique, Palaiseau, France</footnote>

<footnote>d NorthWest Research Associates, Boulder, Colorado</footnote>

<footnote>ᶜMesoscale and Microscale Meteorology, National Center for Atmospheric Research, Boulder, Colorado</footnote>

<footnote>f Department of Meteorology and Geophysics, University of Vienna, Vienna, Austria</footnote>

<footnote>Denotes content that is immediately available upon publica-
tion as open access.</footnote>

<footnote>Corresponding author: B. Kärcher, bernd.kaercher@dlr.de</footnote>
inertial subrange of turbulence. The inner scale represents the size of the smallest eddies and may be identified with the Kolmogorov microscale $\eta$ (Landau and Lifshitz 1987). At scales below $\eta$ [$C(1)$ cm in the UT], in the viscous dissipation range, mechanical flow energy transforms into heat and molecular diffusion removes sharp scalar field gradients. At scales above the outer scale, $L_o$ [$C(1)$ cm in the UT], fluid motions are impacted by buoyancy effects in a stably stratified atmosphere. Decreasing the grid size in cloud-resolving models to resolve the inner scale is not feasible in the foreseeable future due to a corresponding massive increase in computational demand.

Turbulence drives entrainment and subsequent mixing. These processes are known to exert significant effects on the lifetime and coverage of liquid-phase clouds (Hoffmann 2023), and aircraft measurements suggest that they may be relevant for cirrus clouds as well (Atlas and Bretherton 2023). While cloud model studies addressed the impact of turbulence or vertical winds on the development of cirrus and contrail cirrus clouds, often without direct communication between turbulence parameterization and ice microphysics (Gu and Liou 2000; Dobbie and Jonas 2001; Spichtinger and Gierens 2009; Dinh et al. 2010; Jensen et al. 2012; Lewellen et al. 2014; Paoli and Shariff 2016; Gasparini et al. 2019; Ohno et al. 2019; Nugent et al. 2022; Turbeville et al. 2022; Atlas et al. 2024), microscale variability in ice supersaturation in the cloud-free UT, let alone its potential impact on aerosol-mediated ice crystal nucleation, has not yet been addressed. The interaction with turbulence may also matter for the dissipation of cirrus.

Satellite observations illustrate the many scales of variability of relative humidity in the atmosphere that are difficult to reproduce in models (Gettelman et al. 2006). Large-scale ice-supersaturated regions, where cirrus and contrail cirrus clouds predominantly form and evolve, were inferred from a combination of space-borne remote sensing and aircraft measurements (Lamquin et al. 2012). Variability in UT ice supersaturation on the mesoscale due to GWs was studied with a process-oriented model constrained by in situ measurements, both within cirrus and in cloud-free conditions (Kärcher et al. 2023).

Previous research into atmospheric supersaturation variability generated by turbulence employed stochastic models with a focus on moist convection and droplet growth by condensation in liquid-phase clouds (Politovich and Cooper 1988; Jensen and Baker 1989; Kulmala et al. 1997; Krueger et al. 1997; Su et al. 1998; Lanotte et al. 2009; Paoli and Shariff 2009; Siebert and Shaw 2017; Abade et al. 2018; Chandrakar et al. 2020, 2022; Lim and Hoffmann 2023). The present study closes in on ice supersaturation fluctuations in the cloud-free UT. We investigate variability in ice supersaturation due to inertial range turbulence and entrainment in stably stratified conditions. An important goal of our study is to identify physical and numerical parameters affecting vertical ice supersaturation distributions most strongly. Our work is the first step toward developing a model to explore the role of microscale supersaturation fluctuations in ice nucleation in the homogeneous droplet freezing region.

We apply the linear eddy model (LEM) (Kerstein 1988), in which all turbulence length scales down to $\eta$ are resolved in a one-dimensional framework. The LEM is a stochastic mixing

<figure><img src="image_2.png" /><figcaption>FIG. 1. EPM setup. The vertically oriented, Lagrangian numerical LEM domain inside the parcel follows the blue temperature profile. The effect of the triplet map on a segment of an initially linear scalar profile within the domain is also shown (inset). Molecular diffusion subsequently smoothens the sharp edges of the mapped profile.</figcaption></figure>

model to which turbulence parameters are prescribed and thus allows us to address parametric sensitivities. Finite-rate turbulent mixing and molecular diffusion of both, H₂O and T, are treated distinctly. An important goal of our study is to identify physical and numerical parameters related to UT turbulence, which affect ice supersaturation distributions most strongly.

To this end, we investigate variability in $s$ brought about by fluctuations in moisture and temperature fields in the LEM that is coupled to an entraining parcel model (EPM). Vertical motion is forced by a mesoscale updraft. In a one-dimensional representation of the evolution of individual flow realizations within the parcel, we examine the impact of stochastic turbulent temperature fluctuations on the vertical distribution of $s$ down to the Kolmogorov scale. The explicit treatment of entrainment allows us to study how turbulent mixing of environmental air affects $s$ in the rising air column.

We present the EPM setup in section 2 and the LEM in section 3. Molecular diffusion processes are outlined in section 4, and simulation scenarios are described in section 5. We present and discuss our results in section 6 and conclude the paper in section 7. A list of symbols and acronyms used here is provided in Table C1 in appendix C.

## 2. Entraining parcel model setup

We seek the simplest possible model setup to explore UT ice supersaturation variability due to turbulence. For convenience of discussion, we define the ice saturation ratio, $S = s + 1$, equivalent to ice supersaturation. The basic EPM framework is shown schematically in Fig. 1.

In the EPM, an air parcel is temporarily subject to a constant updraft speed $w$ induced by small-scale fluid motions (i.e., in Kelvin–Helmholtz billows) or by GWs. In the context of our study, the air parcel is assumed to represent cloud-free UT conditions but may actually contain ice crystals few enough so that they do not affect the evolution of ice supersaturation on time scales of a few minutes relevant for our
study (section 3c). Supersaturation quenching times within cirrus clouds typically exceed 0.5 h (Kärcher et al. 2023).

Initially, temperature $T_0$ and pressure $p_0$ in the parcel and the environment are uniform. The rising parcel may experience a number of entrainment events $n_b$ at different altitudes, whereby a portion of the parcel air is replaced by one or more blobs of environmental air. In the baseline configuration, $n_b = 1$.

The environmental air pressure follows from hydrostatic balance for a constant lapse rate $\gamma$ and constrains the air pressure in the parcel: $p_e = p_0(T_e/T_0)^\kappa$, $\kappa = g/(\gamma R_g)$, with the acceleration due to gravity $g$ and the specific gas constant for dry air $R_g$. Atmospheric pressure differences are assumed to be instantaneously damped by the emission of acoustic waves, which have a typical phase speed of about 300 m s$^{-1}$. The associated time scale for a parcel dimension of, e.g., 10 m, is therefore about 30 ms. The air parcel temperature adjusts rapidly to the dry-adiabatic lapse rate, $\Gamma$: $T(h) = T_0 - \Gamma h$, as the time scale of the imposed vertical motion is much shorter than radiative relaxation times. The environmental temperature follows from $T_e(h) = T_0 - \gamma h \ge T(h)$. This implies that entrainment causes a larger temperature perturbation (warming) when it occurs at a higher altitude.

The initial value $S_0$ defines the H$_2$O mixing ratio in the parcel, $q_0 = S_0 p_s(T_0)/p_0$, which is conserved during adiabatic ascent. Thus, the mean ice saturation increases due to the uplift. We prescribe an initial ice saturation ratio in the environment $S_e$, which we set to $S_0$ in the baseline configuration, i.e., identical H$_2$O mixing ratios in the parcel and its environment. Otherwise, the latter is given by $q_e = (S_e/S_0)q_0$. In both cases, we keep $q_e$ constant, so that $S_e$ increases with altitude because $p_s(T_e)$ decreases.

Entrainment of warmer environmental air decreases the average supersaturation in the parcel to a degree determined by the total mixing fraction (ratio of entrained vs parcel air volume). Entrainment can only lead to an increase in the parcel supersaturation if $q_e$ is significantly larger than $q_0$.

The EPM is coupled to the LEM, which is discussed in more detail in section 3. This framework differs from previous applications of the LEM. Krueger (1993) used the LEM to examine the effects of mixing on buoyancy of entrained air in stratus clouds. Krueger et al. (1997) combined the LEM and entrainment in an ascending parcel model allowing condensation using saturation adjustment but no droplet microphysics. Su et al. (1998) extended the latter study to include the condensational growth of Lagrangian aerosol particles, including their activation into water droplets.

In the EPM, the computational domain (LEM grid) in the air parcel is oriented vertically and has a length $L_+$. Associated altitude levels $z \in [0, L_+]$ are defined as a set of discrete stations with constant grid spacing, $\delta z = L_- / 6$ (section 3d). Here, $L_-$ and $L_+$ represent the smallest and largest model eddy sizes (the size of an eddy represents the distance over which the fluid velocity varies substantially). Placing the domain in the vertical allows us to represent stochastic cooling and warming events due to turbulent temperature fluctuations (TTFs), caused by randomly ascending or descending eddies

of various sizes. TTFs tend to establish a dry-adiabatic lapse rate across the LEM domain.

We opt to set the entire LEM to the parcel values $q_0$ and $T_0$ in view of future model applications to cirrus ice formation, except in a control simulation designed to examine the general mixing behavior of our model in section 6a. This will ultimately allow us to compare results with those obtained based on the traditional parcel approach used to estimate nucleated ice numbers in cirrus (Lin et al. 2002). We stop the simulations at the point where homogeneous freezing of aqueous supercooled aerosol droplets would commence in the real atmosphere, i.e., around $S_{hom} = 1.5$ at 220 K (Kärcher and Lohmann 2002). Ice crystals would form and begin to deplete the supersaturated vapor, causing $S$ to relax toward ice saturation (freezing relaxation).

In the absence of H$_2$O losses and spatial gradients, $S$ evolves over the parcel time $t$ as $S \approx S_0 \exp(awt)$, with the initial value, $S(t=0) = S_0$, and a thermodynamic factor $a$ related to adiabatic expansion and contraction (Korolev and Mazin 2003). Assuming $a \approx 10^{-3}$ m$^{-1}$ to be constant in a small interval below $T_0$ defines the simulation time after which $S_{max}$ is reached:

$$t_{max} = \ln(S_{hom}/S_0)/(aw), \quad (1)$$

and the air parcel's maximum lifting distance, $H = wt_{max}$. To minimize the impact of larger scale processes, we choose $S_0 \le S_{hom}$, i.e., $S_0 = 1.45$. As $a$ is weakly T-dependent, Eq. (1) is only approximate and the maximum value, $S_{max} = 1.504$, is slightly larger than $S \approx S_{hom}$.

Altitudes of entrainment events, $h \in [0, H]$, are sampled randomly from a uniform distribution. Additional stochastic processes are treated in the LEM, which is implemented as a Monte Carlo simulation of many individual flow field realizations.

### 3. Linear eddy model

The LEM constitutes a model for the transport and mixing of diffusive tracers by turbulence. As such, it does not predict turbulence and parameters describing the turbulence are inputs that do not evolve. It has been employed in a wide range of research fields including turbulent mixing (Kerstein 1991) and reactive flow combustion studies (Menon and Kerstein 2011), cloud physics (Su et al. 1998; Hoffmann and Feingold 2019), and supernovae explosion simulations (Woosley et al. 2009), as well as in studies of gas phase chemistry and aerosol dynamics in jet aircraft exhaust plumes (Menon and Wu 1998; Wu and Menon 2001).

The LEM models turbulent mixing due to eddy motions in the inertial range by stochastic, mass- and energy-conserving scalar stirring events. The representation of a three-dimensional scalar field is interpreted as a line segment aligned locally with the scalar gradient. In other words, the LEM grid depicted in Fig. 1 can be viewed as a particular vertical line in the flow field, explicitly representing spatial variability in one dimension due to entrainment and turbulent mixing down to the smallest turbulence scale.
a. Turbulence parameters

Scalar flow constituents (such as moisture or heat) are affected by turbulent mixing and molecular diffusion. Irreversible mixing only happens at scales near the Kolmogorov scale, which is defined as

$$\eta = (\nu^3/\epsilon)^{1/4}, \quad (2)$$

and is determined by the dissipation rate of turbulence kinetic energy (TKE) per unit mass of air $\epsilon$, and the kinematic viscosity of air $\nu$.

The variance of an initial scalar gradient with characteristic length scale $\ell$ imposed on a homogeneous turbulent flow field is significantly reduced after $t_\ell$, the turbulent mixing time for an eddy of size $\ell$ (Landau and Lifshitz 1987):

$$t_\ell = (\ell^2/\epsilon)^{1/3}. \quad (3)$$

The eddy diffusivity $D_\ell$ characterizes the rate of turbulent mixing of a tracer from the scale of the largest eddies down to $\eta$. Following previous studies (e.g., Podglajen et al. 2017), we use the relationship proposed by Lilly et al. (1974):

$$D_\ell = \epsilon/(3N^2), \quad (4)$$

where $N$ is the Brunt–Väisälä (buoyancy) frequency:

$$N^2 = g(\Gamma - \gamma)/T. \quad (5)$$

Equation (4) assumes steady-state, isotropic turbulence and is an appropriate choice in stably stratified conditions ($N^2 > 0$). The buoyancy Reynolds number is defined by $Re = D_\ell/\nu$; active turbulence only appears for $Re \gg 1$.

As a measure of turbulence intensity, $\epsilon$ and $N$ are key parameters driving the LEM, since they combine to determine $D_\ell$. The UT is mostly stably stratified; local episodes of strong (active) turbulence do not occur frequently (Quante and Starr 2002; Podglajen et al. 2017; Dörnbrack et al. 2022). Vertical shear of the horizontal wind field and breaking of gravity waves dominate turbulence production. We examine observational data sources in appendix A to constrain $\epsilon$ in our simulations. Then, we prescribe a value for $N$ to estimate $D_\ell$.

It is possible to have strong shear overcoming stratification and generating turbulence despite moderate or even large values for $N$. For very weak stratification, little opposes vertical motions of air and the occurrence of stronger turbulence appears more likely. However, a universal relationship between $\epsilon$ and $N$ does not exist. As a measure of stability of air to vertical displacements, we therefore treat $N$ as an independent parameter. Finally, we derive the environmental lapse rate, $\gamma \le \Gamma$, from $N$ using Eq. (5). We recall that $\gamma$ determines the perturbation of $s$ with regard to $T$ in entrainment events (section 2).

b. Entrainment

The turbulent mixing process can be affected by discrete events of entrainment of environmental air into the LEM domain (Fig. 1). Entrainment events may be characterized by

simulations of the hydrodynamic flow field at high resolution. For example, direct numerical simulations (DNS) of stratified shear turbulence revealed highly intermittent entrainment zones at the edges of shear layers (Werne and Fritts 1999). However, harvesting results from such simulations to constrain entraining motions is not straightforward. For this reason, we treat mixing and entrainment distinctly in our study and seek the simplest possible representation of entrainment, which is assumed to leave the imposed wind speed unaffected.

When coupled to the EPM, the LEM applies a discrete set of entrained blobs of air at randomly selected parcel altitudes (section 2), in line with Krueger et al. (1997). The fraction of the LEM domain of length $L_+$ (Fig. 1) occupied by a single blob is $\beta \le 1$. In these new segments of parcel air, $q$ and $T$ are replaced by their corresponding values in the environment at the respective altitude. The location of each entrained blob within the LEM domain is also sampled randomly from a uniform distribution.

It is generally difficult to constrain precisely parameters such as $n_b$ and $\beta$. Observations, when available, may not be able to constrain them along with the domain size. For stratosphere clouds, Krueger et al. (1997) discuss how they used available observations to assign entrainment parameters. For cirrus clouds, no observations are available to determine $n_b$ and there are no fixed rules to determine $\beta$. Variations in both $\beta$ and $n_b$ change the average temperature and air mixing fraction in the parcel. Our default choices are $n_b = 1$ due to the short parcel ascent time and $\beta = 0.2$.

We do not vary $\beta$ but investigate the sensitivity of our results to the entrainment process by increasing $n_b$. For $n_b > 1$, we scale the fractional blob size, $\beta/n_b$, leaving the total fractional size of inserted blobs unchanged. We note that due to mixing and diffusion, it is possible that a new entrainment event replaces air that was already affected by a previous blob. Therefore, the total mixing fraction may be smaller than $\beta$ in the case of multiple blobs. The time scale for the turbulent breakdown of a blob, $\beta^{2/3}t_{L_+}$, follows from Eq. (3) with $\ell = \beta L_+$, where $t_{L_+}$ is the mixing time for the largest eddy represented in the LEM.

c. Model turbulent length scales

The $L_-$ and $L_+$ denote the model inner and outer turbulence length scales used in the LEM, which can differ from values for $\eta$ and $L_o$, respectively.

1) INNER SCALE

Air viscosity acts as momentum diffusivity. Thus, the associated Kolmogorov length scale quantifies the inner scale for velocity fluctuations. The inner scale for fluctuations in concentrations of $H_2O$ or $T$ may differ from $\eta$. The Schmidt number, $Sc = \nu/D_m$, relates $\nu$ to molecular diffusivity, $D_m$. For $Sc > 1$, velocity fluctuations will be damped out by viscosity at a scale larger than for the scalar fluctuations (Batchelor scale). For $Sc < 1$, the smallest scalar length scale appears in the inertial range (Obukhov–Corrsin scale). Typically, $Sc \approx 1$ in air. Similar arguments hold for thermal diffusivity.
In principle, we could scale $L_-$ = $c_1\eta$ with appropriate factors $c_1$ determined separately for H₂O and T by comparison with data. However, due to the lack of constraints on $c_1$, we set $L_-$ = $\eta$, except in sensitivity studies investigating the model resolution dependence of our results.

## 2) OUTER SCALE

We choose a value for the model outer length scale $L_+$. The true (physical) outer scale $L_o$ characterizes the size of the largest eddies in the inertial range. It is related to the Ozmidov length scale,

$$L_O = \sqrt{\epsilon/N^3}, \qquad (6)$$

which represents the largest scale at which turbulent eddies can overturn against the stratification (Riley and Lindborg 2008). At larger scales, the flow is anisotropic. The proper choice of the factor $c_2$ when setting $L_o = c_2L_O$ is largely empirical. Individual $c_2$ values range from 1 to ≈10 (Weinstock 1978; Waite 2011; Dörnbeck et al. 2022); however, $c_2$ may not even be a constant. Resulting estimates for $L_O$ are ⊕(1) m based on $N = 0.015$ s⁻¹ and typical UT values for $\epsilon$ (appendix A). Other outer length scales have been suggested, such as the buoyancy and Thorpe length scales (Wilson 2004).

Equating the eddy turnover time scale from Eq. (3) up to a constant factor with the time scale for turbulent diffusion from dimensional analysis, $t_t = \ell^2/D_t$, yields the inertial range (Kolmogorov) scaling: $D_t \propto \ell^{4/3} \epsilon^{1/3}$. Inserting $\ell = L_o \propto L_O$ gives $D_t \propto \epsilon/N^2$, consistent with the physical diffusivity, Eq. (4). Contributions of larger scales to turbulent mixing are suppressed by stable buoyancy. The scaling $D_t \propto \ell^{4/3}$ lends itself to the definition of a model Reynolds number, $Re_{LEM} = (L_+/L_-)^{4/3}$.

Given the uncertainty in determining $c_2$ and in view of future applications to cirrus ice formation, we opt to identify $L_+$ with the vertical thickness (depth) of homogeneous freezing layers $L_f$. The corresponding parcel ascent time scales $\propto 1/w$ and measures the typical duration of an entire freezing event, from first-to-last nucleation. Below we give a practical method for determining the layer depth based on microphysical modeling of homogeneous freezing events.

Values for $L_f$ observed in cirrus range from about 20 to 100 m according to aircraft measurements in the TTL (Jensen et al. 2022) but might be as low as a few meters (Jensen et al. 2013). However, homogeneous freezing layer depths cannot be well constrained by such measurements due to unknown times of observations past nucleation allowing for diffusive broadening of the layers and due to the fast horizontal aircraft motions. These factors make it hard to disentangle horizontal and vertical cloud structures, primarily affecting inferred larger $L_f$ values.

In numerical parcel simulations of homogeneous freezing of supercooled aerosol droplets (Kärcher 2017), we assume, for given $w$, that the calculated layer depth measures the altitude difference from the point where 1% of all ice crystals have formed to the point where the nucleated ice crystal number mixing ratio increases by less than 0.1%. This pragmatically captures the full width of the freezing pulse without letting its wings have a strong influence. At $p_0 = 230$ hPa and

$T_0 = 220$ K, the simulations predict layer depths 11–15 m within the range $w = 0.01–1$ m s⁻¹. (For much colder conditions in the TTL, the simulations yield a similar range, $L_f = 10–22$ m.) In stronger updrafts, air parcels reach higher levels faster, and at the same time, freezing proceeds more rapidly. Both effects approximately offset each other so that the layer depth should depend only weakly on $w$.

Taken together, we use $L_+ = L_f = 15$ m as a baseline value in the LEM. Coincidentally, this value is in line with the above estimates of $L_o$. We employ different values for $L_+$ in sensitivity studies. With these choices, all LEM parameters are defined.

### d. Turbulent mixing

Turbulent mixing is represented by stirring events, in which spatial positions of eddies along the LEM grid are altered without intermixing of their content while obeying Kolmogorov scaling laws (Kerstein 1988). Interfacial filaments created by the stirring eventually become very thin, and molecular diffusion (section 4) removes the resulting large scalar gradients (Menon and Kerstein 2011). Stirring events are not associated with an intrinsic time scale but occur at a certain frequency. While eddies take a finite time to transport scalars, stirring is implemented as a specific sequence of instantaneous permutations of fluid elements (triplet maps), which conserve fluid properties without causing discontinuities and significant error in the statistically mean fields (Kerstein 1991; Wu and Menon 2001).

In the triplet map, three copies of the segment to be stirred are made, each is compressed by a factor of 3, and the middle segment is flipped (Fig. 1). The latter ensures that the map imitates the rotational folding effect of real turbulent motions. Applying the map to a function $g(z)$ (taken to be $q$ and $T$) at any time over the interval $[z_0, z_0 + \ell]$ leads to (Kerstein 1991)

$$\hat{g}(z) = \begin{cases} g(3z - 2z_0) & z_0 \le z \le z_0 + \ell/3 \\ g(-3z + 4z_0 + 2\ell) & z_0 + \ell/3 \le z \le z_0 + 2\ell/3 \\ g(3z - 2z_0 - 2\ell) & z_0 + 2\ell/3 \le z \le z_0 + \ell, \\ g(z) & \text{else} \end{cases} \qquad (7)$$

where $z_0$ is the lower endpoint of the mapping interval and $\ell$ represents the size of an eddy. In this algorithm, eddies of size 3 (in terms of the LEM grid spacing $\delta z$) do not perturb a scalar profile, as their triplet segments collapse to points. Therefore, an eddy must be at least of size 6 and a multiple of 3. Thus, the size of LEM grid cells is $\delta z = L_-/6$. (We neglect small changes in $\delta z$ during entrainment that are brought about by differences in air density between the parcel and the environment.) The largest eddy with $z_0 = 0$ covers the entire LEM domain. The triplet map is not applied if $z_0 + \ell > L_+$. If an air parcel is displaced by $m$ grid cells, it experiences a dry-adiabatic temperature change (i.e., a TTF) of size $\pm m\Gamma\delta z$.

The boundary condition means that fewer mixing events are realized near the lower and upper domain boundaries. The lapse rate in the LEM domain approaches neutral stability in the presence of TTFs for sufficiently strong vertical mixing, implying that the dry lapse rate develops later near the
boundaries than in the interior of the domain. Physically, we are thus placing the LEM in a fluid-dynamical environment characterized by larger static stability, which acts to dampen vertical mixing at the edges of the domain.

In the next step, we estimate the frequency of stirring (rear-rangement) events based on a probability distribution $f(\ell)$ of eddy sizes. The size of a rearrangement event is equal to a selected eddy size. The turbulence diffusivity associated with a 1D random walk of a marker particle (fluid element) is given by $\mathcal{D}_t = \Lambda\ell(\ell^2)/2$, where $\langle\ell^2\rangle$ is the mean square displacement and $\Lambda\ell$ is the event frequency ($\Lambda$ is the stirring frequency per unit eddy length). If an eddy displaces a marker particle, the center of the event must lie within $\ell/2$ of that particle and its frequency is given by $\Lambda\ell$. Within the triplet map, the particle's mean square displacement is given by $\langle\ell^2\rangle = 4\ell^2/27$ (Kerstein 1991); hence, $\mathcal{D}_t = 2\Lambda\ell^3/27$.

Accounting for a probability distribution of eddy sizes, $f(\ell)$, yields the total eddy diffusivity

$$D_t = \frac{2}{27} \Lambda \int_{L_{-}}^{L_{+}} \ell^2 f(\ell) d\ell. \quad (8)$$

In Eq. (8), we have replaced $\mathcal{D}_t$ by $D_t$ from Eq. (4), because stirring by eddies of all sizes should recover the total diffusivity stochastically. Note that $D_t$ is viewed as the average turbulence diffusivity for the entire LEM domain. The eddy sizes are obtained from the normalized distribution

$$f(\ell) = \frac{5 \ell^{-8/3}}{3 L_{-}^{5/3} - L_{+}^{5/3}}, \quad \int_{L_{-}}^{L_{+}} f(\ell) d\ell = 1. \quad (9)$$

The power-law coefficient $-8/3$ is chosen such that Kolmogorov scaling, in particular $\mathcal{D}_t \propto \ell^{4/3}$, is satisfied (section 3c). Using this distribution in Eq. (8) allows us to solve for the stirring frequency per unit eddy length:

$$\Lambda = \frac{54}{5} D_t \frac{L_{-}^{5/3} - L_{+}^{5/3}}{L_{+}^{4/3} - L_{-}^{4/3}}. \quad (10)$$

The average frequency of stirring events is given by

$$v_s = \Lambda L_{+}. \quad (11)$$

Inasmuch as $D_t$ represents mixing due to all eddies up to size $L_+$, it can be shown that in the limit $L_+ \gg L_-$, $\Lambda$ becomes independent of $L_+$. In the same limit, the average eddy size in stirring events is $\int_{L_{-}}^{L_{+}} \ell f(\ell) d\ell = (5/2)L_-$. Thus, with $L_- = \eta \approx 10$ mm, the average TTF magnitude is $\Delta T = (5/2)\Gamma\eta \approx 0.25$ mK. However, when occasionally large eddies are sampled from Eq. (9), $\Delta T$ can be significantly larger.

To implement the triplet map, we first select a random location for the stirring event within the LEM domain from a uniform distribution. Then, an eddy size is chosen from $f$ by calculating the cumulative distribution,

$$F(\ell) = \int_{L_{-}}^{\ell} f(\hat{\ell}) d\hat{\ell} = \frac{\ell^{-5/3} - L_{-}^{-5/3}}{L_{+}^{-5/3} - L_{-}^{-5/3}}. \quad (12)$$

After sampling a uniformly distributed, random number $0 \le R \le 1$, the desired value of $\ell$ follows from solving $R = F(\ell)$ for $\ell$.

Stirring occurs at randomly determined times during parcel ascent. Assuming statistical independence of individual stirring events, we sample the number of stirring events after a waiting time from a Poisson probability distribution with the expectation value $v_s\Delta t$, where $\Delta t$ is the time elapsed since the last stirring event (waiting time). The minimum waiting time is set by the numerical model time step. For low turbulence intensity, the waiting time for a new stirring event to occur may cover many time steps. For strong turbulence, many stirring events may occur in a single time step.

We expect our results to be relatively insensitive to $L_+$, because the largest eddies do not contribute much to turbulent mixing. Furthermore, in the limit $L_+ \gg L_-$, we have shown that the mean size of eddies is on the order of $L_-$ and the stirring frequency depends only weakly on $L_+$. Finally, molecular diffusion is an inherently small-scale process unaffected by entrainment until scalar length scales have been reduced to near $L_-$.

## 4. Molecular diffusion

Diffusion coefficients for vapor and heat characterize scalar mixing at the dissipation scale. Pruppacher and Klett (1997) provide a relationship for the diffusivity of H₂O in air $D_v$ as a function of air pressure and temperature. Thermal diffusivity quantifies how fast heat diffuses through a material. Since both diffusion coefficients are approximately equal for the conditions of this study, we use a constant value $D_m$ for both vapor and heat diffusion, corresponding to $Sc = 0.7$ (section 3c).

While turbulent convection is represented as a random sequence of eddy motions, viscous and conductive motions are incorporated as deterministic processes governed by evolution equations defined on the LEM grid, i.e., two parabolic, linear diffusion equations: $\partial_t y = D_m \partial_z^2 y$, with $y = q$, $T$. These equations are cast into Crank–Nicholson difference schemes (Crank and Nicolson 1947). Implicit solutions are obtained by solving the resulting systems of tridiagonal matrix equations at every time step with zero gradient (no flux) boundary conditions, i.e., $\partial_z y = 0$, across the LEM domain boundaries. This is consistent with the noncyclic boundary condition in the LEM (section 3d). A numerical time step that yields sufficiently accurate solutions is given by $\delta t = 0.5(\delta z)^2/D_m$, with $\delta z = L_-/6$ (section 3d), setting the overall time step in the EPM.

Coarse LEM grid sizes, obtained by setting $L_- > \eta$, may be used to perform sensitivity studies addressing how poorer spatial resolution affects our results (section 3c). In such cases, $D_m$ is replaced by the diffusivity due to eddies with sizes up to $L_-$:

$$D_m \rightarrow D_t(L_-/L_+)^{4/3}, \quad (13)$$

with the diffusion coefficient $D_t$ from Eq. (4).

## 5. Simulation scenarios

We reiterate that input parameters to the EPM are $w$, $p_0$, $T_0$, $S_0$, and $S_e$. The LEM is driven by the parameters $\epsilon$, $D_t$, and $n_b$, where $N$ is additionally used to estimate $\gamma$ and $D_t$. Moreover, $L_- = \eta$ and $L_+ = L_f$. The environmental
parameters $q_e$ and $T_e$ are determined by $\gamma$, $T_0$, and $S_e$, respectively. The environmental air pressure follows from hydrostatic balance.

We describe parameters characterizing the simulation scenarios discussed in section 6. The maximum parcel altitude in all scenarios is $H = 36.2 \text{ m}$ (section 2). TTFs are accounted for in all scenarios, except in a discussion of the general mixing behavior of our model in section 6a. All random numbers in the EPM are derived consecutively from the Mersenne Twister algorithm (Matsumoto and Nishimura 1998).

As for the time sequence of processes in the EPM, we first determine the new parcel altitude, forcing the mean supersaturation in the parcel, and then update the environmental temperature and H₂O mixing ratio. One or more entrainment events occur according to their randomly determined times (but never at the same time), followed by a series of stirring events (eddy size selection and triplet maps) determined stochastically based on their average rate $v_s$ (section 3d). Finally, we solve the predictive equations for molecular diffusion (section 4).

The result of the kth individual simulation is the ice supersaturation distribution (SD) across the LEM grid (i.e., vertical profile), with indices $j = 1, ..., J$, at $t_{max}$: $s^{(k)}(z_j)$. For each scenario, we generate $K = 50$ different statistical representations (section 6a). We evaluate ensemble-mean SDs, $s_m(z_j) = (1/K) \sum_k s^{(k)}(z_j)$, at each LEM grid point. We also compute the corresponding domain-average supersaturation, AVG, variance, VAR, and associated standard deviation, $SDEV = \sqrt{VAR}$, and dispersion, $\delta_s = SDEV/AVG$ of the ensemble mean SD. The latter is useful as a relative measure of statistical spread as AVG increases over time due to cooling (section 2).

The baseline (Base) scenario assumes $w = 0.1 \text{ m s}^{-1}$, $p_0 = 230 \text{ hPa}$, $T_0 = 220 \text{ K}$, $S_e = S_0 = 1.45$, $\epsilon = 10^{-5} \text{ m}^2 \text{ s}^{-3}$, $N = 0.015 \text{ s}^{-1}$, $L_+ = 15 \text{ m}$, $n_b = 1$, and $\beta = 0.2$. This gives $D_t = 0.015 \text{ m}^2 \text{ s}^{-1}$, $\gamma = 4.72 \text{ km}^{-1}$, and $Re = 377$. The inner model length scale is $L_- = \eta = 8.8 \text{ mm}$, the LEM grid size is $\delta z = 1.47 \text{ mm}$, and the time step is $\delta t = 19 \text{ ms}$. The LEM domain extends over 10206 grid cells. The maximum parcel altitude is reached after $t_{max} \approx 5.6 \text{ min}$. The large-eddy mixing time in this scenario is $t_L = 4.7 \text{ min}$, somewhat shorter than $t_{max}$. The total number of stirring events is estimated from $N_s = v_s t_{max} = 57961$.

Most $D_t$ values in our simulations are smaller than the measurement estimates reported by Podglajen et al. (2017). Representing averages over an entire aircraft campaign, the latter are dominated by relatively rare, large magnitude events.

The following sensitivity studies are defined relative to scenario Base, changing one parameter at a time.

The updraft speed is varied in scenarios Wind-s ($w = 0.02 \text{ m s}^{-1}$, slow updraft) and Wind-f ($w = 0.5 \text{ m s}^{-1}$, fast updraft), encompassing a typical range of updraft speeds due to gravity waves (Kärcher and Podglajen 2019). Such updrafts are approximately constant over a fraction of a buoyancy period (a few minutes). These scenarios explore the dependence of our results on the large-scale driver of supersaturation. The simulations stop at $t_{max} = 28 \text{ min}$ (Wind-s) and $t_{max} = 1.1 \text{ min}$ (Wind-f); the total number of stirring events changes accordingly.

Scenarios Turb-l ($\epsilon = 10^{-6} \text{ m}^2 \text{ s}^{-3}$) and Turb-h ($\epsilon = 10^{-4} \text{ m}^2 \text{ s}^{-3}$) assume lower and higher values and allow us to study the impact of turbulence intensity. They encompass a range of commonly observed values in the UT (appendix A). We enhance and lower the atmospheric stability by decreasing and increasing the buoyancy frequency in scenario Stab-l ($N = 0.01 \text{ s}^{-1}$) and Stab-h ($N = 0.02 \text{ s}^{-1}$), respectively.

Scenarios Blob-0 ($n_b = 0$) and Blob-3 ($n_b = 3$) allow us to investigate how pure TTFs (no entrainment) and multiple entrainment events affect our results. We study the sensitivity with regard to moisture entrainment by assuming a slightly drier or moister environment in scenarios Env-d ($S_e = 1.4$) and Env-m ($S_e = 1.5$), respectively.

Scenarios Inner-0.1 (Inner-1) assume larger inner model length scales, hence, coarser grid spacing, $L_- = 0.1 \text{ m}$ (1 m), to explore the effect of a poorer resolution of eddy sizes. Molecular diffusion is in these scenarios represented by the diffusion due to eddies up to size $L_- > \eta$ [Eq. (13)]. Sensitivity studies Outer-5 (Outer-25) assume a smaller (larger) outer model length scale, $L_+ = 5 \text{ m}$ (25 m), to assess the sensitivity of our results on the assumed depth of a homogeneous freezing layer.

Parameters that characterize all scenarios are summarized in Table 1. These scenarios may be compared to a reference case without turbulence and entrainment. This case represents an adiabatic air parcel, wherein a spatially uniform supersaturation field evolves up to the final value $s_{max} = 0.504$ (section 2). We refer to the values of AVG and SDEV also listed in Table 1 in section 6.

# 6. Results and discussion

We present and discuss the results from the scenarios described in section 5, focusing in greater detail on scenario Base (section 6a). We then discuss scenarios dealing with variations in the vertical wind forcing, turbulence intensity, atmospheric stability, the entrainment process, and scenarios in which we vary the model length scales of turbulence (section 6b).

## a. Scenario Base

We first examine the general mixing behavior of our model by modifying scenario Base with respect to the treatment of entrainment and the consideration of TTFs. To this end, we inspect the temporal evolution of the domain-average supersaturation variance, VAR, after an initial perturbation. As for the latter, we insert one or more blobs of environmental air at the start of each simulation. All blobs have the same H₂O mixing ratio, but a slightly different temperature, 0.1 K, relative to the parcel temperature. This mimics the perturbation in T that occurs during entrainment happening around a parcel altitude $h = 20 \text{ m}$ in scenario Base.

Scalar variance is reduced by molecular diffusion. For Péclet numbers $Pe \gg 1$ ($\ll 1$), advective (diffusive) processes dominate mass transport (appendix B). In scenario Base, $Pe \gg 1$, which means that advection (i.e., turbulent mixing) accelerates the decay of H₂O variance. A similar argument can be made for heat transport. Therefore, it is meaningful to analyze ice supersaturation fluctuation variance as a function
<table><caption>TABLE 1. Model parameters for the baseline and sensitivity scenarios. Those include updraft speed $w$, environmental ice saturation ratio $S_e$, number of entrainment events $n_b$, TKE dissipation rate $\epsilon$, buoyancy frequency $N$, and the inner and outer model length scales ($L_-$ and $L_+$). We also list values for the turbulent diffusivity $D_t$, the environmental lapse rate $\gamma$, the eddy diffusion time scale $t_e$, and the expected total number of stirring events $N_s$. All scenarios assume $p_0 = 230$ hPa, $T_0 = 220$ K, $S_0 = 1.45$, and $\beta = 0.2$. Model input parameter changes relative to scenario Base are in bold. Model variables that change according to altered input parameters are italicized. The last two entries are results of our simulations: mean ice supersaturation (AVG) and associated standard deviation (SDEV) of each ensemble-mean vertical supersaturation profile taken at the final parcel altitude.</caption><thead><tr><th>Scenario</th><th>w (m s<sup>-1</sup>)</th><th>S<sub>e</sub></th><th>n<sub>b</sub></th><th>&epsilon; (m<sup>2</sup> s<sup>-3</sup>)</th><th>N (s<sup>-1</sup>)</th><th>L<sub>-</sub> (mm)</th><th>L<sub>+</sub> (m)</th><th>D<sub>t</sub> (m<sup>2</sup> s<sup>-1</sup>)</th><th>&gamma; (K km<sup>-1</sup>)</th><th>t<sub>e</sub> (min)</th><th>N<sub>s</sub></th><th>AVG</th><th>SDEV</th></tr></thead><tbody><tr><td>Base</td><td>0.1</td><td>1.45</td><td>1</td><td>10<sup>-5</sup></td><td>0.015</td><td>8.8</td><td>15</td><td>0.015</td><td>4.72</td><td>253</td><td>57 961</td><td>0.501</td><td>0.0023</td></tr><tr><td>Wind-s</td><td><b>0.02</b></td><td>1.45</td><td>1</td><td>10<sup>-5</sup></td><td>0.015</td><td>8.8</td><td>15</td><td>0.015</td><td>4.72</td><td>253</td><td>57 961</td><td>0.501</td><td>0.0046</td></tr><tr><td>Wind-f</td><td><b>0.5</b></td><td>1.45</td><td>1</td><td>10<sup>-5</sup></td><td>0.015</td><td>8.8</td><td>15</td><td>0.015</td><td>4.72</td><td>253</td><td>57 961</td><td>0.501</td><td>0.0012</td></tr><tr><td>Turb-l</td><td>0.1</td><td>1.45</td><td>1</td><td><b>10<sup>-6</sup></b></td><td>0.015</td><td>15.7</td><td>15</td><td>0.0015</td><td>4.72</td><td>2531</td><td>2220</td><td>0.501</td><td>0.0011</td></tr><tr><td>Turb-h</td><td>0.1</td><td>1.45</td><td>1</td><td><b>10<sup>-4</sup></b></td><td>0.015</td><td>5</td><td>15</td><td>0.15</td><td>4.72</td><td>25.3</td><td>1 512 851</td><td>0.501</td><td>0.0063</td></tr><tr><td>Stab-l</td><td>0.1</td><td>1.45</td><td>1</td><td>10<sup>-5</sup></td><td><b>0.01</b></td><td>8.8</td><td>15</td><td>0.033</td><td>7.52</td><td>112.5</td><td>130 413</td><td>0.503</td><td>0.0032</td></tr><tr><td>Stab-h</td><td>0.1</td><td>1.45</td><td>1</td><td>10<sup>-5</sup></td><td><b>0.02</b></td><td>8.8</td><td>15</td><td>0.0083</td><td>0.79</td><td>450</td><td>32 306</td><td>0.499</td><td>0.0027</td></tr><tr><td>Blob-0</td><td>0.1</td><td>1.45</td><td><b>0</b></td><td>10<sup>-5</sup></td><td>0.015</td><td>8.8</td><td>15</td><td>0.015</td><td>4.72</td><td>253</td><td>57 961</td><td>0.504</td><td>0.0021</td></tr><tr><td>Blob-3</td><td>0.1</td><td>1.45</td><td><b>3</b></td><td>10<sup>-5</sup></td><td>0.015</td><td>8.8</td><td>15</td><td>0.015</td><td>4.72</td><td>253</td><td>57 961</td><td>0.501</td><td>0.0020</td></tr><tr><td>Env-d</td><td>0.1</td><td><b>1.4</b></td><td>1</td><td>10<sup>-5</sup></td><td>0.015</td><td>8.8</td><td>15</td><td>0.015</td><td>4.72</td><td>253</td><td>57 961</td><td>0.490</td><td>0.0049</td></tr><tr><td>Env-m</td><td>0.1</td><td><b>1.5</b></td><td>1</td><td>10<sup>-5</sup></td><td>0.015</td><td>8.8</td><td>15</td><td>0.015</td><td>4.72</td><td>253</td><td>57 961</td><td>0.511</td><td>0.0036</td></tr><tr><td>Inner-0.1</td><td>0.1</td><td>1.45</td><td>1</td><td>10<sup>-5</sup></td><td>0.015</td><td>100</td><td>15</td><td>0.015</td><td>4.72</td><td>253</td><td>1013</td><td>0.501</td><td>0.0019</td></tr><tr><td>Inner-1</td><td>0.1</td><td>1.45</td><td>1</td><td>10<sup>-5</sup></td><td>0.015</td><td>1000</td><td>15</td><td>0.015</td><td>4.72</td><td>253</td><td>22</td><td>0.501</td><td>0.0019</td></tr><tr><td>Outer-5</td><td>0.1</td><td>1.45</td><td>1</td><td>10<sup>-5</sup></td><td>0.015</td><td>8.8</td><td><b>5</b></td><td>0.015</td><td>4.72</td><td>28.1</td><td>83 606</td><td>0.501</td><td>0.0020</td></tr><tr><td>Outer-25</td><td>0.1</td><td>1.45</td><td>1</td><td>10<sup>-5</sup></td><td>0.015</td><td>8.8</td><td><b>25</b></td><td>0.015</td><td>4.72</td><td>703</td><td>48 885</td><td>0.501</td><td>0.0013</td></tr></tbody></table>

of time scaled by the large-eddy mixing time, $\tau = t/t_{L+}$, as $t_{L+}$ is the shortest time scale governing mass transport.

To illustrate how VAR evolves over the scaled time in the EPM, we have performed LEM stand-alone Base simulations (no updraft, $w = 0$), without and with TTFs and with $n_b = 1$ and $n_b = 3$. Figure 2 shows 50 individual realizations and the resulting mean VAR evolutions. Without TTFs, VAR decays monotonically, as $s$ is treated as a completely passive scalar, i.e., $T$ fluctuates due to stirring, but the resulting adiabatic changes are set to zero. When TTFs are accounted for, VAR no longer decays monotonically. A dry-adiabatic lapse rate tends to establish across the LEM domain after sufficient turbulent mixing (section 2), which leads to a temperature gradient and, hence, maintains the temperature variance. Once neutral stability is achieved, triplet maps will have no further

impact on the temperature profile. We note that some of the larger TTFs create $VAR > 10^{-4}$ (SDEV > 0.01). Regardless of TTFs, VAR decays faster for $n_b = 3$, because more blobs offer a greater scalar surface area, which enhances mixing and molecular diffusion.

We now turn to a discussion of the unmodified scenario Base (section 5), which differs from the above (Fig. 2) in the treatment of entrainment and includes TTFs. No blobs are present initially, but one blob is inserted randomly during the simulation. Vertical profiles of supersaturation along the LEM domain, $s^{(k)}(z_j)$, taken at $h = H$ ($t = t_{max}$) are displayed in Fig. 3, both individually and as the mean over 50 statistical realizations. Note that the altitude grid points in the LEM are scaled, $\zeta = z/L_f$. The results show that at specific altitudes, supersaturation can vary up to +0.02 and -0.03 around the

<figure><img src="image_3.png" /><figcaption>FIG. 2. Variance of domain-average ice supersaturation fluctuations vs scaled time after initial blob insertion. The results (left) without TTFs and (right) with TTFs; TTFs produce spikes in the supersaturation variance as the temperature changes when a triplet map displaces an air parcel (note the altered VAR scale). The thin gray curves show 50 individual statistical realizations for $n_b = 1$. The thick black (green) curves show the averages over all realizations for $n_b = 1$ (3). Note that the VAR axes have different scales.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3. (left) Vertical profiles of ice supersaturation vs scaled altitude grid in air parcels at the final lifting distance for scenario Base, which includes TTFs and one entrainment event at a randomly selected parcel altitude. The thin gray curves show 50 individual statistical realizations taken at the final parcel altitude, the thick black curve is the distribution averaged over all 50 realizations, the colored curves highlight three individual cases, and the thick black line illustrates the solution for a dry-adiabatic lapse rate. (right) The dispersion of the distributions averaged over *K* up to *K* = 100 realizations. The results (circles) are connected with dotted lines to guide the eye.</figcaption></figure>

domain-average value. Such variations are in line with results from cloud-resolving simulations employing turbulence closure schemes to model turbulent fluctuations in stably stratified situations with a grid spacing of 100 m (Gu and Liou 2000). Absolute changes in *s* up to 0.03 occur rarely and represent the largest excursions seen in the simulations. Due to the limited sample size, we may not capture all statistical outliers present in LEM property fields. While such large excursions are visible, e.g., in the red distribution, the mean SD averaged over 50 realizations shows only small scatter. Figure 3 also shows the solution for a dry-adiabatic lapse rate, representing the limiting case of strong vertical mixing. We recall that mixing is incomplete as $t_{Ls} < t_{max}$ (section 5).

Domain-averaged mean value and standard variation of the associated mean SD are AVG = 0.501 and SDEV = 2.3 × 10⁻³, respectively. The term AVG is slightly lower than the final mean supersaturation, $s_{max}$ = 0.504, in the reference case without turbulence and entrainment (section 5), due to the single entrainment event in scenario Base. The evolution of the distribution spread in terms of dispersion (relative SDEV) as a function of total number of statistical realizations *K* shows that $\delta_s$ decreases

notably after $K \approx 30$ and achieves lower values, $\delta_s \approx 4.5 \times 10^{-3}$, for $K > 45$, and then $\delta_s \approx 4 \times 10^{-3}$, for $K > 80$.

We recall that our goal is to identify the most important parameters affecting SDs. Our conclusions do not rely on the evaluation of statistical significance. We therefore base our analysis in section 6b on an ensemble size of *K* = 50.

To complete the discussion of scenario Base, we discuss with the help of Fig. 4 the temporal evolution of the SD dispersion and two probability distributions of supersaturation fluctuations (around the evolving mean value), $s'$, for one arbitrarily selected statistical representation. As we start the simulation with a homogeneous SD, the initial dispersion is zero and begins to increase due to the action of TTFs. The associated probability distribution (initially a delta function at $s' = 0$, not shown) decays approximately symmetrically.

For this ensemble member, an entrainment event occurs about halfway through the ascent, replacing a fraction of the parcel air with slightly warmer (by 0.076 K) air from the environment. This inhomogeneity causes the SD dispersion to increase abruptly. The perturbation is subsequently broken down by turbulent mixing, which gradually decreases the

<figure><img src="image_5.png" /><figcaption>FIG. 4. Dispersion of the SD vs time (left) for one ensemble member of scenario Base. The arrow marks an entrainment event, and the colored circles mark the times when probability distributions of supersaturation fluctuations are taken as shown in the right panel.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 5. Ensemble-mean vertical profiles of ice supersaturation vs scaled altitude grid in air parcels at the final lifting distance for all sensitivity scenarios listed in Table 1. For ease of comparison, the brown curve in each panel repeats the Base distribution from Fig. 3. The dashed lines in scenario Blob-0 mark the SD without turbulence (black) and the average value for scenario Base (brown).</figcaption></figure>

dispersion over time. This results in a final probability distribution that is broad, yet still peaked at negative fluctuations of about −0.003, reminiscent of the warming event. The entrainment mixing causes supersaturation fluctuations of either sign about an order of magnitude larger than those induced by TTFs prior to entrainment.

b. *Sensitivity studies*

The domain-averaged (mean) final SDs for all scenarios are shown in Fig. 5. Similar to Base, most sensitivity scenarios have mean *s* values around AVG = 0.501 within a typical standard deviation of about 2 × 10⁻³ (Table 1), except scenarios Blob and Env. We recall from section 2 that the mean *s* value at *t*max without turbulence is 0.504 and that we reduce the fractional blob size in entrainment events by the number of inserted blobs (section 3b).

In the Wind scenarios, the simulations terminate earlier or later than in Base without affecting turbulence parameters (Table 1). This basically means more (Wind-s) or less (Wind-f) time available for turbulence to act on the temperature field. The mean SD is more tilted in the case of the slow updraft, because turbulence had more time to induce a dry lapse rate; hence, ∂_z s > 0 across a significant portion of the domain, ζ = 0.3–0.7, corresponding to a 6-m-thick layer. No notable temperature inversion occurs in the fast updraft scenario. The standard deviations of the profiles are thus considerably different, with SDEV = 4.6 × 10⁻³ (Wind-s) and SDEV = 1.2 × 10⁻³ (Wind-f), i.e., about twice and half the value in scenario Base, respectively.

The Turb scenarios show similarities to the Wind scenarios. In Turb-l, the SD remains oriented vertically with weak perturbations due to the strongly reduced stirring event number relative to Base. The distribution spread, SDEV = 1.1 × 10⁻³, is
the smallest across all sensitivity scenarios and comparable to that of Wind-f. In Turb-h, stirring occurs much more frequently. As a result, a dry lapse rate establishes basically across the full domain, i.e., over a region larger than that in scenario Wind-s. This results in a large spread, $SDEV = 6.3 \times 10^{-3}$, in the final SD.

In scenario Wind-s, $\epsilon$ is too low to overcome the dampening effect of the boundary conditions (section 3d), which manifests in the inverted S-shape of the vertical SD. Note that the Wind-s simulation stops shortly after one large-eddy mixing time, $t_{L-} = 4.7$ min, preventing full vertical mixing. By contrast, the distribution in Turb-h is much closer to the dry-adiabatic solution due to stronger and faster vertical mixing (cp. Fig. 3) with $t_{L+} = 2.2$ min.

Decreasing atmospheric stability increases the number of stirring events (eddy diffusivity) and decreases the temperature perturbation in entrainment events (via $\gamma$). In scenario Stab-l, where the stability is weaker than in Base (Table 1), the SD thus shows a larger gradient across the LEM domain with an increased spread, $SDEV = 3.2 \times 10^{-3}$. In scenario Stab-h, with increased stability, mixing and diffusion are slower and temperature perturbations are larger. Accordingly, the final SD shows more microscale variability relative to Stab-l.

In the Blob scenarios, standard deviations, $SDEV \approx 2 \times 10^{-3}$, are similar to Base (with one inserted blob). The main difference in the case without entrainment (Blob-0) relative to Base is the slightly increased average, AVG = 0.504, that coincides with the unperturbed mean value, $s_{max}$, basically shifting the Blob-0 SD slightly toward higher $s$. This is expected, as entrainment increases the domain-average temperature and thus lowers AVG. The Blob-3 SD with an increased total mixing fraction has an average that is coincidentally similar to scenario Base.

The Env scenarios result in the largest changes in the spatial distributions across all scenarios. Dispersion is large in Env-d and Env-m, since both temperature and water vapor mixing ratio is perturbed in the single entrainment event. The water vapor perturbation amounts to $\pm 3.59$ ppm by mass, or about 3.5% relative to Base. A slightly drier (moister) environment in Env-d (Env-m) shifts the SD to the left (right) of the Base profile. The SDs have average values AVG = 0.490 (0.511) that substantially differ from each other and from the Base distribution.

We finally discuss scenarios in which we vary the inner and outer length scales that bound the inertial subrange of turbulence in the LEM. Scenario Inner increases $L_-$ more than tenfold to 0.1 m (Inner-0.1) and more than 100-fold to 1 m (Inner-1) to assess the effect of degrading spatial resolution. Scenario Outer changes $L_+$ from its default value, 15 m, to a smaller (5 m, Outer-5) and a larger value (25 m, Outer-25).

Setting $L_- > \eta$ decreases the stirring frequency and increases the mean size of eddies in stirring events and thus the magnitudes of TTFs (section 3d). The final SDs are not sensitive to this parameter: all SDs are basically indistinguishable from Base. This may be expected, since the molecular diffusivity is enhanced to include the turbulent mixing contribution from the missing eddies in the size range [$\eta, L_-$] (section 4).

For high Reynolds number flow, the LEM results are expected to depend only weakly on $L_+$ (section 3d). Indeed, we find that the SDs from scenarios Outer deviate only moderately from Base. The occasional sampling of larger eddies with longer breakdown time scales causes more fine structure in the SD from scenario Outer-25 relative to Outer-5, which basically acts to smoothen the Base SD.

## 7. Summary and conclusions

Traditional air parcel models have no internal variability and represent air and cloud properties by assuming instantaneous mixing. In this way, no information on the internal microstructure of clouds can be obtained. Dynamical forcing of parcel motions take into account large-scale motions, at best down to mesoscale motions such as GWs, but usually disregard turbulence. Cloud-resolving simulations of cirrus have not yet explicitly treated ice microphysics on the microscale. In the present study, we coupled a parcel model to a well-established turbulent mixing model to overcome this situation.

In a first step toward a better understanding of how small-scale mixing and molecular diffusion of moisture and heat affect UT ice supersaturation, we simulated the temporal evolution of vertical distributions of ice supersaturation due to turbulence in rising air parcels over a vertical domain of length 15 m, our best estimate for the thickness of homogeneous freezing layers. Initial profiles of moisture and temperature, hence ice supersaturation, were uniform. The parcels were lifted with a constant updraft speed that represents a (sub)mesoscale forcing. Turbulence parameters were specified based on upper-tropospheric observations for stable stratification. We chose initial conditions close to the homogeneous freezing limit (high ice supersaturation) in view of future applications to cirrus formation, but our results apply to other magnitudes of supersaturation as well.

We inspected the ability of the model to predict the decay of supersaturation variance after an initial perturbation. Furthermore, we examined and characterized in more detail a baseline scenario and defined a number of sensitivity studies in which we varied one parameter at a time. We analyzed ensemble-average SDs (vertical profiles) taken after a lifting distance of about 36 m based on 50 statistical realizations in each scenario.

Scenario Base, based on typical UT turbulence intensity, stability, and GW-induced updraft speed, pointed to a large variability in individual SDs from the ensemble, with substantial yet rare supersaturation excursions of $\pm 0.01$ or more away from the mean value (0.5). Stochastic turbulent vertical displacements of air parcels cause turbulent temperature fluctuations. TTFs alone led to a standard deviation of the mean final supersaturation of about 0.002, which is enhanced after entrainment events with considerable sensitivity to the humidity difference between parcel and environmental air. Large case-to-case variability in SDs was also found in all other scenarios. Entrainment-mixing may cause supersaturation fluctuations much larger than those induced by TTFs alone.

Moisture and temperature fluctuations led to small-scale variability in supersaturation, small variations in the ensemble-mean supersaturation within $\pm 0.01$ around the mean value,
and standard deviations in the range 0.001–0.006. Entrainment of warmer air into the parcel decreased the mean supersaturation. In slow updrafts, turbulence acts longer on SDs and makes meter-scale layers more neutrally stratified.

In the presence of enhanced stratification, turbulence produced SDs with more spatial variability. SD standard deviations substantially increased in scenarios with high turbulence intensity and entrainment of drier or moister air. The former scenario generates a dry-adiabatic lapse rate over the full computational domain, resulting in a significant vertical gradient of supersaturation of about $10^{-3}$ m$^{-1}$. The two latter scenarios changed the mean supersaturation relative to Base most strongly. These findings call for a better characterization of the number of, and the properties of ambient air in, entrainment processes in ice-supersaturated UT air parcels, including variability in the moisture field on scales of order $L_o$.

Minimum and maximum sizes of turbulent eddies represented in our model may not be identical to the physical inner and outer turbulence length scales. We generally found a low sensitivity of average SD properties upon variation in the model length scales.

Ice supersaturation development in nature is more complex than illustrated here. In contrast to observations, further high-resolution modeling may provide information for the entire range of scales of interest. For example, future work may employ stochastic entrainment modeling in the context of a parcel model (Abade et al. 2018). Alternatively, our results may be augmented by applying the one-dimensional turbulence model (Kerstein and Wunsch 2006), combining the LEM with a dynamical method for choosing eddy sizes and locations. Furthermore, DNS may help judge the simplified treatment of entrainment in our model and pin down the influence of developing turbulence on ice supersaturation distributions (Chandrakar et al. 2022).

Our study is a starting point to better understand turbulence–ice nucleation interactions during cirrus formation. Homogeneous freezing of aqueous solution droplets and heterogeneous activation of ice-nucleating particles (INPs) depend on temperature and ice supersaturation. A typical homogeneous freezing event occurs within a small range of ice supersaturation, within about 0.03 depending on the updraft regime. Many INPs activate across a larger range of supersaturation, within about 0.1.

Given the level of variations in supersaturation and inhomogeneous mixing features discussed here, it is conceivable that cirrus formation is affected by microscale turbulence. To which degree nucleated ice crystal numbers change is hard to estimate without detailed numerical simulations in view of the nonlinear dependence of the ice nucleation process on temperature and supersaturation. In subsequent work, we will extend our model to include Lagrangian aerosol and ice microphysics and study the effect of turbulence on cirrus formation.

*Acknowledgments.* Fruitful discussions with Joan Alexander, Martina Bramberger, Andreas Dörnbrack, and Joe Werne are gratefully acknowledged. BK thanks Suresh Menon for helpful

exchange on implementing the linear eddy model. FH acknowledges the support from the Emmy Noether Program of the German Research Foundation (DFG) under Grant HO 6588/1-1. This work has also received support from École Polytechnique, from the ANR project BOOST3R (ANR-17-CE01-0016-01) and the French-American project Stratéole 2 (CNES). BG acknowledges funding received from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement 101025473.

*Data availability statement.* Datasets generated by our numerical simulations are publicly accessible (Kärcher 2024).

## APPENDIX A

### TKE Dissipation Rate: Atmospheric Measurements

In cloud-free UT regions at northern midlatitudes, $\epsilon$ can be as small as $10^{-8}-2 \times 10^{-7}$ m$^2$ s$^{-3}$ (Schumann et al. 1995; Gultepe et al. 1995; Mühlbauer et al. 2014). Another set of aircraft measurements in the UT over the southern tip of South America found (3–9) $\times 10^{-5}$ m$^2$ s$^{-3}$ (Dörnbrack et al. 2022). High vertical-resolution (5 m) radiosonde data revealed a long-term (12 years) mean value of about $1.8 \times 10^{-4}$ m$^2$ s$^{-3}$ in the UT over the United States (Ko et al. 2019). Median values of $\epsilon$ derived from multiyear (2010–13) radar measurements between 5- and 15-km altitude over Northern Norway were about (1.6–5.2) $\times 10^{-4}$ m$^2$ s$^{-3}$ (Li et al. 2016). Breaking of a mountain wave near the tropopause resulted in about $10^{-3}$ m$^2$ s$^{-3}$ (Whiteway et al. 2003). Aircraft data suggested a wide range of turbulence conditions also in the TTL, $\epsilon = 10^{-8}-10^{-2}$ m$^2$ s$^{-3}$ (Podglajen et al. 2017; Atlas and Bretherton 2023). Lindborg (1999) reports a typical range $10^{-6}-10^{-5}$ m$^2$ s$^{-3}$ for tropopause conditions. Turbulence intensities in tropical anvil cirrus in response to radiative heating reached values up to about $\epsilon = 5 \times 10^{-3}$ m$^2$ s$^{-3}$ in cloud-resolving simulations (Hartmann et al. 2018), possibly contributing to the longevity of these clouds by promoting ice nucleation.

In sum, very low levels of turbulence in the cloud-free UT are characterized by values up to $\epsilon \approx 10^{-6}$ m$^2$ s$^{-3}$. Such values are probably characteristic for the UT over oceans or flat land surfaces away from sources of mesoscale variability and for regions of strong stratification and weak wind shear. The highest values exceed $10^{-3}$ m$^2$ s$^{-3}$ and represent very rare episodes of intense active turbulence. While it is arguably difficult to characterize typical UT conditions in terms of a single value for $\epsilon$ in cloud-free conditions, a range of values, $10^{-6}-10^{-4}$ m$^2$ s$^{-3}$, are supported by observations. In-cloud values can be significantly larger.

## APPENDIX B

### Scalar Transport and Length-Scale Reduction

The Péclet number Pe, as a measure of scalar transport, is most generally defined by the ratio of advective to diffusive scalar transport rates, $Pe = uL/D$. Here, $u$ is a flow velocity, $L$ is a turbulence length scale, and $D$ is a mass
diffusion coefficient. We introduce inertial range scaling, i.e., $\epsilon = (\text{TKE})^{3/2}/L$ and $u = \sqrt{\text{TKE}}$.

In our case, $u = w'$, $L = L_+$, and $D = D_t$. Using $t_t = L_+^2/D_t$, we get $Pe = t_t w'/L_+$. Applying the above scalings shows that the advective time scale, $L_+/w'$, is identical to Eq. (3) for $\ell = L_+$; hence, $Pe = t_t w'/L_+$. Alternatively, for stratified turbulence, we use $D_t = \epsilon/N^2$ and $L_+ = \sqrt{\text{TKE}}/N$ to evaluate the “advective diffusivity,” $w'/L_+$. This leads to $Pe = D_t/D_m$. Note that in these order of magnitude scalings, $D_t$ differs from Eq. (4) by a factor of 3.

As $t_t$ represents the time needed for a size-$\ell$ eddy to traverse the full inertial range, an approximate equation for the eddy size evolution is given by $d\ell/dt = -\ell/t_t \propto (\epsilon\ell)^{1/3}$ (Jensen and

Baker 1989). With the initial condition $\ell(t=0) = L_+$ and the scaled time $\tau = t/t_{L_+}$, the solution is $\ell(\tau) = L_+[1-(2\tau/3)]^{3/2}$. The rate $d\ell/dt$ is viewed as the average rate of scale reduction in the LEM, and the solution shows that this reduction occurs over roughly one large-eddy mixing time $t_{L_+}$.

In sum, for $Pe \gg 1$, $t_{L_+}$ is the shortest transport time scale and scalar features are broken down by turbulence after $\tau \gg 1$.

## APPENDIX C

### List of Symbols and Acronyms

Table C1 lists the symbols and acronyms used in this study.

<table><caption>TABLE C1. Symbols and acronyms used in this study.</caption><thead><tr><th>Symbol</th><th>Description</th><th>Units</th></tr></thead><tbody><tr><td>a</td><td>Dry-adiabatic coefficient</td><td>m<sup>-1</sup></td></tr><tr><td>g</td><td>Acceleration due to gravity</td><td>m s<sup>-2</sup></td></tr><tr><td>h</td><td>Lagrangian air parcel altitude</td><td>m</td></tr><tr><td>f</td><td>Probability distribution of eddy sizes</td><td>m<sup>-1</sup></td></tr><tr><td>ℓ</td><td>Size of turbulent eddy</td><td>m</td></tr><tr><td>n<sub>b</sub></td><td>Number of blobs inserted during entrainment</td><td>—</td></tr><tr><td>p</td><td>Air pressure</td><td>Pa</td></tr><tr><td>p<sub>s</sub></td><td>Ice saturation vapor pressure</td><td>Pa</td></tr><tr><td>p<sub>v</sub></td><td>Water vapor partial pressure</td><td>Pa</td></tr><tr><td>q</td><td>Parcel water vapor mixing ratio</td><td>—</td></tr><tr><td>q<sub>e</sub></td><td>Environmental water vapor mixing ratio</td><td>—</td></tr><tr><td>s</td><td>Ice supersaturation</td><td>—</td></tr><tr><td>s<sub>m</sub></td><td>Gridpoint ensemble-mean s average</td><td>—</td></tr><tr><td>t</td><td>Lagrangian air parcel time</td><td>s</td></tr><tr><td>t<sub>K</sub></td><td>Kolmogorov time scale</td><td>s</td></tr><tr><td>t<sub>ℓ</sub></td><td>Turbulent mixing time for size-ℓ eddy</td><td>s</td></tr><tr><td>t<sub>max</sub></td><td>Maximum simulation time</td><td>s</td></tr><tr><td>t<sub>m</sub></td><td>Time scale for molecular diffusion</td><td>s</td></tr><tr><td>t<sub>t</sub></td><td>Time scale for turbulent diffusion</td><td>s</td></tr><tr><td>t<sub>L<sub>+</sub></sub></td><td>Large eddy mixing time scale</td><td>s</td></tr><tr><td>w</td><td>Updraft speed of air parcel</td><td>m s<sup>-1</sup></td></tr><tr><td>z</td><td>LEM grid altitude in the Lagrangian air parcel</td><td>m</td></tr><tr><td>D<sub>t</sub></td><td>Turbulent diffusivity</td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>D<sub>m</sub></td><td>Molecular diffusivity</td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>F</td><td>Cumulative distribution of eddy sizes</td><td>—</td></tr><tr><td>H</td><td>Maximum parcel altitude w<sub>t,max</sub></td><td>m</td></tr><tr><td>J</td><td>Number of LEM grid cells</td><td>—</td></tr><tr><td>K</td><td>Ensemble size</td><td>—</td></tr><tr><td>L<sub>-</sub> (L<sub>+</sub>)</td><td>Inner (outer) LEM length scale</td><td>m</td></tr><tr><td>L<sub>f</sub></td><td>Thickness of a homogeneous freezing layer</td><td>m</td></tr><tr><td>L<sub>o</sub></td><td>Outer length scale in turbulence</td><td>m</td></tr><tr><td>L<sub>O</sub></td><td>Ozmidov length scale</td><td>m</td></tr><tr><td>N</td><td>Brunt-Väisälä frequency</td><td>s<sup>-1</sup></td></tr><tr><td>N<sub>s</sub></td><td>Total number of stirring events v<sub>s</sub>l<sub>max</sub></td><td>—</td></tr><tr><td>Pe</td><td>Péclet number for mass transport</td><td>—</td></tr><tr><td>R</td><td>Random number [0, 1] from uniform distribution</td><td>—</td></tr><tr><td>R<sub>g</sub></td><td>Specific gas constant for dry air</td><td>J (kg K)<sup>-1</sup></td></tr><tr><td>Re</td><td>Buoyancy Reynolds number</td><td>—</td></tr><tr><td>Re<sub>LEM</sub></td><td>Model Reynolds number</td><td>—</td></tr><tr><td>S</td><td>Ice saturation ratio s + 1</td><td>—</td></tr><tr><td>Sc</td><td>Molecular Schmidt number</td><td>—</td></tr><tr><td>S<sub>e</sub></td><td>Environmental ice saturation ratio</td><td>—</td></tr><tr><td>S<sub>hom</sub></td><td>Homogeneous freezing threshold</td><td>—</td></tr></tbody></table>
<table><caption>TABLE C1. (Continued)</caption><thead><tr><th>Symbol</th><th>Description</th><th>Units</th></tr></thead><tbody><tr><td>S<sub>max</sub></td><td>Maximum ice saturation ratio</td><td>—</td></tr><tr><td>T</td><td>Air parcel temperature</td><td>K</td></tr><tr><td>T<sub>e</sub></td><td>Environmental air temperature</td><td>K</td></tr><tr><td>β</td><td>Size of entrained blob relative to L<sub>+</sub></td><td>—</td></tr><tr><td>δ<sub>s</sub></td><td>Domain-average dispersion of SD</td><td>—</td></tr><tr><td>δt</td><td>Model time step</td><td>s</td></tr><tr><td>Δt</td><td>Waiting time since the last stirring event</td><td>s</td></tr><tr><td>ΔT</td><td>Magnitude of a TTF</td><td>K</td></tr><tr><td>δz</td><td>Size of one grid cell (LEM resolution)</td><td>m</td></tr><tr><td>ε</td><td>Dissipation rate of TKE</td><td>m<sup>2</sup> s<sup>-3</sup></td></tr><tr><td>η</td><td>Kolmogorov microscale</td><td>m</td></tr><tr><td>γ</td><td>Environmental lapse rate -dT<sub>e</sub>/dz</td><td>K m<sup>-1</sup></td></tr><tr><td>Γ</td><td>Dry-adiabatic (parcel) lapse rate -dT/dz</td><td>K m<sup>-1</sup></td></tr><tr><td>λ</td><td>Model length scale ratio L<sub>+</sub>/L<sub>-</sub></td><td>—</td></tr><tr><td>Λ</td><td>Stirring frequency per unit eddy length</td><td>(m s)<sup>-1</sup></td></tr><tr><td>ν</td><td>Kinematic viscosity of air</td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>ν<sub>s</sub></td><td>Average stirring frequency ΛL<sub>+</sub></td><td>s<sup>-1</sup></td></tr><tr><td>τ</td><td>Air parcel time scaled by large-eddy mixing time</td><td>—</td></tr><tr><td>ζ</td><td>Scaled altitude z/L<sub>+</sub></td><td>—</td></tr><tr><td>AVG</td><td>Mean value of SD (domain-average)</td><td>—</td></tr><tr><td>DNS</td><td>Direct numerical simulations</td><td>—</td></tr><tr><td>EPM</td><td>Entraining parcel model</td><td>—</td></tr><tr><td>GW</td><td>Gravity wave</td><td>—</td></tr><tr><td>H<sub>2</sub>O</td><td>Water vapor</td><td>—</td></tr><tr><td>LEM</td><td>Linear eddy model</td><td>—</td></tr><tr><td>TKE</td><td>Turbulence kinetic energy per unit mass of air</td><td>m<sup>2</sup> s<sup>-2</sup></td></tr><tr><td>TTF</td><td>Turbulent temperature fluctuation</td><td>—</td></tr><tr><td>TTL</td><td>Tropical tropopause layer</td><td>—</td></tr><tr><td>UT</td><td>Upper troposphere</td><td>—</td></tr><tr><td>SD</td><td>Supersaturation distribution (vertical profile)</td><td>—</td></tr><tr><td>SDEV</td><td>Standard deviation of SD</td><td>—</td></tr><tr><td>VAR</td><td>Variance of SD</td><td>—</td></tr></tbody></table>

REFERENCES

Abade, G. C., W. W. Grabowski, and H. Pawlowska, 2018: Broadening of cloud droplet spectra through eddy hopping: Turbulent entraining parcel simulations. *J. Atmos. Sci.*, **75**, 3365–3379, https://doi.org/10.1175/JAS-D-18-0078.1.

Atlas, R., and C. S. Bretherton, 2023: Aircraft observations of gravity wave activity and turbulence in the tropical tropopause layer: Prevalence, influence on cirrus clouds, and comparison with global storm-resolving models. *Atmos. Chem. Phys.*, **23**, 4009–4030, https://doi.org/10.5194/acp-23-4009-2023.

Atlas, R. L., C. S. Bretherton, A. B. Sokol, P. N. Blossey, and M. F. Khairoutdinov, 2024: Tropical cirrus are highly sensitive to ice microphysics within a nudged global storm-resolving model. *Geophys. Res. Lett.*, **51**, e2023GL105868, https://doi.org/10.1029/2023GL105868.

Chandrakar, K. K., W. Cantrell, S. Krueger, R. A. Shaw, and S. Wunsch, 2020: Supersaturation fluctuations in moist turbulent Rayleigh-Bénard convection: A two-scalar transport problem. *J. Fluid Mech.*, **884**, A19, https://doi.org/10.1017/jfm. 2019.895.

—, H. Morrison, W. W. Grabowski, G. H. Bryan, and R. A. Shaw, 2022: Supersaturation variability from scalar mixing:

Evaluation of a new subgrid-scale model using direct numerical simulations of turbulent Rayleigh-Bénard convection. *J. Atmos. Sci.*, **79**, 1191–1210, https://doi.org/10.1175/JAS-D-21-0250.1.

Corcos, M., A. Hertzog, R. Plougonven, and A. Podglajen, 2021: Observation of gravity waves at the tropical tropopause using superpressure balloons. *J. Geophys. Res. Atmos.*, **126**, e2021JD035165, https://doi.org/10.1029/2021JD035165.

Crank, J., and P. Nicolson, 1947: A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. *Math. Proc. Cambridge Philos. Soc.*, **43**, 50–67, https://doi.org/10.1017/S0305004100023197.

Dinh, T. P., D. R. Durran, and T. P. Ackerman, 2010: Maintenance of tropical tropopause layer cirrus. *J. Geophys. Res.*, **115**, D02104, https://doi.org/10.1029/2009jd012735.

Dobbie, S., and P. Jonas, 2001: Radiative influences on the structure and lifetime of cirrus clouds. *Quart. J. Roy. Meteor. Soc.*, **127**, 2663–2682, https://doi.org/10.1002/qj.49712757808.

Dörnbrack, A., P. Bechtold, and U. Schumann, 2022: High-resolution aircraft observations of turbulence and waves in the free atmosphere and comparison with global model predictions. *J. Geophys. Res. Atmos.*, **127**, e2022JD036654, https://doi.org/10.1029/2022JD036654.
Fritts, D. C., and M. J. Alexander, 2003: Gravity wave dynamics and effects in the middle atmosphere. *Rev. Geophys.*, **41**, 1003, https://doi.org/10.1029/2001RG000106.

Gasparini, B., P. N. Blossey, D. L. Hartmann, G. Lin, and J. Fan, 2019: What drives the life cycle of tropical anvil clouds?. *J. Adv. Model. Earth Syst.*, **11**, 2586–2605, https://doi.org/10.1029/2019MS001736.

Gettelman, A., W. D. Collins, E. J. Fetzer, A. Eldering, F. W. Irion, P. B. Duffy, and G. Bala, 2006: Climatology of upper-tropospheric relative humidity from the atmospheric infrared sounder and implications for climate. *J. Climate*, **19**, 6104–6121, https://doi.org/10.1175/JCLI3956.1.

Gu, Y., and K. N. Liou, 2000: Interactions of radiation, microphysics, and turbulence in the evolution of cirrus clouds. *J. Atmos. Sci.*, **57**, 2463–2479, https://doi.org/10.1175/1520-0469(2000)057%3C2463:IORMAT%3E2.0.CO;2.

Gultepe, I., D. O'G. Starr, A. J. Heymsfield, T. Uttal, T. P. Ackerman, and D. L. WestPhal, 1995: Dynamical characteristics of cirrus clouds from aircraft and Radar observations in micro and meso-γ scales. *J. Atmos. Sci.*, **52**, 4060–4078, https://doi.org/10.1175/1520-0469(1995)052%3C4060:DCOCFF%3E2.0.CO;2.

Hartmann, D. L., B. Gasparini, S. E. Berry, and P. N. Blossey, 2018: The life cycle and net radiative effect of tropical anvil clouds. *J. Adv. Model. Earth Syst.*, **10**, 3012–3029, https://doi.org/10.1029/2018MS001484.

Hoffmann, F., 2023: The small-scale mixing of clouds with their environment: Impacts on micro- and macroscale cloud properties. *Clouds and Their Climatic Impacts: Radiation, Circulation, and Precipitation, Geophys. Monogr.*, Vol. 281, Amer. Geophys. Union, 255–270, https://doi.org/10.1002/9781119700357.ch12.

—, and G. Feingold, 2019: Entrainment and mixing in stratocumulus: Effects of a new explicit subgrid-scale scheme for large-eddy simulations with particle-based microphysics. *J. Atmos. Sci.*, **76**, 1955–1973, https://doi.org/10.1175/JAS-D-180318.1.

Jensen, E. J., L. Pfister, and T. P. Bui, 2012: Physical processes controlling ice concentrations in cold cirrus near the tropical tropopause. *J. Geophys. Res.*, **117**, D11205, https://doi.org/10.1029/2011JD017319.

—, and Coauthors, 2013: Ice nucleation and dehydration in the tropical tropopause layer. *Proc. Natl. Acad. Sci. USA*, **110**, 2041–2046, https://doi.org/10.1073/pnas.1217104110.

—, G. S. Diskin, J. DiGangi, S. Woods, R. P. Lawson, and T. V. Bui, 2022: Homogeneous freezing events sampled in the tropical tropopause layer. *J. Geophys. Res. Atmos.*, **127**, e2022JD036535, https://doi.org/10.1029/2022JD036535.

Jensen, J. B., and M. B. Baker, 1989: A simple model of droplet spectral evolution during turbulent mixing. *J. Atmos. Sci.*, **46**, 2812–2829, https://doi.org/10.1175/1520-0469(1989)046%3C2812:ASMODS%3E2.0.CO;2.

Kärcher, B., 2017: Homogeneous ice formation in convective cloud outflow regions. *Quart. J. Roy. Meteor. Soc.*, **143**, 2093–2103, https://doi.org/10.1002/qj.3069.

—, 2024: Effects of turbulence on upper tropospheric ice supersaturation. Zenodo, accessed 17 August 2024, https://doi.org/10.5281/zenodo.10876514.

—, and U. Lohmann, 2002: A parameterization of cirrus cloud formation: Homogeneous freezing of supercooled aerosols. *J. Geophys. Res.*, **107**, 4010, https://doi.org/10.1029/2001JD000470.

—, and A. Podglajen, 2019: A stochastic representation of temperature fluctuations induced by mesoscale gravity waves.

*J. Geophys. Res. Atmos.*, **124**, 11506–11529, https://doi.org/10.1029/2019JD030680.

—, E. J. Jensen, G. F. Pokrifka, and J. Y. Harrington, 2023: Ice supersaturation variability in cirrus clouds: Role of vertical wind speeds and deposition coefficients. *J. Geophys. Res. Atmos.*, **128**, e2023JD039324, https://doi.org/10.1029/2023JD039324.

Kerstein, A. R., 1988: A Linear-eddy model of turbulent scalar transport and mixing. *Combust. Sci. Technol.*, **60**, 391–421, https://doi.org/10.1080/00102208808923995.

—, 1991: Linear eddy modeling of turbulent transport: Part 6. Microstructure of diffusive scalar mixing fields. *J. Fluid Mech.*, **231**, 361–394, https://doi.org/10.1017/S0022112091003439.

—, and S. Wunsch, 2006: Simulation of a Stably Stratified Atmospheric Boundary Layer using one-dimensional turbulence. *Bound.-Layer Meteor.*, **118**, 325–356, https://doi.org/10.1007/s10546-005-9004-x.

Ko, H.-C., H.-Y. Chun, R. Wilson, and M. A. Geller, 2019: Characteristics of atmospheric turbulence retrieved from high vertical-resolution radiosonde data in the United States. *J. Geophys. Res. Atmos.*, **124**, 7553–7579, https://doi.org/10.1029/2019JD030287.

Korolev, A. V., and I. P. Mazin, 2003: Supersaturation of water vapor in clouds. *J. Atmos. Sci.*, **60**, 2957–2974, https://doi.org/10.1175/1520-0469(2003)060%3C2957:SOWVIC%3E2.0.CO;2.

Krueger, S. K., 1993: Linear eddy modeling of entrainment and mixing in stratus clouds. *J. Atmos. Sci.*, **50**, 3078–3090, https://doi.org/10.1175/1520-0469(1993)050%3C3078:LEMOEA%3E2.0.CO;2.

—, C.-W. Su, and P. A. McMurtry, 1997: Modeling entrainment and finescale mixing in cumulus clouds. *J. Atmos. Sci.*, **54**, 2697–2712, https://doi.org/10.1175/1520-0469(1997)054%3C2697:MEAFMI%3E2.0.CO;2.

Kulmala, M., U. Rannik, E. L. Zapadinsky, and C. F. Clement, 1997: The effect of saturation fluctuations on droplet growth. *J. Aerosol Sci.*, **28**, 1395–1409, https://doi.org/10.1016/S0021-8502(97)00015-3.

Lamquin, N., C. J., Stubenrauch, K. Gierens, U. Burkhardt, and H. Smit, 2012: A global climatology of upper-tropospheric ice supersaturation occurrence inferred from the Atmospheric Infrared Sounder calibrated by MOZAIC. *Atmos. Chem. Phys.*, **12**, 381–405, https://doi.org/10.5194/acp-12-381-2012.

Landau, L., and E. M. Lifshitz, 1987: *Fluid Mechanics, Course of Theoretical Physics*. Vol. 6, Pergamon Press, 539 pp.

Lanotte, A. S., A. Seminara, and F. Toschi, 2009: Cloud droplet growth by condensation in homogeneous isotropic turbulence. *J. Atmos. Sci.*, **66**, 1685–1697, https://doi.org/10.1175/2008JAS2864.1.

Lewellen, D. C., O. Meza, and W. W. Huebsch, 2014: Persistent contrails and contrail cirrus. Part I: Large-eddy simulations from inception to demise. *J. Atmos. Sci.*, **71**, 4399–4419, https://doi.org/10.1175/JAS-D-13-0316.1.

Li, Q., M. Rapp, A. Schrön, A. Schneider, and G. Stober, 2016: Derivation of turbulent energy dissipation rate with the Middle Atmosphere Alomar Radar System (MAARSY) and radiosondes at Andøya, Norway. *Ann. Geophys.*, **34**, 1209–1229, https://doi.org/10.5194/angeo-34-1209-2016.

Lilly, D. K., D. E. Waco, and S. I. Aldefang, 1974: Stratospheric mixing estimated from high-altitude turbulence measurements. *J. Appl. Meteor.*, **13**, 488–493, https://doi.org/10.1175/1520-0450(1974)013%3C0488:SMEFHA%3E2.0.CO;2.
Lim, J.-S., and F. Hoffmann, 2023: Between broadening and narrowing: How mixing affects the width of the droplet size distribution. *J. Geophys. Res. Atmos.*, **128**, e2022JD037900, https://doi.org/10.1029/2022JD037900.

Lin, R.-F., D. O. Starr, P. J. DeMott, R. Cotton, K. Sassen, E. Jensen, B. Kärcher, and X. Liu, 2002: Cirrus Parcel Model Comparison Project. Phase 1: The critical components to simulate cirrus initiation explicitly. *J. Atmos. Sci.*, **59**, 2305–2329, https://doi.org/10.1175/1520-0469(2002)059%3C2305:CPMCPP%3E2.0.CO;2.

Lindborg, E., 1999: Can the atmospheric kinetic energy spectrum be explained by two-dimensional turbulence? *J. Fluid Mech.*, **388**, 259–288, https://doi.org/10.1017/S0022112099004851.

Matsumoto, M., and T. Nishimura, 1998: Mersenne Twister: A 623-dimensionally equidistributed uniform pseudo-random number generator. *ACM Trans. Model. Comput. Simul.*, **8**, 3–30, https://doi.org/10.1145/272991.272995.

Menon, S., and J. Wu, 1998: Effects of micro- and macroscale turbulent mixing on the chemical processes in engine exhaust plumes. *J. Appl. Meteor.*, **37**, 639–654, https://doi.org/10.1175/1520-0450(1998)037%3C0639:EOMAMT%3E2.0.CO;2.

—, and A. R. Kerstein, 2011: The linear-eddy model. *Turbulent Combustion Modeling*, T. Echekki and E. Mastorakos, Eds., Fluid Mechanics and Its Applications, Vol. 95, Springer, 221–247.

Mühlbauer, A., H. Kalesse, and P. Kollias, 2014: Vertical velocities and turbulence in midlatitude anvil cirrus: A comparison between in situ aircraft measurements and ground-based Doppler cloud radar retrievals. *Geophys. Res. Lett.*, **41**, 7814–7821, https://doi.org/10.1002/2014GL062279.

Murphy, D. M., and T. Koop, 2005: Review of the vapour pressures of ice and supercooled water for atmospheric applications. *Quart. J. Roy. Meteor. Soc.*, **131**, 1539–1565, https://doi.org/10.1256/qj.04.94.

Nugent, J. M., S. M. Turbeville, C. S. Bretherton, P. N. Blossey, and T. P. Ackerman, 2022: Tropical cirrus in global storm-resolving models: 1. Role of deep convection. *Earth Space Sci.*, **9**, e2021EA001965, https://doi.org/10.1029/2021EA001965.

Ohno, T., M. Satoh, and A. Noda, 2019: Fine vertical resolution radiative-convective equilibrium experiments: Roles of turbulent mixing on the high-cloud response to sea surface temperatures. *J. Adv. Model. Earth Syst.*, **11**, 1637–1654, https://doi.org/10.1029/2019MS001704.

Paoli, R., and K. Shariff, 2009: Turbulent condensation of droplets: Direct simulation and a stochastic model. *J. Atmos. Sci.*, **66**, 723–740, https://doi.org/10.1175/2008JAS2734.1.

—, and —, 2016: Contrail modeling and simulation. *Annu. Rev. Fluid Mech.*, **48**, 393–427, https://doi.org/10.1146/annurev-fluid-010814-013619.

—, O. Thouron, J. Escobar, J. Picot, and D. Cariolle, 2014: High-resolution large-eddy simulations of stably stratified flows: Application to subkilometer-scale turbulence in the upper troposphere-lower stratosphere. *Atmos. Chem. Phys.*, **14**, 5037–5055, https://doi.org/10.5194/acp-14-5037-2014.

Podglajen, A., and Coauthors, 2017: Small-scale wind fluctuations in the tropical tropopause layer from aircraft measurements: Occurrence, nature and impact on vertical mixing. *J. Atmos. Sci.*, **74**, 3847–3869, https://doi.org/10.1175/JAS-D-17-0010.1.

Politovich, M. K., and W. A. Cooper, 1988: Variability of the super-saturation in cumulus clouds. *J. Atmos. Sci.*, **45**, 1651–1664, https://doi.org/10.1175/1520-0469(1988)045<1651:VOTSIC>2.0.CO;2.

Pruppacher, H. R., and J. D. Klett, 1997: *Microphysics of Clouds and Precipitation*. Kluwer Academic Publishers, 954 pp.

Quante, M., and D. O. Starr, 2002: Dynamic processes in cirrus clouds: A review of observational results. *Cirrus*, D. K. Lynch et al., Eds., Oxford University Press, 375–396.

Riley, J. J., and E. Lindborg, 2008: Stratified turbulence: A possible interpretation of some geophysical turbulence measurements. *J. Atmos. Sci.*, **65**, 2416–2424, https://doi.org/10.1175/2007JAS2455.1.

Schumann, U., P. Konopka, R. Baumann, R. Busen, T. Gerz, H. Schlager, P. Schulte, and H. Volkert, 1995: Estimate of diffusion parameters of aircraft exhaust plumes near the tropopause from nitric oxide and turbulence measurements. *J. Geophys. Res.*, **100**, 14147–14162, https://doi.org/10.1029/95JD01277.

Siebert, H., and R. A. Shaw, 2017: Supersaturation fluctuations during the early stage of cumulus formation. *J. Atmos. Sci.*, **74**, 975–988, https://doi.org/10.1175/JAS-D-16-0115.1.

Spichtinger, P., and K. M. Gierens, 2009: Modelling of cirrus clouds – Part 1a: Model description and validation. *Atmos. Chem. Phys.*, **9**, 685–706, https://doi.org/10.5194/acp-9-685-2009.

Su, C.-W., S. K. Krueger, P. A. McMurtry, and P. H. Austin, 1998: Linear eddy modeling of droplet spectral evolution during entrainment and mixing in cumulus clouds. *Atmos. Res.*, **47–48**, 41–58, https://doi.org/10.1016/S0169-8095(98)00039-8.

Turbeville, S. M., J. M. Nugent, T. P. Ackerman, C. S. Bretherton, and P. N. Blossey, 2022: Tropical cirrus in global storm-resolving models: 2. Cirrus life cycle and top-of-atmosphere radiative fluxes. *Earth Space Sci.*, **9**, e2021EA001978, https://doi.org/10.1029/2021EA001978.

Vicari, R., C. C. Stephan, T. P. Lane, and Y. Huang, 2024: Analysis of trapped small-scale internal gravity waves automatically detected in satellite imagery. *J. Geophys. Res. Atmos.*, **129**, e2023JD038956, https://doi.org/10.1029/2023JD038956.

Waite, M. L., 2011: Stratified turbulence at the buoyancy scale. *Phys. Fluids*, **23**, 066602, https://doi.org/10.1063/1.3599699.

Weinstock, J., 1978: Vertical turbulent diffusion in stably stratified fluid. *J. Atmos. Sci.*, **35**, 1022–1027, https://doi.org/10.1175/1520-0469(1978)035<1022:VTDIAS>2.0.CO;2.

Werne, J., and D. C. Fritts, 1999: Stratified shear turbulence: Evolution and statistics. *Geophys. Res. Lett.*, **26**, 439–442, https://doi.org/10.1029/1999GL000022.

Whiteway, J. A., E. G. Pavelin, R. Busen, J. Hacker, and S. Vosper, 2003: Airborne measurements of gravity wave breaking at the tropopause. *Geophys. Res. Lett.*, **30**, 2070, https://doi.org/10.1029/2003GL018207.

Wilson, R., 2004: Turbulent diffusivity in the free atmosphere inferred from MST radar measurements: A review. *Ann. Geophys.*, **22**, 3869–3887, https://doi.org/10.5194/angeo-22-3869-2004.

Woosley, S. E., A. R. Kerstein, V. Sankaran, A. J. Aspden, and F. K. Röpke, 2009: Type Ia supernovae: Calculations of turbulent flames using the linear eddy model. *Astrophys. J.*, **704**, 255, https://doi.org/10.1088/0004-637X/704/1/255.

Wu, J., and S. Menon, 2001: Aerosol dynamics in the near field of engine exhaust plumes. *J. Appl. Meteor.*, **40**, 795–809, https://doi.org/10.1175/1520-0450(2001)040<0795:ADITNF>2.0.CO;2.