Elucidating ice formation pathways in the aerosol–climate model ECHAM6-HAM2

Cloud microphysics schemes in global climate models have long suffered from a lack of reliable satellite observations of cloud ice. At the same time there is a broad consensus that the correct simulation of cloud phase is imperative for a reliable assessment of Earth's climate sensitivity. At the core of this problem is understanding the causes for the inter-model spread of the predicted cloud phase partitioning. This work introduces a new method to build a sound cause-and-effect relation between the microphysical parameterizations employed in our model and the resulting cloud field by analysing ice formation pathways. We find that freezing processes in supercooled liquid clouds only dominate ice formation in roughly 6 % of the simulated clouds, a small fraction compared to roughly 63 % of the clouds governed by freezing in the cirrus temperature regime below −35 ^∘C. This pathway analysis further reveals that even in the mixed-phase temperature regime between −35 and 0 ^∘C, the dominant source of ice is the sedimentation of ice crystals that originated in the cirrus regime. The simulated fraction of ice cloud to total cloud amount in our model is lower than that reported by the CALIPSO-GOCCP satellite product. This is most likely caused by structural differences of the cloud and aerosol fields in our model rather than the microphysical parametrizations employed.

1 Introduction

Clouds are an important modulator for Earth's climate. They exert a net radiative effect of approximately −20 W m⁻² and thus significantly cool the planet (Boucher et al., 2013). Compared to this, the forcing induced by well-mixed greenhouse gases since pre-industrial times is almost 1 order of magnitude smaller, approximately 3 W m⁻² (Myhre et al., 2013), and has the opposite sign associated with a warming. Small changes in the cloud radiative effect can therefore easily offset or strengthen a greenhouse-gas-induced warming. In fact, there is a broad consensus that clouds contribute the largest uncertainty for climate projections in state-of-the-art atmospheric general circulation and Earth system models (Flato et al., 2013).

The models that contributed to the Coupled Model Intercomparison Project Phase 5 (CMIP5) generally agree that cloud adjustments to a warming climate will likely further reinforce the initial warming, a positive feedback loop. Clouds are expected to become fewer, higher and optically thicker in the global mean (Zelinka et al., 2013).

The warming response of cloud ice has been hypothesized to counteract and therefore reduce the global mean, net positive cloud feedback through a transition from optically thin ice clouds in the present-day climate to optically thick liquid clouds in a warmer climate (Mitchell et al., 1989; Storelvmo et al., 2015; Ceppi et al., 2016; Frey and Kay, 2018). The magnitude of this so-called cloud phase feedback strongly depends on the simulation of the present-day cloud phase partitioning in models. It is well established that models tend to underestimate the ratio of supercooled liquid water to ice in present-day conditions (Cesana et al., 2015). This leads to an overestimation of the magnitude of the cloud phase feedback (Li and Le Treut, 1992; Terai et al., 2016; Tan et al., 2016) in a greenhouse-gas-induced, warmer climate in some models. Recent modelling efforts challenge the universality of this finding (Lohmann and Neubauer, 2018; Bodas-Salcedo, 2018), calling for a more comprehensive description of ice formation pathways in general circulation models (GCMs).

Especially sensitive to global warming is ice in mixed-phase clouds, predominantly originating from heterogeneous immersion or contact freezing of mineral dust (Atkinson et al., 2013; DeMott et al., 2015; Kanji et al., 2017). Those clouds are usually found in the temperature range between −38 and 0 ^∘C, which we therefore refer to as the mixed-phase temperature regime. Beyond the theoretical and experimental uncertainties of the freezing mechanisms (Welti et al., 2014; Ickes et al., 2015; Marcolli, 2017), a realistic representation of cloud glaciation is further complicated by uncertainties in the parameterization of subsequent ice growth and sedimentation. Due to the lower saturation water vapour pressure over an ice crystal compared to a liquid droplet surface, cloud ice can grow below water saturation, eventually leading to evaporation of the cloud droplets. This process is called the Wegener–Bergeron–Findeisen (WBF) process and can glaciate a cloud within minutes to hours (Korolev and Isaac, 2003), depending on the number of ice crystals, temperature and vertical velocity of the air parcel.

Due to the complex processes governing mixed-phase clouds, models inhibit a very large spread in the simulation of the phase partitioning (Fan et al., 2011; McCoy et al., 2016). A model intercomparison study of Komurcu et al. (2014) found for the six models (ECHAM6, CAM-IMPACT, CAM-Oslo, CAM5.1 MAM3, CAM5.1 MAM7, SPRINTARS) in their study that even for the same ice nucleation parameterizations, the spread in the simulated phase partitioning among the models was not reduced.

A significant portion of cloud ice is found at temperatures below −38 ^∘C. There, ice crystals can freeze homogeneously from pre-activated cloud droplets (Lohmann et al., 2016) and deliquesced aerosols (Koop et al., 2000) or nucleate directly on an ice nucleating particle (INP). The latter two processes do not require water saturation, leading to fundamentally different types of cirrus clouds (Krämer et al., 2016; Wernli et al., 2016; Gasparini et al., 2018) with significant differences in the respective microphysical and thus optical properties, highlighting the need for a correct representation of ice formation pathways of all clouds containing ice, not just those in the mixed-phase regime.

These complex, sub-grid-scale processes can only be studies through idealized laboratory measurements and inferred from comprehensive field studies which are usually conducted over small areas. Being weakly constrained, ice process parameterizations differ greatly among models and because of that different ice formation pathways can be favoured, leading to the wide spread in phase partitioning among models (Cesana et al., 2015; McCoy et al., 2016). To find the parameterizations that control the cloud phase partitioning, we need to quantify the impact each process has on the simulated cloud fields.

The amount of ice that forms following a nucleation event in either a mixed-phase or cirrus cloud strongly depends on the available humidity and, in the case of mixed-phase clouds, liquid water. A strong nucleation event does not necessarily induce a long-lived or massive ice cloud and vice versa. Furthermore, sedimentation of cloud ice disconnects ice initiation and growth regions, which further complicates diagnosis of the causal link between the parameterized microphysical processes and the observed cloud state. Traditional, spatio-temporally averaged model output usually reports the mean state of the cloud fields and sometimes microphysical tendencies, but does not allow the dominant microphysical formation mechanisms to be traced back due to the reasons above. Unlike the real atmosphere, which does not directly reveal its governing processes, the models do so by definition. Here we harvest this additional information that is present in the models by introducing additional tracers to quantify the microphysical ice formation pathways along with the model integration to assess the importance of different ice initiation and growth mechanisms.

Furthermore, we use the pathway analysis to classify clouds and accumulate statistics on their properties. The benefit of cloud classification is well established and has a long history. A common method is the identification of dynamical cloud regimes to compute conditional statistics (Jakob and Tselioudis, 2003; Williams and Tselioudis, 2007; Williams and Webb, 2009; Tsushima et al., 2016). These methods use high-frequency output and short-term temporal aggregation of cloud top temperature vs. cloud optical thickness histograms together with clustering algorithms, classifying the mean dynamic states over the aggregation period. Our method focuses on the microphysics of cloud formation linking parameterizations and the resulting cloud fields. This link is essential to identify the parametrizations responsible for potential model-to-model or model-to-observation mismatches.

This study is based on an adapted version of the predicted bulk particle properties (P3) scheme of Morrison and Milbrandt (2015) for the use in a GCM. Single-category ice phase microphysics schemes have recently also been implemented in the Community Atmosphere Model Version 5 (CAM5) GCM (Eidhammer et al., 2017; Xi et al., 2017). Dietlicher et al. (2018) (hereafter D18) provides a detailed description of the technical aspects of the implementation of P3 in ECHAM6-HAM2, the GCM we use in this study. This scheme has the decisive advantage of representing cloud ice in a consistent manner, predicting the particle size distribution as well as the mass-to-size relationship. It employs one single, prognostic ice phase category, rendering heuristic conversion rates between in-cloud and precipitation type ice hydrometeors unnecessary. Process rates are computed offline and read back from look-up tables. This allows each process rate to be integrated over the entire particle spectrum. For example, our scheme consistently represents the size dependence of the depositional growth rate using a size-dependent ice crystal capacitance. Simplified, size-independent ice crystal growth rates are hypothesized to significantly contribute to the large spread among models in the simulated phase partitioning (Fan et al., 2011).

This study can be divided into two parts. First, we show how the new model performs on the global scale and then introduce the method to elucidate the ice formation pathways in the second part of this study. In Sect. 2 we highlight the main differences in our model compared to the reference model version ECHAM6.3-HAM2.3 and discuss the new cloud fraction parameterization that has been added since the previous description of the version in D18. We evaluate the simulation of clouds, with a special focus on cloud ice, against a series of satellite observations in Sect. 3. In the second part, we introduce additional prognostic equations that are solved to quantify ice formation pathways in Sect. 4. In Sect. 5 we compile a climatology of cloud types based on these ice formation pathways and demonstrate how they can help us understand cloud phase partitioning.

2 Model description

The improvements to the previous model version (ECHAM6.3-HAM2.3; Neubauer et al., 2014) (hereafter called the reference model) are twofold. In the new model, cloud ice is a fully prognostic quantity, i.e. the advection equation for vertical transport (sedimentation) is solved online. Describing all ice particles with only a single category based on Morrison and Milbrandt (2015) no longer requires the heuristic partitioning between in-cloud ice and precipitating snow. The technical details of the model regarding prognostic treatment of sedimentation and the associated numerical stability restrictions on the time step as well as an inter-comparison of the different ice representation schemes can be found in D18. Here we will focus on the global evaluation of the new model both in terms of the feasibility of the results as well as computational performance.

The (P3) scheme presented in Morrison and Milbrandt (2015) predicts four properties of the ice particle distribution: the total ice mass mixing ratio q_i, the total ice number mixing ratio N_i, the rimed ice mass mixing ratio q_rim and the rimed ice density b_rim. Originally developed for the Weather Research and Forecasting model (WRF), the scheme was not intended for the coarse resolution of GCMs and climate projections. As we have shown in D18 by idealized single-column model (SCM) simulations, GCMs are likely unable to properly represent small-scale weather features like squall lines to produce significant rimed ice formation. We test this hypothesis in the global set-up of ECHAM-HAM in this paper.

The single conceptual deviation from the microphysics scheme described in D18 concerns the sub-grid cloud cover parameterization. While this work focuses on the importance of the representation of cloud ice, the sub-grid cloudiness is a fundamental property of clouds in GCMs, governing both cloud–radiation and cloud–precipitation interactions.

The approach taken in D18 tried to minimize the difference between the new and the reference model by extending the approach of Sundqvist et al. (1989) (hereafter S89) to ice clouds with a smooth transition from liquid water to ice clouds. As we will discuss in Sect. 3, this scheme leads to a strong positive bias of cloud cover for cirrus clouds. Therefore, we used a slightly different cloud cover scheme here.

The cloud cover scheme has also implications for the growth by condensation and deposition. Diagnosing sub-grid cloud cover from the large-scale relative humidity creates a tight link between the cloud fraction and growth by condensation or deposition. The fundamental concept of the S89 scheme is that convergence of humidity within a grid box contributes to an increase in total cloud water by condensation and deposition as well as an increase in the cloud fraction through moistening of the grid box. Here we extend this idea to cold clouds and implement a cloud cover scheme which treats the microphysical structure of ice clouds consistently, requiring mixed-phase clouds to form via the liquid phase and allowing supersaturation in the cirrus regime.

2.1 Sub-grid cloud fraction

The original sub-grid cloud cover scheme of S89 does not consider cloud ice and assumes that a grid box will be fully covered by the cloud if the water vapour mixing ratio q_v surpasses the liquid water saturation mixing ratio q_sw, i.e. the relative humidity with respect to liquid water s_w satisfies $s_{\mathrm{w}}=q_{\mathrm{v}}/q_{\mathrm{sw}}=\mathrm{1}$. Cloud cover b increases with increasing s_w as follows:

$$\begin{array}{}\text{(1)}& b=\mathrm{1}-{\left({\displaystyle \frac{s_{max}-s_{\mathrm{w}}}{s_{max}-s_{min}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

The parameter s_min is the threshold relative humidity for sub-grid cloudiness below grid box mean water saturation and $s_{max}=\mathrm{1}$. The assumptions underlying Eq. (1) are reasonable for synoptic-scale warm clouds. However, this scheme cannot be easily extended to cold clouds where a grid box may or may not contain cloud ice in the mixed-phase regime, rendering the global use of s_w for all temperatures inadequate. At the same time, relative humidity with respect to ice (s_i) can reach values higher than 140 % before ice nucleates from deliquesced aerosols (Koop et al., 2000), which is inconsistent with the use of $s_{max}=\mathrm{1}$ for ice clouds at cold temperatures. A few approaches exist to overcome either of these aspects. In the following we present the reasoning that led to the development of the new method used here. For an overview, see Fig. 1.

Figure 1 Regions where clouds can form (grey shaded areas) as a function of temperature and supersaturation with respect to ice. The coloured lines show liquid water and ice saturation ratios s_w∕i (solid lines) and the critical saturation ratios required for first cloud formation (dashed; s_w∕iK). We choose a constant K=0.8 for this illustration. Panel (a) shows the sub-grid cloud scheme of the reference model ECHAM6.3-HAM2.3 which switches between liquid and ice saturation in the mixed-phase regime, (b) shows the mixed-saturation approach of Dietlicher et al. (2018) (D18) and (c) shows the new model.

As illustrated in Fig. 1a, Lohmann and Roeckner (1996) (the method used in the reference model) replaced s_w in Eq. (1) by s_i if a threshold amount of ice is exceeded within the grid box. This approach has the drawback that clouds artificially expand upon glaciation because of the lower saturation water vapour pressure over ice than over liquid water. A similar approach is taken in Morrison and Gettelman (2008) (and D18; illustrated in Fig. 1b) where instead of the discontinuous transition from liquid water to ice saturation, the saturation liquid and ice water vapour mixing ratios are interpolated between q_sw and q_si as a function of temperature in the mixed-phase regime. The drawback here is that this interpolated value is not directly relatable to a physical quantity as it does not arise from a valid solution of the Clausius–Clapeyron equation. Both of these schemes are not designed to allow for sub-grid clouds in the cirrus regime where s_i needs to be well above 1 for ice nucleation to occur; i.e. all cirrus clouds formed by homogeneous freezing will occupy the entire grid box. Supersaturation with respect to cloud ice has been accounted for in the cloud cover scheme by Gettelman et al. (2010), where individual cloud fractions for liquid water and ice clouds are computed with a functional dependency of the cloud fraction b(s) similar to Eq. (1) but relaxing the relative humidity for complete cloud coverage to $s_{max}=\mathrm{1.1}$ for ice clouds and treating the value of s_max as a model parameter.

The cloud fraction parameterization used in our microphysics scheme is very similar to this last approach, but we do not separate cloud cover for liquid and ice clouds. After all, in the mixed-phase regime we expect ice clouds to originate from liquid clouds (Ansmann et al., 2009; Hoose et al., 2010; Hande and Hoose, 2017). Furthermore, the magnitude of supersaturation for complete cloud cover in Gettelman et al. (2010) is arbitrary. The scheme presented here uses a limit which is consistent with the parameterization of ice nucleation based on the theory of Koop et al. (2000). This scheme is illustrated in Fig. 1c.

One can associate s_min and s_max in Eq. (1) with the minimal relative humidity below which clouds start to evaporate and the maximal relative humidity which can be reached before the cloud covers the entire grid box. These two constraints provide some guidelines on how to extend the S89 scheme to cold clouds.

For clouds that form via liquid water, we assume that s_max is reached at liquid water saturation, e.g. when liquid water is assumed to condense immediately, which is a reasonable assumption for stratiform clouds. At temperatures warmer than −35 ^∘C, any ice must form via liquid water, consistent with the above-mentioned lidar observations and modelling studies. For colder temperatures, it has been shown (Koop et al., 2000) that ice can nucleate from deliquesced aerosols already below liquid water saturation.

The minimal relative humidity s_min must be chosen such that ice does not sublimate above ice saturation. Cloud ice can sediment into and form within the mixed-phase regime. For this reason, and in order to retain sub-grid cloudiness for warm clouds, $s_{min}\le \mathrm{1}$ must hold for all temperatures.

2.2 Growth by condensation and deposition

The water vapour available for condensation is linked to the cloud cover b by the following:

$$\begin{array}{}\text{(2)}& Q=-b(\mathrm{\Delta }{q}_{\mathrm{v}}-\mathrm{\Delta }{q}_{\mathrm{sw}}),\end{array}$$

where Δq_v is the moisture convergence in the grid box by the resolved transport and Δq_sw is the change in liquid water saturation water vapour mixing ratio according to the Clausius–Clapeyron equation. Growth by deposition A is computed as a function of the relative humidity with respect to ice s_i (see e.g. D18):

$$\begin{array}{}\text{(3)}& A=\mathrm{\Delta }t{N}_{\mathrm{i}}{\mathit{\alpha }}_{\mathrm{m}}{f}_{\mathrm{v}}{\displaystyle \frac{\mathrm{4}\mathit{\pi }C(s_{\mathrm{i}}-\mathrm{1})}{F_k^i+F_d^i}},\end{array}$$

where α_m=0.5 is the probability of a water vapour molecule to successfully be incorporated into an ice crystal, Δt is the model time step and C is the mean ice particle capacitance integrated over the particle size distribution. The parameters $F_k^i$ and $F_d^i$ are thermodynamic parameters depending only on temperature.

The maximal amount of water that can deposit in a single time step is given by the following:

$$\begin{array}{}\text{(4)}& Q_{\mathrm{i}}=Q+q_{\mathrm{l}}+(q_{\mathrm{v}}-q_{\mathrm{si}}),\end{array}$$

where q_l is the liquid water mixing ratio and q_si is the saturation water vapour mixing ratio with respect to ice. Equation (4) implicitly represents the WBF process in mixed-phase clouds through the inclusion of q_l. The last term in Eq. (4) allows any water vapour which is supersaturated with respect to cloud ice to deposit onto ice crystals and thus requires to add q_i to the water vapour mixing ratio for the computation of the relative humidity in Eq. (1) to prevent clouds from shrinking due to deposition. To account for ice clouds, we replace the relative humidity with respect to liquid water s_w in Eq. (1) by the more general relative humidity term s defined below.

Combining the considerations for growth by condensation/deposition and the computation of cloud fraction, we obtain the key parameters for the computation of sub-grid cloudiness in Table 1. Two new parameters are introduced: the critical relative humidity with respect to liquid water for warm clouds K, which is equivalent to the parameterization for warm clouds used in the reference model and based on Xu and Krueger (1991), and the relative humidity required for nucleation of ice from deliquesced aerosols s_koop according to the activation theory of Koop et al. (2000) in the cirrus regime.

Table 1 The parameters involved in the cloud cover scheme, Eqs. (1) and (4). The rows show different parameters as a function of the temperature regime T. We discuss two different choices for s and Q_i in the text, $s_{\mathrm{1}}/Q_{\mathrm{i},\mathrm{1}}$ and s₂/Q_i,2. We use a linear weighting function $w\left(T\right)=(T-\mathrm{273.15})/(\mathrm{238.15}-\mathrm{273.15})$ for T in K.

Two versions of the s and Q_i terms are given in Table 1 and represent two different methods to treat depositional growth in the mixed-phase regime. They bracket the main problem associated with cloud ice when diagnosing cloud cover from relative humidity. Cloud ice should be able to grow as long as water vapour is supersaturated with respect to ice. Since we are using the relative humidity with respect to liquid water (which is lower than the relative humidity with respect to ice) as maximal relative humidity s_max in the mixed-phase regime, ice growth below water saturation leads to a decrease in relative humidity and therefore lowers the cloud fraction b. Therefore we implement two versions of the new scheme. The first, default, method (denoted by a subscripted 1) links cloud cover to the ice mass mixing ratio. It allows cloud ice to grow below liquid water saturation and adds cloud ice to water vapour for the computation of s in Eq. (1) to avoid cloud shrinking by vapour deposition. However, this coupling makes the sedimentation sink of cloud ice also a sink for cloud fraction. We test the sensitivity of this coupling with a second method (denoted by a subscripted 2) where ice growth is inhibited below water saturation in the mixed-phase regime.

3 Model evaluation

We evaluate the new model against a series of satellite observations as well as the reference model and test sensitivities to parametrization choices. Namely, we extend the investigation of ice properties presented in D18 to the global set-up of ECHAM-HAM and assess the behaviour of the new cloud cover parametrization. Furthermore, we use an updated version of the cirrus parametrization, including heterogeneous ice nucleation of mineral dust at temperatures below −35 ^∘C at low ice supersaturation based on Kärcher et al. (2006) (implementation courtesy of Steffen Münch, personal communication, 2018).

Results from five different model configurations for 10 years from 2003 to 2012 are shown, one with the reference model (ECHAM6.3-HAM2.3; REF) and four with the new scheme. We assess the influence of rime properties in the full P3 scheme (4M) by comparison with a model configuration where the rime properties are set to zero ($q_{\mathrm{rim}}=b_{\mathrm{rim}}=\mathrm{0}$) (2M). For the reasons presented below, we use the latter as the default configuration of the new model. With 2M as the starting configuration, we then investigate the effect of the artificial sedimentation–cloud-cover feedback introduced by the new cloud cover scheme by limiting ice growth to liquid water saturation (LIM_ICE) and use an updated version of the cirrus cloud parametrization which includes heterogeneous ice nucleation on mineral dust below −35 ^∘C (HET_CIR). The simulations and the employed tuning parameters are summarized in Table 2. The tuning parameters have been adjusted such that the model is in radiative balance at the top of the atmosphere (TOA) for the simulations 2M, HET_CIR and REF.

Table 2 Description of the model configurations shown in this paper. Tuning parameters are as follows: γ_r is the scaling factor for warm rain formation, f_fall is a scaling factor for ice sedimentation speed, γ_s is a scaling parameter for snow formation, e_ii is the collision efficiency of ice crystals, and γ_cpr is the conversion rate from cloud water and ice to precipitation in the convection scheme. Dashes (–) denote that those parameters are no longer needed in the new scheme. The last column shows the CPU time for 1 year of simulation relative to the reference model.

3.1 Computational performance

The new model employs an adaptive time-stepping scheme to achieve numerical stability for the solution of the vertical advection equation for prognostic sedimentation of cloud ice. Sub-stepping the microphysics scheme comes at a cost, quantified by the CPU time column in Table 2. The numbers shown are computed from 1 year of simulation. Numbers vary by a few percent between runs, illustrated by the differences between the 2M, LIM_ICE and HET_CIR simulations, which are comparable in terms of model complexity. Compared to the reference model, the new model runs approximately 25 % slower in the 2M/LIM_ICE/HET_CIR set-ups and 40 % slower in the 4M set-up. This difference is mainly due to the high fall speeds of rimed particles in the 4M simulation, which require more sub-steps. Compared to this, the cost of advecting two additional tracers is small.

3.2 Comparison to the reference model

Two main differences between the new model and the reference model are immediately evident from Table 2. The single category no longer requires heuristic parametrizations for falling ice crystals and snow formation. At the same time, the new cloud cover parametrization allows the scaling factors to be reduced, effectively bringing the parametrizations closer to their conceptual origin.

To illustrate the most prominent differences in the simulated cloud field, we compare cloud water contents (ice and liquid) as well as cloud cover in Fig. 2 of the new model in the 2M configuration to the reference model. It is evident that the new model produces a lot more ice than the reference model. A large part of the difference can be explained by the fact that the single-category scheme includes cloud ice and snow while snow is a diagnostic quantity which is not included in the ice water content in the two-category scheme employed in the reference model. Of course there are other contributing factors like parameter tuning, a switch from the parametrization of homogeneous freezing of cloud droplets from Kärcher and Lohmann (2002) to Jeffery and Austin (1997) for technical reasons and the cloud cover parametrization, which are hard to disentangle.

Figure 2 The 10-year zonal mean cloud ice water content (a, b, c), cloud liquid water content (d, e, f) and cloud cover (g, h, i) for the new model (2M configuration; a, d, g), the reference model (REF; b, e, h) and differences between the two schemes (c, f, i).

In the simulated cloud liquid water there is no structural difference; the 2M simulation has substantially higher values everywhere. As the warm phase is not directly affected by the changes in cloud ice and the cloud cover parametrization, the differences are due to different autoconversion tuning parameters (γ_r in Table 2) and interactions with the new cloud cover and ice phase parametrizations.

The cloud cover differs significantly due to the new cloud cover parametrization in the new model. Differences are most pronounced in the cirrus regime where a cloud now only covers the entire grid box at the high supersaturation needed for homogeneous nucleation of solution droplets (see Sect. 2). Since this leads to a substantial difference in the cirrus cloud cover, we compare the reference model and the 2M configuration of the new model to satellite data in the following. For this purpose, we make use of the Cloud Feedback Model Intercomparison Project (CFMIP) Observation Simulator Package (COSP) (Bodas-Salcedo et al., 2011), which has already been implemented in ECHAM6-HAM2, and compare its output to the appropriate satellite product; the Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO), i.e. the GCM-oriented CALIPSO cloud product (GOCCP) (Cesana and Chepfer, 2013). Using a simulator allows for a consistent comparison of model output and satellite-derived cloud fields by minimizing artificial biases introduced by different cloud diagnostics and accounts for instrument sensitivity, which is particularly important for CALIPSO to account for the attenuation of the lidar signal.

The comparison to the cloud cover climatology from 2008 to 2014 of the CALIPSO-GOCCP product in Fig. 3 reveals that the new scheme fits better to the observed cloud structure (higher Pearson correlation coefficient ρ=0.75 and lower root mean squared error RMSE=4.0 % in the new model as compared to ρ=0.72 and RMSE=6.4 % in the reference model), despite an overall underestimation of cloudiness. It is worth noting that this underestimation is a consequence of the new cloud cover parametrization and not a result of model tuning. The distinct transition from the negative to positive bias between 5 and 10 km in the cloud fraction of the reference simulation is a consequence of the discontinuous description of cloud cover from the mixed-phase temperature regime (0 to −35 ^∘C) to the cirrus regime (colder than −35 ^∘C) that has been removed by the new scheme. Despite this improvement in terms of the vertical cloud structure, we find that the total cloud cover is worse in the new model. Figure 4c shows that the total cloud cover differs less between the two models than their difference in high cloud cover. This is probably due to the vertical overlap of high-level clouds with mid- and low-level clouds that already cover the entire grid box. In areas where there are no mid- and low-level clouds, the large amount of high-level clouds in the reference model can partly offset the overall underestimation of the total cloud cover. As a result, the reference model performs better in terms of RMSE for this variable. This is in agreement with Cesana and Waliser (2016), who find that climate models produce too few cloudy columns but that the clouds cover more vertical levels when present.

Figure 3 Comparison of the simulated cloud fraction for simulations with the 2M configuration of the new model and the reference model (REF) for 10 years (2003–2012) and the CALIPSO-GOCCP satellite product (day and night, 7 years, 2008–2014). Model output is computed using the COSP simulator. The Pearson correlation coefficient (ρ) and root mean square error (RMSE) of the full, three-dimensional cloud cover climatology are computed for both models and shown in the difference plots. The dashed line indicates the area where there is a large (more than 5 %) difference between the 2M and REF simulations as shown in Fig. 2i.

Figure 4 The 10-year zonal mean values for longwave (LW CRE) and shortwave cloud radiative effects (SW CRE) (a, b) and total cloud cover and precipitation (c, d). The coloured lines show simulations of different configurations of the new model in Table 2. Black lines show satellite products. For CRE we show data from CERES EBAF (2000 to 2017). Total cloud cover is shown for CALIPSO-GOCCP (2006 to 2012) and the model data are produced using the COSP-simulator. Precipitation data are from the GPCP (2003 to 2012). The Pearson correlation coefficient and root mean square error are reported by text in matching colours for the different simulations and computed using the full two-dimensional fields.

3.3 Model tuning strategy

Model tuning has been conducted for the reference model (REF) and the two configurations of the new model, 2M and HET_CIR, while the 4M and LIM_ICE configurations use the same tuning as 2M. A summary of global, annual mean quantities is shown in Table 3. The main target of the tuning process has been the global, annual mean shortwave (SW) and longwave (LW) fluxes at TOA as well as the sum of the two. With the tuning parameters summarized in Table 2 we have a direct handle on the global, annual mean cloud radiative effect (CRE) at TOA, defined as the difference of all-sky and clear-sky radiative fluxes. To reach our SW and LW radiation as well as the radiation imbalance tuning targets at TOA, we adjust the CRE. A consequence of this tuning strategy is that all model simulations (including the reference model with a substantially different microphysics scheme) have a net CRE roughly between −26 and −27 W m⁻² (depending on the achieved imbalance), which is more negative than any of the observational estimates and points at a structural problem in the model which is not specific to the microphysics scheme.

Table 3 The 10-year global average values for a selection of key microphysical parameters for the model configurations presented in Table 2. Observational data are used from the compilation of Lohmann and Neubauer (2018) with the exception of ice water contents where we use the dataset of Li et al. (2012) and report their average standard deviation.

We compare the simulated zonal, annual mean precipitation to the Global Precipitation Climatology Product (GPCP) (Adler et al., 2018) for the years 2003 to 2012. Total precipitation is largely governed by the rate of evaporation at the surface. As the sea surface temperature is prescribed and most evaporation occurs over the oceans, cloud microphysics and convection parametrizations must produce similar amounts of total surface precipitation.

The lower cloud cover of the 2M and 4M simulations as compared to the LIM_ICE simulation is a direct consequence of the coupling of cloud ice and cloud cover. There are two main contributors: (1) the LIM_ICE configuration does not link the sedimentation sink of cloud ice to cloud cover while the other configurations (2M, 4M and HET_CIR) do and (2) the LIM_ICE configuration limits depositional growth of ice crystals and therefore the removal of water vapour by ice sedimentation is weaker, leaving more humidity in the atmosphere which in turn leads to more clouds.

For the TOA energy fluxes we use the climatology from 2000 to 2017 from the Clouds and the Earth's Radiant Energy System (CERES) Edition 4.0 Energy Balanced And Filled (EBAF) (Loeb et al., 2018) product for the model evaluation.

Without the overestimation of cirrus cloud cover, we can reduce the convective rain formation rate γ_cpr to values closer to the pure ECHAM6 (without online aerosols from the HAM model and a single-moment cloud microphysics scheme; see Stevens et al., 2013) model with a value of $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹ (Mauritsen et al., 2012). This allows more water vapour to be retained in the updraught cores, thus enhancing the formation of tropical cirrus clouds in the outflow regions of convective anvils, associated with a more pronounced peak in the LW CRE in the tropics (Fig. 4). However, the correlation coefficients with the CERES EBAF dataset are very similar to the reference model while the RMSEs are larger in all configurations of the new model. This is most likely a consequence of the different tuning in the two models. In the reference model both LW and SW CRE are stronger while the TOA fluxes are at the lower end of the tuning target in terms of magnitude; see Table 3.

3.4 Cloud ice

In Fig. 5 we evaluate the representation of cloud ice in the new model against the dataset compiled from CloudSat and CALIPSO retrievals by Li et al. (2012). This dataset is based on three different cloud water content products: 2B-CWC-RO4 (Austin et al., 2009), DARDAR (raDAR/liDAR) (Hogan, 2006; Delanoë and Hogan, 2008, 2010) and 2C-ICE (Deng et al., 2013). This dataset has the advantage that considerable effort has been put into the distinction between in-cloud ice in stratiform clouds (CIWC) and snow from stratiform clouds and convective cores (TIWC). The latter can usually not be diagnosed in models which treat snow formation diagnostically. The new single-category ice scheme used in this study only predicts the total stratiform ice water content, including snow. Therefore TIWC is most meaningful for a comparison with this model in the extra-tropics, where the contribution from convective snow is rather small. The reference model predicts stratiform cloud ice explicitly but diagnoses the snow mass flux. Therefore we can make full use of the distinction by Li et al. (2012) and compare the reference model to CIWC.

Figure 5 Ice water content and ice water path compiled from observations by Li et al. (2012) and different versions of ECHAM6-HAM2 (10-year zonal averages) with the new microphysics scheme and the reference model. Panels (a)–(e) show vertical profiles for ice water content, panel (f) shows ice water path. The observations in black include total ice water content/path (TIWC/P; solid lines) and cloud ice water content/path (CIWC/P; dotted lines). Colours show total (stratiform) ice water content from the new model and cloud ice water content from the reference model. The grey shaded areas represent the average standard deviation as reported by Li et al. (2012). The profiles are averaged over different latitude bands: SH high-lat (80 to 60^∘ S), SH mid-lat (60 to 30^∘ S), tropics (30^∘ S to 30^∘ N), NH mid-lat (30 to 60^∘ N) and NH high-lat (60 to 80^∘ N).

Comparing the different simulations of the new scheme, we see that the 2M and 4M simulations are very similar both in terms of ice water content profiles (Fig. 5a–e) and ice water path (TIWP and CIWP respectively; Fig. 5f). The effect of riming on the ice crystal properties is most evident in the mid-latitudes. There, the 4M configuration produces rimed or partially rimed particles with higher fall speeds, which results in a reduced ice water content in the lower atmosphere (below roughly 700 hPa), as compared to the 2M configuration, where riming does not affect particle properties and thus fall speeds. However, the overall similarities suggest that the increased computational cost of the 4M configuration (see Table 2) and parametrization complexity (four versus two prognostic ice moments) do not significantly improve the overall representation of cloud ice in the GCM.

As explained in Sect. 2, the new cloud cover parametrization couples the cloud fraction to the ice water content. The aggregation parametrization depends on the concentration of ice crystals $S_{\mathrm{agg}}\propto N_{\mathrm{i}}^{\mathrm{2}}$, which in turn is computed from the grid box mean value as $N_{\mathrm{i}}={\stackrel{\mathrm{‾}}{N}}_{\mathrm{i}}{b}^{-\mathrm{1}}$, i.e. $S_{\mathrm{agg}}\propto b^{-\mathrm{2}}$. We denote the parametrized source and sink terms by $S_{\cdot }=({\partial }_{\mathrm{t}}{q}_{\mathrm{i}}{)}_{\cdot }$ i.e. the partial time derivative ∂_t of cloud ice q_i restricted to one process. This leads to an artificial feedback loop: aggregation reduces N_i and increases ice particle size and fall speed and thus the sedimentation sink. Because cloud ice is linked to cloud cover (see the definition of s₁ in Table 1, which replaces s_w in Eq. (1) in the new cloud cover parametrization) the sedimentation of cloud ice decreases the cloud fraction b which in turn increases the in-cloud N_i and the ice crystal aggregation rate S_agg artificially. This feedback is the main contributor to the overall slightly lower ice water contents in 2M simulation as compared to the LIM_ICE simulation, which does not have this feedback, and outweighs the effect from the additional humidity that is available for ice growth in the 2M simulation. We acknowledge the existence of this artificial increase in the in-cloud N_i by sedimentation and the resulting feedback, but do not consider it to be important enough to sacrifice the physical correctness of the method used in the 2M configuration.

The most prominent differences between our simulations and the observations are in the tropics. As can be seen from Table 3, two-thirds of the entire surface precipitation is produced by the convective scheme. Since convective precipitation is a diagnostic quantity which is not included in the total ice water content here, a direct comparison of modelled, stratiform IWC/IWP with the observations is limited. With that in mind, we find that the 2M and 4M configurations of the new model underestimate the observed TIWP by roughly a factor of 3 in the tropics (Fig. 5). This is true for both the new model and the reference model (which has to be compared to CIWP/CIWC). Higher values for TIWP are produced by the LIM_ICE and HET_CIR simulations. While the slightly higher values in the LIM_ICE simulation are due to the cloud-fraction–sedimentation feedback explained above, the IWC values in the HET_CIR simulations are almost twice as large as in the 2M or 4M configurations and are a result of the lower ice supersaturation needed for ice crystals to nucleate heterogeneously. This allows for a more frequent formation of ice clouds at cold temperatures.

In the extra-tropics, the stratiform IWC is more representative of the observed TIWC/TIWP and we find slightly better agreement between the models and observations (Fig. 5). The vertical profile of TIWC in the mid-latitudes and high latitudes is not reproduced by the new model in any configuration, a feature shared with many CMIP5 models (Li et al., 2012), while the reference model performs rather well and only slightly underestimates the ice water path in mid-latitudes.

On the basis of the evaluation above, we choose the 2M configuration as default for the new model for its computational efficiency compared to the 4M configuration and physical correctness compared to the LIM_ICE configuration. Heterogeneous nucleation of ice crystals at low ice supersaturation in the cirrus regime (parametrized in the HET_CIR configuration) can be included optionally for future work with the new model and yields promising results but is not used for the analysis in Sects. 4 and 5 since it is in part based on unpublished work.

3.5 Cloud liquid water

We compare the simulations with the new and the reference model to the Multisensor Advanced Climatology of Liquid Water Path (MAC-LWP) (O'Dell et al., 2008; Elsaesser et al., 2017); see Fig. 6. Microwave sensors cannot reliably distinguish precipitation from cloud water. Therefore we only evaluate regions with low precipitation, i.e. where the liquid water path (LWP) and total water path (TWP) have similar magnitude. We use the threshold $\mathrm{LWP}/\mathrm{TWP}>\mathrm{0.8}$ suggested by Elsaesser et al. (2017). The differences in liquid water are very small for the different configurations of the new model, so we only evaluate the default configuration 2M. We find an overall high bias in LWP in the extra-tropics for both models but a more pronounced effect in the new model. As discussed above, the TOA energy balance constraints in the new model require a lower tuning parameter for the autoconversion of cloud droplets to rain (γ_r) and thus thicker liquid clouds that reflect more SW radiation. We consider this the main reason for the higher liquid water path in the 2M as compared to the REF simulation, as the warm phase parametrizations are the same in both models.

Figure 6 LWP as observed by the Multisensor Advanced Climatology of LWP (MAC-LWP), climatology from 2003 to 2012. Panel (a) shows the observations, (b) shows the difference between the new model (2M) and the observations, and (c) shows the difference between the reference model and the observations. The models are both run for 10 years from 2003 to 2012. The grey shaded areas show regions where the liquid water path is dominated by precipitation ($\mathrm{LWP}/\mathrm{TWP}<\mathrm{0.8}$), i.e. where there is no reliable estimate for in-cloud liquid water path.

In the Southern Ocean, the overestimation of LWP translates into an overestimation of the SW CRE, evident from Fig. 4b. Note that the REF simulation with a smaller positive LWP bias overestimates SW CRE less.

3.6 Cloud phase partitioning

Different climate models produce a wide range of ice-to-total-cloud fractions (McCoy et al., 2016). Here we examine the histograms of the ratio of cloud phase (PR; relative fraction of ice clouds to total clouds) to temperature (T) used in the model intercomparison by Cesana et al. (2015) (their Fig. 10) in Fig. 7, where the wide variety of different cloud distributions in PR–T histograms among models has been highlighted. In this section we discuss the distributions produced by our models and compare the simulated average cloud phase ratio per temperature to the CALIPSO-GOCCP satellite product.

Figure 7 Frequency of cloud occurrence per temperature and phase ratio for histogram bin widths of 3 ^∘C and 0.1, respectively. The average ice fraction per temperature bin is shown by the red line for the new model (2M; a) and the reference model (REF; b). The data have been accumulated over 10 years of simulation from 2003 to 2012.

The different configurations of the new model have very similar PR–T histograms, which is why we only show the one for 2M. For the new model we only count cloudy regions, i.e. require b>0, to exclude falling ice (snow) with PR=1. Comparing the 2M to the REF simulation in Fig. 7 shows that the 2M simulation produces slightly less supercooled liquid cloud at warm temperatures above −15 ^∘C but slightly more at temperatures below −15 ^∘C. In contrast to REF, where any ice above 0 ^∘C melts instantly, the finite melting rate employed by the 2M configuration allows for cloud ice at temperatures warmer than 0 ^∘C.

Most clouds are either liquid or fully glaciated; mixed-phase states are rare in our models as well as the satellite product, assuming that the undefined-phase clouds are considered as mixed-phase clouds (Cesana et al., 2016). Due to the finite melting rate employed in the new model not all ice crystals melt entirely within one model time step, making mixed-phase states at these warm temperatures around 0 ^∘C slightly more frequent as compared to the reference model. Nevertheless, the full PR–T histogram does not contain significantly more information than the average phase ratio per temperature.

In Fig. 8 we compare the simulated cloud phase partitioning to the lidar observations from the CALIPSO satellite (CALIPSO-GOCCP product). To allow for a quantitative comparison, we used the COSP lidar simulator, which minimizes the differences between model output and the satellite product due to different sampling methods and lidar attenuation. As suggested by Cesana et al. (2016) we counted the clouds in the undefined cloud phase category in the satellite product as well as the COSP output as completely glaciated. The cloud phase ratio obtained from the lidar simulator in this manner is consistently smaller than the observations. In the following sections we will provide tools to identify which parametrized processes are most likely controlling the structure of this distribution.

Figure 8 Cloud phase ratio per temperature for histogram bin widths of 3 ^∘C. The average ice fraction per temperature bin is shown for the model data (red, solid line), the CALIPSO-GOCCP product (black, solid line) and the output from the COSP simulator applied to the model (red, dashed line). For both the model and satellite data we assumed that clouds with undefined phase as classified by the lidar simulator or retrieval algorithm are dominated by ice and thus contribute to the ice phase (Cesana et al., 2016)

4 Quantifying the ice formation pathways

The microphysical properties of a cloud are defined by the cloud formation processes, i.e. the cloud is simply the product of its formation history. Traditionally, model output reveals a snapshot of the simulated cloud field at any given time. Due to the finite storage capacity, aggregates of the model states are stored and information on the integration from process rates is lost. Most common are temporal averages (e.g. monthly or yearly mean values) or vertical aggregates (e.g. burden, total cloud cover, TOA, and surface radiative fluxes or precipitation). The aggregated cloud states are then compared to observations to infer information about the microphysical parameterizations leading to the cloud state. Inferring the dominant processes responsible for a model-to-observation mismatch is difficult. Usually it is very time consuming because many sensitivity studies are necessary where processes are turned off individually and even then the interaction and competition between processes cannot be assessed. Ultimately, spatio-temporally averaged model output defeats the purpose of simulating any physical system: finding cause and effect between differential equations and observables.

In this section we tackle the problem the other way around. Instead of inferring the formation history based on the current cloud state, we make use of additional prognostic equations to store and quantify ice formation pathways. This information can then be used to compute conditional probabilities to provide a sound cause-and-effect relation between the microphysical parameterizations and the resulting cloud state. This analysis makes direct inference from observables on cloud microphysical parameterizations tangible.

Our model parameterizes three fundamentally different source terms for cloud ice: heterogeneous contact and immersion freezing of cloud droplets in the mixed-phase regime S_het−frz, homogeneous freezing of cloud droplets S_hom−frz, and nucleation of ice from deliquesced aerosols S_nuc. The latter two processes are restricted to temperatures below −35 ^∘C while heterogeneous freezing dominates at warmer temperatures. Note that we focus on the mixed-phase regime and thus do not further separate heterogeneous from homogeneous ice nucleation in the cirrus regime, even though nucleation of cloud ice on mineral dust is implemented in the HET_CIR configuration.

To distinguish ice that formed by either of the process rates above, we introduce two additional prognostic tracers to keep track of the formation history of the cloud mass at any given time in the simulation.

4.1 Mixed-phase heterogeneous freezing origin mass fraction

To separate cloud ice that formed in the mixed-phase temperature regime we introduce the heterogeneously nucleated ice mass mixing ratio q_i,het governed by the following equation:

$$\begin{array}{}\text{(5)}& {\partial }_{\mathrm{t}}{q}_{\mathrm{i},\mathrm{het}}=S_{\mathrm{het}-\mathrm{frz}}-v_{\mathrm{m}}{\partial }_z{q}_{\mathrm{i},\mathrm{het}}+F_{\mathrm{het}}\left(S_{\mathrm{col}}+S_{\mathrm{dep}}-S_{\mathrm{sub}}\right).\end{array}$$

The growth of this tracer depends on the ice mass source and sink terms for collisions with cloud droplets S_col, vapour deposition S_dep, sublimation S_sub and the fraction of cloud ice that has already formed previously by heterogeneous freezing, defined as $F_{\mathrm{het}}=q_{\mathrm{i},\mathrm{het}}/q_{\mathrm{i}}$. We abbreviate the partial derivative with respect to the vertical dimension z (in m) by ∂_z. The tracer sediments along with the total ice mass with the mass-weighted terminal velocity v_m in metres per second (m s⁻¹). Convectively detrained ice is missing here due to a lack of an explicit, aerosol-dependent freezing parameterization in the convection scheme. Water is only detrained as ice below temperatures of −35 ^∘C where we assume that liquid water freezes homogeneously.

4.2 Liquid origin mass fraction

We further distinguish between cloud ice that forms via homogeneous freezing of cloud droplets and ice that nucleates in situ from deliquesced aerosol. This allows in situ cirrus to be disentangled from liquid origin cirrus. To quantify the ice mass fraction initiated through freezing of liquid water, we implement a tracer $q_{\mathrm{i},\mathrm{liq}-\mathrm{o}}$ governed by the following equation:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\partial }_{\mathrm{t}}{q}_{\mathrm{i},\mathrm{liq}-\mathrm{o}}& =S_{\mathrm{het}-\mathrm{frz}}+S_{\mathrm{ice}-\mathrm{cv}}+S_{\mathrm{hom}-\mathrm{frz}}-v_{\mathrm{m}}{\partial }_z{q}_{\mathrm{i},\mathrm{liq}-\mathrm{o}}\\ & +F_{\mathrm{liq}-\mathrm{o}}\left(S_{\mathrm{col}}+S_{\mathrm{dep}}-S_{\mathrm{sub}}\right).\end{array}\end{array}$$

We sum up all liquid to ice mass source terms, namely homogeneous freezing of cloud droplets S_hom−frz and convective detrainment of cloud ice S_ice−cv together with the processes defined above. The liquid origin mass fraction is given by $F_{\mathrm{liq}-\mathrm{o}}=q_{\mathrm{i},\mathrm{liq}-\mathrm{o}}/q_{\mathrm{i}}$.

4.3 Cloud types based on the cloud formation history

The additional information on the cloud formation history from Eqs. (5) and (6) can be combined with the liquid fraction $F_{\mathrm{liq}}=q_{\mathrm{l}}/(q_{\mathrm{i}}+q_{\mathrm{l}})$, cloud vertical thickness $\mathrm{\Delta }p=p_{\mathrm{base}}-p_{\mathrm{top}}$ (i.e. the pressure difference between cloud base p_base and cloud top p_top) and cloud top temperature T_top to build an inclusive set of cloud types. As is evident from Fig. 7, mixed-phase clouds are unstable due to the WBF process, leading to a strongly bimodal distribution of the frequency of phase ratio occurrence. This bimodality exists but is less pronounced for the ice mass source fractions F_het and F_liq−o shown in Fig. 9 together with the liquid fraction F_liq. Ice source mass fractions between zero and one arise from competing formation mechanisms and mixing of clouds of distinct sources.

Figure 9 Frequency of cloud occurrence per ice source fraction for the liquid origin fraction F_liq and the different ice source fractions F_het and F_liq−o. The histogram bin width is 0.1. The lines are aligned with bin centres. Only grid boxes are sampled where all the fractions are well-defined, i.e. only cloudy grid boxes for F_liq and only grid boxes containing ice for F_het and F_liq−o. Data are sampled every time step for 10 years of simulation with the new model (2M).

We make use of the separation of cloud occurrence frequency in the five-dimensional parameter space (F_liq, Δp, F_het, F_liq−o, T_top) to classify clouds into the types defined in Table 4. The labels are based on the analysis of the temperature regime each cloud primarily occurs in and are discussed in more detail below.

Table 4 Definition of cloud types. We categorize and separate clouds by exceeding thresholds of the five predictors: the heterogeneous freezing origin fraction F_het, cloud top to bottom pressure difference Δp, cloud liquid fraction F_liq, liquid origin fraction F_liq−o and the cloud top temperature T_top. If a predictor is not used for a certain class, it is symbolized by a “–” sign.

The classification is computed online and thus allows statistics to be accumulated per cloud type. Examples of such statistics are the cloud top temperature, liquid fraction and liquid origin fraction distributions per cloud type shown in Fig. 10. These distributions allow physical properties of each cloud type and its associated formation pathway to be identified. Below we use them to justify the cloud type labels in Table 4.

Figure 10 Relative frequency of cloud occurrence (normalized by the cloud type occurrence frequency). This corresponds to conditional probability histograms, $P(X_{\mathrm{•}}\in \mathrm{Bin}|\mathit{X}\in {\mathrm{Type}}_i)$ for the model state X with components X_•: cloud top temperature (a), liquid fraction (b) and liquid origin fraction (c) for each cloud type i. Bin widths are 3 ^∘C for temperature and 0.1 for the dimensionless fractions respectively; lines are aligned with bin centres. Data are sampled every time step for 10 years of simulation with the new model (2M).

Homogeneously nucleated clouds are separated into three classes. They are labelled cirrus if they are thinner than 500 hPa, otherwise we refer to them as deep. This name arises from the fact that the vast majority of such clouds have cloud tops that are colder than −35 ^∘C as can be seen in Fig. 10a. We further differentiate between in situ cirrus, clouds that nucleated from deliquesced aerosols or formed heterogeneously below −35 ^∘C in the case of the HET_CIR simulation (not shown), and liquid origin cirrus, clouds that formed from homogeneous freezing of cloud droplets. While there is a clear peak for clouds that formed without any liquid water precursors ($F_{\mathrm{liq}-\mathrm{o}}=\mathrm{0}$), the peak for complete liquid origin clouds $F_{\mathrm{liq}-\mathrm{o}}=\mathrm{1}$ is less pronounced due to competing primary ice source terms (Fig. 10c).

Heterogeneously nucleated clouds are separated into two classes. We find two distinctly different cloud types where ice formed predominantly from heterogeneous freezing of cloud droplets. As shown in Figs. 7 and 9, a truly mixed state with liquid water and ice is unstable and thus rare. Therefore we divide heterogeneously nucleated clouds into those that are dominated by liquid water and ice respectively; see Fig. 10b. We label them both as mixed-phase (liquid and ice dominated respectively) since the cloud tops of such clouds are predominantly found in the mixed-phase temperature regime between 0 and −35 ^∘C. Liquid clouds are further separated into those with cold cloud tops (T<0 ^∘C) and warm cloud tops (T>0 ^∘C).

In the following sections we will use the cloud types and associated formation pathways described here to gain insights into the simulated cloud fields. It is crucial to keep in mind that the labels provided for the cloud types here are based on formation pathways. This is in strong contrast to the traditional definitions based exclusively on the cloud state and can thus lead to results that seem counterintuitive at first. For example, cirrus clouds are commonly defined as clouds with temperatures colder than −35 ^∘C (or more generally the temperature below which homogeneous freezing becomes efficient, depending on the model). Here we only require cirrus clouds to form by the processes that are typically found at temperatures below −35 ^∘C and subsequently track their evolution without imposing strict temperature constraints. This means, however, that a cirrus cloud is still called cirrus, even if its constituent ice crystals fall far into the mixed-phase regime.

5 The ice formation pathways in ECHAM-HAM

We implemented the additional prognostic tracers needed to quantify ice formation pathways introduced in the last section in the new model. Along with the new microphysics scheme we introduced a new code structure which removed the sequential computation of microphysical tendencies in the reference model and clearly segregates the computation of microphysical tendencies and local model state updates (Dietlicher et al., 2018). This allows additional tracers to be easily implemented. Although we did not implement the same tracers in the reference model as well, we expect that the insights gained from this section can be transferred to the reference model since many aspects of the model are the same or very similar, e.g. the diffusion, convection and cloud cover parametrizations (except in the cirrus regime).

With the additional prognostic equations to identify ice formation pathways and the cloud types derived thereof, we are able find the microphysical causes for the macrophysical cloud states simulated by the new model. A clear causal relationship between the simulated cloud fields and the underlying parametrizations is critical for model development. Improving the representation of sub-grid processes is usually a two-step process: first, observational data are used to evaluate the simulated cloud fields and then model sensitivities are assessed to find the parametrizations responsible for potential model-to-observation mismatches. In the following we demonstrate how the ice formation pathways introduced in the last section can be used for the second step using the example of the phase ratio discussed in Sect. 3 and Fig. 7.

Figure 11 Global relative contribution to the 3-D cloud volume (as measured by the air mass occupied) for the cloud types defined in Table 4 based on their formation pathways. Data are sampled every time step for 10 years of simulation with the new model (2M).

Figure 12 Same as in Fig. 5 but for the absolute frequency of occurrence profile (a–e) and relative contribution to the total cloud cover (f). The relative contribution is computed as the fraction of cloud volume occupied per column and cloud type. Data are sampled every time step for 10 years of simulation with the new model (2M).

5.1 Relative cloud type frequencies

We classify the clouds online according to Table 4 to generate a global climatology of the relative contributions of each formation pathway to the 3-D cloud volume (defined as the air mass covered by the cloud) in Fig. 11. From this pie chart it is evident that the majority of the cloud population consists of homogeneously nucleated ice clouds (in situ cirrus and deep clouds) and liquid clouds. The next smallest class is the liquid origin cirrus clouds, making up 10.5 % of the cloud population or almost one-third of all cirrus clouds. Mixed-phase clouds only contribute roughly 6 %, with the fraction of such clouds that are dominated by ice being smaller than those being dominated by liquid water.

We also show zonal means of the cloud type occurrence frequencies in Fig. 12. This view sheds more light on where each formation pathway dominates. Liquid clouds with warm cloud tops dominate in the tropics while their cold top counterparts are primarily found in mid-latitudes. Due to abundant deep convection in the tropics, large amounts of liquid water are transported to high altitudes where homogeneous freezing of cloud droplets sets in, leading to a frequent occurrence of liquid origin cirrus clouds. Heterogeneous freezing only affects clouds at low altitudes where there is no competition with homogeneous nucleation. Seemingly inexistent are mixed-phase clouds where heterogeneous freezing of cloud droplets alone causes glaciation of the host liquid cloud.

5.2 Decomposing the cloud phase partitioning

By sampling the cloud types according to temperature we find the dominant parametrizations controlling the phase ratio in our model. Figure 13 decomposes the average liquid fraction per temperature bin (compare Fig. 8) into the different cloud types. Keep in mind that the cloud types are separated into mostly glaciated clouds (cirrus, mixed-phase (ice) and deep) and liquid clouds (warm and cold liquid as well as mixed-phase (liquid)) as shown in Fig. 10c. The most frequent cloud types simulated in the mixed-phase regime are deep clouds making up roughly 40 % of clouds throughout this temperature range. The second most frequent clouds are in situ cirrus, followed by liquid origin cirrus. The two heterogeneously nucleated mixed-phase cloud types make up the smallest fraction of clouds. The fact that we find homogeneously nucleated clouds with a vertical extent of less than 500 hPa (i.e. cirrus clouds) down to temperatures of 0 ^∘C is in part due to the choice of the threshold value to separate cirrus clouds from deep clouds but does not affect the main conclusion: the dominant source term for cloud ice in the mixed-phase temperature regime is homogeneous freezing of deliquesced aerosols taking place at temperatures below −35 ^∘C.

Figure 13 Cloud type fraction as a function of grid box temperature (a) and frequency of occurrence of frozen liquid fraction (b). We only sample cloudy grid boxes, i.e. those with a cloud fraction component X_b>0. For the frozen liquid fraction component $X_{F_{\mathrm{liq}-\mathrm{f}}}$ we show the relative frequency of occurrence $P(X_{F_{\mathrm{liq}-\mathrm{f}}}\in \mathrm{Bin}|\mathit{X}\in {\mathrm{Type}}_i)$. The latter is the accumulated ice mass fraction that formed through any liquid-water-to-ice conversion process (including the WBF process). Note that this is not equal to the liquid origin mass that quantifies the origin of cloud ice. Bin widths are 3 ^∘C and 0.1 respectively; lines are aligned with bin centres. Data are sampled every time step for 10 years of simulation with the new model (2M).

Analogous to the prognostic tracers introduced in Sect. 4 and formally defined in Eq. (A1), we trace the state of water from which ice forms. The frozen liquid ice mass mixing ratio tracer accumulates the mass of liquid water that has been converted to ice by freezing, riming or the WBF process and allows the frozen liquid fraction to be computed as $F_{\mathrm{liq}-\mathrm{f}}=q_{\mathrm{i},\mathrm{liq}-\mathrm{f}}/q_{\mathrm{i}}$. For the majority of homogeneously nucleated clouds, the frozen liquid fraction is very small (<10 %), from which we conclude that ice predominantly forms directly from deposition of water vapour (excluding the WBF process) and thus without a liquid water precursor; see Fig. 13b. This implies that ice nucleated in the cirrus regime not only dominates glaciation in mixed-phase clouds but even inhibits the formation of supercooled liquid water clouds. As a result, contact and immersion freezing is limited to clouds with cloud tops warmer than −35 ^∘C where there is no competition with sedimenting ice that formed in the cirrus regime.

Revisiting Fig. 8 we can come up with potential shortcomings of our model that could explain the underestimation of the ice fraction from the COSP simulator compared to observations. A straightforward explanation could be that our models simply underestimate the frequency of ice clouds. From what we learned from Fig. 13, the best candidate to increase the ice fraction at temperatures warmer than −25 ^∘C are the cold liquid clouds. These clouds have cloud tops warmer than the temperature of homogeneous freezing, which makes immersion and contact freezing the only possible pathways to glaciation. This suggests that our model lacks INPs capable of nucleating ice in this temperature range. Another explanation could be that the cloud structure is unrealistic such that there is a substantial amount of cloud ice within the deep clouds which are optically thick due to their large vertical extent. In such clouds we expect the lidar signal to be completely attenuated which makes them partly invisible to the COSP simulator.

6 Conclusions

We presented the global performance of the new cloud microphysics scheme in the ECHAM6-HAM2 GCM introduced in D18. The main difference as compared to its predecessor is the consistent description of cloud ice using a single, prognostic category. Thus, it does not rely on poorly constrained conversion parametrizations between in-cloud ice and precipitating snow categories. Different from D18, we introduced a new approach to extend the sub-grid cloud cover scheme of Sundqvist et al. (1989) to ice clouds. This scheme no longer has the positive cloud cover bias at temperatures below −35 ^∘C present in the reference model.

We assessed ice formation pathways quantitatively by introducing additional prognostic equations for the heterogeneously formed and liquid origin ice mass mixing ratios. We found that in our model the majority of cloud ice forms below −35 ^∘C by either homogeneous freezing of cloud droplets or homogeneous nucleation of deliquesced aerosols. Only about 6 % of clouds form by contact and immersion freezing in mixed-phase temperatures, making homogeneous freezing the main source for cloud ice even at temperatures warmer than −35 ^∘C. The Lagrangian perspective on the modelled cloud fields provided by the formation pathway analysis allowed in situ and liquid origin cirrus clouds to be distinguished. We found that roughly one-third of all cirrus clouds form by homogeneous freezing of cloud droplets.

Furthermore, ice formation pathways provide the causal link between the microphysical parametrizations and the resulting cloud fields. We showed how they can be used to break the simulated phase ratio down into the contributions of different microphysical parametrizations. We find that mixed-phase processes like the WBF process and heterogeneous freezing of supercooled cloud droplets play a minor role in our model because most ice originates in the cirrus regime. As a result, mixed-phase freezing parametrizations cannot control the average phase ratio per temperature. Their effect is superseded by the structure of ice clouds (with deep clouds making up most of the cloud ice in the mixed-phase regime) and the absence of INPs that could glaciate cold liquid clouds with cloud tops warmer than −35 ^∘C.

Appendix A: Frozen liquid mass fraction

We not only diagnose the phase which initiated ice formation to separate liquid origin from vapour origin ice but also keep track of the total amount of liquid water that has been converted to ice with a separate tracer $q_{\mathrm{i},\mathrm{liq}-\mathrm{f}}$. The mass mixing ratio for frozen liquid water defined as follows:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{\partial }_{\mathrm{t}}{q}_{\mathrm{i},\mathrm{liq}-\mathrm{f}}& =S_{\mathrm{het}-\mathrm{frz}}+S_{\mathrm{ice}-\mathrm{cv}}+S_{\mathrm{hom}-\mathrm{frz}}+S_{\mathrm{col}}\\ & -v_{\mathrm{m}}{\partial }_z{q}_{\mathrm{i},\mathrm{liq}-\mathrm{o}}-F_{\mathrm{liq}-\mathrm{f}}{S}_{\mathrm{sub}}.\end{array}\end{array}$$

Analogous to the tracers defined in Sect. 4, we define the frozen liquid fraction as $F_{\mathrm{liq}-\mathrm{f}}=q_{\mathrm{i},\mathrm{liq}-\mathrm{f}}/q_{\mathrm{i}}$. In contrast to the these tracers defined previously, the frozen liquid water mass mixing ratio contains information on ice growth history rather than source processes.