Articles | Volume 17, issue 8
Development and technical paper
19 Apr 2024
Development and technical paper |  | 19 Apr 2024

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

Sabine Doktorowski, Jan Kretzschmar, Johannes Quaas, Marc Salzmann, and Odran Sourdeval

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 Klein2000; 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., Palmer2001; 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 Gettelman2008; Larson and Griffin2013; 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 Hogan2008, 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 (qi) including the different process rate terms (Q) is written as follows:

(1) d q i d t = Q Ti + Q sed + Q dep - Q mli - Q sbi + Q fr - Q saci - Q agg ,

including advection, parameterized turbulent diffusion, and convective detrainment of cloud ice QTi; sedimentation of cloud ice Qsed; deposition and sublimation Qdep; melting of cloud ice Qmli; instantaneous sublimation in the cloud-free part Qsbi; homogeneous and heterogeneous freezing of liquid water Qfr; accretion of cloud ice by snow Qsaci; and the aggregation of cloud ice Qagg (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

(2) Q agg = q i Δ t 1 .

Δt1 is defined as the time needed for ice crystals to grow from the mean radius Rvi to the smallest radius of snow class particles RS0:

(3) Δ t 1 = - 2 c 1 log R vi R S 0 3 ,


(4) c 1 = q i ρ a I E ii X ρ i ρ 0 ρ 1 3 ,

where aI=700 s−1 is an empirical constant, Eii=0.1 the collection efficiency, ρi=500 kg m−3 the density of cloud ice, and ρ0=1.3 kg m−3 the reference density of air. Combing Eqs. (2)–(4) leads to the final aggregation rate equation:

(5) Q agg = C γ ρ q i 2 a I E ii X ρ 0 ρ 1 3 - 2 ρ i log R vi R s 0 3 .

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 Klein2000). An unbiased aggregation rate is calculated by using an integral over a distribution function of subgrid-scale cloud ice mixing ratio (Qagg(qi)).

Figure 1Schematic 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 qs up to qt+Δqt. Within the cloudy part we define a distribution over the hydrometeors qi+ql. The half-width of the PDF of qi is defined as C multiplied with Δqt and the cloud ice fraction (fice=qiqi+ql).


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 qt from qt-Δqt to qt+Δqt 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 qs is calculated from the grid box mean temperature considering the saturation with respect to ice at temperatures below 0 °C and if qi is higher than a threshold value γthr=5×10-7 kg kg−1. The integral over the distribution from the saturation humidity qs up to the maximum of the distribution function qt+Δqt defines the fractional cloud cover C, which depends on the calculation of the distribution width Δqt=γqs. The parameter γ varies with height from low values near the surface to larger ones in the free troposphere with a prescribed profile (Quaas2012). The final equation for C, which is used in the model, is given by

(6) C = 1 - 1 - r - r 0 r sat - r 0 ,

where r is the relative humidity, rsat is the saturation value (rsat=1), and r0 is a function of pressure and depends on two different tuning parameters (rtop=0.8 and rsurf=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 (Δqi) has to be re-scaled with Δqi=CΔqtfice, where fice describes the cloud ice fraction (qiqi+ql). Using the PDF over qt the all-sky cloud ice and cloud liquid water mixing ratio is defined as the integral over qtqs from qs to the maximum of the PDF(qt):

(7) q i + q l = q s q t + Δ q t ( q t - q s ) PDF ( q t ) d q t .

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

(8) q i + q l = C 2 Δ q t .

To get qi, Eq. (8) is multiplied with fice, which leads to

(9) q i = ( q i + q l ) f ice = C 2 Δ q t f ice .

As described above, the in-cloud cloud ice mixing ratio qic is defined as qi divided by C, and in combination with Eq. (8) it follows that

(10) q i c = q i C = C Δ q t f ice = Δ q i .

This implies that the width of the cloud ice distribution corresponds to the in-cloud grid box mean cloud ice mixing ratio qic. 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[0,1]. This yields the following equation:

(11) q i , new c = q i c + Δ q i ( 2 r - 1 ) .

With Eq. (10),

(12) q i , new c = q i c ( 2 r ) .

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 (Qagg(qi)) with the unbiased process rate (Qagg(qi)) we replace for each time step and grid box qi mean in Eq. (5) with the integral over the entire distribution of qi. 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:

(13) σ q i 2 = 0 2 Δ q i ( q i - q i ) 2 PDF ( q i ) d q i .

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

(14) σ q i 2 = 0 2 Δ q i ( 1 - C ) δ ( q i = 0 ) + C 1 2 Δ q i ( q i - q i ) 2 d q i σ q i 2 = ( q i c ) 2 4 3 C - C 2 .

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 2Zonally averaged annual mean ice water content (mg kg−1) from the DARDAR dataset. (a) Total ice mixing ratio (qi,total), which includes cloud ice from any clouds, and precipitating ice. (b) Cloud ice mixing ratio (qi), where precipitating ice and convective cloud ice are removed. (c) Difference between qi,total and qi.


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 qi, 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 (qi,total) and cloud ice mixing ratio (qi) obtained from the DARDAR product. qi,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 qi is much less than qi,total. But we have to keep in mind that the result of qi 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 3Multiyear mean of the cloud ice mixing ratio (kg kg−1) 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 4Multiyear 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 qi, which is comparable with the model data, is not perfect.

Figure 5Same as Fig. 3 but for cloud ice variance (kg2 kg−2).


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 6Zonally averaged multiyear mean cloud ice (mg kg−1). (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−1 h−1, 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 7Joint histograms of aggregation rate (mg kg−1 h−1) and cloud ice (mg kg−1) 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 8Vertical 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) (Wegener1911; Bergeron1935; Findeisen et al.2015) process. This process takes place if the temperature is between 35 and 0 °C and the qi 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−2), longwave (0.195 W m−2), and net radiation (+0.03 W m−2) 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 Michibata2022). 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 9Averaged 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 qi 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 qi 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 qi sections. The main difference between these methods is that AGGsample needs an additional integration over qi 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:

(A1) σ q i 2 = 0 2 Δ q i ( q i - q i ) 2 PDF ( q i ) d q i .

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

(A2) σ q i 2 = 0 2 Δ q i [ ( 1 - C ) δ ( q i = 0 ) cloud-free part + C 1 2 Δ q i cloudy part ] ( q i - q i ) 2 d q i .

Solving the Dirac function of the cloud-free part yields

(A3) σ q i 2 = ( 1 - C ) q i 2 + C 2 Δ q i 0 2 Δ q i ( q i - q i ) 2 d q i .

Solving the integral with qi=CΔqi yields the following (see Eq. 8).

(A4) σ q i 2 = ( 1 - C ) q i 2 + C 6 Δ q i ( 2 Δ q i - q i ) 3 + q i 3 σ q i 2 = ( 1 - C ) ( C Δ q i ) 2 + C 6 Δ q i ( 2 Δ q i - C Δ q i ) 3 + ( C Δ q i ) 3 σ q i 2 = ( 1 - C ) ( C Δ q i ) 2 + C 6 Δ q i 8 Δ q i 3 - 12 C Δ q i 3 + 6 C 2 Δ q i 3 σ q i 2 = ( Δ q i ) 2 4 3 C - C 2

Replacing Δqi with Eq. (10) yields

(A5) σ q i 2 = ( q i c ) 2 4 3 C - C 2 .

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.

Code availability

The source code of this study is available at (Doktorowski2024). Updates of the code can be found at (last access: 16 April 2024). Information about the license and the respective README files are also included. Access to ICON source code can be obtained at (ICON partnership2024).

Data availability

The 2B-CLDCLASS (Sassen and Wang2008) CloudSat data and DARDAR-Nice dataset (, Sourdeval et al.2018) are downloaded from the ICARE Data and Services Center Archive (, last access: 16 April 2024), which requires a user account.

Author contributions

SD and JQ conceived this study. JK provided expertise on how to run the model. OS gave advise about the satellite data. JK, JQ, and MS helped by checking the theoretical approach of this study. SD prepared the article with contributions from all co-authors.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


This work used resources of the German Climate Computing Centre (Deutsches Klimarechenzentrum, DKRZ) granted by its Scientific Steering Committee (WLA) under project ID 1143. We thank the AERIS/ICARE Data and Services Center for providing access to the DARDAR and 2B-CLDCLASS data used in this study. We thank Takuro Michibata and an anonymous reviewer for comments that led to improvements in the analysis and paper.

Financial support

This research received funding from the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) via its priority program “Fusion of Radar Polarimetry and Atmospheric Modelling” (SPP-2115, PROM), specifically the project PARA (FKZ QU 311/21-1). Johannes Quaas further acknowledges funding by the EU Horizon 2020 project “FORCES” (GA 821205).

Review statement

This paper was edited by Sylwester Arabas and reviewed by Takuro Michibata and one anonymous referee.


Bergeron, T.: On the physics of clouds and precipitation, Proc. 5th Assembly UGGI, Lisbon, Portugal, 156–180, (last access: 16 April 2024), 1935. a

Berner, J., Achatz, U., Batté, L., Bengtsson, L., de la Cámara, A., Christensen, H. M., Colangeli, M., Coleman, D. R. B., Crommelin, D., Dolaptchiev, S. I., Franzke, C. L. E., Friederichs, P., Imkeller, P., Järvinen, H., Juricke, S., Kitsios, V., Lott, F., Lucarini, V., Mahajan, S.,Palmer, T. N., Penland, C., Sakradzija, M., von Storch, J.-S., Weisheimer, A., Weniger, M., and Williams, P. D.: Stochastic parameterization: Toward a new view of weather and climate models, B. Am. Math. Soc., 98, 565–588,, 2017. a

Boutle, I. A., Abel, S. J., Hill, P. G., and Morcrette, C. J.: Spatial variability of liquid cloud and rain: observations and microphysical effects, Q. J. Roy. Meteor. Soc., 140, 583–594,, 2014. a, b

Crueger, T., Giorgetta, M. A., Brokopf, R., Esch, M., Fiedler, S., Hohenegger, C., Kornblueh, L., Mauritsen, T., Nam, C., Naumann, A. K., Peters, K., Rast, S., Roeckner, E., Sakradzija, M., Schmidt, H., Vial, J., Vogel, R., and Stevens, B.: ICON-A, The Atmosphere Component of the ICON Earth System Model: II. Model Evaluation, J. Adv. Model Earth Sy., 10, 1638–1662,, 2018. a

Delanoë, J. and Hogan, R. J.: A variational scheme for retrieving ice cloud properties from combined radar, lidar, and infrared radiometer, J. Geophys. Res., 113, D07204,, 2008. a

Delanoë, J. and Hogan, R. J.: Combined CloudSat-CALIPSO-MODIS retrievals of the properties of ice clouds, J. Geophys. Res., 115, D00H29,, 2010. a

Doktorowski, S.: Subgrid-scale variability of cloud ice in the ICON-AES 1.3.00, Zenodo [code],, 2024. a

Findeisen, W., Volken, E., Giesche, A. M., and Brönnimann, S.: Colloidal meteorological processes in the formation of precipitation, Meteorol. Z., 24, 443–454,, 2015. a

Giorgetta, M. A., Roeckner, E., Mauritsen, T., Bader, J., Crueger, T., Esch, M., Rast, S., Kornblueh, L., Schmidt, H., Kinne, S., Hohenegger, C., Möbis, B., Krismer, T., Wieners, K.-H., and Stevens, B.: The atmospheric general circulation model ECHAM6-model description, Berichte zur Erdsystemforschung/Max-Planck-Institut für Meteorologie, 135, 172 pp.,, 2013. a, b

Giorgetta, M. A., Brokopf, R., Crueger, T., Esch, M., Fiedler, S., Helmert, J., Hohenegger, C., Kornblueh, L., Köhler, M., Manzini, E., Mauritsen, T., Nam, C., Raddatz, T., Rast, S., Reinert, D., Sakradzija, M., Schmidt, H., Schneck, R., Schnur, R., Silvers, L., Wan, H., Zängl, G., and Stevens, B.: ICON-A, the Atmosphere Component of the ICON Earth System Model: I. Model Description, J. Adv. Model Earth Sy., 10, 1613–1637,, 2018. a, b, c, d

Golaz, J.-C., Larson, V. E., and Cotton, W. R.: A PDF-based model for boundary layer clouds. Part I: Method and model description, J. Atmos. Sci., 59, 3540–3551,<3540:APBMFB>2.0.CO;2, 2002a. a

Golaz, J.-C., Larson, V. E., and Cotton, W. R.: A PDF-based model for boundary layer clouds. Part II: Model results, J. Atmos. Sci., 59, 3552–3571,<3552:APBMFB>2.0.CO;2, 2002b. a

Hill, P. G., Morcrette, C. J., and Boutle, I. A.: A regime-dependent parametrization of subgrid-scale cloud water content variability, Q. J. Roy. Meteor. Soc., 141, 1975–1986,, 2015. a

ICON partnership (DWD; MPI-M; DKRZ; KIT; C2SM): ICON release 2024.01, World Data Center for Climate (WDCC) at DKRZ [code],, 2024. a

Imura, Y. and Michibata, T.: Too Frequent and Too Light Arctic Snowfall With Incorrect Precipitation Phase Partitioning in the MIROC6 GCM, J. Adv. Model. Earth Sy., 14, e2022MS003046,, 2022. a

Larson, V. E. and Griffin, B. M.: Analytic upscaling of a local microphysics scheme. Part I: Derivation, Q. J. Roy. Meteor. Soc., 139, 46–57,, 2013. a

Larson, V. E., Wood, R., Field, P. R., Golaz, J.-C., Haar, T. H. V., and Cotton, W. R.: Systematic biases in the microphysics and thermodynamics of numerical models that ignore subgrid-scale variability, J. Atmos. Sci., 58, 1117–1128,<1117:SBITMA>2.0.CO;2, 2001. a

Lebsock, M., Morrison, H., and Gettelman, A.: Microphysical implications of cloud-precipitation covariance derived from satellite remote sensing, J. Geophys. Res.-Atmos., 118, 6521–6533,, 2013. a, b

Levkov, L., Rockel, B., Kapitza, H., and Raschke, E.: 3D mesoscale numerical studies of cirrus and stratus clouds by their time and space evolution, Beitr. Phys. Atmosph., 65, 35–58, 1992. a

Li, J.-L., Waliser, D., Chen, W.-T., Guan, B., Kubar, T., Stephens, G., Ma, H.-Y., Deng, M., Donner, L., Seman, C., and Horowitz, L.: An observationally based evaluation of cloud ice water in CMIP3 and CMIP5 GCMs and contemporary reanalyses using contemporary satellite data, J. Geophys. Res., 117, D16105,, 2012. a, b

Li, J.-L. F., Waliser, D., Woods, C., Teixeira, J., Bacmeister, J., Chern, J., Shen, B.-W., Tompkins, A., Tao, W.-K., and Köhler, M.: Comparisons of satellites liquid water estimates to ECMWF and GMAO analyses, 20th century IPCC AR4 climate simulations, and GCM simulations, Geophys. Res. Lett., 35, L19710,, 2008. a, b

Lohmann, U. and Roeckner, E.: Design and performance of a new cloud microphysics scheme developed for the ECHAM general circulation model, Clim. Dynam., 12, 557–572,, 1996. a

Michibata, T., Suzuki, K., and Takemura, T.: Snow-induced buffering in aerosol–cloud interactions, Atmos. Chem. Phys., 20, 13771–13780,, 2020. a

Morrison, H. and Gettelman, A.: A New Two-Moment Bulk Stratiform Cloud Microphysics Scheme in the Community Atmosphere Model, Version 3 (CAM3). Part I: Description and Numerical Tests, J. Climate, 21, 3642–3659,, 2008. a

Mülmenstädt, J., Sourdeval, O., Delanoë, J., and Quaas, J.: Frequency of occurrence of rain from liquid-, mixed-, and ice-phase clouds derived from A-Train satellite retrievals, Geophys. Res. Lett., 42, 6502–6509,, 2015. a

Murakami, M.: Numerical modeling of dynamical and microphysical evolution of an isolated convective cloud, J. Meteorol. Soc. Jpn., 68, 107–128,, 1990. a

Palmer, T. N.: A nonlinear dynamical perspective on model error: A proposal for non-local stochastic-dynamic parametrization in weather and climate prediction models, Q. J. Roy. Meteor. Soc., 127, 279–304,, 2001. a

Pincus, R. and Klein, S. A.: Unresolved spatial variability and microphysical process rates in large-scale models, J. Geophys. Res.-Atmos., 105, 27059–27065,, 2000. a, b

Quaas, J.: Evaluating the “critical relative humidity” as a measure of subgrid-scale variability of humidity in general circulation model cloud cover parameterizations using satellite data, J. Geophys. Res.-Atmos., 117, D09208,, 2012. a

Sassen, K. and Wang, Z.: Classifying clouds around the globe with the CloudSat radar: 1-year of results, Geophys. Res. Lett., 35, L04805,, 2008. a

Sourdeval, O., Gryspeerdt, E., Krämer, M., Goren, T., Delanoë, J., Afchine, A., Hemmer, F., and Quaas, J.: Ice crystal number concentration from satellite lidar-radar observations (DARDAR-Nice), AERIS [data set],, 2018. a, b

Stephens, G. L., Vane, D. G., Boain, R. J., Mace, G. G., Sassen, K., Wang, Z., Illingworth, A. J., O'connor, E. J., Rossow, W. B., Durden, S. L., Miller, S. D., Austin, R. T., Benedetti, A., Mitrescu, C., and the CloudSat Science Team: THE CLOUDSAT MISSION AND THE A-TRAIN: A New Dimension of Space-Based Observations of Clouds and Precipitation, B. Am. Meteorol. Soc., 83, 1771–1790,, 2002. a, b

Stevens, B., Giorgetta, M., Esch, M., Mauritsen, T., Crueger, T., Rast, S., Salzmann, M., Schmidt, H., Bader, J., Block, K., Brokopf, R., Fast, I., Kinne, S., Kornblueh, L., Lohmann, U., Pincus, R., Reichler, R., and Roeckner, E.: Atmospheric component of the MPI-M Earth System Model: ECHAM6, J. Adv. Model Earth Sy., 5, 146–172,, 2013.  a

Sundqvist, H., Berge, E., and Kristjánsson, J. E.: Condensation and cloud parameterization studies with a mesoscale numerical weather prediction model, Mon. Weather. Rev., 117, 1641–1657,<1641:CACPSW>2.0.CO;2, 1989. a, b, c, d

Thayer-Calder, K., Gettelman, A., Craig, C., Goldhaber, S., Bogenschutz, P. A., Chen, C.-C., Morrison, H., Höft, J., Raut, E., Griffin, B. M., Weber, J. K., Larson, V. E., Wyant, M. C., Wang, M., Guo, Z., and Ghan, S. J.: A unified parameterization of clouds and turbulence using CLUBB and subcolumns in the Community Atmosphere Model, Geosci. Model Dev., 8, 3801–3821,, 2015. a

Trömel, S., Simmer, C., Blahak, U., Blanke, A., Doktorowski, S., Ewald, F., Frech, M., Gergely, M., Hagen, M., Janjic, T., Kalesse-Los, H., Kneifel, S., Knote, C., Mendrok, J., Moser, M., Köcher, G., Mühlbauer, K., Myagkov, A., Pejcic, V., Seifert, P., Shrestha, P., Teisseire, A., von Terzi, L., Tetoni, E., Vogl, T., Voigt, C., Zeng, Y., Zinner, T., and Quaas, J.: Overview: Fusion of radar polarimetry and numerical atmospheric modelling towards an improved understanding of cloud and precipitation processes, Atmos. Chem. Phys., 21, 17291–17314,, 2021. a, b

Weber, T. and Quaas, J.: Incorporating the subgrid-scale variability of clouds in the autoconversion parameterization using a PDF-scheme, J. Adv. Model Earth Sy., 4, M11003,, 2012. a

Wegener, A.: Thermodynamik der atmosphäre, JA Barth, (last access: 16 April 2024), 1911. a

Winker, D. M., Pelon, J., Coakley, J. A., Ackerman, S. A., Charlson, R. J., Colarco, P. R., Flamant, P., Fu, Q., Hoff, R. M., Kittaka, C., Kubar, T. L., Treut, H. L., Mccormick, M. P., Mégie, G., Poole, L., Powell, K., Trepte, C., Vaughan, M. A., and Wielicki, B. A.: The CALIPSO Mission: A Global 3D View of Aerosols and Clouds, B. Am. Meteorol. Soc., 91, 1211–1230,, 2010. a, b

Short summary
Especially over the midlatitudes, precipitation is mainly formed via the ice phase. In this study we focus on the initial snow formation process in the ICON-AES, the aggregation process. We use a stochastical approach for the aggregation parameterization and investigate the influence in the ICON-AES. Therefore, a distribution function of cloud ice is created, which is evaluated with satellite data. The new approach leads to cloud ice loss and an improvement in the process rate bias.