Introduction

Cirrus clouds, made of ice crystals and present at low temperatures (below 235 K) and high altitudes, cover approximately 17% of the Earth1. They control the hydrological balance of the upper troposphere and the longwave cloud feedback2, 3. Outside of areas of strong convection, cirrus form by ice nucleation on aerosol particles. Human activities (e.g., transportation and energy production) and natural phenomena (e.g., volcanic eruptions and dust storms) that modify the concentration of atmospheric aerosol influence the properties of cirrus4, impacting the hydrological balance of the atmosphere and ultimately Earth’s climate2, 5,6,7. For example, observational studies suggest that volcanic ash leads to the increase of ice-containing clouds in the upper troposphere6, 7. There is also evidence of aerosol particles affecting the dehydration efficiency of cirrus in the tropical tropopause, hence the concentration of water vapor entering the stratosphere8. Cirrus modification has been suggested as a geoengineering strategy to counterbalance greenhouse warming9. However, modeling studies show discrepancy on the predicted response of cloud properties to the emission of distinct aerosol species5, 10, 11, introducing uncertainty in the reliability of cirrus geoengineering on climate12.

The representation of cirrus in general circulation models (GCMs) remains challenging. Wind velocity, temperature, relative humidity and the physicochemical properties of the aerosol may vary at smaller scales (~1 to 20 km) than the typical resolution used in GCMs (~100–200 km)13, 14. Moreover, the ice-nucleating efficiency of aerosol depends on particle size and chemistry, and is not completely understood15. Ice formation is also realized by two different processes: homogeneous ice nucleation (HOM), i.e., the spontaneous freezing of supercooled liquid droplets (typically sulfuric acid and sulfate solutions), and, heterogeneous ice nucleation (HET), which requires the presence of ice nucleating particles (INP). Although sparse, INP can profoundly alter the evolution of clouds5. Much investigation has been devoted over the last decade to elucidating the nature of INP showing that a small fraction of atmospheric aerosol (typically some organics, mineral dust and black carbon) are capable of nucleating ice15. The HOM and HET mechanisms can occur simultaneously, and in fact they “compete” during cloud formation. HOM nucleation requires higher relative humidity (RH) than HET. Thus, by depleting the increase in RH from cooling driven by vertical motion, ice crystals formed by the HET mechanism may inhibit and even prevent HOM nucleation16. However, compared to HOM, HET nucleation leads to lower ice crystal concentration17, 18.

Vertical motion determines the maximum relative humidity in a cloudy parcel and drives ice nucleation19. 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 cirrus22,23,24. Modeling studies have also shown that low T and high W favors the HOM over the HET mechanism25,26,27,28,29, which is reflected in the global distribution of N i simulated in GCMs5, 11, 30. Field measurements however suggest predominance of HET and lower N i than implied by HOM10, 14, 24, consistent with sustained levels of supersaturation found in cirrus clouds23, 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 observations33,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, \({\rm{\Phi }}(\bar{W},{\sigma }_{w})\), largely determined by its standard deviation, σ w , since the mean resolved vertical transport is slow compared to cloud-scale motions13, 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 study37 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 studies24, 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 \({\rm{\Phi }}(\bar{W},{\sigma }_{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–200739,40,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 Experiments42. The G5NR is capable of resolving mesoscale systems and organized convection within large scale midlatitude cyclones, important for the simulation of cirrus40.

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 \({\rm{\Phi }}(\bar{W},{\sigma }_{w})\) at the 1° (~100 km) resolution for the two-year period of the simulation. Monthly averages were calculated by averaging \({\sigma }_{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 formation22. However, high frequency waves with periods of a few minutes are responsible for the formation of “pockets” of high N i within clouds14, 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, \({\sigma }_{w}^{2}\), is represented as the sum of the resolved and unresolved components of the high resolution run,

$${\sigma }_{w}^{2}={\sigma }_{w\mathrm{,7}{\rm{km}}}^{2}+{\sigma }_{w,{\rm{unres}}}^{2},$$
(1)

where \({\sigma }_{w\mathrm{,7}{\rm{km}}}^{2}\) is the spatial variance in W calculated from the 7 km resolution output, and \({\sigma }_{w,{\rm{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 \({\sigma }_{w,{\rm{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,

$$\frac{{W}_{{r}_{0}}}{{W}_{{r}_{1}}}={\alpha }_{\mathrm{1,0}}={(\frac{1+\frac{{r}_{1}}{{\rm{\Delta }}Z}}{1+\frac{{r}_{0}}{{\rm{\Delta }}Z}})}^{\mathrm{1/2}},$$
(2)

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 cloud14, 29, 44, 45.

For the special case of a normal distribution at resolution r 1,

$${W}_{{r}_{1}}\sim {\rm{\Phi }}\,(\bar{W},{\sigma }_{w\mathrm{,7}{\rm{km}}})$$
(3)

where \(\bar{W}\) is the grid-scale vertical velocity. Combining Eqs (2) and (3) it can be readily shown that,

$${W}_{{r}_{0}}\sim {\rm{\Phi }}\,({\alpha }_{\mathrm{1,0}}\bar{W},{\alpha }_{\mathrm{1,0}}{\sigma }_{w\mathrm{,7}{\rm{km}}}).$$
(4)

Which implies σ w  = α 1,0 σ w,7km. Equation (4) provides a simple argument to scale the resolved variance calculated directly from the G5NR output, \({\sigma }_{w\mathrm{,7}{\rm{km}}}^{2}\), to obtain the total variance driving cloud formation, \({\sigma }_{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 available13. 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 retrievals13 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 dBz13. 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 method13.

Figure 1
figure 1

Vertical velocity (m s−1) distribution within cirrus at the SGP site from cloud Doppler radar measurements13 (red lines) and high resolution simulations (blue bars). Measurements correspond to the period 2000–2010 and model results to 2005–2007. Scaling was applied to the model results as described in Methods. Black lines represent the “unscaled distribution” from the G5NR. Positive W values indicate updraft. The legends correspond the mean value of σ w for each case.

Figure 2
figure 2

Similar to Fig. 1 but for the Manus site.

To account for vertical motion with scales below 7 km, likely resulting from in situ turbulence and high-frequency gravity waves13, 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 \({\rm{\Phi }}(\bar{W},{\sigma }_{w})\). However Figs 1 and 2 suggest that \({\rm{\Phi }}(\bar{W},{\sigma }_{w})\) can be adequately approximated using Gaussian functions, in line with literature reports26, 29, 46,47,48.

Comparison of the measured and simulated \({\rm{\Phi }}(\bar{W},{\sigma }_{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 \({\rm{\Phi }}(\bar{W},{\sigma }_{w})\) compared to the measurements (Fig. 1, black lines). Scaling brings \({\rm{\Phi }}(\bar{W},{\sigma }_{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 \({\rm{\Phi }}(\bar{W},{\sigma }_{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 G5NR40.

Figure 3
figure 3

Monthly mean standard deviation in vertical velocity (m s−1) at the SGP and Manus sites calculated from 7 km global output (model) and from radar retrievals (obs). Error bars have been omitted for clarity however the standard error in the observations13 is about 0.15 m s−1.

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 elsewhere49. Briefly, HOM occurs on sulfate whereas mineral dust, black carbon and glassy organics are considered INP50. The ice nucleation rate is weighted by the distribution of vertical velocity and explicitly integrated. Large scale deposition, ice cloud fraction, and ice nucleation are coupled49. To minimize the effect of model uncertainties from the transport of aerosol and the meteorological conditions, the 100 km simulation is constrained using the horizontal wind velocity, temperature and water vapor from the second version of the Modern-Era Retrospective Analysis for Research and Applications (MERRA-2)39, 51, using a relaxation time scale of 6 h. Monthly-averaged σ w derived from the G5NR was used to drive ice nucleation. To perform ice nucleation studies, the model was run for five years (2005–2010) using a c90 cubed-sphere grid (spatial resolution of about 1°) and 72 vertical levels. Cloud physical and optical properties obtained in similar simulations have been shown to be in agreement with available in situ and satellite observations49.

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 studies17, 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/.

Results

Vertical velocity distribution

Figure 4 shows the annual global distribution of σ w derived from the G5NR. 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 troposphere46, 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 observations53 however differs from previous work where W, instead of σ w , was assumed to decrease with height54, 55. Figure 4 suggests that high values of W do occur at high altitude, but they become less frequent near the tropopause since \({\rm{\Phi }}(\bar{W},{\sigma }_{w})\) becomes narrow.

Figure 4
figure 4

Annual mean standard deviation in vertical velocity (m s−1) at a horizontal resolution of 1°, calculated from 7 km global output. Left panel: σ w at 250 hPa. Right panel: zonal mean σ w . Maps generated using the NCAR Command Language (Version 6.3.0) Software. (2016). Boulder, Colorado: UCAR/NCAR/CISL/TDD. http://dx.doi.org/10.5065/D6WD3XH5.

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 clouds1. 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 sources56. Although likely, they are progressively less frequent as σ w increases. Globally about one in 104 non-convective cirrus formation events are forced by vertical motion with σ w  > 0.8 m s−1. High values of σ w are more likely in the Tropics driven by underlying convection. Removing grid-cells with vigorous convection results in about a factor of five lower probability of finding σ w  > 0.2 m s−1. Values of σ w greater than 0.5 m s−1 seem to be equally likely in SH and NH indicating that they may result from strong orographic uplift in the mountain ranges.

Figure 5
figure 5

Probability distribution of σ w for different regions, and for grid cells with no underlying convection.

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 jets57. 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 \({\rm{\Phi }}(\bar{W},{\sigma }_{w})\) for the same location48, 57. A longer simulation period is however required to further study the interannual variability in \({\rm{\Phi }}(\bar{W},{\sigma }_{w})\) and will be the subject of a future study.

Figure 6
figure 6

Difference between the monthly and annual σ w at 100 km resolution calculated from the G5NR 7 km resolution output at the 250 hPa pressure level. Maps generated using the NCAR Command Language (Version 6.3.0) Software. (2016). Boulder, Colorado: UCAR/NCAR/CISL/TDD. http://dx.doi.org/10.5065/D6WD3XH5.

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 generation52. As similar parameterizations are used in many GCMs, this exercise shows that Fig. 4 may provide a way to evaluate their performance.

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 .

Figure 7
figure 7

Global distribution of the frequency of cirrus events dominated by homogeneous ice nucleation, vertically weighted by cloud fraction (a) and zonal mean (b). Maps generated using the NCAR Command Language (Version 6.3.0) Software. (2016). Boulder, Colorado: UCAR/NCAR/CISL/TDD. http://dx.doi.org/10.5065/D6WD3XH5.

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 NH58. 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.

Figure 8
figure 8

Monthly mean homogeneous ice nucleation frequency for the Tropical (latitude −30° to 30°) and the Northern (NH, latitude 30° to 60°), and Southern (SH, latitude −30° to −60°) extratropical regions.

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 reports58. This is in part by design as the heterogeneous ice nucleation spectrum used in the model is derived from a collection of available field data50. 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 mechanism5, 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 formation10, and, lower N i and higher saturation ratios than implied by HOM8, 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 temperatures59. Results for the NH, where most cirrus field campaigns have taken place23 show similar tendencies (see Supplementary Fig. S3) with slightly lower variability in N i and better agreement with observations for T > 230K.

Figure 9
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 observations23.

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 reports5, 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 temperature60, 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 tropopause14, 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 dynamics29, 33, 45 and an episodic nature of HOM in cirrus14.

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.