Direct estimation of the global distribution of vertical velocity within cirrus clouds

Cirrus clouds determine the radiative balance of the upper troposphere and the transport of water vapor across the tropopause. The representation of vertical wind velocity, W, in atmospheric models constitutes the largest source of uncertainty in the calculation of the cirrus formation rate. Using global atmospheric simulations with a spatial resolution of 7 km we obtain for the first time a direct estimate of the distribution of W at the scale relevant for cirrus formation, validated against long-term observations at two different ground sites. The standard deviation in W, σ w, varies widely over the globe with the highest values resulting from orographic uplift and convection, and the lowest occurring in the Arctic. Globally about 90% of the simulated σ w values are below 0.1 m s−1 and about one in 104 cloud formation events occur in environments with σ w > 0.8 m s−1. Combining our estimate with reanalysis products and an advanced cloud formation scheme results in lower homogeneous ice nucleation frequency than previously reported, and a decreasing average ice crystal concentration with decreasing temperature. These features are in agreement with observations and suggest that the correct parameterization of σ w is critical to simulate realistic cirrus properties.

Vertical motion determines the maximum relative humidity in a cloudy parcel and drives ice nucleation 19 . Early studies showed that introducing vertical velocity perturbations within parcel model simulations resulted in high variability in ice crystal concentration, N i 20,21 , which has been confirmed by aircraft observations within cirrus 22-24 . Modeling studies have also shown that low T and high W favors the HOM over the HET mechanism [25][26][27][28][29] , which is reflected in the global distribution of N i simulated in GCMs 5,11,30 . Field measurements however suggest predominance of HET and lower N i than implied by HOM 10,14,24 , consistent with sustained levels of supersaturation found in cirrus clouds 23,31,32 . This suggests that the frequency of HOM in cirrus is likely overestimated by GCMs. Recent studies point at overestimation in the parameterized cloud-scale W as the likely cause of the discrepancy between models and observations [33][34][35][36] .
Because of the separation between the relevant scale for ice nucleation and the scale resolved by the GCMs, it is likely that several cloud formation events occur within a model grid cell. This unresolved variability is characterized by a "subgrid" distribution of vertical velocity, σ Φ W ( , ) w , largely determined by its standard deviation, σ w , since the mean resolved vertical transport is slow compared to cloud-scale motions 13,14,29 . The latter are typically not resolved by CGMs due to their coarse resolution (~100 km). To represent subgrid W variability, most GCMs rely on either poorly constrained parameterizations or empirical correlations representing particular cloud realizations. A recent study 37 suggested that vertical wind velocity may be responsible for about 90% variation in calculated ice crystal formation rates, although it is not clear whether the same relation applies to other cloud properties. Field campaign analyses and cloud modeling studies 24,29,38 also suggest a strong relation between the effect of aerosol emissions on cloud properties and W. These highlight the importance of improving the representation of subgrid W variability in GCMs.
In this work we develop a method, using ultra high resolution global simulations to directly calculate the global distribution of subgrid vertical velocity affecting cirrus formation. By implementing our estimates in a global model and running experiments constrained by observations we show that the global distribution of σ w largely determines the balance between homogeneous and heterogeneous ice nucleation during the formation of cirrus.

Methods
Vertical velocity distribution. We developed a method to calculate σ Φ W ( , ) w using results from a global climate simulation carried out using the non-hydrostatic version of the NASA Goddard Earth Observing System (GEOS-5) over the period 2005-2007 39-41 . The simulation was completed using a cubed-sphere grid with a nominal spatial resolution of 7 km and 72 vertical levels, extending from the surface to 0.1 mbar. The time step for integration was set to five minutes and three-dimensional instantaneous meteorological fields were saved every hour. The vertical resolution of the model in the upper troposphere is about 0.5 km. This simulation is referred as the GEOS-5 "nature run" (G5NR). The output of the simulation amounts to about four petabytes and is used to perform Observation System Simulation Experiments 42 . The G5NR is capable of resolving mesoscale systems and organized convection within large scale midlatitude cyclones, important for the simulation of cirrus 40 .
The global distribution of vertical velocity was calculated collocating W from the 7 km G5NR output to 1° (~100 km) horizontal resolution (except for the validation studies, for which 0.5° sections were used) so that at least 256 W values from the G5NR were used to calculate σ w for each 1° grid cell. This procedure resulted in a hourly, three-dimensional global characterization of σ Φ W ( , ) w at the 1° (~100 km) resolution for the two-year period of the simulation. Monthly averages were calculated by averaging σ w 2 using hourly output. Sensitivity tests were carried calculating σ w for 200 km global resolution, producing essentially similar results as our 100 km calculation, indicating that very little of the total W variance is resolved at scales greater than 100 km. A second sensitivity test was performed collocating results from a short term 3.5 km resolution run spanning 15 days in May 2006 to the 1° spatial resolution. After scaling was applied (α 1,0 = 1.26, as explained below, Eq. 4), this test produced similar σ w as when the 7 km simulation was used (Supplemental Fig. S1).
The G5NR resolves waves and vertical motion with periods from minutes to a few hours which are the main drivers of cirrus formation 22 . However, high frequency waves with periods of a few minutes are responsible for the formation of "pockets" of high N i within clouds 14,23 . The scale of such waves would likely be smaller than the 7 km resolution of the G5NR. To account for such a unresolved variability a "scaling" method is developed, as follows. The total spatial vertical velocity variance at the 100 km resolution, σ w 2 , is represented as the sum of the resolved and unresolved components of the high resolution run, ,unres 2 where σ w ,7km 2 is the spatial variance in W calculated from the 7 km resolution output, and σ w ,unres 2 the variance resulting from vertical motion with characteristic scales below 7 km. The contribution to the W spatial variability from motion with scales greater than 100 km is neglected. This is justified because motion with scales greater than 100 km is resolved at the low resolution. To estimate σ w ,unres 2 , the approach of Pauluis, et al. 43 is employed. Using a discretized version of the equations of motion of an anelastic fluid, the authors derived a relation for the vertical velocity resolved at two different horizontal resolutions, e.g., r 0 and r 1 , where ΔZ = 6 km corresponds to the resolution at which the increase in kinetic energy from buoyancy is equally distributed between its horizontal and vertical components. Pauluis, et al. 43 showed that Eq. (2) is accurate for horizontal resolutions between 2 and 20 km, hence it is suitable to scale the G5NR, 7 km results, to finer resolutions, hence to estimate the variance unresolved by the G5NR. In this study we assume r 0 = 0.1 km as the horizontal scale below which the cloud properties can be considered uniform. Using r 1 = 7 km (i.e., the resolution of the G5NR) results in α 1,0 = 1.46. For r 1 = 3.5 km, α 1,0 = 1.26. The choice of r 0 is rather arbitrary, however for small enough values has little effect on α 1,0 (e.g., α 1,0 = 1.41 for r 0 = 0.5 km). Here r 0 is selected much smaller than the typical cloud size to account for short-lived fluctuations that may affect the relaxation of supersaturation in the cloud 14,29,44,45 .
For the special case of a normal distribution at resolution r 1 , where W is the grid-scale vertical velocity. Combining Eqs (2) and (3) it can be readily shown that, Which implies σ w = α 1,0 σ w,7km . Equation (4) provides a simple argument to scale the resolved variance calculated directly from the G5NR output, σ w ,7km 2 , to obtain the total variance driving cloud formation, σ w 2 , at the 100 km resolution.
Validation. The vertical velocity fields from the 7 km simulation were validated against ground-based measurements within cirrus for two different sites with diverse topography and synoptic conditions, and for which long term radar retrievals are available 13 . These correspond to the Southern Great Plains of North America (SGP, 36° 36′ 18″ N, 97° 29′ 6″ W) and the Manus island in the tropical western Pacific (Manus, 2° 3′ 39.64″ S, 147° 25′ 31.43″ E). SGP is a mid-latitude continental site with large variability in temperature and cloud occurrence. Manus is an oceanic site off the coast of Australia with frequent tropical convection. Ground-based radar retrievals 13 of vertical velocity at each site over the period (2000-2010) were obtained from the Atmospheric Radiation Measurement program (www.arm.gov/sites). The retrieval algorithm is based on a decomposition of the Doppler vertical velocity. The typical error in vertical velocity is about 15 cm s −1 for a minimum reflectivity of about −40 dBz 13 . The whole data set spans about 10 years for each site. Data obtained at 10 s intervals for each month were averaged over 5 min to match the time step of the G5NR simulation.
To generate G5NR vertical velocity distributions at the SGP ( Fig. 1) and Manus ( Fig. 2) sites, instantaneous W values over a 0.5° × 0.5° area centered at each site were obtained from the 7 km simulation at three hour intervals. This corresponds to about 15360 values for each monthly distribution at each site. Notice that this includes the spatial and temporal components of the variance since it is difficult to separate them in the observations. The model results were restricted to ice mixing ratios above 5 × 10 −5 kg kg −1 and corresponding ice water content of about 50 mg m −3 , selected to match the maximum sensitivity of the retrieval method 13 . To account for vertical motion with scales below 7 km, likely resulting from in situ turbulence and high-frequency gravity waves 13 , the distributions were scaled to 0.1 km horizontal resolution using the method outlined in the previous section. Equation (2) was applied directly to the simulated distributions, i.e., without invoking the assumption of a normal W ( , ) . However Figs 1 and 2 suggest that σ Φ W ( , ) w can be adequately approximated using Gaussian functions, in line with literature reports 26,29,[46][47][48] .
Comparison of the measured and simulated W ( , ) w σ Φ (Figs 1 and 2) shows that the G5NR is capable of generating realistic vertical velocity distributions, and in reasonable agreement with observations. At the SGP site it is evident that using only the raw G5NR data would result in a too-narrow σ Φ W ( , ) w compared to the measurements (Fig. 1, black lines). Scaling brings σ Φ W ( , ) w closer to the observations suggesting that a significant fraction of the W variance lies in the sub-7 km range. On the other hand, for the Manus site the raw G5NR distribution (Fig. 2, black lines) approximates the observed distribution (Fig. 2, red lines), indicating that most of the W variance is resolved at the 7 km resolution. In the latter, scaling may lead to overestimation in σ w since Eq. (2) implicitly assumes that a significant fraction of the W variance lies at the small scales. Thus, σ w can be considered a upper limit to the vertical velocity variance.
Seasonal variation in the large scale environment may lead to differences between the simulated and observed distributions. Such deviations are typically within the margin of error of the observations. However they can also signal systematic errors in the simulated W ( , ) . To study the latter, σ w was calculated for each month over the whole observational record (resulting in ten data points per month) and from the G5NR (which are available for two separate years). The results are plotted in Fig. 3. At both sites, the simulated and observed σ w show little interannual variability, i.e., for the same month the spread in σ w between years is typically below 0.05 m s −1 . σ w at the SGP site shows a strong annual cycle, whereas at the Manus site it is relatively constant over the year. This suggests that location plays an important role in determining σ w . The distinctive behavior of σ w at each site is well represented by the simulation. However the G5NR tends to predict a stronger annual cycle of σ w at the SGP site than implied by observations, with the maximum σ w occurring too early during the year and an underestimation in σ w between August and October, likely related to the low frequency of convective events during the Fall season in the G5NR 40 .
GCM implementation. The calculated σ w was used to drive ice nucleation in GEOS-5 and study the impact on the balance between HOM and HET processes. Due to the high computational expense of the 7 km simulations, aerosol-cloud interactions are analyzed in the lower resolution, 100 km, simulation, but using σ w obtained from the 7 km run. The goal of the GCM implementation is to analyze the statistics of the ice crystal concentration when the G5NR-generated distribution of σ w is used. The main premises behind this approach are that the interannual variability in σ w is small and that σ w is mostly influenced by local features and convection. This is supported by the low interannual variability in σ w found in our 7 km simulation, and in the 10-year-long time series of radar retrievals at the SGP and Manus sites (Fig. 3).
The model's cloud microphysics and ice nucleation schemes are described elsewhere 49 . Briefly, HOM occurs on sulfate whereas mineral dust, black carbon and glassy organics are considered INP 50 . The ice nucleation rate is A cirrus formation event was considered dominated by homogeneous freezing when at least 80% of the ice crystals were produced by the HOM mechanism. Modeling studies 17,18 show that cloud formation becomes HET dominated within a relative narrow range of INP concentration, typically less than a factor of two. Thus the 80% limit represents the INP concentration at which N i becomes strongly affected by HET nucleation. Results using different thresholds, i.e., 20%, 50% and 90% are discussed below. Bivariate N i vs. T distributions were calculated by counting the number of data points falling within a particular N i and T combination, using 80 logarithmic bins for N i and 80 linear bins for T. The frequencies were then normalized by the most frequent N i vs. T combination within the entire domain. Using this method the expected value of N i is located around the maximum frequency at each temperature. Data Availability. The GEOS-5 source code is available under the NASA Open Source Agreement http:// opensource.gsfc.nasa.gov/projects/GEOS-5/. Figure 4 shows the annual global distribution of σ w derived from the G5NR.

Vertical velocity distribution.
High values of σ w are found around the Inter-Tropical Convergence Zone (ITCZ) and over continental mountain ranges, confirming the notion that underlying convection and orography are the main drivers of dynamic variability in the upper troposphere 46,49,52 . For the same reason σ w is typically higher in the Subtropical Northern hemisphere (NH) than in the Southern hemisphere (SH), except over the Andes mountains where σ w tends to be high due to orographic uplift. Within the troposphere σ w decreases slightly with height, and becomes small above the tropopause due to gravity wave breaking. This is in agreement with observations 53 however differs from previous work where W, instead of σ w , was assumed to decrease with height 54,55 . Figure 4 suggests that high values of W do occur at high altitude, but they become less frequent near the tropopause since W ( , ) Aside from the poles, subtropical regions, usually associated with persistent stratocumulus decks and large scale subsidence (i.e., the coasts of Peru, California and Namibia), display the smallest σ w over the globe. These areas are also associated with a low frequency of cirrus clouds 1 . Low σ w results in low N i and large ice settling velocities, which may contribute to a more efficient cloud removal and explain in part the low cirrus frequency of these regions. At constant pressure σ w tends to decrease poleward from the Tropics in both hemispheres due in part to the breaking of gravity waves at higher pressure levels. In NH the lowest σ w values are found at the northernmost latitudes (above 80°). In SH the lowest values of σ w are found around −60°, however σ w increases south of −60° revealing the orographic effect of the Antarctic continent on W, which may impact the formation of polar stratospheric clouds. Figure 5 shows the global PDF of σ w , P(σ w ), within cirrus, i.e., positive ice water content. As expected, P(σ w ) decreases exponentially with σ w . However several features of P(σ w ) stand out. P(σ w ) peaks at around 0.02 m s −1 and decreases rapidly with increasing σ w so that about 90% of the values of σ w are below 0.1 m s −1 . This suggests that in most cases cirrus formation is driven by large scale vertical transport and inertial gravity waves. Higher values of σ w (up to 0.8 m s −1 ) are associated with high frequency gravity waves from convective and orographic sources 56  There is a slight seasonal variation in the global distribution of σ w . Over the Tropics the highest σ w coincides with the displacement of the ITCZ (Fig. 6). Over the continents and away from orographic features, σ w is influenced by large scale meteorological patterns, i.e., cold fronts, midlatitude cyclones, and subtropical jets 57 . In the subtropics σ w tends to be larger during the summer months, particularly over land. The two years of the G5NR simulation show similar patterns: highest in the ITCZ and in orographic regions and lowest in the high latitudes with little variation in the Subtropical regions. Location seems to be the main factor determining σ w . This is supported by the data at the Manus and SGP sites (Fig. 3), and it is in line with literature reports showing low interannual variability in σ Φ W ( , ) w for the same location 48,57 . A longer simulation period is however required to further study the interannual variability in σ Φ W ( , ) w and will be the subject of a future study. The results shown in Fig. 4 can be used to evaluate the parameterization of σ w currently used in GEOS-5 (Barahona et al. 49 , c.f. Fig. 4). Compared to the parameterized results, σ w is lower in the marine midlatitudes and higher in the Tropics. This parameterization used the orographic stress and turbulence to estimate σ w . Over the ocean the total surface stress was used, leading to an overestimation of σ w in the midlatitudes. On the other hand, the parameterization did not explicitly account for gravity waves generated by convection leading to underestimation of σ w in the Tropical region. Over land the agreement is better since orographic features drive gravity wave   Implications for ice cloud formation. The annual mean frequency of homogeneous freezing events during cirrus formation calculated with GEOS-5 (forcing ice nucleation with σ w calculated from the G5NR output) is shown in Fig. 7. High σ w and low T tend to lead to high annual HOM frequency. This is because under such conditions a larger INP concentration is required to offset the increase in RH from expansion cooling, and is evidenced by the correspondence in the spatial patterns seen in Figs 4 and 7. In the Tropical region where high values of σ w are more likely, HOM dominates about 50-60% of the cloud formation events, and up to 80% in the coldest regions of the Tropical tropopause, over the Central Pacific, where a lack of INP also contributes to diminish the frequency of HET nucleation. Outside the Tropical and Subtropical regions, the frequency of cloud events dominated by HOM nucleation is generally low (below 30%), particularly in the NH. In fact, in the Arctic clouds tend to form almost exclusively by HET (HOM frequency is below 10%), mostly resulting from the low σ w which precludes efficient HOM nucleation, even though the INP availability is low. HOM-dominated cirrus events are also less frequent over North Africa due to the presence of dust and black carbon, and to the absence of orographic features and convection that would produce high σ w .
While the global mean HOM-dominated frequency is relatively constant over the year at ~37%, in both the SH and NH it peaks in the winter months, since low temperature favors HOM (Fig. 8). Maximum HOM occurrence in the SH is 45% during the winter months, while it is only 30% in the NH winter. Even though σ w is higher on average in the NH than in the SH, HOM is more prevalent in the latter during most of the year since INP like mineral dust are abundant in the NH 58 . The seasonal differences in the SH are more pronounced than in the NH due to the larger temperature and aerosol variability, and lower σ w in the former. In the Tropics, a seasonal cycle is also present reflecting the strength and position of the ITCZ influencing σ w and the annual variation in black carbon and dust concentration.
The balance between HOM and HET during cloud formation is significantly influenced by the concentration of ice nucleating particles, N INP . Supplementary Figure S2 shows the frequency of N INP calculated at the maximum RH during cloud formation, RH max . N INP increases steeply around RH max = 15% to a maximum global average value of about 3 L −1 at RH max = 30% (~5 L −1 for NH). The apparent decrease in N INP at high RH max is caused by competition between HOM and HET nucleation. For RH max < 40% where HOM does not occur, the simulated N INP shows similar characteristics as available reports 58 . This is in part by design as the heterogeneous ice nucleation spectrum used in the model is derived from a collection of available field data 50 . Figure S2 shows that GEOS-5 simulates realistic N INP statistics.
Globally, about 70% of the time cirrus form in situations where HOM and HET occur simultaneously, with HET being more prevalent. This finding is at odds with previous modeling studies where HOM is the predominant cirrus formation mechanism 5,49,54 . Our results are however in better agreement with field campaign data suggesting a significant role of dust and other INP species in cirrus formation 10 , and, lower N i and higher saturation ratios than implied by HOM 8,14 . This suggests that the parameterization of σ w may be responsible for the higher HOM frequency typically found in modeling studies. Notably, forcing ice nucleation with our estimate of σ w also results in good agreement of N i with field campaign data at very low temperature (Fig. 9, T < 200 K), where most modeling studies show high overestimation of N i 5, 36 . This suggests that the overestimation in many models may be a result of poorly constrained σ w . GEOS-5 however seems to overestimate the frequency of low N i for T > 230 K. It must be noticed that ice shattering on in-situ probes, leading to overestimation in the in-situ N i , are a likely artifact of the measurements at such temperatures 59 . Results for the NH, where most cirrus field campaigns have taken place 23 show similar tendencies (see Supplementary Fig. S3) with slightly lower variability in N i and better agreement with observations for T > 230K. The effect of selecting different thresholds to define a HOM-dominated cirrus is shown in the Supplementary  Fig. S4. Changing the definition of a HOM-dominated cloud from 90% to 50% of N i produced by HOM, increases the HOM frequency from 33% to 43%. The latter corresponds to the minimum threshold where a cloud could be referred as HOM-dominated since for any value below 50% most of the ice crystals would in fact be produced by HET nucleation. Even at the 50% threshold, the HOM frequency is still much lower than reported values showing that our conclusions are not significantly influenced by the selected threshold. Using a more strict definition, where HOM still occurs but is not dominant, i.e., 20% of ice crystals produced by HOM (actually a HET-dominated cloud) leads to global HOM frequency of about 48% which is still low compared to current reports 5, 49 .

Discussion and Conclusions
This work reports for the first time the direct estimation of the subgrid spatial variance in vertical velocity at the scale relevant for cirrus formation. Figure 4 shows that σ w has considerable spatial and seasonal variability. It is important for GCM parameterizations of σ w to reproduce such a spatial variability. This indicates that parameterizations based on individual field campaigns should be used with caution when applied to the global scale. Large-scale dynamics, turbulence, orographic uplift, and underlying convection impact σ w . Even State-of-the-Art parameterizations of σ w 49, 52 neglect the effect of convection generating gravity waves that can impact cirrus formation, which may result in a too-low frequency of high σ w as indicated by Fig. 5.
Using the direct estimate of σ w to drive the GEOS-5 ice nucleation scheme results in a lower predominance of homogeneous ice nucleation than previously simulated, in better agreement with field measurements. Such lower frequency of HOM events also results in a better agreement of simulated N i with observations, particularly at low T (Fig. 9). This is significant as cloud formation theory typically predicts high N i for low temperature 60, 61 , a trend also found in GCM studies. One way to reconcile theory and observations is to assume very low W at low T 30,36,54 , which however conflicts with observations of high W (>0.5 m s −1 ) near the tropopause 14, 23, 62 . Our results provide a way to reconcile these two views. High values of W do occur at low T, however the low σ w at high altitude limits the frequency of cloud formation events driven by high W, and on average leads to lower N i . This is agreement with previous work suggesting the structuring of cirrus by dynamics 29,33,45 and an episodic nature of HOM in cirrus 14 .
Our new estimate provides a way to validate new parameterizations of σ w at a global scale. Several uncertainties however remain in the modeling of cirrus clouds, the most significant being the concentration of INP in the upper troposphere. Even though our results are in relative good agreement with available reports (Supplementary Fig. S2), still few studies provide observations of N INP with a wide range of aerosol concentrations, temperature and relative humidity and enough spatial coverage that would allow comprehensive validation of GCM predictions. The characterization of W at low T also requires better techniques with smaller errors, since in some cases they can be as high as σ w 13 . In turn an improved representation of σ w in GCMs may help reducing the uncertainty surrounding the estimation of the impact of aerosol emissions on cirrus, and eventually lead to a better prediction of future climate. Figure 9. Global frequency distribution of in-cloud ice crystal number concentration as a function of temperature from GEOS-5 output over a 2-year subset (2005)(2006). Solid lines represent the 25% and 75% quantiles from a compilation of field campaign observations 23 .