Effects of thermodynamics, dynamics and aerosols on cirrus clouds based on in situ observations and NCAR CAM6

Cirrus cloud radiative effects are largely affected by ice microphysical properties, including ice water content (IWC), ice crystal number concentration (N_i) and mean diameter (D_i). These characteristics vary significantly due to thermodynamic, dynamical and aerosol conditions. In this work, a global-scale observation dataset is used to examine regional variations of cirrus cloud microphysical properties, as well as several key controlling factors, i.e., temperature, relative humidity with respect to ice (RHi), vertical velocity (w) and aerosol number concentrations (N_a). Results are compared with simulations from the National Center for Atmospheric Research (NCAR) Community Atmosphere Model version 6 (CAM6). Observed and simulated ice mass and number concentrations are constrained to ≥62.5 µm to reduce potential uncertainty from shattered ice in data collection. The differences between simulations and observations are found to vary with latitude and temperature. Comparing with averaged observations at ∼100 km horizontal scale, simulations are found to underestimate (overestimate) IWC by a factor of 3–10 in the Northern (Southern) Hemisphere. Simulated N_i is overestimated in most regions except the Northern Hemisphere midlatitudes. Simulated D_i is underestimated by a factor of 2, especially for warmer conditions (−50 to −40 ^∘C), possibly due to misrepresentation of ice particle growth/sedimentation. For RHi effects, the frequency and magnitude of ice supersaturation are underestimated in simulations for clear-sky conditions. The simulated IWC and N_i show bimodal distributions with maximum values at 100 % and 80 % RHi, differing from the unimodal distributions that peak at 100 % in the observations. For w effects, both observations and simulations show variances of w (σ_w) decreasing from the tropics to polar regions, but simulations show much higher σ_w for the in-cloud condition than the clear-sky condition. Compared with observations, simulations show weaker aerosol indirect effects with a smaller increase of IWC and D_i at higher N_a. These findings provide an observation-based guideline for improving simulated ice microphysical properties and their relationships with key controlling factors at various geographical locations.

1 Introduction

Cirrus clouds represent one of the most ubiquitous cloud types with an estimated global coverage of approximately 20 % to 40 % (Mace and Wrenn, 2013; Sassen et al., 2008). According to the fifth assessment report of the United Nations Intergovernmental Panel on Climate Change (IPCC) (Boucher et al., 2013), the largest uncertainty in estimating future climate change stems from clouds and aerosols. Unlike most other cloud types, cirrus clouds may produce a net positive or negative radiative forcing depending on their microphysical properties (Stephens and Webster, 1981; Zhang et al., 1999), which are affected by meteorological conditions and aerosol distributions. Tan et al. (2016) showed that the radiative effects of misrepresenting the prerequisite condition of cirrus clouds – ice supersaturation (ISS, where relative humidity with respect to ice (RHi) >100 %) – can lead to an average bias of +2.49 W m⁻² at the top of the atmosphere. Other modeling studies found large differences in the net cloud radiative forcing depending on the fraction of activated ice-nucleating particles (INPs) and the nucleation mechanisms (i.e., homogeneous and heterogeneous nucleation) through which the clouds form (Liu et al., 2012; Storelvmo and Herger, 2014). The large uncertainties in cirrus cloud radiative forcing illustrate the need for further study of cirrus cloud microphysical properties as well as their controlling factors in various geographical locations.

Ideally, a comprehensive quantification of cirrus cloud microphysical properties globally based on high-resolution, in situ observations would mitigate many uncertainties. However, challenges remain in field measurements to achieve such spatial coverage. Previously, efforts have been made to understand cirrus cloud properties based on their geographical locations. Diao et al. (2014b) performed a hemispheric comparison of in situ cirrus evolution and found little difference in the clear-sky ISS frequency as well as the proportion of each evolution phase between the Northern Hemisphere and Southern Hemisphere (NH and SH, respectively). In situ observations of tropical, midlatitude and polar cirrus clouds have shown that ice water content (IWC) can vary orders of magnitude depending on the geographical locations (Heymsfield, 1977; Heymsfield et al., 2005, 2017; Mcfarquhar and Heymsfield, 1997; Schiller et al., 2008). Wolf et al. (2018) used balloon-based in situ observations to analyze microphysical properties of Arctic ice clouds and found differences in particle size distributions (PSDs) depending on the cloud origin. Krämer et al. (2016, 2020) developed a cirrus cloud climatology, focusing on tropical and midlatitude cirrus clouds, and showed that cloud thickness is larger at lower altitudes and thus produces a more negative radiative forcing. Moving from north to south using lidar-based observations from two research cruises starting from Leipzig, Germany, one to Punta Arenas, Chile, and the other one to Stellenbosch, South Africa, Kanitz et al. (2011) observed a decrease in the efficiency of heterogeneous nucleation in the SH, which could be a result of fewer INPs. This hemispheric difference in aerosol indirect effects is consistent with significantly higher aerosol number concentrations in the NH (Minikin et al., 2003). Using satellite observations from the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO), Mitchell et al. (2018) showed the dependence of ice particle effective diameter on temperature, latitude, season and topography. Thorsen et al. (2013) used CALIPSO data to examine cloud fraction of tropical cirrus clouds and showed dependence on altitude and diurnal cycle. Tseng and Fu (2017) used CALIPSO and Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) data and found that the tropical cold-point tropopause temperature is a controlling factor of cirrus cloud fraction in the tropical tropopause layer.

Regional and hemispheric variations of cirrus microphysical properties are produced by various controlling factors, such as thermodynamics (i.e., temperature and RHi), dynamics (e.g., vertical velocity) and aerosols (e.g., number concentration and composition). The effects of temperature have been extensively studied from in situ observations (Heymsfield et al., 2017; Luebke et al., 2013, 2016; Schiller et al., 2008), showing an increase of IWC towards warmer temperatures. A number of studies focused on distributions of RHi have found that in-cloud RHi occurs most frequently at or near 100 % (Jensen et al., 2001; Krämer et al., 2009). Another study by Diao et al. (2017) found that using different RHi thresholds (e.g., 108 % to 130 %) for ice nucleation in simulations can influence IWC and ice crystal number concentrations (N_i) in convective cirrus. In addition, the spatial scales of ice-supersaturated regions can vary from the micro- to mesoscales, largely depending on the spatial variability of water vapor (Diao et al., 2014a). The distributions of vertical velocity have been investigated in different types of cirrus clouds, such as in ridge-crest cirrus, frontal cirrus and anvil cirrus (Muhlbauer et al., 2014a, b). Stronger updrafts are found to be associated with higher occurrence frequency of ISS inside anvil and convective cirrus (D'Alessandro et al., 2017). Regarding the effects of aerosols, Cziczo et al. (2013) and Cziczo and Froyd (2014) investigated ice crystal residuals from in situ observations and discovered that the majority of midlatitude cirrus clouds form via heterogeneous nucleation on mineral dust and metallic particles. Anthropogenic aerosols, such as secondary organic aerosols, were found to be less effective INPs compared with mineral dust (Prenni et al., 2009). Based on remote sensing data, Zhao et al. (2018, 2019) showed that the correlations between ice crystal sizes and aerosol optical depth can be either positive or negative depending on the meteorological conditions in convective clouds. Chylek et al. (2006) showed an increase in ice crystal size during the more polluted winter months compared with cleaner summer months over the eastern Indian Ocean, which the authors speculate to be due to heterogeneous nucleation occurring at lower ice supersaturation compared with homogeneous nucleation, therefore reducing the ambient ice supersaturation magnitude and making homogeneous nucleation a more difficult pathway. Using a global-scale dataset of multiple flight campaigns, Patnaude and Diao (2020) isolated individual effects on cirrus clouds from temperature, RHi, vertical velocity (w) and aerosol number concentrations (N_a). They found that when N_a is 3–10 times higher than average conditions, it shows strong positive correlations with cirrus microphysical properties such as IWC, N_i and number-weighted mean diameter (D_i). These aerosol indirect effects are also susceptible to whether or not thermodynamic and dynamical conditions are controlled, demonstrating the importance of conducting a comprehensive analysis of various key controlling factors altogether.

More recently, in situ observations have been used to evaluate and improve cirrus cloud parameterizations in global climate models (GCMs). Two types of simulations have been frequently used for model evaluation, i.e., free-running (Eidhammer et al., 2014, 2017; Wang and Penner, 2010; Zhang et al., 2013) and nudged (D'Alessandro et al., 2019; Kooperman et al., 2012; Wu et al., 2017) simulations. For free-running simulations, a comparison on statistical distributions of ice microphysical properties is often used for model validation (e.g., Penner et al., 2009). The nudged simulation would nudge certain meteorological conditions towards reanalysis data, such as horizontal wind and temperature (e.g., D'Alessandro et al., 2019; Wu et al., 2017). These nudged simulations can also be output to a similar location and time to those of the aircraft observations. Given the importance and limited understanding of how aerosols interact with cirrus clouds, much attention has been dedicated to the parameterization of aerosol indirect effects (Kärcher and Lohmann, 2002, 2003; Kuebbeler et al., 2014; Wang et al., 2014a). Shi et al. (2015) added the effects of pre-existing ice into the Community Atmosphere Model version 5 (CAM5) and found a decrease in N_i with increasing aerosol concentration due to the reduction of homogeneous nucleation frequency. Other studies also investigated the effect of updraft velocity on simulated N_i and aerosol indirect effects (Zhou et al., 2016; Penner et al., 2018).

This study aims to bridge the knowledge gap on how cirrus clouds vary depending on geographical locations and environmental conditions by using a comprehensive in situ observation dataset that includes seven US National Science Foundation (NSF) flight campaigns. Observations were collected aboard the NSF/National Center for Atmospheric Research (NCAR) Gulfstream V (GV) research aircraft. Descriptions of the seven flight campaigns, instrumentations, model configurations of the NCAR Community Atmosphere Model version 6 (CAM6) are provided in Sect. 2. Both observations and simulations are used to examine the regional variations in the statistical distributions of cirrus microphysical properties, including IWC, N_i and D_i (Sect. 3). Impacts of several key controlling factors, i.e., temperature, RHi, w and N_a, are examined in Sect. 4. Discussions on observation-based findings and model evaluation results are included in Sect. 5.

Table 1 Descriptions of seven NSF flight campaigns conducted with the NSF/NCAR GV research aircraft.

* N, Northern Hemisphere; S, Southern Hemisphere; T, tropics; M, midlatitudes; P, polar regions; A, all regions.
^a,b UCAR/NCAR (2009, 2019a). ^c−f UCAR/NCAR (2019b–e). ^g UCAR/NCAR (2019f). ^h UCAR/NCAR (2018a). ⁱ UCAR/NCAR (2019g). ^j UCAR/NCAR (2018b). ^k UCAR/NCAR (2018c). Full citations of each dataset are included in the reference list.

2 Data and methods

2.1 In situ observations and instrumentations

In this study, in situ airborne observations at 1 Hz are provided by instruments aboard the NSF High-Performance Instrumented Airborne Platform for Environmental Research (HIAPER) GV research aircraft. A comprehensive global dataset is compiled based on seven major flight campaigns funded by the NSF, including START08 (Pan et al., 2010), HIPPO deployments 2–5 (Wofsy, 2011), PREDICT (Montgomery et al., 2012), TORERO (Volkamer et al., 2015), DC3 (Barth et al., 2015), CONTRAST (Pan et al., 2017) and ORCAS (Stephens et al., 2018). Table 1 provides a detailed summary of the seven flight campaigns, including location, duration of flights, total flight hours of all temperatures and flight hours for in-cloud and clear-sky conditions at temperatures $\le -\mathrm{40}$ ^∘C only. Maps comparing the flight tracks of in situ observations and the collocated CAM6 nudged simulations (hereafter named “CAM6-nudg” data) are shown in Fig. 1.

Figure 1 Global maps of flight tracks representing the seven campaigns in this study for (a) in situ observations and (b) CAM6-nudg. Colors denote different campaigns.

For this study, ice particle measurements are provided by the Fast 2-Dimensional Cloud particle imaging probe (Fast-2DC) with a 64-diode laser array for a range of 25–1600 µm. Larger particles can be reconstructed up to 3200 µm. The mass–dimensional relationship of Brown and Francis (1995) is used to calculate IWC for the Fast-2DC probe, which was previously used in other studies of the Fast-2DC probe aboard the NSF GV aircraft (Diao et al., 2014a, b, 2015). Number-weighted mean diameter (D_i) is calculated by summing up the size of particles in each bin using the bin center and then dividing it by the total number of particles. In order to mitigate the shattering effect, particles with diameters <62.5 µm (i.e., first two bins) are excluded in the Fast-2DC measurements when calculating IWC, N_i and D_i. The Rosemount temperature probe was used for temperature measurements, which has an accuracy and precision of $\sim \pm \mathrm{0.3}$ and 0.01 K, respectively. All analyses are restricted to temperatures $\le -\mathrm{40}$ ^∘C, in order to exclude the presence of supercooled liquid droplets in this study. Laboratory calibrated and quality-controlled water vapor data were collected using the vertical-cavity surface-emitting laser (VCSEL) hygrometer (Zondlo et al., 2010), with an accuracy of ∼6 % and precision of ≤1 %. Both temperature and water vapor are used at 1 Hz resolution for this analysis. Aerosol measurements were collected from the Ultra-High-Sensitivity Aerosol Spectrometer (UHSAS), which uses 100 logarithmically spaced bins ranging from 0.06–1 µm. RHi is calculated using saturation vapor pressure with respect to ice from Murphy and Koop (2005). The combined RHi uncertainties from the measurements of temperature and water vapor range from 6.9 % at −40 ^∘C to 7.8 % at −78 ^∘C. Measurements are separated by cloud condition where in-cloud condition is defined by the presence of at least one ice crystal from the Fast 2-DC probe (N_i>0 L⁻¹). The same in-cloud definition has been used by several previous studies (D'Alessandro et al., 2017; Diao et al., 2014a, b, 2015, 2017; Tan et al., 2016), and all other samples are defined as clear sky. For regional variation analysis, data are binned by six latitudinal regions in the two hemispheres, that is, NH polar (60–90^∘ N), SH polar (60–90^∘ S), NH midlatitudes (30–60^∘ N), SH midlatitudes (30–60^∘ S), NH tropics (0–30^∘ N) and SH tropics (0–30^∘ S). The majority of observations in the SH midlatitude and tropical regions are located over the oceans, while the observations of NH midlatitude and polar regions are predominantly over land.

Figure 2 (a) Vertical profiles of temperature, (b) potential temperature vs. temperature and (c) vertical profiles of potential temperature based on in situ observations at temperatures $\le -\mathrm{40}$ ^∘C. Number of samples (N) for 1 Hz observations is shown in the figure legend. Colors denote six latitudinal regions.

The vertical profiles of observed in-cloud temperature, clear-sky potential temperature (Θ) and their correlations are shown in Fig. 2. The observations sampled temperatures from −78 to −40 ^∘C and altitudes from 5–15 km, while a previous study of Krämer et al. (2020) sampled −91 to −30 ^∘C and 5–19 km (their Fig. 2). The lowest temperatures are found in the tropical regions and at the highest altitudes, whereas polar regions show more observations at lower altitudes that satisfy temperature $\le -\mathrm{40}$ ^∘C. Distributions of cirrus cloud properties (i.e., IWC, N_i, D_i), in-cloud and clear-sky RHi, and clear-sky water vapor mixing ratio for the observation dataset are shown in Fig. 3. D_i increases with decreasing altitudes, IWC slightly increases with decreasing altitudes, and N_i is almost independent of altitudes. Clear-sky RHi and water vapor mixing ratio both increase with decreasing altitudes, while in-cloud RHi is centered around 100 % and shows smaller dependency on altitudes. Compared with Fig. 3 in Krämer et al. (2020), 48 % of their ice particle samples have D_i<40 µm, which is below the size cut-off used in this study. The higher D_i in this study also leads to lower range of N_i (0.01–1000 L⁻¹) and higher range of IWC (10⁻⁵–10 g m⁻³) compared with that previous study (i.e., N_i from 0.1–10⁵ L⁻¹ and IWC from 10⁻⁷–1 g m⁻³), representing the sampling bias towards larger particles in this study. The relationships of IWC with respect to meteorological conditions (i.e., temperature and RHi) and other microphysical properties (i.e., N_i and D_i) are shown in Fig. S1 in the Supplement. The distributions of IWC samples are relatively uniform at various temperature and RHi, while more IWC samples are correlated with D_i between 100 and 300 µm.

Figure 3 Latitude and altitude distributions of (a) IWC, (b) N_i, (c) D_i, (d) in-cloud RHi, (e) clear-sky RHi and (f) clear-sky water vapor volume mixing ratio at temperatures $\le -\mathrm{40}$ ^∘C. Total measurement hours and number of samples for given intervals are shown for each variable. Note that the measurement ranges shown in the upper right corner are not the full ranges (see Table 2 for the full ranges).

2.2 Climate model description and experiment design

This study uses model simulations based on the NCAR CAM6 model. Compared with its previous version (the CAM5 model), CAM6 implemented a new scheme, the Cloud Layers Unified by Binomials (CLUBB) for representations of boundary layer turbulence, shallow convection and cloud macrophysics (Bogenschutz et al., 2013). CLUBB is a higher-order turbulence closure scheme that calculates prognostic higher moments based on joint probability density functions (PDFs) for vertical velocity, temperature and moisture (Golaz et al., 2002). An improved bulk two-moment cloud microphysics scheme has been implemented (Gettelman and Morrison, 2015) that replaces diagnostic treatment of rain and snow with prognostic treatment of all hydrometeors (i.e., rain and snow). This is coupled with a four-mode aerosol model (MAM4) (Liu et al., 2016) for simulations of aerosols and aerosol–cloud interactions. It allows ice crystals to form via homogeneous freezing of sulfate aerosols and heterogeneous nucleation of dust particles (Liu et al., 2007; Liu and Penner, 2005). The model uses Wang et al. (2014b) for ice nucleation, which implemented and improved Hoose et al. (2010) by considering the probability density function of contact angles for the classical nucleation theory. The model also uses Shi et al. (2015) for modifications of pre-existing ice. Finally, the deep convection scheme (Zhang and McFarlane, 1995) has been tuned to include sensitivity to convection inhibition.

Results from in situ observations are compared with two types of CAM6 simulations: nudged and free-running simulations. Simulations are based on a finite-volume dynamical core (Lin, 2004) with a horizontal resolution of 0.9^∘ × 1.25^∘ and 32 vertical levels. All simulations are conducted using prescribed sea-surface temperature and present-day aerosol emissions and include a 6-month spin-up time. CAM6 nudged simulations are nudged spatially and temporally with meteorological data (i.e., 2-D horizontal wind and temperature) from the Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA2) (Gelaro et al., 2017) and collocated with aircraft flight tracks in space and time. A nudged simulation was conducted for each campaign independently and was combined into one dataset to compare with observations. One free-running simulation was conducted for the duration of all flight campaigns from July 2008 to February 2016. To reduce the size of model output when comparing with observations, a total of 24 instantaneous output data from the free-running simulation are combined into one dataset (“CAM6-free” hereafter), which includes 00:00 and 12:00 UTC for the first day of each month in 2010. Additional sensitivity tests on different model output from the free-running simulation show very minor differences in the statistical distributions of cirrus microphysical properties and the correlations with their controlling factors when selecting different years, seasons and days in a month.

In order to examine observations and simulations on more comparable scales, a running average of 430 s was calculated for meteorological parameters (i.e., temperature and RHi) and microphysical properties (i.e., IWC, N_i and D_i), which translates to ∼100 km horizontal scales since the mean true air speed below −40 ^∘C for all campaigns was 230 m s⁻¹ (Fig. S2). When applying the running average, both in-cloud and clear-sky conditions (i.e., where N_i and IWC values are zero) are included in the averages. Grid-mean quantities from model output are used in comparisons with observations, including “IWC”, “NUMICE”, “QSNOW” and “NSNOW”, which are mass and number concentrations of ice particles and snow, respectively. Another type of comparison between 1 Hz observations and in-cloud quantities from model output is shown in the Supplement. Both methods have been previously used in model evaluation, such as D'Alessandro et al. (2019), who compared 200 s averaged aircraft observations with simulated grid-mean quantities, and Righi et al. (2020), who compared 1 Hz aircraft observations with simulated in-cloud quantities.

Table 2 Ranges of meteorological conditions and ice microphysical properties for in situ 1 Hz observations, 430 s averaged observations, CAM6-nudg and CAM6-free data used in this study.

A summary of the ranges of meteorological conditions and ice microphysical properties for in situ observations and simulations is shown in Table 2. Simulated RHi is calculated from simulated specific humidity and temperature, and the calculation of saturation vapor pressure with respect ice is based on the equation from Murphy and Koop (2005). Simulated ice and snow are restricted to ≥62.5 µm based on the size cut-off of the Fast-2DC probe by applying methods from Eidhammer et al. (2014). Based on their equations 1 to 5, we followed their assumption that the shape parameter μ equals 0 when calculating the slope parameter λ. Mass and number concentrations of ice and snow are further calculated based on integrals of incomplete gamma functions from 62.5 µm to infinity. The simulated values of IWC, N_i and D_i are calculated based on the combined ice and snow population after applying the size restriction. In-cloud conditions in simulations are defined by concurring conditions of IWC $>{\mathrm{10}}^{-\mathrm{7}}$ g m⁻³ and $N_{\mathrm{i}}>{\mathrm{10}}^{-\mathrm{4}}$ L⁻¹ based on size-restricted grid-mean quantities. These thresholds are the lower limits from observations after calculating the 430 s averages. Note that due to the ice crystal size constraint, some thin cirrus may not be detected. In addition, analysis of simulated cirrus clouds is restricted to similar pressure ranges as those measured in the seven campaigns. An additional constraint on cloud fraction $>{\mathrm{10}}^{-\mathrm{5}}$ was applied to both nudged and free-running simulations to exclude extremely low values.

Figure 4 (a) Observed size distribution (black line) and reconstructed size distributions from simulated ice (blue) and snow (cyan). Size truncations to diameters <62.5 µm (dashed lines) are shown for simulated hydrometeors, while the remaining particles (≥62.5 µm) (solid lines) are used for comparisons with observations. Size distributions for combined ice and snow in the simulations (purple) are also shown before and after the size restriction. Panels (b) and (c) show PDFs of N_i and IWC in the simulation before and after size truncation.

To visualize the impact of the size truncation on simulated data, we employed methods similar to Gettelman et al. (2020) and reconstructed the simulated particle size distributions for snow and ice in Fig. 4a using gamma functions from Morrison and Gettelman (2008). Compared with the observations, the number density for combined ice and snow is overestimated for smaller particles (<400 µm) and underestimated for larger particles (>1000 µm). After applying the size restriction, the PDFs of simulated N_i and IWC show increasing probability of small N_i and decreasing probability of small IWC due to the removal of small particles (Fig. 4b and c).

Finally, simulated aerosol number concentrations are further categorized by diameters >500 and >100 nm (i.e., N_a,500 and N_a,100, respectively), by summing the size-restricted concentrations of the Aitken, accumulation and coarse aerosol modes. Previously, field experiments found that N_a,500 correlates well with INP number concentrations (DeMott et al., 2010). Even though that correlation was only determined based on observations warmer than −36 ^∘C, the separation of N_a,500 and N_a,100 can help to examine the effects of larger and smaller aerosols in this work.

Figure 5 Geometric means of (a–c) IWC and (d–f) N_i, as well as (g–i) linear averages of D_i at 5 ^∘C temperature intervals between −80 and −40 ^∘C, compared between 1 Hz in situ observations (blue lines), 430 s averaged observations (black lines) and CAM6-nudg (red lines). Observed and simulated microphysical properties are binned by six latitudinal regions, where NH is denoted by solid lines, and SH is denoted by dashed lines. The number of samples for 1 Hz observations at temperatures $\le -\mathrm{40}$ ^∘C in the Northern (Southern) Hemisphere tropical, midlatitude and polar regions are 173 930 (15 569), 100 615 (3809) and 6704 (2606), respectively. The number of samples for 430 s averaged observations in these regions are 355 082 (40 683), 233 546 (26 850) and 24 083 (10 252), respectively. The number of samples for CAM6-nudg data in these regions are 3 241 592 (653 110), 2 052 353 (590 503) and 478 844 (209 662), respectively.

3 Regional variations of cirrus cloud characteristics

3.1 Cirrus cloud microphysical properties with respect to temperature

Three cirrus cloud microphysical properties (IWC, N_i and D_i) are examined in relation to temperature at six latitudinal regions (Fig. 5). The standard deviations of the IWC, N_i and D_i in each temperature bin are shown in Fig. S3. The 1 Hz observations of IWC and N_i in the NH indicate clear latitudinal differences with the highest values occurring in the midlatitudes, followed by tropics, then polar regions for temperatures between −40 and −60 ^∘C, while for colder temperatures the NH tropical region shows the highest IWC. In the SH, the highest IWC and N_i occur in the tropics, followed by the polar regions and midlatitudes. Comparing the two hemispheres, IWC and N_i show significant reductions by ∼1 order of magnitude from NH midlatitudes to SH midlatitudes (Fig. 5b, e). These hemispheric differences in midlatitudes may be due to air-mass differences between NH (more continental) and SH (more oceanic) and/or more anthropogenic emissions in the NH. The IWC, N_i and D_i are relatively similar between NH and SH tropical regions, while IWC and D_i are higher in the SH polar region than in NH polar regions.

The simulations are further compared with averaged observations at a similar horizontal scale of ∼100 km. After applying 430 s running averages for observations, the average IWC and N_i values decrease by 0.5–1.5 orders of magnitude compared with 1 Hz observations depending on temperature and geographical region. Hemispheric differences are mostly consistent between 1 and 430 s averaged observations except for polar regions. CAM6-nudg data show a similar trend of average IWC, N_i and D_i with respect to temperature as seen in observations; that is, the average IWC increases with increasing temperature, consistent with previous observational studies (Krämer et al., 2016; Luebke et al., 2013; Schiller et al., 2008), average N_i shows no clear trend with temperature, and average D_i increases with increasing temperature. Differing from observations, CAM6 produces the highest IWC and N_i in the tropical regions, followed by midlatitudes, then polar regions for both hemispheres. The simulated IWC, N_i and D_i also show smaller differences between hemispheres and latitudes. The CAM6-nudg data underestimate and overestimate IWC in the NH and SH by 0.5–1 orders of magnitude, respectively, with the largest discrepancies in the midlatitudes. The simulations overestimate N_i in the tropics and polar regions in both hemispheres by 0.5–1 orders of magnitude and overestimate N_i in the SH midlatitudes by 1–2 orders of magnitude. The simulated D_i is about half of the observed values in most regions except polar regions. This result indicates “too many” and “too small” simulated ice in most regions. The low bias of simulated D_i indicates possible misrepresentation of ice particle growth and sedimentation in the model parameterization.

A sensitivity test is conducted by comparing 1 Hz observations with in-cloud quantities from model output (Fig. S4). Larger differences are seen between simulated and observed IWC and N_i in Fig. S4 compared with Fig. 5. The directions (i.e., positive or negative) of model biases of IWC, N_i and D_i are generally consistent in both comparisons.

A previous study by Righi et al. (2020) evaluated the ice microphysical properties in the EMAC-MADE3 aerosol–climate model (i.e., ECHAM/MESSy Atmospheric Chemistry – Modal Aerosol Dynamics model for Europe adapted for global applications, third generation) by comparing in-cloud quantities from model output with 1 Hz in situ observations of multiple aircraft field campaigns from 75^∘ N to 25^∘ S (Krämer et al., 2009, 2016, 2020). Although that study included a greater number of smaller ice particles (3–1280 µm) compared with this study, they still showed low biases of simulated D_i at 190–243 K and low biases of simulated IWC at 205–235 K, as well as high biases of simulated N_i above 225 K, which are generally in the same direction as the biases we found in CAM6 model. Note that Righi et al. (2020) implemented different cloud microphysics parameterizations compared with the CAM6 model, including a two-moment cloud microphysics scheme of Kuebbeler et al. (2014) and the ice nucleation parameterization for cirrus clouds (T<238.15 K) from Kärcher et al. (2006), which account for both homogeneous and heterogeneous nucleation and the competition between the two mechanisms. Additional future intercomparison studies of these models are warranted to examine the reasons behind the similar biases.

Figure 6 Distributions of RHi at various temperatures and geographical locations from in situ observations under (two left columns) clear-sky and (two right columns) in-cloud conditions. Solid and dashed black lines represent ice and liquid saturation, calculated based on saturation vapor pressure with respect to ice and liquid from Murphy and Koop (2005), respectively. Dash-dotted line denotes the homogeneous freezing threshold for 0.5 µm aerosols based on Koop et al. (2000).

Figure 7 Similar to Fig. 6 but for 430 s scale.

Figure 8 Similar to Fig. 6 but for CAM6-nudg data. RHi values for simulations are calculated using simulated specific humidity and temperature, based on the equation of saturation vapor pressure with respect to ice from Murphy and Koop (2005).

3.2 RHi and σ_w distributions for in-cloud and clear-sky conditions

Regional distributions of RHi for clear-sky and in-cloud conditions are shown for 1 Hz observations (Fig. 6), 430 s averaged observations (Fig. 7) and simulations (Fig. 8). The 1 Hz observations show RHi magnitudes ranging from <5 % up to ∼180 % in both clear-sky and in-cloud conditions, and are mostly located below the homogeneous freezing line except for the NH tropical region. A few samples exceed the liquid saturation line but are within the measurement uncertainties of RHi. This result agrees with the RHi distributions based on previous midlatitudinal observations (Cziczo et al., 2013). Differing from 1 Hz observations, 430 s averaged observations show much lower RHi magnitudes for both clear-sky and in-cloud conditions, ranging from <5 % to 120 %–140 %. For clear-sky conditions, the majority of the observed and simulated RHi values are below 100 %, while the CAM6-nudg data show fewer RHi exceeding ice saturation. For in-cloud conditions, both 1 Hz observations and simulations show that RHi frequently occur within ∼20 % of ice saturation, consistent with previous observation and modeling studies (Diao et al., 2014a, 2017; D'Alessandro et al., 2017, 2019; Krämer et al., 2009), while almost no simulated RHi data exceed the homogeneous freezing threshold. The higher RHi observed in the NH tropical region was also observed by Krämer et al. (2009). Such features can be explained by the competition between higher updrafts seen in the tropics and the depletion of water vapor from newly nucleated ice particles as discussed in Kärcher and Lohmann (2002). For the polar regions, in-cloud RHi is skewed towards ISS in both observations and simulations, indicating less effective water vapor depletion likely due to lower N_i values (Fig. 5e). Note that the simulation samples in the tropical regions show peak frequencies at certain temperatures due to larger bin sizes of pressure levels in the lower latitudes.

Figure 9 Distributions of σ_w calculated for every 40 s using the 1 Hz observations under (two left columns) clear-sky and (two right columns) in-cloud conditions.

Figure 10 Similar to Fig. 9 but for 430 s averaged observations.

Figure 11 Similar to Fig. 9 but for the CAM6-nudg data.

Regional distributions of the variance of w (σ_w) for 1 Hz observations at 40 and 430 s scales and CAM6 nudged simulations are shown in Figs. 9, 10 and 11, respectively. σ_w in the observations is calculated as the variance of w within each 40 and 430 s of data, which correspond to a horizontal scale of ∼10 and 100 km, respectively. The σ_w in simulations is based on the “w_sub” variable, which is calculated from the square root of turbulent kinetic energy (TKE) (Gettelman et al., 2010). Observed σ_w shows the highest values in the tropical and midlatitude regions, reaching up to ∼3 m s⁻¹, while the polar regions show updrafts up to ∼1 m s⁻¹. A similar decreasing trend of maximum σ_w is seen in the simulations from the lower to higher latitudes. The observations show similar σ_w maximum values between clear-sky and in-cloud conditions, while the simulations show much higher maximum σ_w for in-cloud conditions in the tropics (1 m s⁻¹), midlatitude (1 m s⁻¹) and polar regions (0.5 m s⁻¹), compared with those values in clear sky (i.e., 0.5, 0.25 and 0.1 m s⁻¹, respectively). This result suggests that the model has a stronger dependence on higher σ_w for cirrus cloud formation compared with observations. We further examine the potential impact of convection in simulations and observations. Figure S5 shows the locations where w>1 m s⁻¹ is seen in the observations as well as where w_sub>0.5 m s⁻¹ is seen in the CAM6-nudg data for in-cloud conditions. Since w_sub in CAM6 is based on the turbulent scheme, higher w_sub values indicate that the convection scheme may be active and produce detrained ice in convective outflows. The majority of observed and simulated in-cloud samples do not appear to have high w or w_sub, indicating that detrained ice from the convection is unlikely a significant contribution. More future investigation is needed to track cirrus cloud origins and quantify impacts from convection.

Figure 12 Probability density functions (PDFs) for (a, d, g, j) temperature, (b, e, h, k) RHi and (c, f, i, l) σ_w, compared among (a–c) 1 Hz observations, (d–f) 430 s averaged observations, (g–i) CAM6-nudg and (j–l) CAM6-free data. Note that σ_w in panel (c) is calculated for every 40 s.

4 Individual impacts of key controlling factors on cirrus clouds

4.1 Probability density functions of temperature, RHi and σ_w

PDFs of temperature, RHi and σ_w are shown in Fig. 12. The PDFs are normalized by the total number of samples of both clear-sky and in-cloud conditions. The observations are located mostly around −68 to −40 ^∘C, and the simulations show similar temperature distributions. For the PDFs of RHi, the observations and simulations all show a peak position around 100 % for in-cloud condition. However, a secondary peak is shown in simulations at 80 % RHi, which is likely due to the parameter of RHi_min for ice cloud fraction calculation being set at 80 % for representing variance of humidity in a grid box (more details on RHi_min are described in Gettelman et al., 2010). In addition, the maximum RHi values for in-cloud conditions are 170 %, 154 %, 160 % and 257 % for 1 Hz observations, 430 s averaged observations, CAM6-nudg and CAM6-free, respectively. The maximum RHi values for clear-sky conditions are 175 %, 135 %, 108 % and 181 %, respectively. The CAM6-free data show higher maximum RHi values than CAM6-nudg data, likely due to additional data from tropical regions at temperatures below −70 ^∘C (Fig. 12j). When using a lower size cut-off (1 µm) of ice particles for the simulation data, the number of in-cloud samples increases by 4 % (Fig. S6). However, negligible differences are seen in the PDFs of temperature, RHi and σ_w for the two simulations between Figs. 12 and S6.

PDFs of σ_w show consistent results with Figs. 9–11, with simulations showing much higher maximum σ_w for in-cloud conditions than clear-sky conditions, while observations show similar maximum σ_w in both conditions. The lower maximum values of σ_w in simulations are most likely a result of the model missing representations of gravity waves from topography, fronts and convection, and only including σ_w from turbulence.

Figure 13 Correlations between RHi and in-cloud IWC, N_i and D_i (columns 1–3, respectively), compared among (a–c) in situ 1 Hz observations, (d–f) 430 s averaged observations, (g–i) CAM6-nudg and (j–l) CAM6-free data. Black lines and whiskers denote geometric means and standard deviations, respectively.

4.2 Effects of RHi and σ_w on ice microphysics

The relationships between ice microphysical properties and RHi are examined in Fig. 13. For the 1 Hz observations, the maximum IWC and N_i occur slightly above ice saturation at 110 % RHi, while the maximum D_i occur at 130 % RHi. The average IWC and N_i increase 1.5 orders of magnitude from 40 % to 110 % RHi and decrease 0.5 orders of magnitude (i.e., a factor of 3) from 110 % to 130 % RHi. The maximum IWC and N_i do not occur at the highest RHi most likely due to the consumption of water vapor by ice deposition. High D_i values at lower RHi (∼30 %) are likely a result of sedimenting large ice crystals, which has been previously observed by Diao et al. (2013) when investigating the evolutionary phases of cirrus clouds. For 430 s averaged observations, the peak IWC, N_i and D_i occur at 100 %, 100 % and 115 % RHi. The maximum IWC and N_i values are nearly the same between 1 and 430 s averaged observations (i.e., 0.04 g m⁻³ and 10 L⁻¹, respectively) near saturation, while the 430 s averaged observations show lower minimum IWC and N_i at very low RHi (<10 %), which are 1.5 orders of magnitude lower than 1 Hz observations. This feature is due to in-cloud segments being longer around saturation compared with subsaturated conditions as shown in Diao et al. (2013), which means that fewer clear-sky conditions are being included in the 430 s averages around saturation and therefore show little reduction of the IWC and N_i due to spatial averaging.

In contrast to observations, both CAM6-nudg and CAM6-free simulations show bimodal distributions of IWC and N_i with the primary peak at 100 % RHi and the secondary peak at 80 % RHi. The secondary peak at RHi 80 % is likely produced by the RHi_min parameter reflecting subgrid-scale RHi variance as mentioned above (Gettelman et al., 2010), which was set at the default value (80 % RHi) for both simulations. The primary peak at 100 % RHi is likely a result of the minimum threshold for heterogeneous ice nucleation being set at 120 % as well as a subgrid variability scaling factor of 1.2 being considered (Wang et al., 2014a). Similar to 430 s averaged observations, IWC and N_i show steep increase (i.e., 3–4 orders of magnitude) from 40 % to 100 %. Increases of average IWC and N_i are seen in the simulations as RHi increases from 120 % to 160 %, differing from the decreasing trend seen in the observations. These higher values of IWC and N_i near 160 % are possibly due to RHi reaching the homogeneous nucleation thresholds, where ice nucleation becomes more dependent upon temperature and updraft speed (Liu and Penner, 2005). Note that at this same point as IWC and N_i increase, there is a decrease in D_i, which also suggests homogeneous nucleation in the model. For D_i–RHi correlations, both simulations show similar results to the observations, with the maximum D_i around 130 % RHi and some large ice particles in the subsaturated conditions. The large variability of observed ice microphysical properties is also significantly underestimated in the model for ISS conditions. Standard deviations are 0.5–1 orders of magnitude lower for IWC and N_i and a factor of 2 lower for D_i compared with observations.

Figure 14 Similar to Fig. 13 but for correlations with σ_w. Note that σ_w values of observations are calculated for panels (a–c) every 40 s and (d–f) every 430 s using the 1 Hz observations.

Comparing the correlations with σ_w (Fig. 14), the simulations show increasing IWC and N_i with higher σ_w, which agree with observations, although the increases of IWC and N_i are smaller in the simulations than the 430 s observations. The simulated D_i is relatively constant with increasing σ_w, which differs from the slight positive correlation between D_i and σ_w in the observations. This slight positive D_i–σ_w correlation is likely due to the growth of ice particles as cirrus clouds evolve with continuous updrafts that supply excess water vapor above ice saturation, which was previously discussed in a cirrus cloud evolution analysis (Diao et al., 2013). The simulations may overlook this positive correlation for several reasons, such as the lack of temporal resolution to resolve cirrus evolution in the growth phase, the lack of vertical velocity subgrid variabilities (as discussed in Zhou et al., 2016) and a dry bias (i.e., lower RHi) in the model (as discussed in Wu et al., 2017).

Comparing the performance of two types of simulations, both CAM6-nudg and CAM6-free show bimodal distributions for IWC–RHi and N_i–RHi correlations, and they both show positive correlations for IWC–σ_w and N_i–σ_w. This result indicates that the general trends in these correlations are statistically robust and less affected by sampling sizes and geographical locations. For correlations with RHi, the maximum IWC value in CAM6-nudg and CAM6-free is lower than the 430 s averaged observations by a factors of 25 and 100, respectively. The maximum N_i value in CAM6-nudg is similar to the 430 s averaged observations, while that value in CAM6-free is lower by a factor of 3. For correlations with σ_w, there are no significant differences for the maximum IWC between the two simulation types. The maximum N_i value in CAM6-nudg and CAM6-free is higher than the 430 s averaged observations by factors of 3 and 10, respectively. These results show that CAM6-nudg data, which are collocated with flight tracks, produce IWC and N_i values closer to the 430 s averaged observations than CAM6-free, possibly due to the variabilities of IWC and N_i in different geographical locations as shown in Fig. 5.

Figure 15 Aerosol indirect effects from logarithmic-scale N_a,500 on (a–c) IWC, (d–f) N_i, (g–i) D_i and (j–l) cloud fraction, compared among 1 Hz observations (left column), 430 s averaged observations (middle column) and CAM6-nudg data (right column). Number of samples of each bin is shown in the bottom row (m–o). Cloud fraction is calculated as the number of in-cloud samples over the total number of samples for a given temperature and N_a bin.

Figure 16 Similar to Fig. 15 but examined for log₁₀(N_a,100).

4.3 Aerosol indirect effects

The effects of larger and smaller aerosols (i.e., N_a,500 and N_a,100) on ice microphysical properties are further examined for observations and CAM6-nudg data (Figs. 15 and 16). Cloud fraction is calculated in each temperature–N_a bin by normalizing the number of in-cloud samples with the total number of samples in that bin for both observations and simulations. For three cirrus microphysical properties (i.e., IWC, N_i and D_i), positive correlations are seen in 1 Hz observations with respect to N_a,500 and N_a,100. In addition, higher N_a,500 (>10 cm⁻³) and N_a,100 (>100 cm⁻³ values are associated with significant increases in cloud fraction. At −70 to −60 ^∘C, higher IWC, N_i and cloud fraction are seen when N_a,500 is observed, with positive correlations of IWC and N_i with respect to N_a,500. This finding indicates that larger aerosols provide an effective pathway of ice particle formation for colder conditions. The higher IWC and N_i are only shown in much higher N_a,100 (>100 cm⁻³) between −70 and −60 ^∘C, demonstrating that larger aerosols facilitate ice formation more effectively than smaller aerosols at this temperature range, possibly due to the activation of larger aerosols as INPs for heterogeneous nucleation. Compared with 1 Hz observations, 430 s averaged observations show weaker correlations of IWC, N_i and D_i with respect to N_a. However, they do still show higher IWC and N_i between −70 and −40 ^∘C associated with higher N_a (i.e., $N_{\mathrm{a},\mathrm{500}}>\mathrm{1}$ cm⁻³ and $N_{\mathrm{a},\mathrm{100}}>\mathrm{30}$ cm⁻³).

The CAM6-nudg simulation shows increasing average IWC, average N_i and cloud fraction with increasing N_a,500, consistent with the observations. But at temperatures below −60 ^∘C, simulated IWC and N_i do not show a sudden increase with higher N_a,500 as shown in the observations. The simulated D_i slightly decreases with increasing N_a,500, differing from the increasing trend seen in observations. For aerosol indirect effect analysis based on N_a,100, the comparison results are similar to N_a,500; that is, the CAM6-nudg simulation is able to represent positive correlations of IWC, N_i and cloud fraction with respect to N_a,100. However, the CAM6-nudg simulation shows smaller (larger) increase of IWC (N_i) at very high N_a (i.e., $N_{\mathrm{a},\mathrm{500}}>\mathrm{1.6}$ cm⁻³ and $N_{\mathrm{a},\mathrm{100}}>\mathrm{30}$ cm⁻³) compared with the 430 s averaged observations. The model also misses positive correlations between D_i and N_a seen in both 1 Hz and 430 s averaged observations.

5 Discussion and conclusions

In this study, we investigate the statistical distributions of cirrus cloud microphysical properties (i.e., IWC, N_i and D_i) as well as several key controlling factors (i.e., temperature, RHi, σ_w and N_a) using a comprehensive in situ observational dataset and GCM simulations. Regional variations of cirrus cloud microphysical properties are examined for six latitudinal regions in two hemispheres. Two types of CAM6 simulations are evaluated, i.e., nudged and free-running simulations.

Regarding the regional variations in 1 Hz observations, the highest and lowest IWC values were observed in NH midlatitudes and SH midlatitudes, respectively, while the polar regions show the lowest N_i and highest D_i at warmer conditions (i.e., −55 to −40 ^∘C) (Fig. 5). The hemispheric differences between NH and SH midlatitudes indicate a possible role of anthropogenic aerosols and/or land–sea contrast in controlling ice microphysical properties. Thermodynamic and dynamic conditions can also affect nucleation mechanisms. For example, the tropical regions show the highest IWC and N_i at temperatures below −55 ^∘C possibly due to convection anvils with the droplet freezing from down below or homogeneous nucleation in gravity waves generated by convection. This feature is corroborated by the fact that tropical regions show the highest RHi values for both clear-sky and in-cloud conditions (Fig. 6), while the midlatitude and polar regions show fewer samples exceeding the homogeneous nucleation threshold. The higher RHi values in tropics are likely contributed by higher updrafts (indicated by higher σ_w in Fig. 9). These results demonstrate the important roles of these controlling factors on cirrus clouds at different latitudinal and temperature ranges.

Evaluating the model simulations of cirrus microphysical properties, different model performance results are seen in different regions. For example, simulations underestimated the IWC in NH (Fig. 5), possibly due to model dry bias to form ice clouds (as discussed in Wu et al., 2017) and/or smaller aerosol indirect effects on IWC in the simulations (Figs. 15 and 16). Differences in the particle size distribution, such as lower number density of larger particles (>1000 µm) in the simulation (Fig. 4a), may also contribute to the underestimation of IWC by the simulation. All the comparison results on IWC, N_i and D_i are only applicable to the size range being evaluated (≥62.5 µm). For RHi distributions, both simulations represent a similar peak position at ice saturation for in-cloud RHi PDFs compared with observations but CAM6-nudg underestimates the frequency and magnitude of ISS for clear-sky condition. For σ_w distributions, simulations represent similar regional variations of σ_w compared with observations, with σ_w decreasing from lower to higher latitudes. The model performs well for representing the effects of RHi and σ_w on ice microphysical properties, specifically for showing the maximum IWC and N_i at 100 % RHi and the positive correlations with σ_w. Some differences include the simulated average IWC and N_i showing a secondary peak position at 80 % RHi, likely due to the minimum RHi threshold used in the model parameterization. Both simulation types show similar correlation trends of ice microphysical properties with respect to RHi and σ_w. CAM6-nudg performs better for representing IWC and N_i magnitudes than CAM6-free, possibly due to better collocation between CAM6-nudg and observations.

For aerosol indirect effects, the simulations underestimate IWC, N_i, D_i as well as cloud fraction at colder conditions ($<-\mathrm{60}$ ^∘C) when larger aerosols exist, indicating that the effectiveness of larger aerosols is underestimated at the colder conditions. The observations also show higher D_i than simulations by a factor of 3–4 at warmer temperatures (−50 to −40 ^∘C), indicating misrepresentation of ice particle growth and/or sedimentation in the simulations. In addition, the IWC, N_i and D_i in 430 s averaged observations show an increase at higher N_a,500 (>1 cm⁻³) and N_a,100 (>30 cm⁻³), while simulations only show a significant increase of N_i. This result indicates that aerosol indirect effects may be underestimated especially for higher N_a values. It is possible that small ice crystals <62.5 µm may have formed under high N_a but are excluded due to the size constraint. Additionally, because INP activation is highly dependent upon temperature, we acknowledge the limitation of using N_a,500 to indicate INP concentrations. The assumption of ice mass and dimension relationship from Brown and Francis (1995) may also lead to uncertainties due to various ice habits. These caveats call for more investigation of small ice measurements, INP measurements at temperature $\le -\mathrm{40}$ ^∘C and measurements of various ice habits.

Overall, the global-scale observational dataset used in this study provides statistically robust distributions of cirrus cloud microphysical properties, which can be used to evaluate the effects of thermodynamics, dynamics and aerosols on cirrus clouds in a global climate model. Extending from previous studies that investigated climate model sensitivity to individual cirrus cloud controlling factors, i.e., w (Shi and Liu, 2016), RHi (D'Alessandro et al., 2019), water vapor (Wu et al., 2017) and aerosols (Wang et al., 2014a), this study provides an analysis of all factors. In addition, further attention was given towards evaluating these factors in the simulations based on geographical locations. For both observations and simulations, higher ice supersaturations and stronger vertical motions are shown in tropical and midlatitude regions, which possibly lead to increased homogeneous nucleation and convection-generated cirrus, consistent with higher IWC and N_i and lower D_i in these regions compared with polar regions. In addition, underestimating aerosol indirect effects in the simulations likely contributes to the underestimation of IWC in the NH. Even though small ice particles (<62.5 µm) are excluded in this study, correlations between ice microphysical properties and these key controlling factors are still clearly seen in the observation dataset. In addition, using two methods that compare observations on the horizontal scales of 230 m and 100 km with simulations, both methods show similar signs for model biases of IWC, N_i and D_i, while smaller model biases are seen when comparing against the coarser resolution observations. This study underscores the importance of correctly representing the thermodynamic, dynamic and aerosol conditions in climate models at various regions, as well as accurately simulating their correlations with ice microphysical properties. Failing to do so may result in biases of cirrus cloud microphysical properties depending on different regions and temperatures, leading to biases in cirrus cloud radiative effects on a global scale.