Subgrid-scale variability of cloud ice in the ICON-AES 1.3.00

This paper presents a stochastic approach for the aggregation process rate in the ICOsahedral Nonhydrostatic general circulation model (ICON-AES), which takes subgrid-scale variability into account. This method creates a stochastic parameterization of the process rate by choosing a new specific cloud ice mass at random from a uniform distribution function. This distribution, which is consistent with the model's cloud cover scheme, is evaluated in terms of cloud ice mass variance with a combined satellite retrieval product (DARDAR) from the satellite cloud radar CloudSat and the Cloud–Aerosol Lidar and Infrared Pathfinder Observations (CALIPSO). The global patterns of simulated and observed cloud ice mixing ratio variance are in a good agreement, despite an underestimation in the tropical regions, especially at lower altitudes, and an overestimation in higher latitudes from the modeled variance. Due to this stochastic approach the yearly mean of cloud ice shows an overall decrease. As a result of the nonlinear nature of the aggregation process, the yearly mean of the process rates increases when taking subgrid-scale variability into account. An increased process rate leads to a stronger transformation of cloud ice into snow and therefore to a cloud ice loss. The yearly averaged global mean aggregation rate is more than 20 % higher at selected pressure levels due to the stochastic approach. A strong interaction of aggregation and accretion, however, lowers the effect of cloud ice loss due to a higher aggregation rate. The new stochastic method presented lowers the bias of the aggregation rate.

1 Introduction

A correct representation of clouds and cloud-related microphysical process rates, which describes the time-dependent source and sink terms of cloud ice or liquid water, is one of the central challenges in global climate modeling. Since global climate models typically run at a rather coarse resolution (of the order of 100 km), it is important to look into the unresolved microphysical process rates. Most climate models use a cloud cover parameterization. Microphysical process rates are computed based on grid box mean in-cloud ice or liquid water mixing ratios. The in-cloud cloud ice mass mixing ratio is the cloud ice mass per cloudy area. Considering subgrid-scale variability of in-cloud variables reduces the biases of the nonlinear microphysical process rates (Pincus and Klein, 2000; Larson et al., 2001). For example, Weber and Quaas (2012) numerically integrated the process rate over the probability density function (PDF), which is a very accurate method but needs additional computational time. Another method, which works with no additional computational time, is a stochastic approach for the process rates by taking a randomly chosen value per time step and grid box (e.g., Palmer, 2001; Berner et al., 2017). With this method a randomly disturbed process rate is created in order to give a better representation of the state of the atmosphere. Another method is the Cloud Layers Unified By Binormals (CLUBB) (Golaz et al., 2002a, b), which works with a set of PDFs for all cloud types to avoid difficulties in coupling of stratocumulus and shallow convection parameterizations. Since CLUBB does not parameterize subgrid-scale variability of cloud ice, Thayer-Calder et al. (2015) include cloud ice in the CLUBB PDFs. The PDFs are sampled by the subgrid importance Latin hypercube sampler (SILHS), which connects CLUBB with the microphysics for stratiform and convective clouds. They found improvements, e.g., of liquid water path (LWP), precipitable water, and shortwave cloud forcing, but also a degradation in precipitation. Including the subgrid-scale effect in the autoconversion and accretion rate in warm clouds reduces the bias significantly and leads to an enhancement of the process rate (Boutle et al., 2014; Lebsock et al., 2013). Since previous studies mainly focus on warm-rain formation processes (e.g., Morrison and Gettelman, 2008; Larson and Griffin, 2013; Boutle et al., 2014; Lebsock et al., 2013), it is also important to concentrate on snow formation effects by taking subgrid-scale effects into account.

Mülmenstädt et al. (2015) emphasized that most global rain is produced via ice-phase processes. Thus, we focus on cloud-ice-related process rates, especially on the precipitation processes initiated via the ice phase. In the ICOsahedral Nonhydrostatic general circulation model (ICON-AES) snow is formed by the aggregation process. Aggregation describes the process where cloud ice particles grow to snowflake sizes by sticking together. Therefore, in this study we implement a stochastic aggregation parameterization into the ICON-AES (Giorgetta et al., 2018) by taking subgrid-scale variability of cloud ice into account. Subgrid-scale variability of the total water mixing ratio (sum of cloud liquid water, cloud ice, and water vapor) is already used for determining the cloud cover (Sundqvist et al., 1989). Here, we use this uniform distribution approach to create a distribution of cloud ice within the cloudy part of the grid box. Instead of taking a grid box mean in-cloud ice mixing ratio for the nonlinear aggregation parameterization, we feed the process rate with a cloud ice mass randomly chosen from the distribution of cloud ice mass assumed in the cloud scheme.

To evaluate the uniform distribution of cloud ice at a global scale, large-scale observations of cloud ice are necessary. The combined dataset of spaceborne radar from the CloudSat satellite (Stephens et al., 2002) and the lidar from the Cloud–Aerosol Lidar and Infrared Pathfinder Observations satellite (CALIPSO, Winker et al., 2010) allows retrieving a global dataset of cloud ice. Here, we use the lidar–radar (DARDAR; Delanoë and Hogan, 2008, 2010) dataset. Comparing observed cloud ice with modeled cloud ice is a challenge. Since the cloud ice from the ICON-AES does not include snow (precipitating cloud ice) and convective cloud ice, it is necessary to remove falling and convective cloud ice from the DARDAR dataset in order to give a meaningful comparison between model and observations. A flag method from Li et al. (2008, 2012) is used to estimate the precipitating and convective cloud ice part per grid box, which is removed from the dataset.

In this study, we investigate an important process rate which transforms ice to snow, the aggregation rate, and how it is treated in the ICON-AES. We include the stochastic approach in the aggregation in order to quantify the influence of taking subgrid-scale variability into account. The selected distribution function of cloud ice is evaluated with the DARDAR dataset and the effect of the stochastic approach in the ICON-AES is investigated. As an additional evaluation, we compare an unbiased process rate with the stochastic approach in order to investigate how well this simple method performs.

2 Methods

For all simulations the ICOsahedral Nonhydrostatic general circulation model (ICON-AES-1.3.00, Giorgetta et al., 2018) in its global version is used. It includes the Max Planck Institute physics package based on the ECHAM6 physics (Stevens et al., 2013). All runs were performed for 6 years with prescribed sea surface temperature and sea ice boundary conditions for a period from 2004 to 2009 and with instantaneous diagnostics output every 6 h by using a model time step of 10 min, a horizontal resolution of 160 km, and 47 vertical hybrid sigma levels up to 80 km height (Crueger et al., 2018; Giorgetta et al., 2018). To avoid any effect of the model spin-up, we ignored the first year from the model results. Therefore all multiyear averages were done for the time period 2005–2009. Afterwards, the ICON-AES data are interpolated to selected pressure levels. The ICON-AES contains prognostic equations for water vapor and for cloud liquid water and cloud ice in stratiform clouds. Stratiform cloud cover is computed by using a diagnostic cloud cover scheme (Sundqvist et al., 1989). Stratiform cloud microphysics is parameterized following Lohmann and Roeckner (1996). The prognostic equation for the grid box mean cloud ice mixing ratio (${\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}$) including the different process rate terms (Q) is written as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}}{\mathrm{d}t}}& =Q_{\mathrm{Ti}}+Q_{\mathrm{sed}}+Q_{\mathrm{dep}}-Q_{\mathrm{mli}}-Q_{\mathrm{sbi}}\\ & +Q_{\mathrm{fr}}-Q_{\mathrm{saci}}-Q_{\mathrm{agg}},\end{array}\end{array}$$

including advection, parameterized turbulent diffusion, and convective detrainment of cloud ice Q_Ti; sedimentation of cloud ice Q_sed; deposition and sublimation Q_dep; melting of cloud ice Q_mli; instantaneous sublimation in the cloud-free part Q_sbi; homogeneous and heterogeneous freezing of liquid water Q_fr; accretion of cloud ice by snow Q_saci; and the aggregation of cloud ice Q_agg (Giorgetta et al., 2013). Microphysical processes are determined from in-cloud cloud ice and liquid water, respectively. The in-cloud values are calculated by dividing the grid box mean cloud ice mixing ratio and cloud liquid water mixing ratio by the fractional cloud cover (C). C is calculated by a diagnostic cloud cover scheme by Sundqvist et al. (1989) (more details in Sect. 2.2). Cloud liquid water and cloud ice are considered to be well-mixed if they coexist. Therefore, there are no separate liquid or ice parts in the cloud. The ECHAM6 physics includes diagnostic rain and snow profiles in the columns. It is not transported by advection (Giorgetta et al., 2018).

2.1 Aggregation parameterization in the ICON-AES

In this study, we focus on the aggregation parameterization. In the current version of the ICON-AES, with ECHAM6 physics, the conversion rate of ice to snow by aggregation as a nonlinear process is given by Levkov et al. (1992) based on the work of Murakami (1990), with

$$\begin{array}{}\text{(2)}& Q_{\mathrm{agg}}={\displaystyle \frac{q_{\mathrm{i}}}{\mathrm{\Delta }{t}_{\mathrm{1}}}}.\end{array}$$

Δt₁ is defined as the time needed for ice crystals to grow from the mean radius ${\stackrel{\mathrm{‾}}{R}}_{\mathrm{vi}}$ to the smallest radius of snow class particles R_S0:

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }{t}_{\mathrm{1}}=-{\displaystyle \frac{\mathrm{2}}{c_{\mathrm{1}}}}\log {\left({\displaystyle \frac{R_{\mathrm{vi}}}{R_{\mathrm{S}\mathrm{0}}}}\right)}^{\mathrm{3}},\end{array}$$

with

$$\begin{array}{}\text{(4)}& c_{\mathrm{1}}={\displaystyle \frac{q_{\mathrm{i}}\mathit{\rho }{a}_I{E}_{\mathrm{ii}}X}{{\mathit{\rho }}_{\mathrm{i}}}}{\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{0}}}{\mathit{\rho }}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}},\end{array}$$

where a_I=700 s⁻¹ is an empirical constant, E_ii=0.1 the collection efficiency, ρ_i=500 kg m⁻³ the density of cloud ice, and ρ₀=1.3 kg m⁻³ the reference density of air. Combing Eqs. (2)–(4) leads to the final aggregation rate equation:

$$\begin{array}{}\text{(5)}& Q_{\mathrm{agg}}=C\mathit{\gamma }{\displaystyle \frac{\mathit{\rho }{\stackrel{\mathrm{‾}}{q_{\mathrm{i}}}}^{\mathrm{2}}{a}_I{E}_{\mathrm{ii}}X{\left(\frac{{\mathit{\rho }}_{\mathrm{0}}}{\mathit{\rho }}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}}{-\mathrm{2}{\mathit{\rho }}_{\mathrm{i}}\log {\left(\frac{R_{\mathrm{vi}}}{R_{\mathrm{s}\mathrm{0}}}\right)}^{\mathrm{3}}}}.\end{array}$$

The process may be tuned with a tuning parameter γ, which is currently set to γ=95. The parameterization uses the grid box mean in-cloud ice mixing ratio to calculate the process rate of aggregation from ice to snow, which introduces biases in the process rates, since the aggregation is a nonlinear process (Pincus and Klein, 2000). An unbiased aggregation rate is calculated by using an integral over a distribution function of subgrid-scale cloud ice mixing ratio ($\stackrel{\mathrm{‾}}{Q_{\mathrm{agg}}\left(q_{\mathrm{i}}\right)}$).

Figure 1 Schematic overview of cloud ice variability as part of the cloud cover parameterization: the cloud cover is defined as the part of the PDF which exceeds the saturation humidity q_s up to ${\stackrel{\mathrm{‾}}{q}}_{\mathrm{t}}+\mathrm{\Delta }{q}_{\mathrm{t}}$. Within the cloudy part we define a distribution over the hydrometeors q_i+q_l. The half-width of the PDF of ${\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}$ is defined as C multiplied with Δq_t and the cloud ice fraction ($f_{\mathrm{ice}}=\frac{{\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}}{{\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}+{\stackrel{\mathrm{‾}}{q}}_{\mathrm{l}}}$).

2.2 Subgrid-scale variability of cloud ice in the aggregation process

The ICON-AES determines the fractional cloud cover according to Sundqvist et al. (1989). A uniform distribution function of the total water mixing ratio q_t from ${\stackrel{\mathrm{‾}}{q}}_{\mathrm{t}}-\mathrm{\Delta }{q}_{\mathrm{t}}$ to ${\stackrel{\mathrm{‾}}{q}}_{\mathrm{t}}+\mathrm{\Delta }{q}_{\mathrm{t}}$ is considered (Fig. 1). The total water mixing ratio describes the sum of water vapor, cloud liquid water, and cloud ice mixing ratio. The saturation specific humidity q_s is calculated from the grid box mean temperature considering the saturation with respect to ice at temperatures below 0 °C and if q_i is higher than a threshold value ${\mathit{\gamma }}_{\mathrm{thr}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}$ kg kg⁻¹. The integral over the distribution from the saturation humidity q_s up to the maximum of the distribution function ${\stackrel{\mathrm{‾}}{q}}_{\mathrm{t}}+\mathrm{\Delta }{q}_{\mathrm{t}}$ defines the fractional cloud cover C, which depends on the calculation of the distribution width Δq_t=γq_s. The parameter γ varies with height from low values near the surface to larger ones in the free troposphere with a prescribed profile (Quaas, 2012). The final equation for C, which is used in the model, is given by

$$\begin{array}{}\text{(6)}& C=\mathrm{1}-\sqrt{\mathrm{1}-{\displaystyle \frac{r-r_{\mathrm{0}}}{r_{\mathrm{sat}}-r_{\mathrm{0}}}}},\end{array}$$

where r is the relative humidity, r_sat is the saturation value (r_sat=1), and r₀ is a function of pressure and depends on two different tuning parameters (r_top=0.8 and r_surf=0.968), which defines the condensation threshold.

To define a new subgrid-scale cloud ice mass, we use this fractional cloud cover approach (Fig. 1). The half-width of the cloud ice PDF (Δq_i) has to be re-scaled with Δq_i=CΔq_tf_ice, where f_ice describes the cloud ice fraction ($\frac{{\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}}{{\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}+{\stackrel{\mathrm{‾}}{q}}_{\mathrm{l}}}$). Using the PDF over q_t the all-sky cloud ice and cloud liquid water mixing ratio is defined as the integral over q_t−q_s from q_s to the maximum of the PDF(q_t):

$$\begin{array}{}\text{(7)}& {\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}+{\stackrel{\mathrm{‾}}{q}}_{\mathrm{l}}=\underset{q_{\mathrm{s}}}{\overset{{\stackrel{\mathrm{‾}}{q}}_{\mathrm{t}}+\mathrm{\Delta }{q}_{\mathrm{t}}}{\int }}(q_{\mathrm{t}}-q_{\mathrm{s}})\mathrm{PDF}\left(q_{\mathrm{t}}\right)\mathrm{d}{q}_{\mathrm{t}}.\end{array}$$

Here, q_i and q_l are considered to be well-mixed within the cloudy part of the grid box, so q_i and q_l occur in the same volume. Solving the Eq. (7) yields

$$\begin{array}{}\text{(8)}& {\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}+{\stackrel{\mathrm{‾}}{q}}_{\mathrm{l}}=C^{\mathrm{2}}\mathrm{\Delta }{q}_{\mathrm{t}}.\end{array}$$

To get ${\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}$, Eq. (8) is multiplied with f_ice, which leads to

$$\begin{array}{}\text{(9)}& {\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}=({\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}+{\stackrel{\mathrm{‾}}{q}}_{\mathrm{l}})\cdot f_{\mathrm{ice}}=C^{\mathrm{2}}\mathrm{\Delta }{q}_{\mathrm{t}}\cdot f_{\mathrm{ice}}.\end{array}$$

As described above, the in-cloud cloud ice mixing ratio $q_{\mathrm{i}}^{\mathrm{c}}$ is defined as q_i divided by C, and in combination with Eq. (8) it follows that

$$\begin{array}{}\text{(10)}& {\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}^{\mathrm{c}}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}}{C}}=C\mathrm{\Delta }{q}_{\mathrm{t}}\cdot f_{\mathrm{ice}}=\mathrm{\Delta }{q}_{\mathrm{i}}.\end{array}$$

This implies that the width of the cloud ice distribution corresponds to the in-cloud grid box mean cloud ice mixing ratio ${\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}^c$. To obtain the subgrid-scale cloud ice mixing ratio, we use a Monte Carlo approach. Here, we choose an element from the cumulative distribution function (CDF) with the help of a random number $r\in [\mathrm{0},\mathrm{1}]$. This yields the following equation:

$$\begin{array}{}\text{(11)}& q_{\mathrm{i},\mathrm{new}}^c={\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}^c+\mathrm{\Delta }{q}_{\mathrm{i}}(\mathrm{2}r-\mathrm{1}).\end{array}$$

With Eq. (10),

$$\begin{array}{}\text{(12)}& q_{\mathrm{i},\mathrm{new}}^c={\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}^c\left(\mathrm{2}r\right).\end{array}$$

Finally, the subgrid-scale cloud ice mixing ratio only depends on the grid box mean cloud ice mass mixing ratio and the choice of the random number. This new specific cloud ice mass replaces the grid box mean cloud ice in Eq. (5). To compare the current, biased aggregation rate ($Q_{\mathrm{agg}}\left({\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}\right)$) with the unbiased process rate ($\stackrel{\mathrm{‾}}{Q_{\mathrm{agg}}\left(q_{\mathrm{i}}\right)}$) we replace for each time step and grid box ${\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}$ mean in Eq. (5) with the integral over the entire distribution of q_i. This comparison will be shown in the last part of the Results section.

To evaluate the distribution function, the cloud ice variance is calculated for all-sky conditions because all output cloud variables are also all-sky variables. The variance can be written as follows:

$$\begin{array}{}\text{(13)}& {\mathit{\sigma }}_{q_{\mathrm{i}}}^{\mathrm{2}}=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathrm{\Delta }{q}_{\mathrm{i}}}{\int }}(q_{\mathrm{i}}-{\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}{)}^{\mathrm{2}}\text{PDF}(q_{\mathrm{i}})\mathrm{d}{q}_{\mathrm{i}}.\end{array}$$

As we focus on all-sky variance, this equation can be separated into two parts: the clear-sky part of the distribution (q_i=0) and the cloudy part (q_i>0). The cloud-free part can be described with Dirac's delta function, which yields the following equation (derivation is shown in the Appendix):

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{\mathit{\sigma }}_{q_{\mathrm{i}}}^{\mathrm{2}}& =\underset{\mathrm{0}}{\overset{\mathrm{2}\mathrm{\Delta }{q}_{\mathrm{i}}}{\int }}\left[(\mathrm{1}-C)\mathit{\delta }(q_{\mathrm{i}}=\mathrm{0})+C{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }{q}_{\mathrm{i}}}}\right](q_{\mathrm{i}}-{\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}{)}^{\mathrm{2}}\mathrm{d}{q}_{\mathrm{i}}\\ {\mathit{\sigma }}_{q_{\mathrm{i}}}^{\mathrm{2}}& =(q_{\mathrm{i}}^{\mathrm{c}}{)}^{\mathrm{2}}\cdot \left({\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}C-C^{\mathrm{2}}\right).\end{array}\end{array}$$

2.3 Satellite retrievals of cloud ice water content

As mentioned above, we will focus on ice-phase processes and their impact on the global ice water content. To evaluate the results, a combined global ice water product of the CloudSat (Stephens et al., 2002) and Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) (Winker et al., 2010), the DARDAR-Nice dataset (Sourdeval et al., 2018), is used. Both satellites are part of the A-train constellation and are flying with a time interval of only 15 s between them. The W-band (94 GHz) cloud-profiling radar (CPR) on CloudSat provides vertical cloud profiles with a minimum detectable reflectivity of −28 dBz. CALIPSO contains the Cloud–Aerosol Lidar with Orthogonal Polarization (CALIOP), which measures co-aligned pulses of 532 and 1064 nm wavelength. Both satellite instruments provide global retrievals of clouds and their properties. A combined dataset of lidar and radar measurements provides global information about larger as well as smaller ice particles and therefore a more detailed overview of the cloud. Outside the temperature range of −40 to 0 °C, supercooled liquid water does not exist. The time period used from the dataset is January 2007 until December 2010, since the datasets are available for the entire years.

Figure 2 Zonally averaged annual mean ice water content (mg kg⁻¹) from the DARDAR dataset. (a) Total ice mixing ratio (q_i,total), which includes cloud ice from any clouds, and precipitating ice. (b) Cloud ice mixing ratio (q_i), where precipitating ice and convective cloud ice are removed. (c) Difference between q_i,total and q_i.

For an accurate comparison between satellite data and model data the two-dimensional DARDAR data have to be interpolated to the same horizontal grid as the ICON-AES. For each grid box, the mean of all variables is calculated and vertically interpolated to the same pressure levels as for the ICON-AES. The modeled q_i, which is used for the aggregation parameterization, does not include precipitating and convective cloud ice. Therefore, an estimate of the convective and precipitating cloud ice has to be removed from the satellite dataset. To provide a dataset that is comparable with the ICON-AES, a flag method after Li et al. (2008, 2012) is used. All columns that are flagged as “precipitating at the surface” are not considered to ensure that larger falling ice particles are removed from the dataset. This is because the model distinguishes between precipitating ice and cloud ice in two bulk classes. In contrast, it is challenging for the satellite retrievals to separate falling ice and cloud ice. Additionally, cloud ice, which is flagged as originating from “deep convection” or “cumulus” (2B-CLDCLASS dataset), is also removed. In order to compare modeled cloud ice with observations, alternative methods are possible (e.g., removing ICON nonzero surface precipitation points). Since we want to evaluate the cloud ice distribution that is used for the aggregation, we had to adjust the DARDAR data in order to find the most consistent way.

Cloud variables simulated by the ICON-AES are output for grid box mean all-sky conditions. Therefore, grid box mean cloud ice and cloud ice variance are calculated for all-sky conditions in the observational data as well.

Figure 2 shows a multiyear mean of zonally averaged total ice mixing ratio (q_i,total) and cloud ice mixing ratio (q_i) obtained from the DARDAR product. q_i,total shows the highest values over the Equator at lower pressure fields and over the midlatitudes at higher pressure fields. Most of the removed cloud ice is detected due to the precipitation flag. As a result, the remaining cloud ice water content (Fig. 2b) is much lower, especially in the midlatitudes at lower regions and in the tropics. The maxima of cloud ice are shifted upward, since most of the precipitating cloud ice is below or in the lower regions of the cloud. The plot highlights the importance of removing precipitation from the DARDAR dataset in order to give a more reliable comparison with the ICON-AES. This is due to the fact that the remaining q_i is much less than q_i,total. But we have to keep in mind that the result of q_i strongly depends on the accuracy of the different flags because it selects which cloud ice is detected as convective or precipitating cloud ice and therefore which cloud ice is removed from the dataset. Additionally, the initial satellite data are measured on a 1D curtain, while the model uses 3D grid boxes. Hill et al. (2015) calculated a measure of the difference in standard deviation considering a 2D cloud field compared to a 1D cross-section. They estimated a 30 % larger standard deviation in 2D fields compared to the 1D track. Therefore, we should consider that this also has an effect on our the cloud ice variance calculation. However, there are limitations in availability of satellite data. Therefore we tried to be as consistent as possible in the comparison between simulations and observations.

3 Results

3.1 Employing a distribution of cloud ice in the aggregation parameterization of the ICON-AES

As described above, for the evaluation of the cloud ice distribution the all-sky cloud ice variance is calculated for each grid box. Equation (14) was used for the modeled cloud ice variance, while for the DARDAR data the spatial variance within each GCM grid box was calculated from all footprints within a grid box. But before comparing the cloud ice variance, we have to make sure that the cloud ice mixing ratios from DARDAR data and ICON-AES show the same order of magnitude. Figure 3 shows the annually averaged global distribution of cloud ice from DARDAR and ICON as well as the differences between them at different pressure levels. ICON-AES shows higher values in the midlatitudes at higher pressure fields. At 400 hPa the main difference is around the Equator, where DARDAR shows higher values than ICON-AES, ICON-AES still has maximum values in the midlatitudes. However, there is good agreement between the pattern of cloud ice mixing ratio from ICON and DARDAR.

Figure 3 Multiyear mean of the cloud ice mixing ratio (kg kg⁻¹) at four different pressure levels calculated for the DARDAR data (a, d, g, j), the ICON data (b, e, h, k), and the difference between ICON and DARDAR (c, f, i, j).

Figure 4 Multiyear cloud ice water path (CIWP) for the DARDAR data and ICON data as well as the difference between DARDAR and ICON.

The cloud ice water path (CIWP), which is the column-integrated ice water content, is given in Fig. 4. ICON overestimates the CIWP in the middle and higher latitudes, while it underestimates the CIWP in the tropics. As is already visible in Fig. 3, DARDAR shows larger cloud ice values down the lower altitudes over the tropics, which leads to an increased CIWP compared to the modeled CIWP. Especially in the midlatitudes at higher pressure fields, the model tends to overestimate the cloud ice. One should keep in mind that the way the DARDAR data are filtered to get the CIWP or q_i, which is comparable with the model data, is not perfect.

Figure 5 Same as Fig. 3 but for cloud ice variance (kg² kg⁻²).

Figure 5 shows the comparison of cloud ice variance between ICON-AES and DARDAR data for the same pressure levels chosen in Fig. 3. In general, the pattern of cloud ice variance follows the pattern of the global cloud ice distribution. Higher values are visible over the tropics at higher altitude and over the storm tracks in the midlatitudes at lower altitudes. Besides the good agreement in the distribution pattern, there are some regional differences. At 300 hPa, the maxima over the northern part of North America and central Africa are more intense in the DARDAR data than in the ICON data. At 400 hPa ICON underestimates the cloud ice variance compared to DARDAR. At lower altitudes the cloud ice variance is higher in the ICON data in most of the midlatitudes compared to the DARDAR data. ICON underestimates the cloud ice variance in the tropics, especially in lower regions, since ICON shows less cloud ice in the tropical region at the same pressure levels. However, we should keep in mind that the measured variance contains uncertainties regarding the method of filtering precipitation and convection. Additionally the modeled cloud ice variance makes an assumption of a distribution, while the DARDAR data show the variance from discrete measurements. Despite these discrepancies, we conclude that the simple uniform distribution of cloud ice as written in Eq. (14) basically captures the measured distribution of cloud ice in each grid box. Therefore it is suitable for use in the aggregation parameterization.

Figure 6 Zonally averaged multiyear mean cloud ice (mg kg⁻¹). (a) Control run of the ICON-AES. (b) Difference between the new stochastic aggregation parameterization and the control run (adapted from Trömel et al., 2021).

As described in Sect. 2.2, with the help of this cloud ice distribution a new specific cloud ice mass is used for the aggregation parameterization. Figure 6 shows its influence on the zonally averaged cloud ice adapted from Trömel et al. (2021), where the same tuning parameters were used. The control run of the ICON model (CTRL) gives a maximum of cloud ice at levels of lower pressure over the tropics. There are also two maxima over the midlatitudes between 600 and 400 hPa, with a more pronounced maximum in the Southern Hemisphere than in the Northern Hemisphere. An overall reduction of cloud ice is visible due to the stochastic aggregation scheme (AGGstoch).

Since aggregation is a nonlinear process, it becomes stronger due to taking a distribution of cloud ice into account. A stronger aggregation process leads to a higher conversion rate from cloud ice to snow. Therefore, more cloud ice is removed due to the process. Figure 7 shows the effect of using the stochastic approach instead of the default parameterization directly for the aggregation rate in a joint histogram and the corresponding multiyear zonally averaged aggregation rate. While the CTRL simulation produces an intense maximum aggregation rate around ca. 30 mg kg⁻¹ h⁻¹, the AGGstoch run reveals a larger spread around a less intense maximum. The higher spread of the maximum stems from the randomly picked cloud ice mass for the aggregation process. Therefore, there are values which are much higher and much lower than the original grid box mean cloud ice mass. Instead of having this intense maximum aggregation rate, higher and lower aggregation rate cases are visible. Therefore, the annual average over all of these aggregation rate cases yields a higher value for the AGGstoch run compared to CTRL run, caused by the nonlinearity of the process, which is also visible in Fig. 7f.

Figure 7 Joint histograms of aggregation rate (mg kg⁻¹ h⁻¹) and cloud ice (mg kg⁻¹) for CTRL (a), AGGstoch (b), and the difference between the two (c). Zonally averaged aggregation rate of the multiyear mean for CTRL (d), AGGstoch (e), and the difference (f).

Figure 8 Vertical profiles of the global mean process rates which are related to cloud ice microphysical processes for CTRL (solid lines) and AGGstoch (dashed lines): (Aci), freezing rate (Frz), deposition rate (Dep), melting rate (Mlt), evaporation rate (Evp), sedimentation rate (Sed), aggregation rate (Agg), and accretion rate (transport terms are not included).

The microphysical processes in clouds are strongly interlinked with each other and react whenever one process rate changes. Figure 8 shows how the different cloud-ice-related process rates change using the stochastic aggregation parameterization. The different vertical profiles of the global mean process rates are depicted. Negative process rates indicate a cloud ice loss, while positive process rates lead to a cloud ice gain. So, a more negative process rate is linked with a stronger cloud ice loss. It is visible that the deposition rate is the most effective cloud-ice-gaining process. It is triggered by the Wegener–Bergeron–Findeisen (WBF) (Wegener, 1911; Bergeron, 1935; Findeisen et al., 2015) process. This process takes place if the temperature is between −35 and 0 °C and the q_i exceeds a specific threshold (Giorgetta et al., 2013). The strongest negative process rate is the accretion rate. The accretion rate (Aci) describes the growth of snow by collecting the surrounding ice crystals. Therefore, the accretion rate is strongly linked with aggregation rate, since the aggregation process produces the snow first. Aggregation and accretion describe the formation and growing of snow by coagulation of cloud ice, and thus they lead to a cloud ice reduction. The riming process (not shown here) is described like the linear accretion process and leads to an increase in snowfall, but with the difference that the equation depends on cloud water. However, the accretion rate is the more efficient process compared to aggregation. Figure 8 confirms the higher cloud ice loss in the AGGstoch compared to the CTRL simulation due to the more negative aggregation rate. Moreover, due to the stronger aggregation rate at levels of high pressure, less cloud ice may be transported into lower regions. Hence, the sedimentation process slightly decreases. Due to the change in the aggregation process rate, the accretion rate becomes less intense and leads to more cloud ice. The reduced accretion rate may be explained by the cloud ice loss due the stronger aggregation rate. Less cloud ice is collectable for growing snowflakes. No significant changes are visible for freezing, deposition, melting, and evaporation rate. Due to the stochastic method, the global mean aggregation rate is intensified by more than 20 % in the middle and upper troposphere. In contrast, the maximum change in global mean accretion rate is less than 8.9 %. Since the accretion rate is much higher, this smaller change leads to a compensation of the aggregation rate increase as described above. However, the increase in aggregation rate by more than 20 % is significant. Due to the small change in the microphysical properties, no important change in multiyear global mean shortwave (+0.161 W m⁻²), longwave (−0.195 W m⁻²), and net radiation (+0.03 W m⁻²) at the top of the atmosphere is visible (not shown here). Overall, we found a cloud ice loss in the AGGstoch run of up to 5 %, but the reduction of cloud ice is compensated for by a less intense accretion rate. The effect on radiation could be increased if this stochastic approach was implemented in other processes, since there has to be a change of a factor of 2 or more in microphysical process rates to see an effect on radiation (Michibata et al., 2020; Imura and Michibata, 2022). Using the new approach just for one single process makes it easier in the beginning to see the effect of changing one process, since all processes are connected. Besides that, the aggregation process is the only nonlinear cloud ice process rate in the ICON-AES. Since we focused on cloud-ice-related processes we just implemented the subgrid-scale approach in the aggregation parameterization. In future studies one can think about additionally including the subgrid-scale approach in cloud-water-related processes in order to see a stronger effect on radiation, cloudiness, and precipitation as well.

Figure 9 Averaged aggregation rate at different cloud ice mixing ratio bins and the corresponding 2D histograms for CTRL, AGGstoch, and AGGsample.

As already mentioned, using a stochastic approach in nonlinear process rates lowers the bias of the process rate. To create an unbiased aggregation rate (AGGsample), as described above, we make use of the entire cloud ice distribution function. Figure 9 shows a comparison of the averaged aggregation process rate at different q_i bins and the 2D histograms of CTRL, AGGstoch, and AGGsample. AGGstoch shows good agreement with AGGsample, while CTRL produces much lower process rates. Both AGGstoch and AGGsample produce a higher aggregation rate compared to CTRL, especially for higher q_i values. Figure 9b gives the corresponding joint histogram of CTRL and (c)–(d) the difference histogram to CTRL of cloud ice and the aggregation rate. AGGstoch shows a higher spread in the distribution, as already mentioned, while in AGGsample the distribution is shifted towards higher aggregation rates. However, both methods produce higher aggregation rates on average for different q_i sections. The main difference between these methods is that AGGsample needs an additional integration over q_i for the calculation of the aggregation rate, while AGGstoch needs no extra computational time. Therefore, the described stochastic method, which is used in AGGstoch, allows an improvement of the process rate bias against CTRL.

4 Summary

We introduce a stochastic approach of the aggregation parameterization in the ICON-AES by considering, consistently with the models' cloud scheme, a uniform distribution function of cloud ice and randomly choose a cloud ice water content from this distribution function for each grid box and time step. The chosen distribution function of cloud ice is evaluated with the help of a combined lidar–radar dataset (DARDAR). An estimate of the precipitating and convective cloud ice mass is removed from the dataset in order to allow a more consistent comparison to the cloud ice in the model. The comparison of the simulated and the observed cloud ice variance shows good agreement in pattern, with just a few regional differences.

Overall we show that the uniform distribution function of cloud ice is usable for the microphysical process rates, e.g., the aggregation. From this uniform distribution function a randomly chosen cloud ice value is implemented into the nonlinear aggregation rate in order to represent subgrid-scale variability. Due to this stochastic method, the aggregation rate is intensified on average, since aggregation is a nonlinear process. As a result, more cloud ice is transformed to snow, which leads to a cloud ice loss. However, the decrease in the accretion rate that results from using the stochastic aggregation scheme acts against the more intense aggregation rate. Therefore, the cloud ice loss is not as strong as expected. This indicates that changing only one of the two process rates of snow formation does not lead to a change as large as one might have initially expected because aggregation and accretion interact very strongly with each other. However, the effect of taking subgrid-scale variability into account for process rates has a significant impact on the microphysical process (20 % stronger global averaged process rate).

An unbiased process rate is calculated by integrating the aggregation rate over the entire cloud ice distribution function. The new stochastic method shows better agreement of the aggregation rates with the unbiased method compared to the control run. It follows that the new method lowers the bias of the aggregation rate, which does not need additional computational time. Therefore, this study suggests that using a stochastic approach for microphysical process rates helps to improve the representation of clouds and precipitation processes in global climate models.

Appendix A: Variance of cloud ice

Here, the calculation of cloud ice variance, which is used in the ICON-AES, is given step by step. The general equation is written as follows:

$$\begin{array}{}\text{(A1)}& {\mathit{\sigma }}_{q_{\mathrm{i}}}^{\mathrm{2}}=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathrm{\Delta }{q}_{\mathrm{i}}}{\int }}(q_{\mathrm{i}}-{\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}{)}^{\mathrm{2}}\text{PDF}(q_{\mathrm{i}})\mathrm{d}{q}_{\mathrm{i}}.\end{array}$$

Separating the integral in the cloud-free part, which is expressed as a Dirac function, and the cloudy parts yields

$$\begin{array}{}\text{(A2)}& {\mathit{\sigma }}_{q_{\mathrm{i}}}^{\mathrm{2}}=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathrm{\Delta }{q}_{\mathrm{i}}}{\int }}[\underset{\text{cloud-free part}}{\underbrace{(\mathrm{1}-C)\mathit{\delta }(q_{\mathrm{i}}=\mathrm{0})}}+\underset{\text{cloudy part}}{\underbrace{C{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }{q}_{\mathrm{i}}}}}}](q_{\mathrm{i}}-{\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}{)}^{\mathrm{2}}\mathrm{d}{q}_{\mathrm{i}}.\end{array}$$

Solving the Dirac function of the cloud-free part yields

$$\begin{array}{}\text{(A3)}& {\mathit{\sigma }}_{q_{\mathrm{i}}}^{\mathrm{2}}=(\mathrm{1}-C){\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}^{\mathrm{2}}+{\displaystyle \frac{C}{\mathrm{2}\mathrm{\Delta }{q}_{\mathrm{i}}}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathrm{\Delta }{q}_{\mathrm{i}}}{\int }}(q_{\mathrm{i}}-{\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}{)}^{\mathrm{2}}\mathrm{d}{q}_{\mathrm{i}}.\end{array}$$

Solving the integral with ${\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}=C\mathrm{\Delta }{q}_{\mathrm{i}}$ yields the following (see Eq. 8).

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}{\mathit{\sigma }}_{q_{\mathrm{i}}}^{\mathrm{2}}& =(\mathrm{1}-C){\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}^{\mathrm{2}}+{\displaystyle \frac{C}{\mathrm{6}\mathrm{\Delta }{q}_{\mathrm{i}}}}\left[(\mathrm{2}\mathrm{\Delta }{q}_{\mathrm{i}}-{\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}{)}^{\mathrm{3}}+{\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}^{\mathrm{3}}\right]\\ {\mathit{\sigma }}_{q_{\mathrm{i}}}^{\mathrm{2}}& =(\mathrm{1}-C)(C\mathrm{\Delta }{q}_{\mathrm{i}}{)}^{\mathrm{2}}\\ & +{\displaystyle \frac{C}{\mathrm{6}\mathrm{\Delta }{q}_{\mathrm{i}}}}\left[(\mathrm{2}\mathrm{\Delta }{q}_{\mathrm{i}}-\mathrm{C}\mathrm{\Delta }{q}_{\mathrm{i}}{)}^{\mathrm{3}}+(C\mathrm{\Delta }{q}_{\mathrm{i}}{)}^{\mathrm{3}}\right]\\ {\mathit{\sigma }}_{q_{\mathrm{i}}}^{\mathrm{2}}& =(\mathrm{1}-C)(C\mathrm{\Delta }{q}_{\mathrm{i}}{)}^{\mathrm{2}}\\ & +{\displaystyle \frac{C}{\mathrm{6}\mathrm{\Delta }{q}_{\mathrm{i}}}}\left[\mathrm{8}\mathrm{\Delta }{q}_{\mathrm{i}}^{\mathrm{3}}-\mathrm{12}C\mathrm{\Delta }{q}_{\mathrm{i}}^{\mathrm{3}}+\mathrm{6}{C}^{\mathrm{2}}\mathrm{\Delta }{q}_{\mathrm{i}}^{\mathrm{3}}\right]\\ {\mathit{\sigma }}_{q_{\mathrm{i}}}^{\mathrm{2}}& =(\mathrm{\Delta }{q}_{\mathrm{i}}{)}^{\mathrm{2}}\cdot \left({\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}C-C^{\mathrm{2}}\right)\end{array}\end{array}$$

Replacing Δq_i with Eq. (10) yields

$$\begin{array}{}\text{(A5)}& {\mathit{\sigma }}_{{\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}}^{\mathrm{2}}=({\stackrel{\mathrm{‾}}{q}}_{\mathrm{i}}^{\mathrm{c}}{)}^{\mathrm{2}}\cdot \left({\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}C-C^{\mathrm{2}}\right).\end{array}$$

In the end, the variance of cloud ice just depends on the cloud cover and grid box mean of in-cloud ice. So, the variance, which is calculated for the aggregation process, can be checked directly for the specific amount of cloud ice available for the aggregation.