A Stochastic Representation of Temperature Fluctuations Induced by Mesoscale Gravity Waves

Ubiquitous mesoscale gravity waves generate high cooling rates important for cirrus formation. We make use of long‐duration balloon observations to devise a probabilistic model describing mesoscale temperature uctuations (MTF) away from strong wave sources. We define background conditions based on observed probability distributions of temperature and underlying vertical wind speed fluctuations. We show theoretically that MTF are subject to damping at a rate near the Coriolis frequency when the vertical wind speed fluctuations are autocorrelated over a fraction of a Brunt‐Väisälä period. We find that for background wave activity, a decrease in temperature of 1K translates into cooling rate standard deviations and mean updraft speeds of 4–8Kh−1 and ≈ 10–20 cms−1, respectively, depending on latitude and stratification. We introduce an effective Coriolis frequency to generate cooling rates in equatorial regions consistent with balloon data. Above ice saturation, MTF are large enough to affect ice crystal nucleation. Our results help constrain uncertainty in aerosol‐cirrus interactions, provide insights to better meet challenges in comparing measurement data with model simulations, and support the development of cutting‐edge ice cloud schemes in global models.


Introduction
Mesoscale temperature fluctuations (MTF) and the associated cooling and heating rates (CHR) exert a crucial influence on ice crystal number concentrations in cirrus and polar stratospheric clouds by perturbing mean thermodynamic states of Lagrangian air parcels Carslaw et al., 1999;Dinh et al., 2016;Jensen et al., 2016;Murphy, 2003;Spichtinger & Krämer, 2013). The current inability to constrain cloud models and parameterizations regarding MTF limits the fidelity with which changes in cirrus clouds due to anthropogenic activities can be simulated (Kärcher, 2017).
The need to consider effects of MTF on supersaturation and ice nucleation in microphysical cloud models and cloud schemes of global models has clearly been recognized. Considering only synoptic temperature variations leads to dramatic underestimations of cooling rates and therefore to substantial errors in predictions of microphysical properties of ice clouds (Jensen & Pfister, 2004;Jensen & Toon, 1994;Hoyle et al., 2005;Kärcher & Ström, 2003;Kienast-Sjögren et al., 2015;Meilinger et al., 1995;Murphy & Gary, 1995;Shi & Liu, 2016).
Gravity waves (GW) are an important factor influencing meteorological processes in the stably stratified upper troposphere (UT) and lower stratosphere (LS) on the mesoscale (Craig & Selz, 2018). Balloon, aircraft, and radar measurements link vertical air motion variability to GW activity (Bacmeister et al., 1996;Ecklund et al., 1986;Hoyle et al., 2005;Li et al., 2018;Podglajen et al., 2016;Murphy & Gary, 1995;Schoeberl et al., 2017). Measurements show that MTF are widespread globally and occur frequently, even in conditions not directly affected by convective and mountain wave activity or by the jet stream, introducing the notion of background conditions (Gary, 2006). Low frequency GW associated with long wavelengths (>100-1,000 km) control where cirrus occur by producing large ice supersaturated areas (Kim et al., 2016) and influence their evolution through wind and temperature perturbations (Podglajen et al., 2018). High frequency GW associated with large cooling rates strongly influence the number concentration of nucleated ice crystals.
In models without explicit representation of the full mesoscale GW spectrum from ≈1 km to several 100 km of horizontal extension, vertical wind speeds must be parameterized if processes sensitive to cooling rates such as cloud ice crystal formation are to be estimated. In global weather and climate models, wave-induced updraft speeds and equivalent cooling rates are notoriously under-resolved, and physically based parameterizations of their impact on cirrus formation are lacking. The goal of the present study is to devise a model with only few parameters capable of explaining those aspects of observations that are arguably most crucial for studies of ice crystal nucleation in cirrus. The improved understanding and representation of Lagrangian cooling rates-those measured along air parcel trajectories-will ultimately allow cloud modelers to reduce uncertainty in cirrus simulations and advance parameterizations of cirrus formation.
Previous approaches representing MTF include expansion of wave perturbations into harmonic functions Hoyle et al., 2005;Jensen & Pfister, 2004), multifractal random noise with shortterm autocorrelation (Murphy, 2014), an autoregressive model constrained by long-duration, quasi-Lagrangian balloon observations , regional weather forecast model data at high temporal and spatial resolution (Kienast-Sjögren et al., 2015), and global model simulations (Barahona et al., 2017;Penner et al., 2018). The wide spectrum of GW and the multitude of tropospheric GW sources motivates the development of a probabilistic approach-applicable in background conditions in the global UTLS-which allies basic stochastic methods with mesoscale GW.
Section 2 recalls main features of GW and MTF in the UTLS and reviews observational data sets. Section 3 offers an operational definition of background conditions and presents an analytical model connecting MTF and CHR amplitudes in such conditions. In section 4, this model is compared to and interpreted with a probabilistic approach to predict MTF and its numerical solution. Section 5 suggests how updraft speeds used as drivers for ice formation parameterizations in global models might be improved and explores how MTF affect ice supersaturation fluctuations. Section 6 summarizes results of this study and proposes topics for further investigation.

Mesoscale Temperature Fluctuations
We recall key properties of mesoscale GW and emphasize the significant difference between quantifying MTF and CHR variations in Lagrangian and Eulerian frames of reference. Basic information is summarized here and later combined with model development. We assort experimental results to indicate the range of GW-induced fluctuation amplitudes measured with a number of observation techniques in various locales, seasons, and meteorological situations. Finally, we comment on issues associated with datamodel intercomparisons.

Gravity Waves
Any process that displaces air parcels can trigger atmospheric GW. GW sources are multifold and occur across disparate temporal and spatial scales (Nappo, 2013). Sources of mesoscale GW include flow over land and elevated topographic barriers; convection and storm systems; thermal generation by latent heat release in clouds or evaporation of rain; and strong wind shears and geostrophic adjustment near jet streams or frontal systems. GW are capable of carrying momentum and energy into areas far away from their source regions. They travel horizontally and vertically affected by variable wind and density profiles; in their fields, air parcels are displaced along slant paths. Breaking GW generate turbulence and mixing. They break due to convective or shear instability or when they approach critical levels, that is, locations where the wave phase speed matches the background wind speed. Individual GW tend to form localized wave packets containing a spectrum of frequencies.
We study effects of mesoscale GW influenced by Earth's rotation in an atmosphere that is not only statically but also inertially stable. Intrinsic wave frequencies, ω-those measured relative to the mean air motionare constrained by a dispersion relation. In its simplest form for an inviscid atmosphere, it is given by (Holton, 2004): where N is the Brunt-Väisälä frequency, f is the Coriolis frequency, and m is the magnitude of the vertical and k,l are the magnitudes of the horizontal wavenumbers. This relationship holds for sufficiently small vertical wavelengths, Λ (m=2π/Λ≫1/(2H)), that is, when non-hydrostatic effects involving the scale height of the unperturbed air mass density field ρ ð Þ; H ¼ −½∂lnðρÞ=∂z −1 ≈7 km, can be neglected; z is the altitude. The frequency range of GW with real wave numbers is constrained to lie within |f| and N. Ignoring the effect of rotation leads to the dispersion relation for pure buoyancy waves. Waves with frequencies exceeding N are evanescent and can no longer propagate. Waves with frequencies close to f dissipate due to dynamical instability or radiative damping.
Considering GW activity, the mesoscale region is in terms of characteristic time scales constrained by N and f. The Brunt-Väisälä frequency is given by N ¼ ½g ∂lnðθÞ=∂z 1=2 , where g is the acceleration due to gravity and θ is the potential temperature. The buoyancy oscillation period is 2π/N≃9 min for average tropospheric conditions (N=0.012 rad/s) and about half as long in the LS. The Coriolis frequency is defined by f ¼ 2ΩsinðϕÞ, where Ω is the Earth's rotation rate and ϕ is the geographical latitude. Oscillation periods are 2π/f≃16 hr at midlatitudes (ϕ=50°, f≃10 −4 rad/s), 12 hr at the poles, and >1 day in the tropics. Vertical wavelengths are typically few kilometers. Horizontal wavelengths of GW range from few kilometers to thousands of kilometers. At times, we use the term small scale to denote the sub-kilometer domain encompassing the transition region to inertial turbulence (Fritts et al., 2009;Paoli et al., 2014).
Intrinsic frequencies and vertical wave numbers may vary strongly as GW propagate vertically through background wind shear and temperature gradients (Fritts & Alexander, 2003). Wave group velocities describe energy transport. Low frequency GW (ω≪N) have slower vertical than horizontal group velocities than high frequency GW (ω≫f). Purely vertical displacements occur only when GW travel horizontally, and air parcels affected by them oscillate at their highest possible frequency. Purely horizontal displacements with circular trajectories occur only when GW travel vertically, in which case air parcels undergo lowfrequency oscillations. At intermediate frequencies, wave motions in the plane perpendicular to the wave vector are elliptical. Since horizontal temperature gradients are weaker than vertical variations, it is the vertical displacement component due to GW that gives rise to the rapid CHR in air parcels we are interested in.
The above implies that low frequency GW are capable of affecting large horizontal areas, creating much of the CHR variability that we associate with background conditions in section 3.1. By contrast, high frequency GW, predominantly generated by deep convection and mountains, influence mainly the regions above which they are generated and are capable of triggering rather localized and short-term large-amplitude perturbations (Podglajen, 2017).
Tropical GW behave distinctly differently from those in the extratropics, as their properties become increasingly affected by the dependence of the Coriolis frequency on latitude. While f is assumed to be constant in deriving equation (1), latitudinal variations in f cannot be neglected when analyzing properties of GW in the tropics. The low-frequency range is the domain of inertia GW, quasi-geostrophic Rossby waves, Kelvin waves, and mixed Rossby-GW. Free equatorial GW are shallow waves that move predominantly zonally with restricted north-south extension. Kelvin waves have relatively high intrinsic frequencies when their wavelengths are sufficiently long. Past studies established a tight connection between these wave types and deep convection (Kiladis et al., 2009).
The magnitude of temperature variations, Oð1Þ K, their spectral properties, and their ubiquity in the UTLS strongly suggest a key role of mesoscale GW as the primary cause of MTF. Vertical scales associated with small fluctuation amplitudes are Oð100Þ m, using the dry adiabatic lapse rate Γ=g/c p ≃10 K/km, where c p is the specific heat of dry air at constant pressure. As shown in this study, this translates into cooling rates >1 K/hr, exceeding the magnitudes of radiative and synoptic temperature tendencies at least tenfold.

Fluctuations in Lagrangian Versus Eulerian Reference Frames
The distinction between fluctuation amplitudes evaluated in the two reference frames is important for interpreting observations and relating them to modeled MTF and associated CHR. To quantify the different magnitudes, we assume that only vertical stratification and pressure gradients in the unperturbed air are relevant. Due to the short time scales (up to few hours) considered here, radiative processes can be ignored. This means that in cloud-free conditions, the dry potential temperature, θ, is conserved when following air parcels: where D·/Dt is the material time derivative and R is the gas constant of dry air. The variables x={θ,p.T,ρ} are decomposed into small fluctuations around a dominant, vertically varying background: A similar decomposition is made for the vertical wind speed field, w, assuming that w ′ ≫w . The fluctuations vary on shorter time scales than the mean values. Without loss of generality, we ignore horizontal wind components.
Lagrangian fluctuations may be viewed as the thermodynamical response (compression and expansion) of air parcels as they move vertically. The Lagrangian CHR follows immediately from equation (2): using the ideal gas law, p=ρRT, and assuming a hydrostatic pressure gradient, ∂p=∂z ¼ −ρ g. While the definition of κ embraces both cooling and heating-induced temperature changes, we refer at times to κ>0 as absolute values of cooling rates generated in updrafts (w ′ >0) in the remainder of this study.
In the Eulerian framework, the fluctuations arise from the net flow of air at a given location (advective response). Since θ is assumed to behave like a conserved tracer, motion along a vertical potential temperature gradient results in θ-advection producing ∂θ=∂t≃−w ′ ∂θ=∂z. This translates into the local Eulerian Ttendency (Fritts & Alexander, 2003): with N 2 ¼ g∂lnðθÞ=∂z.
For a given vertical wind fluctuation, Lagrangian and Eulerian temperature tendencies may have different magnitudes. This can be best illustrated by considering vertical motion in a neutral atmosphere; for ∂θ=∂z ¼ 0, Eulerian temperature fluctuations vanish. On the contrary, thermodynamic effects do still occur due to the hydrostatic pressure gradient leading to non-vanishing Lagrangian CHR.
Fluctuations evaluated in the two reference frames differ not only by their magnitude but also by their vertical and temporal variations. Figure 1 shows a time mean temperature profile along with the Eulerian variance, ðσ 2 T Þ EUL , and the Lagrangian variance, ðσ 2 T Þ LAG , derived from radiosonde observations. Assuming conservative vertical wave propagation (no wave dissipation), altitude variations of the Lagrangian fluctuation variance follow from ðσ 2 T Þ LAG ∝1=ðρN Þ , whereas ðσ 2 T Þ EUL ∝T 2 N 3 =ρ (Fritts & Alexander, 2003). In Figure 1, ðσ 2 T Þ EUL is estimated directly from the observations, while ðσ 2 T Þ LAG is deduced from a combination of equations (3) and (4): While ðσ 2 T Þ EUL increases strongly with altitude across the tropopause, ðσ 2 T Þ LAG decreases with altitude. The opposite evolution of MTF with respect to the tropopause in the Eulerian and Lagrangian frames of reference is largely caused by their different response to changes in the static stability, N . For comparison, Figure 1 also shows the profiles that would be expected if the fluctuations were caused solely by conservative wave propagation with the observed variance taken at the approximate mean cold point tropopause height, 17.5 km. The overestimation at higher altitudes and underestimation at lower altitudes observed in both cases likely result from loss of wave energy when propagating from lower to higher levels, as well as air motions unrelated to GW being more prevalent at the bottom of the profiles.

Data Sets on Temperature Fluctuations
Measurements on floating superpressure balloons (SPBs) provide direct evidence of GW as the key mechanism responsible for mesoscale vertical air motion variability. Importantly, they provide Lagrangian MTF properties suitable for cloud studies. By contrast, aircraft, radar, and sounding balloons provide data on horizontal or vertical transects and thus get only limited information that is often difficult to interpret. Data from various sources discussed below are summarized in Table 1. Moreover, we briefly comment on how our approach relates to previous attempts to model MTF.

Extensive SPB and Aircraft Data Sets
A detailed characterization of MTF properties from SPB measurements was provided for conditions above the tropical tropopause layer (TTL) and in the LS over Antarctica, including probability density functions (PDFs) of vertical wind speed and temperature fluctuations . During the long-term campaigns, research balloons floated for several months on isopycnic surfaces sampling a wide variety of forcings, encompassing day/night, convective/quiescent, cloudy/clear-sky, and land/ocean conditions. The measurements provided fluctuation amplitudes, T ′ and κ, as standard deviations; both variables were not measured independently. Standard deviations of corresponding w ′ -values for the long-term equatorial (polar) SPB measurements in the LS were σ w =0.11 (0.13) m/s. Analyses of wind and pressure fluctuations measured on more than 1,500 SPBs floating mostly between 18 and 21 km in the Southern Hemisphere LS are described in Schoeberl et al. (2017). These measurements extend the above SPB campaigns to include the Southern Hemisphere jet region. The data confirm the power law dependence of MTF and kinetic energy with frequency, with spectral slopes near −2. The slopes were found to be almost independent of latitude and altitude; MTF amplitudes were higher over land than over ocean. The mesoscale variance (power) of variable x={T,κ,w}, σ 2 x , is given as a frequency integral over the corresponding power spectral density (PSD), S x (ω), between f and N, To connect S T with S κ , we note that the latter follows from squaring the time derivative of Fourier components of temperature, _ T ðωÞ∝ωexp ð−iωtÞ: The observed near −2-slope of logðS T Þ versus logðωÞ implies that the CHR (or equivalently, the vertical wind speed) spectrum is approximately flat and non-zero within [f,N]. The dependencies S κ ∝ω 0 and S T ∝ω −2 also follow from dimensional analysis of properties of non-breaking GW in the midfrequency range (f≪ω≪N), that is, neglecting Coriolis and non-hydrostatic effects (Dewan, 1997). According to equation (3), σ w =σ κ /Γ.
Mesoscale temperature variations of UTLS air parcels, caused by adiabatic vertical displacements of isentrope surfaces, were derived from a large number of aircraft measurements in the northern and southern hemisphere (Gary, 2006;2008). MTF amplitudes were found to vary in systematic ways with altitude, latitude, season, and topography. The author proposed an empirical parameterization of MTF amplitudes scaled to different altitudes by an air pressure-dependent correction factor. The author notes that the proposed parameterization of MTF should not be applied in situations in proximity to, and therefore directly affected by, strong GW sources. Bacmeister et al. (1999) analyzed aircraft measurements in the midlatitude LS at flight altitudes 18-20 km, mostly over open ocean. The authors reported a mean MTF amplitude 0.89 K and a peak cooling rate 3.6 K/hr at horizontal scales 3-25 km. Aircraft measurements in the UT show w-fluctuation amplitudes in the range 0.05-0.2 m/s (Ström et al., 1997). Kärcher and Ström (2003) inferred mean w ′ =0.26 (0.23) m/s from aircraft measurements in the eastern North Atlantic (Southern Ocean). Hoyle et al. (2005) carried out an analysis of aircraft data for conditions in the UT over the continental United States and removed flights directly affected by lee waves from the data set prior to analysis. They report a distribution of MTF variances peaking within 0.5-2 K 2 around a mean value of 1.24 K 2 . Jensen et al. (2013) analyzed aircraft measurements Note. Observations include aircraft (AC), ground-based Doppler radar (DR), superpressure balloon (SPB) measurements, and sounding balloons (SB). Fluctuation amplitudes are denoted by T ′ (temperature), κ (cooling/heating rate), or w ′ (vertical wind speed), making no distinction between standard deviations and mean values. Only two data sources provide the fluctuation amplitude ratio, ν, determining the cooling rate associated with a given temperature fluctuation.

Limited Aircraft, Radar, and Sounding Balloon Data Sets
All data sets represent background conditions as defined here only to a degree and are not necessarily intercomparable; only SPB data represent Lagrangian fluctuations. If κ and w ′ were not reported simultaneously, we used the dry adiabatic lapse rate for data conversion. LS = lower stratosphere; NH = northern hemi-  2015) inferred a statistic of mean vertical wind speed variance from ascent data of 71 sounding balloons. The statistic is representative of the region in the vicinity of Zurich and encompasses data taken over a period of 5 years, confirming that mean variances in the same area can be highly variable. Ground-based Doppler radar retrievals of vertical wind speed in cirrus (Kalesse & Kollias, 2013) are reported in two locations over the central United States (Southern Great Plains) and the tropical Western Pacific (Manus), yielding mean w ′ =0.1−0.2 m/s in both cases (Barahona et al., 2017).
In agreement with SPB measurements, flat vertical wind speed PSDs are inferred from aircraft measurements on horizontal scales >6 km mostly over oceans (Bacmeister et al., 1996), and in radar measurements in low background wind conditions, minimizing Doppler shift effects that distort intrinsic wave frequencies (Ecklund et al., 1986). Measurements performed with the Middle Atmosphere Alomar Radar System installed in northern Norway showed mildly sloped S w frequency dependencies in low wind conditions (Li et al., 2018), probably due to mountain wave activity.

Model Simulations
Lagrangian MTF time series have been computed based on superposition of harmonics with independent random phases and assumed temperature PSDs Hoyle et al., 2005;Jensen & Pfister, 2004). In microphysical applications, the fluctuations are then superimposed on large-scale temperature trajectories. Here, we improve on these efforts to synthesize wave-induced MTF by specifying PDFs of the underlying vertical wind speed fluctuations and developing a physical model connecting both variables. Kienast-Sjögren et al. (2015) present a comprehensive analysis of temperature and vertical wind velocity PSDs for areas over Switzerland from numerical weather prediction model-derived trajectories at variable horizontal and temporal resolution. The trajectories are used in conjunction with a Lagrangian cloud model constrained with measured vertical profiles of cirrus optical extinction.
Global fields of vertical wind speeds were estimated using 2 years of high (7 km horizontal) resolution global model simulations (Barahona et al., 2017). The simulations resolved waves with horizontal wavelengths of several tens of kilometers and employed a theoretical scaling of unresolved w-variance. This scaling has been developed in the context of improving updraft speeds associated with deep convection in coarse resolution models. It is unclear how well these simulations capture unresolved variance due to high frequency GW.
The above modeling studies emphasize the need to account for unresolved high frequency temperature fluctuations to accurately simulate or parameterize cirrus formation and evolution. Here, we develop a versatile model-constrained by quasi-Lagrangian measurements-capable of predicting MTF time series that cover the full spectral range of GW-induced perturbations.

A Cautionary Note on Data Intercomparison
Intercomparison of fluctuation amplitudes between different data sets is hindered, since not all measurements represent similar conditions as defined in section 3.1. Aircraft, radar, and sounding balloon data do not provide MTF or CHR in a Lagrangian reference frame; approximate intrinsic spectral information is only available in low background wind conditions. Uncertainties in w ′ from aircraft measurements can approach values similar to those due to GW depending on temporal data resolution (Mallaun et al., 2015). Individual measurements of T ′ in the (Gary, 2006;2008) data sets scatter by ±0.4 K-a sizeable amount of GWinduced MTF.
With regard to our model, we interpret amplitudes of fluctuating variables as standard deviations, σ\, of their respective probability distributions. This notion is consistent with the SPB measurements. The associated variances are viewed as averages of a great number of air parcels, that is, ensemble averages of different statistical realizations of fluctuation time series. When using measured fluctuation amplitudes as model input, it must be known whether the measured data represent mean, mode, median values, absolute or standard deviations of the underlying PDF and whether data sampling or instrumental limitations cause biases in the determination of these variables. These issues also introduce uncertainty when comparing results from different observation methods.

Definition of Background Conditions
The probability of encountering intermittency-the fact that high CHR are significantly more abundant than in a Gaussian distribution-increases with proximity to strong, high frequency GW source regions. SPB measurements reveal exponentially distributed vertical wind speed fluctuations and thereby underline the importance of intermittency . By contrast, associated temperature PDFs may or may not be Gaussian. To further illustrate this point, Figure 2 depicts SPB time series of MTF and associated CHR taken over the course of one season in different regions. The data show that intermittent events are absent in equatorial MTF samples and occur very rarely in polar MTF samples (top panel). However, intermittency is observed in both, equatorial and polar time series of CHR, as indicated by the relatively more abundant large-amplitude GW events (bottom panel).
Aircraft data that sample a range of meteorological conditions revealed that strong, infrequent temperature perturbations created by large-amplitude mountain waves produce non-Gaussian wings in PDF(T ′ ), especially at horizontal scales below about 10 km . Only when such events are removed from data sets during analysis or are absent in the first place, MTF are approximately normally distributed.
These observations motivate us to introduce two PDFs describing GW-induced Lagrangian fluctuations, w ′ and T ′ , with zero means. A normal distribution for the T ′ -statistic with variance σ 2 T (Gary, 2006): and an underlying w ′ -statistic approximated by a Laplacian with variance σ 2 w : where μ w is the mean value of the updraft speed distribution.
The kurtosis K-defined as the fourth PDF moment normalized by the square of its variance-measures the flatness of a PDF around its core relative to the steepness of its wings; data within ±σ w contribute very little  . The data were collected during boreal spring (equatorial flights) and during austral spring (polar flights). They have been filtered to retain only frequencies between 2π/1 day and N/2 to remove contributions from low frequency planetary waves and avoid complications with the interpretation of balloon motions at high frequencies.
to K. The kurtosis serves as an indicator for the likelihood of occurrence of outlier events with values far exceeding the mean. The excess kurtosis, K−3, is used to describe tail extremities of a given PDF relative to that of a Gaussian that has K=3. Equation (8) yields K=6; thus, PDF(w ′ ) is leptokurtic (K>3), with comparatively fat tails and a slender core.
The ubiquity of GW-induced MTF in the UTLS suggested by all types of observations motivates the definition of background conditions as UTLS regions not directly influenced by high frequency (mainly orographic and convective) GW sources. We propose conditions be called "background," if PDF (w ′ )=L(σ w ) and PDF(T ′ )=N(σ T ). Non-background situations are realized when PDF(w ′ ) is non-Gaussian with K>6 or when PDF(T ′ ) exhibits non-Gaussian behavior.
The requirements expressed in equations (7) and (8) rest on an idealization that in nature may be realized only to a degree. However, Figure 2 along with the associated PDF(T ′ ) and PDF(κ\) presented by Podglajen et al. (2016) provides evidence that this background concept has some practical applicability. Therefore, the notion of background conditions is meaningful in analyzing and categorizing atmospheric observations.

Cooling Rate Model
We use PSDs to derive a relationship between MTF and CHR. We recall from section 2 that the SPB measurements link MTF directly to mesoscale GW and that the PSD of CHR fluctuations follows a top hat distribution: Combining equation (9) with equation (6) yields: In the above, all frequencies are non-zero and positive, and the variances are ensemble averages over many air parcel trajectories. After taking the limit f≪N, we find the fluctuation amplitude ratio: this result relates σ T to GW-induced background cooling rates, σ κ . As a geometric mean, ν lies midway between f and N. Viewed as an inverse time scale, ν determines the average magnitude of the cooling rate in response to a decrease in temperature due to GW.
Given σ T , the CHR variance, σ 2 κ , increases with N, as greater static stability leads to faster buoyancy oscillations of air parcels and therefore faster temperature changes. CHR decrease as ν approaches f, that is, σ 2 κ decreases with latitude. The f-dependence is generally in line with latitudinal variations of GW-induced temperature variance (Alexander et al., 2002). For f→0, a correction to equation (11) is required to obtain cooling rates consistent with the balloon data in the equatorial region. We propose such a correction below.
An open issue is the observation of an increase in CHR variance as ω→N (Ecklund et al., 1986;Podglajen et al., 2016), modifying the top hat distribution equation (9). Quantification of this effect, potentially caused by wave reflection and resonance, is hampered in the case of SPB measurements due to the non-isopycnic balloon response at frequencies near N/2 (Podglajen, 2017).
Two SPB data sets are available providing concomitant Lagrangian measurements of σ T and σ κ in polar (Antarctica) and equatorial (above TTL) LS conditions . Values N=0.022 (0.026) rad/s were obtained from a global numerical weather prediction model, and f=2π/12 hr (<2π/1 day−2π/ 2 days) represent long-term averages for polar (equatorial) conditions in the LS. However, close comparison is hampered because (a) the polar data do not represent pure background conditions (probably due to mountain wave activity over the Antarctic peninsula), as PDF(T ′ ) was found to be non-Gaussian violating equation (7), and (b) equation (11) predicts vanishing cooling rates at the equator.
ν=8.4 hr −1 , only by about 30%; at this point, it remains unclear whether this slight discrepancy is a matter of data-model comparability (section 2.4). Regarding (b), we use the data to constrain a lower-limit f-value. Given an oscillation period of 2 days, f≃Ω/2=3.65×10 −5 s −1 , we get ν=3.5 hr −1 , a value that also lies within 30% of the measurement, ν=5 hr −1 . We redefine f via f ¼ 2Ω sinðϕÞ : ϕ≥ϕ eq 2Ωα : ϕ<ϕ eq ( : The imposed lower-limit frequency, Ω/2, corresponds to α=1/4 and a cutoff latitude of ϕ eq ¼ arcsinðαÞ of about 15°. We display the corrected amplitude ratio in Figure 3 as a function of latitude for two representative N-values. As a first attempt to theoretically derive a lower-limit f-value, we relate α>0 to a dispersion relation of equatorial GW in Appendix A. The order of magnitude of the cooling rates listed in Table 1 is reasonably well captured by equation (11).
The two SPB measurements are in line with the predicted dependence of ν on latitude ϕ, in part because the balloon time series were filtered to keep only frequencies exceeding f. More latitudinal coverage is needed to better judge the quality of the background cooling rate model.

Probabilistic Description of MTF
This section presents a minimal approach fusing important aspects of observations of GW-induced mesoscale variability in vertical wind speeds and temperature within a stochastic model framework, explores its theoretical underpinning based on statistical physics, discusses numerical solutions, and comments on how to apply it in cirrus models.

Damped-Driven Temperature Fluctuations
We regard Lagrangian air parcels as dynamical systems, whose temperature is adiabatically perturbed by a large number of random vertical velocity fluctuations. Vertical wind speeds, w, determine the slowly varying or constant mean conditions of air parcel ensembles. Here, we set w ¼ 0. The fluctuations, w ′ , due to GW, account for the non-Gaussian shape of its statistics. We introduce a Newtonian relaxation of the resulting MTF, T ′ . While T ′ and w ′ are induced by GW, we regard T ′ as being driven (forced) by w ′ and use the term damping to refer to dynamic instabilities leading to dissipation of low frequency waves. We employ κ instead of w ′ in the following discussion.
Temperature fluctuations around a mean value, T ≫T ′ (either constant or changing on a long time scale depending on w), are described by with the MTF damping rate λ>0; ψ κ defines the random CHR. In addition to representing κ=−Γw ′ according to equation (3), Ψ κ indicates that κ-values (a) evolve around a zero mean; (b) are drawn from equation (8) with the (ensemble-averaged) total variance, σ 2 κ ; and (c) are generally a function of time altitude, latitude, season, and topography. According to the central limit theorem of probability theory, when a variable is forced by a large number of random, statistically independent influences, its probability distribution becomes Gaussian. This holds regardless of the nature of the forcing requiring only that it has a well-defined variance. It is the fundamental reason why MTF obtained from this model follow equation (7). Here, we interpret random influences as the local action of many GW packets inducing random displacements of air parcels in the UTLS.
After these general remarks, we revisit the relationship between the stationary MTF amplitude, σ T , and the associated CHR fluctuation amplitude, σ κ (section 3.2). We link their ratio to the fundamental statistical physics underlying equation (13) and consider observed PSD properties as embodied in equation (11) as an additional constraint.

Fluctuation Amplitude Ratio
Equation (13) is a stochastic differential equation with linear, additive forcing. It describes a random walk in T-space (or equivalently, air parcel altitudes) and is formally similar to Brownian motion. This analogy becomes readily apparent when air parcels are replaced by microscopic particles immersed in a fluid; T ′ is replaced by fluctuations of the particle's velocity; −λT ′ is interpreted as the friction force due to fluid viscosity limiting the particle's velocity; and Ψ κ is regarded as the fluctuating particle acceleration due to thermal agitation by surrounding molecules. In our case, damping ensures that the temperature variance averaged across an ensemble of air parcels becomes stationary for t≫1/λ. In the absence of damping, σ 2 T grows linearly with time and without limit, reminiscent of a purely diffusive process.
We derive the relationship between σ 2 T and the variance of the forcing, σ 2 κ . We define changes in T ′ within a time interval τ due to random CHR fluctuations: where τ is chosen small enough for the resulting T ′ -values to stay small, but long enough for the fluctuations to alter their sign many times. In the case of GW, we demand that τ be of the order of a buoyancy period, 2π/ N. The solution of equation (14) is nontrivial as ψ κ is only defined at discrete times, but it is not necessary to solve the stochastic integral explicitly.
On the one hand, we know that the stationary limit of the temperature fluctuations is related to the timeintegrated temperature forcing, equation (14), via: Equation (15) is obtained by integrating equation (13)  On the other hand, < G 2 τ > is equivalent to a double integral over the ensemble-averaged correlation function of the CHR fluctuations: which we solve as follows. By design, ψ κ (t) oscillates in a highly irregular manner around a zero mean value. Therefore, the averaged correlation function, Ψ, will be zero, too, for sufficiently large values of |t ′ −t| required to be smaller than a short, characteristic coherence time, t c <τ ¼ Oð2π=NÞ . Therefore, we can assume that Ψ depends only on the absolute value of time differences: In line with the Markovian character of equation (13), κ-values arising from ψ κ are sharply correlated only within t c . The correlation function equation (17) embodies this behavior.
We finally transform equations (16) and (17) to the variables u=t ′ +t and v=t ′ −t and note that Ψ is independent of u. With dt dt ′ =du dv/2, this leads to Combining equations (15) and (18) yields the fluctuation amplitude ratio relating the power, σ 2 κ , of the small-amplitude CHR fluctuations (forcing) to the magnitude, λ, of the linear dissipative response (damping) of MTF at long times.
Our final task is to constrain t c and λ. Comparing equation (19) with (11) implies that these variables cannot be chosen independently. We already have the constraint t c <τ and explicitly consider the PSD of the temperature fluctuations to constrain λ. According to equation (13), damped MTF are correlated via (Becker, 1978). According to the Wiener-Khinchin theorem (Gardiner, 2004), S T (ω) is given by the Fourier transform of the autocorrelation function equation (20): so that σ 2 T represents the MTF variance contained within [f,N]. This implies λ ¼ Oðf Þ in order to dampen frequencies 0<ω<f (since f is the lowest possible intrinsic GW frequency) and to render S T ∝1/ω 2 for ω≫f consistent with observations (section 2.3). Consistent with equations (11) and (19), we set Equation (21) shows that λ prevents S T from diverging as ω→0. Therefore, σ T becomes asymptotically stationary, and S T approaches a constant scaling ∝σ 2 T =f for low frequencies. In analyzing the SPB measurements, Podglajen et al. (2016) used a filter to remove low frequency (ω<f) motions. In our approach, no such filtering is necessary due to the inclusion of damping.
A physical theory consists of a mathematical framework and the interpretation of its predictions. While we have exploited the mathematical similarity between the random walk statistics of particle movement forced by molecular collisions to that of air parcels forced by GW-induced, vertical wind speed perturbations, we acknowledge limits in the interpretation of this analogy.
First, while in the case of Brownian motion, t c -set by the time between molecular collisions-is very much shorter than any associated macroscopic time scale, values t c ¼ 1=N ¼ Oð1 minÞ should not be neglected in microphysical applications of equation (13). The short autocorrelation time is responsible for the high frequency variability in CHR influencing ice crystal formation (Kärcher et al., 2019). Essentially, non-zero t cvalues from equation (22) produce a short-term MTF autocorrelation that results in an additional cutoff of S T (ω) at high frequencies. We revisit this issue in section 4.3.
Second, in the case of Brownian motion, fluctuations and dissipation have the same origin (molecular collisions). Their rates are linked by the fluctuation-dissipation theorem, as the system of molecules and particles is constrained by thermodynamic equilibrium: The stationary velocity variance averaged over many particles is given by the mean thermal speed. No such constraint exists in our case; instead, we have used PSDs of MTF and CHR to determine the fluctuation amplitude ratio.
Our model replicates aspects of atmospheric observations that are highly relevant in future microphysical applications, but the compatibility with equations (2) and (3) is not immediately obvious. We consider equation (13) as a minimal approach based on few parameters and assumptions without providing tighter links to GW physics. For example, we specify the noise source term, ψ κ , empirically based on observed PDFs and PSDs of GW-induced fluctuations, and we ignore any possible scale or frequency dependence of λ that might tie the damping to actual wave dissipation mechanisms. Future studies might improve on this situation as needed.

Numerical Solutions
Equation (13) is amenable to Monte Carlo simulation. We solve it recursively over a time step, Δt, with an implicit damping term (λ=f) in the difference form with the time step-integrated temperature forcing Discretizing the stochastic integral equation (24) over time to turn it into a regular (Riemann) integral introduces an autocorrelation time of magnitude δt (set equal to 1/N to capture each fluctuation in w).
We explore whether numerical solutions obtained in this way reproduce analytical solutions. For illustration, we prescribe ϕ=50°(f=1.12×10 −4 rad/s), N=0.012 rad/s, and σ κ =3.53 K/hr. We run the model across an ensemble of 5,000 air parcels for up to three damping time scales, 3/f≈7.5 hr, to ensure that the temperature fluctuations evolved into a quasi steady state. Figure 4 shows the temporal evolution of the MTF variance, σ 2 T ðtÞ, with and without (λ=0) damping. Without damping, σ 2 T ∝t, consistent with the diffusive nature of the random walk process. The simulations also confirm that damping causes σ 2 T to approach a stationary value: Averaged over the time period 2/f−3/f, we find that the simulated σ T -value agrees with the expected value, σ κ /ν≃0.85 K, within 1%. This means that fluctuations more than a few σ T around the mean are highly unlikely to occur, as expected for mesoscale variations. The relationship between MTF and clear-sky ice supersaturation variance is explored in section 5.2. Figure 5 depicts the PDFs of w ′ and T ′ taken at a time (7.5 hr) when the temperature fluctuations are stationary. The solid curves display (a) the vertical wind speed fluctuations calculated according to equation (8) used as input to the simulations and (b) the resulting temperature fluctuations. The dashed curves follow from equations (7) + (8) using the simulated values σ T and σ κ , respectively. Simulated and expected PDFs agree very well. In particular, the numerical results confirm that MTF evolve into Gaussians as expected from the central limit theorem. We also note that the ensemble mean of the simulated w ′ perturbations is practically zero (<w ′ >≈0.01 cm/s).
Three examples of MTF time series, T ′ (t), are displayed in Figure 6. Those might be superimposed, for example, onto a synoptic trajectory, T ðtÞ. A range of mesoscale frequencies can be identified visually, the highest frequencies dominating cooling rates locally within short periods of time on the order of 1/N=1.39 min, with

Journal of Geophysical Research: Atmospheres
important ramifications for ice crystal nucleation in cirrus (Kärcher et al., 2019). According to the statistical law of low leads, only a small subset of MTF paths behave like the black curve oscillating closely around zero net temperature excursion. Along most MTF paths, air parcels experience multiple consecutive cooling and heating events; over time, the resulting temperature excursions can become unrealistically large, if damping were neglected (cf. Figure 4). The colored curves represent the most extreme cooling and heating cases that occurred in this set of simulations. For example, the blue curve leads to an overall decrease in temperature of 1.8 K within the first 3.6 hr, corresponding to an average cooling rate of 0.5 K/hr. Inspection of the vertical wind speeds driving the MTF along this trajectory shows that significantly faster cooling rates (e.g., 18 K/hr) occur locally due to high frequency contributions (ω≈N) in w ′ .
We have also investigated the temperature fluctuation PSD, S T (ω), of the air parcel ensemble by deriving the PSD via Fourier analysis of all 5,000 MTF time series and averaged the results (not shown). As expected, the spectrum follows a red noise behavior in the midfrequency range leveling off at low frequencies due to damping as predicted by equation (21). Furthermore, S T rolls off at high frequencies as explained in section 4.3, because the forcing loses power due to autocorrelation of unresolved vertical wind speed fluctuations.

Implementation in Cirrus Models
We provide with equation (23) a numerical tool to simulate GW-induced background cooling rates in microphysical cloud models. To illustrate its application in a parcel model framework, we summarize the most important quantities and equations.
Selecting the Coriolis frequency, f, fixes the damping rate, λ, and selecting the Brunt-Väisälä frequency, N, fixes the autocorrelation time, t c =1/N, according to equation (22). This determines the time step δt=t c , which may be larger than the time step, Δt, used in a microphysical model to accurately resolve ice crystal nucleation events. Lagrangian fluctuations according to equation (3) relate CHR, κ, to vertical wind speed perturbations, w ′ . Discrete w ′ -values are sampled from the Laplacian distribution equation (8) at integer multiples of t c and interpolated as stair steps to arbitrary times, t. MTF time series, T ′ (t), then follow from applying the stochastic equations (23) + (24) with a prescribed variance, σ 2 w . In practice, is added to temperature tendencies from mean vertical motion and latent heat exchange before numerical integration. As a result, MTF are superimposed onto temperature trajectory data that lack variability below times 1/f. If the trajectories are taken from a host model that resolves GW activity at shorter time scales, 1/ω ⋆ with ω ⋆ >f, then f must be replaced by ω ⋆ to avoid double counting of MTF variance.
We emphasize that σ 2 w must be selected with care to enable meaningful comparisons with observations. When using equation (11), the vertical wind speed variance can only be determined if σ T is known. Alternatively, σ 2 w may directly be obtained from estimates of GW momentum fluxes. Improvements of ice nucleation parameterizations in cloud schemes of global models are discussed in section 5.1.

Atmospheric Implications
This section explores two variables that must be known with confidence to improve cirrus modeling and better constrain uncertainties associated with aerosol-cloud interactions. First, we indicate how cirrus parameterizations driving ice nucleation parameterizations in global models might be advanced by a better representation of unresolved wave-induced updraft speeds. Second, we indicate how advection- condensation simulations based on large-scale model trajectories might be improved by superimposing unresolved MTF.

Cloud-Scale Updraft Speeds
Mesoscale variability in vertical wind speeds has a strong impact on homogeneous aerosol freezing and alters effects of heterogeneous ice-nucleating particles in cirrus formation. Lohmann and Kärcher (2002) proposed to use, for first exploration, a GCM's turbulence scheme to estimate unresolved updraft speeds driving the model's ice nucleation parameterization. The latter were associated with turbulence kinetic energy (TKE) and added to the resolved (large-scale) updraft speeds, which are typically much smaller. This dynamical forcing of ice nucleation provided reasonable agreement of interactively simulated cirrus properties and UTLS ice supersaturation statistics, despite the fact that the mesoscale w-fluctuations in the UTLS are controlled by GW and therefore TKE is not the most relevant forcing process.
Nonetheless, the TKE approach is still used in climate models; some models use a single updraft speed to calculate nucleated ice crystal numbers. Two models take TKE-induced variations of w into account by employing a Gaussian probability distribution of w ′ with constant lower and upper limit values (Salzmann et al., 2010;Bacer et al., 2018). Accounting for variability in w in this way does not bring the representation of cirrus clouds closer to the actual processes governing mesoscale cooling rates in the UTLS. Recently, Penner et al. (2018) used equation (8) to sample updraft speeds, taking into account the non-Gaussian nature of PDF(w ′ ). The underlying mean updraft speed fluctuation was taken from the equatorial SPB measurements discussed above along with a static stability scaling of the autocorrelation time after which the values from stochastic w ′ -time series were updated. Moreover, seasonal and geographical variations were introduced based on the σ T -scalings provided by Gary (2006). Global models parameterize tropospheric GW sources to represent momentum flux generation and deposition in the atmosphere. Since high frequency GW carry significant momentum flux vertically, a major focus lies on examining GW effects on the middle atmosphere circulation. For consistency between the model representation of different GW-dependent processes, a sensible approach to represent GW-induced CHR in cloud schemes would be to employ a GW momentum flux parameterization. For a monochromatic wave packet, the vertical wind speed variance associated with the vertical flux of horizontal momentum due to GW is given by (Ern et al., 2004), where ||F p || is the norm of the momentum flux vector. This allows total σ 2 w to be estimated from the superposition of variances from different wave sources, capturing the pronounced variability in observed mean σ 2 w -values (Table 1). The precise implementation of equation (26) will depend on the type of parameterization, for example, by ray tracing that relies on geometric optics to track paths and properties (amplitudes, wavelengths) of wave packets propagating away from source regions at their local group velocities.
As an intermediate solution to parameterize σ 2 w due to background (non-orographic) GW, equation (11) may be employed with an effective Coriolis frequency in equatorial regions, equation (12). Following Penner et al. (2018), σ T may be taken from (Gary, 2006;2008) and validated against other data as far as possible. Care must be taken to avoid double counting of w-variance when this approach is used in combination with an orographic GW parameterization. As models for global weather and climate prediction move toward finer resolution, they might resolve waves with intrinsic frequencies higher than f. In such cases, it may be necessary to replace f with a higher, resolution-dependent wave frequency (section 4.4). This issue is relevant globally, but particularly challenging in the deep tropics. As tropical waves are coupled to deep convection, the interplay between GW-induced vertical wind speed fluctuations and deep convection warrants closer scrutiny.
To advance cirrus parameterizations in global models, it is not sufficient to alter σ w . First, detailed microphysical ensemble simulations are needed to assess how mesoscale CHR variability affects the ice formation process-for example, in terms of expectation values for nucleated ice crystal numbers-that will differ from those obtained with a purely deterministic approach. Kärcher et al. (2019) report a first attempt to tackle this issue. Second, small-scale turbulence is associated with entrainment and mixing of air parcels. In the presence of strong turbulence, the outcome of ice nucleation events may differ from adiabatic simulations (Kärcher & Jensen, 2017). Therefore, coupled turbulence-microphysics studies resolving the influence of inhomogeneous mixing of air and molecular diffusion of water vapor on ice nucleation at the dissipation scale (<1 cm) are needed to characterize and quantify the TKE-related contribution to cirrus formation.

Clear-Sky Ice Supersaturation Fluctuations
The stochastic processes discussed above operate along air parcel trajectories. Studies have used trajectories based on global model wind and temperature fields to estimate the effect of advection and condensation on the evolution of water vapor (Fueglistaler et al., 2014;Pierrehumbert et al., 2007;Schoeberl et al., 2019). In such studies, fluctuations of water vapor concentration and temperature at scales smaller than those resolved in a large-scale model must be parameterized. While they may not substitute cloud-resolving models, such kinematic studies might also be used to estimate conditions conducive to ice nucleation by tracking ice supersaturation prior to ice cloud formation. We present a diagnostic technique to compute the magnitude of Lagrangian supersaturation fluctuations induced in the UTLS by mesoscale GW.
We relate MTF to instantaneous fluctuations in ice supersaturation, defined as the ratio between the water vapor partial pressure, ê, and the saturation vapor pressure over ice, e s , above ice saturation: s ¼ ðê=e s Þ−1; here, we use e s ¼ e ⋆ exp½−L s =ðR w TÞ, with the gas constant of water vapor, R w , the latent heat of sublimation, L s /R w =6132.9 K, and e ⋆ =3.4452×10 10 hPa (Marti & Mauersberger, 1993). The term clear sky includes the presence of ice clouds with number concentrations of ice crystals low and sizes small enough to render deposition and sublimation processes inside cloud inefficient, which would otherwise lead to strong microphysical damping of supersaturation fluctuations (Kärcher et al., 2014).
We view s as a function of the stochastic variables ê and T and derive a relationship connecting their variances including adiabatic correlation. Expanding sðê; TÞ in a Taylor series to second order in the fluctuations of the independent variables yields for the supersaturation variance, σ 2 s ¼ ðs−sÞ 2 : where the overbar denotes a mean value, ρ c is the linear correlation coefficient, and L ¼ L s =ðR w T Þ.
Fluctuations in ê may be thought of being composed of two additive, uncorrelated contributions,ê ¼ Δe ad þ e. The first arises from adiabatic temperature change, and the second shall be connected with non-adiabatic processes (e.g., turbulent mixing). In Appendix C, we derive the associated correlation coefficient: where η=7/2 is the adiabatic index and the dispersion (relative spread) of a fluctuating variable, x, is denoted by δ Values δ e =0.1−0.2 have been derived from aircraft measurements near ice saturation (Gierens et al., 2007). However, these results represent horizontal variability in water vapor as captured during aircraft sampling. Spatial variability has been studied in the same way by Diao et al. (2014). As such, this information is not directly applicable in the context of changes due to GW. In air parcels undergoing small vertical displacements, events in which air parcels entrain and eventually mix with ambient air (leading to δ e >0) are rare in background conditions, for example, for typical turbulence energy dissipation levels in the TTL . Therefore, the only way by which ê is affected in air parcels undergoing stable oscillations is due to adiabatic change (i.e., ê ¼ Δe ad and therefore δ e =0 and ρ c =1), in which case equation (27) reverts to Equation (29) implies that changes in the spread of the supersaturation statistic scale in proportion to σ T and adiabatic motion reduce the contribution of temperature fluctuations, Lσ T =T, to the supersaturation dispersion by the factor ð1−η=LÞ. The relevance of Lagrangian MTF in determining the clear-sky supersaturation variance increases with decreasing T and increasing s.
Evaluating equation (29) with σ T =1 K and s ¼ 0 yields σ s ≃0.11 for T ¼ 220 K (L≃28). Such MTF levels generate non-Gaussian supersaturation statistics with significant positive skewness (Kärcher & Haag, 2004). For increasing s>0, this causes non-equilibrium states with levels of ice supersaturation (Kärcher & Burkhardt, 2008) large enough to influence the competition of different types of ice-nucleating particles for available water vapor during ice formation events. The combination of unresolved ice supersaturation variance and better estimates of σ 2 w as outlined above enables improved studies of cirrus occurrence within the trajectory-based advection-condensation modeling framework in future studies.

Summary
A series of long-term, quasi-Lagrangian observations suggests that MTF due to GW could be represented as a first-order autoregressive process . Here, we generalize this approach to the global UT and LS and devise a minimal probabilistic model describing the Lagrangian evolution of temperature fluctuations driven by a random wave spectrum and damped by a low frequency dissipation mechanism. The model addresses the cloud physics community and may serve both theoretical investigations and numerical simulations.
A brief literature review substantiates the notion that fluctuations due to GW are commonplace in cirrus levels. We bring the issue one step further by developing a stochastic model for background fluctuations away from clearly identified, high frequency wave sources, consistent with observations. Properties include the probability of occurrence of fluctuation amplitudes and their magnitude and frequency components. We stress that Lagrangian fluctuations are to be used to drive cloud models and make use of radiosonde observations to illustrate and quantify the potentially large differences in fluctuation amplitudes evaluated in Eulerian and Lagrangian reference frames. We define background conditions based on the shapes of observed, exponential probability distributions of vertical wind speed and associated normally distributed temperature fluctuations. We constrain the ratio of cooling rate and temperature fluctuation amplitudes by their observed power spectra. This ratio determines the average cooling or heating rate associated with a temperature fluctuation due to mesoscale GW in background conditions. In the stochastic model, temperature fluctuations are driven by random, autocorrelated vertical wind speeds and are subject to damping. We draw a mathematical analogy to Brownian motion, allowing us to identify the damping rate with the Coriolis frequency and the autocorrelation time with the inverse Brunt-Väisälä frequency. Low frequency damping ensures that the temperature fluctuation variance remains finite for all frequencies. The non-Gaussian wind speed fluctuations-illustrated by SPB measurement time series -are autocorrelated over a fraction of a buoyancy period. The model is consistent with the above fluctuation amplitude ratio, containing dependencies on atmospheric stratification and geographical latitude thereby enabling applications in the tropics, at midlatitudes, and in the polar regions. In the tropics, effects of mesoscale and larger scale waves are difficult to disentangle. At this point, we suggest to use an effective equatorial Coriolis frequency to generate cooling rates in the tropics that agree with balloon data. More research is needed to improve this ad hoc solution. This model may be used as driver for ensemble simulations of ice crystal nucleation leading to cirrus cloud formation.
We discuss implications of our findings for the representation of updraft speeds in cirrus ice crystal nucleation parameterizations used in the cloud schemes of global models. We suggest to modify the representation of cirrus formation in low-resolution models to account for unresolved updraft speeds that represent widespread and ever-present background, and more localized and sporadic orographic and convective GW activity and represent their effect on nucleated ice crystal number concentrations. Such a framework may be based on parameterizations of GW momentum fluxes. To aid kinematic studies of processes influencing the humidity distribution along large-scale air parcel trajectories, we estimate the magnitude of unresolved Lagrangian fluctuations of ice supersaturation. Improved cooling rate and ice supersaturation fluctuation amplitudes might be combined to better predict ice crystal formation conditions in an advectioncondensation model framework.

Outlook
Future studies might provide tighter links between our stochastic model and the physics of GW generation, propagation, and dissipation. While the model results discussed here are generally in line with observations, data from different measurement platforms that may serve as constraints are not easily intercomparable for a number of reasons. The predicted dependence of the fluctuation amplitude ratio on Coriolis and Brunt-Väisälä frequencies should be systematically investigated as functions of altitude, latitude, season, and topography in future work, as current measurements are sparse. We also suggest to test the proposed definition of background conditions by concomitant measurements of probability distributions of mesoscale temperature and cooling/heating rate fluctuations. Long-duration Lagrangian measurements are ideally suited for model validation. On the theoretical side, nonlinear coupling of random perturbations is responsible for the non-Gaussian shape of probability distributions of vertical wind speeds. More advanced analyses should attempt to formulate stochastic equations with nonlinear, multiplicative noise terms that additionally allow for non-Gaussian temperature fluctuation statistics caused by intermittent, large-amplitude waves.
Our findings are relevant for advancing research into the physics of ice clouds. This applies to processoriented simulations, the improvement of parameterization schemes, and the interpretation of field measurements. With regard to process simulations, proper assessment of the interaction between different aerosol particle types during ice crystal nucleation requires coupling of detailed microphysical models with wave-induced vertical wind speed and temperature variability. Concerning efforts to parameterize ice clouds in models with coarse resolution, it will be necessary to establish a predictive capability for waveinduced vertical wind speed variances and link them to synoptic situations. Relating field measurements of cirrus formation to process-oriented models will benefit from advancing ice crystal nucleation simulations. In view of the localized, random nature of ice formation events, this requires measurement planning and analysis strategies explicitly designed for this purpose. Ultimately, progress along these lines will advance fully interactive global model simulations of cirrus.
It is important to predict cirrus formation with confidence on the process level before parameterizing aerosol-cloud interactions. Therefore, progress in characterizing the mesoscale dynamical environment of cloud ice formation in concert with dedicated observations will greatly improve our ability to carry out climate projections. While we recognize the significant progress in measuring number concentrations and chemical composition of ice-nucleating aerosol particles (DeMott et al., 2011), we can only lament the scarcity of in situ measurements characterizing the local dynamical environment during cirrus formation events.

Appendix A: Effective Coriolis Frequency
To derive an effective Coriolis frequency that leads to cooling rates in the equatorial region consistent with observations, we recall equation (1) and note that the dispersion relation for fixed-f hydrostatic GW with long horizontal wavelengths (k 2 +l 2 ≪m 2 ) takes the form: This may be compared with the dispersion relation for equatorial inertia GW (k≫l; Matsuno, 1966):

10.1029/2019JD030680
Journal of Geophysical Research: Atmospheres where β=Ω/R E is the Rossby parameter, R E =6,371 km is the radius of the Earth, and q=0,1,2,… is the wave mode number. This comparison lends itself to the definition of an effective Coriolis frequency describing conditions in the equatorial area: or, applying equation (12) and m=2π/Λ, the sine of the equivalent cutoff latitude ϕ eq : A shorter vertical wavelength, Λ, confines the gray-shaded zone depicted in Figure 3 to a narrower latitude band. Inserting N=0.022 rad/s and α≃1/4 obtained for the equatorial SPB measurements (section 3.2) leads to the constraint Λ≃33 km/(2q+1) for vertically propagating waves, generating the observed fluctuation amplitude ratio of ν=5 hr −1 . For equatorial latitudes, ν ¼ ffiffiffiffiffiffiffiffiffi ffi f eq N q ∝N 3=4 Λ 1=4 . equation of state connecting e and T. We associate Δe ad and ΔT with small but macroscopic changes from the mean states, ē and T , and set the macroscopic changes equal to standard deviations, σ: We derive the correlation coefficient, ρ c , between temperature and H 2 O partial pressure due to adiabatic coupling by evaluating the covariance of the joint water-temperature distribution. To this end, we split the H 2 O partial pressure into a part caused by adiabatic temperature changes, Δe ad , and an additional part, e, that is assumed to statistically independent from (orthogonal to) the T-fluctuations: ê ¼ e þ Δe ad ; ê ¼ ē.
The variances σ 2 e and σ 2 ead of the uncorrelated probability distributions, F, are additive, σ 2 ê ¼ σ 2 e þ σ 2 ead , and the joint distribution, G, factorizes when formulated in terms of the variables e and Δe ad : Gðe; e ad Þ ¼ F ead ðΔe ad Þ·F e ðê−Δe ad Þ: The linear correlation coefficient between two variables, z and y, is defined by ρ c ðz; yÞ ¼ zy−zy σ z σ y : We evaluate zy ¼ ∫∫zyGðz; yÞ dzdy ¼ ∫∫zy F y ðyÞF x ðz−yÞ dzdy ¼ ∫y F y ðyÞ ∫ðx þ yÞF x ðxÞ dx h i dy ¼ ∫yðx þ yÞ F y ðyÞ dy ¼ y 2 þ xy ; and therefore arrive at ρ c ðz; yÞ ¼ y 2 þ x y−ðx þ yÞy σ z σ y ¼ σ y σ z : This general result does not resort on a specific functional form of the joint distribution. Setting ê ¼ z , Δe ad =y, and e=x yields ρ c ¼ σ ead ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi σ 2 e þ σ 2 ead q : The limiting cases are ρ c →0 for σ ead →0 (equivalent to σ T →0) and ρ c →1 for σ e →0. B. K. thanks Henrike Wilms for performing Fourier analyses of temperature fluctuation time series and providing many valuable comments on the manuscript, and Joan Alexander, Thomas Birner, and Eric Jensen for helpful discussions. A. P. thanks Riwal Plougonven, Albert Hertzog, and Leonhard Pfister for fruitful discussions. Both authors are grateful to Thomas Peter and an anonymous reviewer for their careful reviews. The superpressure balloon data were collected in the frame of the Concordiasi project. Concordiasi was built by an international scientific group and is currently supported by the following agencies: Météo-France, CNES, IPEV, PNRA, CNRS/INSU, NSF, NCAR, the Concordia consortium, the University of Wyoming, and Purdue University. ECMWF also contributes to the project through computer resources and support, as well as scientific expertise. The two operational polar agencies, PNRA and IPEV, are thanked for their support at Concordia station. Concordiasi is part of the THORPEX-IPY cluster within the International Polar Year effort. The radiosonde data are available online (http://weather. uwyo.edu/upperair/sounding.html).