Evaluation of biases in mid-to-high-latitude surface snowfall and cloud phase in ERA5 and CMIP6 using satellite observations

Supercooled liquid-containing clouds (sLCCs) play a significant role in Earth's radiative budget and the hydrological cycle, especially through surface snowfall production. Evaluating state-of-the-art climate models with respect to their ability to simulate the frequency of occurrence of sLCCs and the frequency with which they produce snow is, therefore, critically important. Here, we compare these quantities as derived from satellite observations, reanalysis datasets, and Earth system models from Phase 6 of the Coupled Model Intercomparison Project (CMIP6) and find significant discrepancies between the datasets for mid- and high latitudes in both hemispheres. Specifically, we find that the ERA5 reanalysis and 10 CMIP6 models consistently overestimate the frequency of sLCCs and snowfall frequencies from sLCCs compared to CloudSat–CALIPSO satellite observations. The biases are very similar for ERA5 and the CMIP6 models, which indicates that the discrepancies in cloud phase and snowfall stem from differences in the representation of cloud microphysics rather than the representation of meteorological conditions. This, in turn, highlights the need for refinements in the models’ parameterizations of cloud microphysics in order for them to represent cloud phase and snowfall accurately. The thermodynamic phase of clouds and precipitation has a strong influence on simulated climate feedbacks and, thus, projections of future climate. Understanding the origin(s) of the biases identified here is, therefore, crucial for improving the overall reliability of climate models.

1 Introduction

Snowfall and snow cover have a significant impact on the Earth's energy budget and the hydrological cycle, especially in mid- and high-latitudes and thereby strongly impact ecosystems and human societies. Snow cover influences vegetation growth, animal populations, and ecosystem processes, while also impacting economic activities, infrastructure, and health (Callaghan et al., 2011; Bokhorst et al., 2016). Indeed, snowfall is also an important water resource globally (Barnett et al., 2005), and snow cover increases the surface albedo, reducing the absorption of incoming solar energy (Zhang, 2005). However, heavy snowfall events have the potential to negatively impact local communities and infrastructure (Eisenberg and Warner, 2005; Scott et al., 2008; Fox et al., 2023).

Cloud phase determines the formation process of snow. Snow can originate from pure ice clouds, as well as supercooled liquid-containing clouds (sLCCs), where both supercooled water droplets and ice crystals coexist at temperatures between −40 and 0 °C (Jiusto and Weickmann, 1973; Korolev et al., 2017). These sLCCs dominate at latitudes higher than 45° and cover up to 20 %–30 % of the Earth, depending on the season (Warren et al., 1988; Matus and L'Ecuyer, 2017). The cloud phase affects not only the rate and intensity of snowfall but also the microphysical properties of the snowflakes that reach the ground (Jiusto and Weickmann, 1973; Liu, 2008). Clouds reflect solar radiation (cooling the surface) and trap terrestrial radiation (warming the surface). The balance between the two effects is partly influenced by the phase composition of the clouds due to distinct scattering properties of liquid and ice (Shupe and Intrieri, 2004; Ehrlich et al., 2009; Cesana et al., 2012). Despite having similar properties otherwise, a cloud with more liquid reflects more solar radiation back to space. Therefore, representing the cloud phase correctly in Earth system models (ESMs) has substantial implications for the simulation of cloud radiative properties (Matus and L'Ecuyer, 2017), especially in future climate projections with a warmer climate. This represents an important component of the cloud feedback, i.e., how clouds respond to changes in surface air temperature and how these changes, in turn, influence temperature. The cloud feedback is, at present, the most uncertain physical climate feedback, significantly contributing to the spread of climate sensitivity in ESMs (Zelinka et al., 2020). Among the most poorly understood cloud feedbacks is the one associated with cloud-phase changes (Bjordal et al., 2020; Zelinka et al., 2020).

Atmospheric warming can lead to cloud-phase changes through two different pathways: less cloud ice due to higher temperatures or more ice due to potentially more ice-nucleating particles (INPs). INPs are microscopic particles present in the air that serve as the starting point for ice crystal formation in clouds (DeMott et al., 2010). Several studies, including Murray et al. (2021), showed that increased temperatures, especially in polar regions, have caused a shift in the mixed-phase clouds towards higher latitudes and altitudes due to the ice reduction in the atmosphere in these regions. The shift in cloud phase towards more liquid and less ice leads to a reduction in the fraction of precipitation falling as snow in previously snowy areas in the Northern Hemisphere (NH), but snowfall events will happen further north in the future (Chen et al., 2020). This shift would result in an expected decrease in snowfall events and duration of the snowfall season for most regions in the NH (Danco et al., 2016; Chen et al., 2020). However, potentially partly counteracting this, the reduction in Arctic sea ice may also facilitate local emission of INPs (Carlsen and David, 2022). Furthermore, the well-established general increase in total precipitation with warming may also lead to increased snowfall in some regions (Douville et al., 2023). Hence, cloud phase, snowfall amounts, and snowfall frequency of occurrence will likely change with warming in ways that are both complex and currently poorly understood (Danco et al., 2016; Chen et al., 2020; Quante et al., 2021).

There are several processes governing snow formation within sLCCs. The transition from a fully liquid to a completely frozen cloud can follow various pathways (Costa et al., 2017) driven by differences in the saturation vapor pressures of ice and liquid at temperatures below 0 °C (Wegener, 1912; Bergeron, 1928; Findeisen, 1938), where ice crystals will grow at the expense of supercooled liquid water (Korolev and Mazin, 2003; Korolev, 2007; Storelvmo and Tan, 2015; Korolev et al., 2017). This phenomenon is called the Wegener–Bergeron–Findeisen (WBF) process and can lead to the rapid growth of ice crystals, which eventually fall out as snow. Increasing INP concentrations can cause cloud glaciation, amplified by the WBF process, increasing the precipitation at the surface and consequently shortening the cloud lifetime (Lohmann and Diehl, 2006; Rosenfeld et al., 2011; Storelvmo and Tan, 2015). Furthermore, ice crystals can be introduced into this sLCC by either falling into it from above (e.g., Proske et al., 2021) or by turbulent mixing.

The microphysical processes described above occur on scales smaller than the grid resolution used in ESMs, and their influence on cloud macroscopic properties must therefore be parameterized. Cloud microphysical parameterizations in ESMs are known to often be overly crude and simplistic. For example, ESMs generally assume that liquid and ice are uniformly mixed within a grid box, while observations have shown that sLCCs typically consist of cloud pockets exclusively composed of liquid or ice (Korolev and Milbrandt, 2022). Additionally, the treatment of primary ice production has an influence on the amount of supercooled liquid water in climate models (Vergara-Temprado et al., 2018). Only a small number of ESMs include parameterizations of primary ice production that depend on the presence of aerosols with the ability to act as INPs, while the majority include simpler parameterizations that rely only on temperature. Cesana et al. (2015) found that models with more complex microphysics (e.g., prognostic ice and liquid water content, heterogeneous freezing, riming, accretion, and the WBF process) tend to provide a more accurate representation of the ice phase. Similarly, Komurcu et al. (2014) noted that the variability in the cloud phase among models was influenced by the specific ice nucleation scheme used and the representation of other microphysical processes associated with ice.

ESMs have previously been shown to not represent the cloud phase accurately by often underestimating liquid and overestimating ice compared to satellite measurements. This has particularly been the case for high-latitude regions (Komurcu et al., 2014; Cesana et al., 2015; Tan and Storelvmo, 2016; Kay et al., 2016; McIlhattan et al., 2017; Bruno et al., 2021; Shaw et al., 2022). However, most studies to date were conducted based on the previous generation of ESMs (Coupled Model Intercomparison Project Phase 5 (CMIP5); Taylor et al., 2012), and a wide range of different cloud-phase metrics have been applied in the past, often without a clear strategy that allows like-for-like comparison with satellite observations. In other words, uncertainties remain in understanding the complex processes governing cloud phase and their representation in ESMs (Komurcu et al., 2014). Because cloud phase and snowfall are tightly linked through the processes outlined earlier, any inaccuracies in representing cloud phase could lead to biases in the simulation of snow growth, formation, and the precipitation reaching the ground in solid or liquid form (Mülmenstädt et al., 2021; Stanford et al., 2023). It is important to note, however, that while such biases in cloud-phase representation might exist, other compensating model biases could nevertheless lead to accurate precipitation simulations.

CloudSat (Stephens et al., 2002) and the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO; Winker et al., 2010), flying in the afternoon train constellation, have provided global estimates of cloud properties and snowfall since 2006. CloudSat is equipped with a cloud-profiling radar (CPR) that detects large cloud and precipitation particles, while CALIPSO, with its Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP; Winker et al., 2007), can, among other things, determine the phase of cloud layers. The CloudSat–CALIPSO constellation overpasses the same regions approximately every 16 d, providing long-term periodic monitoring of cloud and snowfall characteristics across the globe. Previous studies have contributed to a better understanding of the uncertainties associated with satellite measurements of clouds and precipitation (e.g., Stephens and Kummerow, 2007; Hiley et al., 2011). For example, Stephens and Kummerow (2007) identified two primary sources of uncertainty in retrieval methods, namely errors in distinguishing between cloudy and clear-sky scenes and between precipitating and non-precipitating clouds. Furthermore, the forward models used are highly sensitive to their input parameters, particularly the radiative transfer and atmospheric models. Hiley et al. (2011) demonstrated that snowfall retrievals are also influenced by retrieval assumptions and the use of different ice particle models, which can significantly affect the estimated snowfall rates. Despite limitations, the value of satellite data like CloudSat and CALIPSO for validating reanalysis and ESM data, which both rely on parameterized microphysics (Forbes and Ahlgrimm, 2014; McIlhattan et al., 2017; Milani et al., 2018; Edel et al., 2020; Daloz et al., 2020), is well established. Previous studies suggest that there is a notable difference in snowfall estimates between ESMs and satellite observations. For example, Heymsfield et al. (2020) found that the Met Office Unified Model and the Community Atmosphere Model 6 produce double the amount of snowfall relative to satellite observations. Additionally, Roussel et al. (2020) discovered that the ensemble median of CMIP5 experiments tends to show a positive bias in snowfall rates compared to the CloudSat average. Of particular relevance for the analysis presented here is the study by McIlhattan et al. (2017), who investigated the causes and impacts of liquid-containing cloud (LCC) biases in the Arctic region. They found that the Community Earth System Model Large Ensemble (CESM-LE) underestimates the frequency of Arctic LCCs compared to observations. Moreover, the CESM-LE overestimated the snowfall frequency from these clouds, possibly indicating an overactive WBF process leading to snowfall that is too frequent and cloud lifetimes that are too short (McIlhattan et al., 2017). In a comparable study of the Southern Hemisphere (SH), Roussel et al. (2020) showed a weak annual cycle of snowfall rates in the Antarctic plateau regions and found that the CMIP5 and CMIP6 (Eyring et al., 2016) simulations tended to overestimate the average precipitation rate in Antarctica. Despite the improvements in surface temperature representation from CMIP5 to CMIP6, there is no corresponding improvement in the representation of large-scale mean precipitation rate and seasonality in the region (Roussel et al., 2020). Precipitation errors in ESMs are not primarily driven by first-order physical links but by atmospheric circulation and cloud microphysics (Roussel et al., 2020). Antarctica and the surrounding ocean are also known to be a region in which ESMs have large biases, which are thought to be attributable to cloud-phase biases (e.g., Vergara-Temprado et al., 2018). These persistent discrepancies between ESMs and observations further highlight the importance of understanding the link between cloud phase and solid precipitation on different scales and improving the representation of cloud phase in ESMs.

Reanalysis products employ numerical weather prediction models that assimilate various observations to generate a continuous spatial and temporal dataset. The European Center for Medium-Range Weather Forecasts (ECMWF) reanalysis datasets use data collected from satellites, weather stations, and ocean buoys to combine with a numerical weather prediction model to provide a comprehensive view of the global atmospheric climate. Nonetheless, the precision in depicting cloud and surface snowfall relies on the underlying model and the assimilated dataset (Boisvert et al., 2018; Daloz et al., 2020; Boisvert et al., 2020). Challenges may arise due to limited spatial resolution and sparse observations, particularly in remote and complex topographical regions (Boisvert et al., 2018; Daloz et al., 2020; Boisvert et al., 2020). Wang et al. (2019) studied snowfall within two ECMWF reanalysis datasets, namely ERA-Interim (Simmons et al., 2007; Dee et al., 2011) and ERA5 (Hersbach et al., 2020). The latter includes higher horizontal resolution, improvements to the numerical model, improved data assimilation, and different cloud physics schemes (Forbes and Ahlgrimm, 2014). Wang et al. (2019) showed that ERA5 led to a better representation of snowfall than ERA-Interim. Nevertheless, Milani et al. (2018) concluded that ERA-Interim reanalysis data produce mean annual snowfall patterns similar in magnitude to CloudSat over Antarctica. In the Arctic, the ERA-Interim data qualitatively represented the interannual snowfall rates and seasonal cycle well but underestimated high snowfall rates significantly during summer and overestimated weak snowfall rates over open water compared to CloudSat (Edel et al., 2020). Seasonal biases tended to be higher in colder months when heavier snowfall occurred (Edel et al., 2020). Boisvert et al. (2020) compared different reanalysis products, including ERA-Interim and ERA5, and focused on snowfall in the Southern Ocean. They found similar spatial and interannual snowfall patterns, where ERA-Interim produced the least snowfall compared to ERA5 and the other reanalysis datasets. Sea ice representation, atmospheric moisture content, temperature, cloud microphysics schemes, and data assimilation used in the reanalysis were all identified to contribute to the disagreement regarding snowfall magnitude at lower latitudes (north of the Antarctic continent; Boisvert et al., 2020). They also showed that the difference between subsequent reanalysis iterations is due to these factors and model resolution changes. However, it is challenging to attribute the differences further to specific factors.

This study adds to the previous literature discussed above by connecting cloud phase and surface snowfall in ERA5 reanalysis, CMIP6 models, and CloudSat–CALIPSO observations through the analyses of the frequency of occurrence of sLCCs and associated surface snowfall in these datasets. The frequency of sLCC occurrence (f_sLCC) represents the fraction of time (in percent) that sLCCs are observed, while the frequency of snowfall from sLCC (f_snow) denotes the percentage of time that sLCCs are snowing. This information is crucial for understanding when and where sLCCs produce snowfall and is vital for predicting and modeling precipitation patterns in high-altitude and latitudinal areas. By studying the f_snow, valuable insights into the occurrence of snowfall can be gained, further enhancing our understanding of its role in the water cycle, as well as the representation of the WBF process in reanalysis and ESMs.

Another novel contribution is the assessment of cloud-phase biases based on a relatively new cloud-phase metric that aims for a like-for-like comparison with satellite retrievals to the extent possible and the relation of these biases to associated snowfall biases in reanalysis and ESMs. The study area spans the latitudinal ranges of 45–82° N and 45–82° S. We select these regions because they cover the mid-latitudes and polar areas in which mixed-phase clouds and surface snowfall are predominantly observed (Warren et al., 1988; Komurcu et al., 2014; Korolev et al., 2017; Matus and L'Ecuyer, 2017; Chen et al., 2020). The central research questions of this study are as follows: what are the cloud-phase biases in the NH and SH mid-to-high latitudes in ERA5 and “historical” simulations from CMIP6 with respect to the relatively new cloud-phase metrics, and how do these biases relate to snowfall biases? By addressing these research questions, we aim to improve our understanding of cloud-phase representation and its connection to snowfall and ultimately contribute to the advancement of climate modeling and prediction.

In the subsequent sections, we will investigate the frequency and geographic distribution of sLCCs and associated snowfall in order to distinguish their contribution to the overall snowfall patterns. Section 2 of this paper includes a description of the datasets and methodologies utilized in this study. Section 3 presents the results for the f_sLCC and f_snow values for the NH and SH mid-to-high latitudes between 2007 and 2010. Section 4 is a discussion on CloudSat–CALIPSO and the use of different time sampling for the various datasets. In Sect. 5, the connection between cloud-phase and snowfall biases is discussed, and possible next steps for future studies are given.

2 Data and methods

Our primary emphasis is on sLCCs in the mid-to-high latitudes and understanding how frequently and where they produce snowfall. We utilize CloudSat–CALIPSO satellite observations, ERA5 reanalysis, and CMIP6 simulations to assess biases in the spatial distribution and frequency of occurrence of sLCCs, as well as their respective snowfall frequency. This section introduces the different datasets and presents the statistical methods used for the cloud-phase and snowfall analysis. Table 1 lists the units and data properties for the variables used from CloudSat–CALIPSO, ERA5, and CMIP6, respectively.

Table 1 Overview of variable and original data field name, units, and data properties for CloudSat–CALIPSO, ERA5, and CMIP6 datasets. CMIP6 model variables represent variable short names.

2.1 CloudSat–CALIPSO satellite retrievals

In our investigation of the relationship between sLCCs and surface snowfall biases, we rely on satellite observations from CloudSat and CALIPSO for the comparison with ERA5 and CMIP6 models. Specifically, we use Release 5 (R05) versions of the 2B-CLDCLASS-LIDAR product (Sassen et al., 2008) for cloud-phase determination, the 2C-SNOW-PROFILE product (Wood and L'Ecuyer, 2018) for the surface snowfall rate estimation, and 2 m temperatures from the ECMWF-AUX dataset that contains ancillary ECMWF state variable data interpolated to each vertical radar bin (see Table 1).

The 2B-CLDCLASS-LIDAR product utilizes the CALIOP lidar and the CPR on CloudSat and their different sensitivities to liquid droplets and ice crystals in order to retrieve the phase of a cloudy layer. The algorithm uses a temperature-dependent radar reflectivity threshold (Zhang et al., 2010), the integrated attenuated lidar backscattering coefficient, and cloud-base and cloud-top temperatures from atmospheric reanalysis data to discriminate between ice, mixed, or liquid water clouds (Wang and Sassen, 2019). In accordance with McIlhattan et al. (2017), we focus on the phase of the bottom cloud layer within the atmospheric column and neglect the cloud-phase quality flag, irrespective of the confidence level indicated. The rationale for the latter is the lidar's ability to robustly detect any liquid (liquid or mixed), and a low-confidence value for the phase determination could stem from the uncertainty in the distinction between purely liquid clouds and mixed-phase clouds. This is particularly true for the frequent structure of polar mixed-phase clouds with the liquid layer at cloud top (McIlhattan et al., 2017).

To estimate the surface snowfall rate, we use the 2C-SNOW-PROFILE product (Wood and L'Ecuyer, 2018), which is based on an optimal estimation algorithm to retrieve profiles of parameters of the snow size distribution. The optimal estimation uses the radar reflectivity profile of the snow layer, ancillary meteorological information, and assumptions about snow microphysical and scattering properties. As the CPR cannot reliably measure near-surface reflectivities due to ground clutter, the surface snowfall rate is estimated based on the lowest clutter-free radar bin (Wood et al., 2014). The truncation height of the snow profile due to ground clutter is surface-dependent (two range bins above surface and over ocean (∼500 m); four range bins over land and sea ice (∼1000 m); Wood et al., 2014).

We focus on the years 2007 to 2010 before CloudSat switched to daytime-only operations due to a battery malfunction. Furthermore, we exclude September 2008 and December 2009 due to insufficient CALIOP data and the CloudSat battery failure (Keys, 2010). We aggregate the profile-by-profile data to a horizontal grid with a resolution of 3.75°×1.9° for each month. Subsequently, we calculate the monthly f_sLCC and f_snow, as described in Sect. 2.4.

The sea ice concentration (SIC) data used in this study are obtained from the Institute of Environmental Physics, University of Bremen. The data are derived using the ARTIST Sea Ice (ASI) algorithm developed by Spreen et al. (2008). The ASI retrieval method utilizes microwave radiometer data from the Advanced Microwave Scanning Radiometer for the Earth Observing System (EOS) instrument on the Aqua satellite and the Advanced Microwave Scanning Radiometer 2 instrument on the Global Change Observation Mission–Water (GCOM-W1) satellite. Both datasets were reprocessed in 2018 with the same parameters. To present the 20 % SIC, a seasonal average of the SIC is calculated for the years 2007 to 2010 and then gridded onto the 3.75°×1.9° grid using bilinear interpolation.

2.2 ERA5 reanalysis

The ERA5 reanalysis data employ the 4D-Var data assimilation of the ECMWF Integrated Forecast System. ERA5 incorporates an improved stratiform cloud and precipitation scheme, enhancing the representation of mixed-phase clouds compared to ERA-Interim (Hersbach et al., 2020) and includes prognostic variables for both rain and snow (Forbes et al., 2011; Forbes and Tompkins, 2011; Forbes and Ahlgrimm, 2014). The meteorological values are output at a resolution of 0.25×0.25°.

In order to detect the presence of sLCCs in ERA5 data, specific criteria involving the 2 m temperature (T2m) threshold and liquid water path (LWP) threshold must be met (see Sect. 2.4). However, it is worth noting that ERA5 does not readily provide LWP values. Therefore, we calculated the LWP using the total column cloud liquid water (tclw) and the total column rainwater (tcrw; expressed in kg m⁻²; Table 1). We incorporate tcrw with the T2m threshold below 0 °C as this threshold is used to exclude any rainwater below the melting layer. Using daily mean values of tclw and tcrw and applying the daily man temperature threshold, we can analyze the role of supercooled liquid water within clouds and the contribution of liquid water to the snowfall precipitation process in ERA5.

To calculate the f_snow, we rely on the mean snowfall rate variable (msr, Table 1). The mean snowfall rate (expressed in units of $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) provides a combined measurement of both large-scale and convective snowfall. It is calculated as an average over 24 h. Here we use the ERA5 Climate Data Store tool to average all of the output variables into daily mean values based on hourly data (Cucchi et al., 2021; Copernicus Climate Change Service, 2018). It is important to note that 2t, tclw, and tcrw (Table 1) are reported as instantaneous values every hour, meaning that when they are averaged to obtain daily mean values, they may not accurately represent the actual daily mean values. Nevertheless, a non-zero model daily mean LWP means that there was an sLCC in the grid box at some point during the day. In contrast, the mean snowfall rate is reported as a temporal average and accurately represents the daily mean snowfall rate.

2.3 CMIP6 models

Our intermodel comparison relies on data from CMIP6, involving simulations from various modeling centers. Specifically, we focus on the “historical” simulations (1850–near-present; Eyring et al., 2016), selecting 10 models with daily mean outputs from different ensemble members, as outlined in Table A1. In most of the figures shown in this study, we present the CMIP6 multimodel mean.

In our analysis, we consider the following variables (Table 1): the mass mixing ratio of cloud liquid water (clw) encompassing liquid water in both large-scale and convective clouds provided at each model level (in units of kg kg⁻¹). Precipitating hydrometeors within the clw variable are included only if they affect the model's radiative transfer calculations. Near-surface air temperature (tas), reflecting the temperature at a height of 2 m above the surface (measured in K), and snowfall flux (prsn), representing precipitation in all forms of frozen water (expressed in $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) are also considered.

As none of the models listed in Table A1 directly provides LWP as a daily output variable, we derive LWP based on clw. To calculate the LWP, we interpolate the CMIP6 hybrid sigma pressure to isobaric pressure (Appendix A). The UKESM1-0-LL and HadGEM3-GC31-LL models provide values on orographic vertical coordinates. Consequently, for these models, we use the CMIP6 variable phalf (Table 1) instead of calculating isobaric pressure levels. Subsequently, we utilize the hydrostatic equation to calculate the liquid water content for each vertical grid box. The liquid water content is then summed up over the vertical to obtain the LWP for each horizontal grid box (see Appendix A).

CMIP6 data are gathered from 2006 to 2009 as these years represent the closest available overlapping time range for all CMIP6 models (some models ended the simulations by December 2009 and are marked with an asterisk in Table A1). The slight mismatch in the time range is of limited relevance as CMIP6 model simulations are not designed to reproduce the exact temporal evolution of past weather but instead generate their own internal variability (e.g., El Niño–Southern Oscillation (ENSO) cycles), so they should only be viewed as generally representative of the time period in question.

2.4 Calculation of f_sLCC and f_snow

The instrument sensitivities of the CPR and CALIOP, as well as their different temporal and spatial sampling in contrast to the gridded ERA5 and CMIP6 data, require different metrics that ensure a like-for-like comparison to the extent possible. In the following, we describe the calculation of these metrics (f_sLCC and f_snow) and call attention to differences between the satellite observations and ERA5 or CMIP6. For ERA5 and for CMIP6 models, we remap all models to the coarsest resolution of 3.75°×1.9°, which corresponds to the grid of IPSL-CM5A2-INCA (Table A1), with a nominal resolution of approximately 500 km. We use area-weighted averages for both the NH and SH to calculate spatial means.

Table 2 Calculation of f_sLCC and f_snow and overview of applied thresholds. The first row corresponds to CMIP6/ERA5 metrics, while rows denoted with “CC” describe CloudSat–CALIPSO metrics.

2.4.1 Frequency of sLCC

We defined the occurrence of an sLCC in a cloudy CloudSat–CALIPSO profile when the lidar detected any liquid in the lowermost cloud layer (phase flag “liquid” or “mixed”) and when the surface temperature was below freezing (Temperature_2m≤0 °C). The temperature threshold was applied to ensure that the clouds were supercooled. In ERA5 and CMIP6, we applied the same temperature threshold and defined an sLCC when the LWP is above 5 g m⁻². This is in accordance with McIlhattan et al. (2017) and based on their sensitivity estimate when comparing CloudSat with ground-based microwave radiometer observations at Summit, Greenland. This definition of sLCCs is likely a very conservative estimate of the actual frequency of clouds containing supercooled liquid in the atmosphere because we require temperatures below freezing at the surface, while, in reality, there could still be supercooled liquid even if the surface is above freezing. However, as we are investigating the frequency of snowing sLCCs in the second part of this study, we want to ensure that any precipitation from the sLCCs is primarily in the form of snow.

In ERA5 and CMIP6, we calculate f_sLCC as the number of days with sLCCs present (N_sLCC) divided by the total number of days (N_all) over a specific time period (see Table 2). However, we cannot utilize daily statistics on a global scale from CloudSat–CALIPSO due to their transect-based sampling in combination with a revisit time of 16 d and, consequently, insufficient horizontal coverage (e.g., Kotarba, 2022; von Lerber et al., 2022). Thus, after aggregating the CloudSat–CALIPSO profiles to the same grid as ERA5 and CMIP6, we calculate f_sLCC from satellite data by dividing the number of profiles in that month with an sLCC by the total number of profiles (see Table 2).

2.4.2 Frequency of snowfall in sLCCs

For the occurrence of snowfall from sLCCs, we define a snowfall threshold of $\mathrm{0.01}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$. Again, this is based on McIlhattan et al. (2017) and aims to mitigate biases due to instrument sensitivities. The f_snow in ERA5 and CMIP6 is then calculated as the number of days with snowfall (N_sLCC_sf) divided by the total number of days with sLCCs present (N_sLCC; see Table 2). The ERA5 and CMIP6 snowfall rates and T2m are given as daily mean values. It is important to note that as f_sLCC is calculated using daily mean values, the simulated precipitation can, in principle, occur as rain (or supercooled rain), depending on the temperature at the time of the precipitation. However, as no additional information is available, we here assume that all precipitation from sLCCs is in the form of snow. To ensure comparability among the three datasets (CloudSat–CALIPSO, ERA5, and CMIP6), we perform a unit transformation, as outlined in Appendix B.

To account for the differences in sampling in space and time, we calculate f_snow from CloudSat–CALIPSO as the number of observations with snowing sLCCs divided by the number of observations with sLCCs (Table 2).

3 Results

The following subsections examine and intercompare f_sLCC and f_snow, respectively, from the three datasets examined (CloudSat–CALIPSO, ERA5, and the CMIP6 models).

3.1 Frequency of sLCCs

Figure 1a–d present the seasonal variability in f_sLCC observed by CloudSat–CALIPSO in the mid-to-high latitudes of the NH. Evident from the figure is a clear seasonal progression in the spatial distribution of f_sLCC. In summer (June–August, JJA), f_sLCC reaches its minimum spatial average, with non-negligible values primarily within the Arctic basin. In NH autumn (September–November, SON), f_sLCC values increase significantly in the Arctic basin but also spread poleward to land areas at about 60° N. In boreal winter (December–February, DJF), maximum values shift away from the Arctic basin and are instead found over the ocean close to the average seasonal sea ice edge (Fig. 1a; red line).

Figure 1 The NH mid-to-high-latitude seasonal averages of f_sLCC. Panels (a)–(d) show CloudSat–CALIPSO, panels (e)–(h) display ERA5 data, and panels (i)–(l) show the CMIP6 model mean. Each map includes an area-weighted average for the study area (lower-left corner). These averages are calculated for areas in which CloudSat–CALIPSO have valid observations (between 45–82° N) and exclude the dotted area (e–l). The dashed black line represents the seasonal mean 2 m temperature 0 °C isotherm for each individual product. The red line (a–d) shows the average sea ice edge of 20 % sea ice concentration (SIC) between 2007 and 2010 for the given season.

The Greenland ice sheet (GRIS) is a region of particular interest for this study, as sLCCs have been proposed to accelerate GRIS surface melt (Bennartz et al., 2013; Hofer et al., 2019), while snowfall represents the primary ice sheet growth mechanism. The GRIS, which is characterized by its high elevation and snow-covered surface, experiences numerous sLCCs during boreal summer and autumn due to sufficient moisture and temperatures below 0 °C in the CloudSat–CALIPSO observations (Fig. 1c, d). However, during boreal winter and spring (March–May, MAM), the f_sLCC is low (< 15 %) as the region becomes too dry and cold (Shupe et al., 2013) to support the formation of sLCCs. The seasonal cycle of the f_sLCC as observed by CloudSat–CALIPSO is largely what would be expected from the combined seasonal influence of atmospheric circulation, moisture availability, and temperature conditions in the mid-to-high-latitudes in the NH.

Comparing the CloudSat–CALIPSO observations with ERA5 and the CMIP6 model mean, we generally observe that ERA5 and the CMIP6 models show a similar seasonal progression of f_sLCC to that seen in the satellite observations. However, both ERA5 and the CMIP6 mean overestimate f_sLCC to various extents during all seasons. The most significant discrepancy is observed during boreal spring, with area-weighted-averaged differences between CloudSat–CALIPSO and ERA5 and the CMIP6 model mean of −11 % and −14 %, respectively (Fig. C1f, j). Nevertheless, the spatial extent of areas with f_sLCC amounts larger than 0 % in CloudSat–CALIPSO is replicated well in ERA5 and the CMIP6 model mean (Fig. 1). Although ERA5 and the CMIP6 models have remarkably similar spatial patterns of f_sLCC, the CMIP6 model mean generally exhibits slightly larger area-weighted biases in the NH (6 %–14 %, Fig. C1i–l) than ERA5 (1 %–11 %, Fig. C1e–h) when compared to CloudSat–CALIPSO.

The fact that these discrepancies are about equally present in both ERA5 and the CMIP6 models provides valuable insights into their potential causes. They are most likely linked to the microphysical parameterizations of cloud processes that govern the cloud phase in ERA5 and CMIP6 models. It is reasonable to assume that the temperature in the ECMWF-AUX product used in CloudSat–CALIPSO is quite similar to the ERA5 daily mean as indicated by the seasonal mean 0 °C isotherm in Figs. 1 and 2. However, ERA5 shows a slight variation in the 0 °C isotherm line over central Europe during DJF compared to ECMWF-AUX (Fig. 1e). Furthermore, a comparison of the 2 m temperature between ECMWF-AUX and ERA5 shows a latitudinal average difference of 0.24 K ± 0.22 K (Fig. D1). Atmospheric circulation and overall cloud cover should be well constrained by the observations used in the ERA5 reanalysis. However, the cloud phase is not as evident in the identified biases. This finding will be discussed in greater depth in Sect. 5.1.

Figure 2 The SH mid-to-high-latitude seasonal averages of f_sLCC. Panels (a)–(d) show CloudSat–CALIPSO, panels (e)–(h) display ERA5 data, and panels (i)–(l) show the CMIP6 model mean. Each map includes an area-weighted average for the study area (lower-left corner). These averages are calculated for areas in which CloudSat–CALIPSO have valid observations (between 45–82° S) and exclude the dotted area (e–l). The dashed black line represents the seasonal mean 2 m temperature 0 °C isotherm for each individual product. The red line (a–d) shows the average sea ice edge of 20 % sea ice concentration (SIC) between 2007 and 2010 for the given season.

Unlike the NH mid-to-high latitudes, the SH does not experience significant seasonal variations in f_sLCC, according to the CloudSat–CALIPSO observations (Fig. 2a–d). The f_sLCC remains relatively constant across all seasons and is the highest in a band bounded by the Antarctic continent in the south and approximately the 60° parallel in the north, with values between 17 %–24 %. However, the area-weighted averages of f_sLCC in the SH are generally higher than in the NH. The relatively constant f_sLCC indicates that a more persistent cloud regime exists in the SH, with extensive sLCC cover present all year round. This is consistent with the previous literature reporting that the Southern Ocean is the region of the world with the most extensive mixed-phase cloud cover (Matus and L'Ecuyer, 2017).

Like the GRIS, the Antarctic ice sheet has a low f_sLCC amount, with the lowest frequency in austral winter (JJA; ≤15 %). With its high elevations, the East Antarctic Ice Sheet exhibits a lower f_sLCC amount than the flatter West Antarctic Ice Sheet during austral summer (DJF), autumn (MAM), and spring (SON). Just like the NH, the SH shows a reduced area of the f_sLCC in austral summer (Fig. 2c) but with a larger extent towards lower latitudes (NH: > 70° N; SH: > 60° S). In all seasons, f_sLCC activities are observed primarily southward of the Antarctic Circumpolar Current, which may be due to the presence of warm water that results in surface temperatures above the sLCC threshold (Fig. 2a–d).

In ERA5 and the CMIP6 model mean, there is a prominent gradient of sLCCs between the Southern Ocean and the Antarctic ice sheet, with f_sLCC of up to 100 % over the Southern Ocean and < 15 % over the Antarctic ice sheet. For the SH spatial mean, ERA5 and the CMIP6 model mean moderately overestimate f_sLCC (1 %–8 %), with the CMIP6 model mean, in general, performing slightly better than ERA5 in comparison to CloudSat–CALIPSO (Fig. C2).

Figure 3 The annual cycle of monthly f_sLCC for (a) NH (45–82° N) and (b) SH (45 and 82° S). Shown are area-weighted averages for CloudSat–CALIPSO (black), ERA5 (pink), and the CMIP6 model mean (green). Each box represents the interquartile range (IQR) from the 25th (Q1) to the 75th (Q3) percentile. The whiskers extend to the minimum and maximum values, defined as $\mathrm{Q}\mathrm{1}-\mathrm{1.5}\times \mathrm{IQR}$ and $\mathrm{Q}\mathrm{3}+\mathrm{1.5}\times \mathrm{IQR}$, respectively. Any points falling outside these whiskers are considered outliers and are marked with crosses. The green dots represent the minimum and maximum CMIP6 model values over all years.

The biases in ERA5 and the CMIP6 models regarding f_sLCC, relative to CloudSat–CALIPSO, are further demonstrated through the area-weighted averages for each month in both hemispheres (Fig. 3). For the area-weighted averages, we consider only the locations for which CloudSat–CALIPSO provided valid data (values between |45°–82°|). CMIP6 models generate a large spread of f_sLCC, with both CloudSat–CALIPSO and ERA5 falling within the range of the model maximum and minimum values (green dots in Fig. 3). Again, it is evident that ERA5 and the CMIP6 model mean have higher f_sLCC than CloudSat–CALIPSO in the NH (Fig. 3a). The SH area-weighted averages show a better agreement between CloudSat–CALIPSO, ERA5, and the CMIP6 model mean than in the NH (Fig. 3b). The CMIP6 model mean underestimates f_sLCC in comparison to CloudSat–CALIPSO during austral summer and autumn but overestimates during winter and spring, while ERA5 overestimates f_sLCC to various degrees all year round in the SH.

Figure 4 Magnitude of seasonal area-weighted averages (between |45–82°|) of the f_sLCC for CloudSat–CALIPSO, ERA5, and the CMIP6 models (numbers). The heatmap colors correspond to the differences in the area-weighted averages of f_sLCC between CloudSat–CALIPSO and ERA5 and CMIP6 models. Green (pink) values indicate the underestimation (overestimation) of the individual model with respect to CloudSat–CALIPSO. Panel (a) shows the mid-to-high-latitude NH, and panel (b) shows the mid-to-high-latitude SH. Per season, the smallest (largest) seasonal and spatial differences are outlined with blue (red) lines.

Finally, to assess the model performance relative to CloudSat–CALIPSO for ERA5 and the individual CMIP6 models for f_sLCC, we analyze the difference in the seasonal and spatial averages for each of the 11 models (Fig. 4) for the NH and SH mid-to-high latitudes. The most significant discrepancy is seen in boreal winter and spring for AWI-ESM-1-1-LR, MPI-ESM1-2-LR, and IPSL-CM6A-LR with a negative difference > 15 % (Fig. 4a). In the NH, MPI-ESM1-2-LR is generally the model with the larger overestimation (> 15 %), despite a difference of < 10 % in boreal summer. IPSL-CM5A2-INCA has a low difference, < 5 %, in boreal summer, autumn, and winter. As previously highlighted, the difference values for f_sLCC are lower in the SH (Fig. 4b), and none of the models stands out in terms of overall performance in either direction.

3.2 Frequency of snowfall from sLCC

To assess the ability of ERA5 reanalysis and the CMIP6 models to represent snowfall processes, we compare the f_snow between them and the CloudSat–CALIPSO observations. Based on CloudSat–CALIPSO observations, NH sLCCs frequently produce snowfall, as can be seen by the non-zero values of f_snow in Fig. 5a–d. However, during boreal summer, the majority of sLCCs rarely produces snowfall, as shown by the relative increase in the areal extent of zero f_snow values (Fig. 5c). Notably, during boreal winter, spring, and autumn, the seasonal average f_snow is highest over open-ocean regions, such as the Greenland, Norwegian, and Barents seas (30 %–65 %), the Labrador Sea (40 %–80 %), and the Bering Sea (40 %–100 %; Fig. 5a, b, d). This is likely due to plentiful available moisture for snowfall to occur and is possibly linked to cold-air outbreaks (CAOs; Young et al., 2017). CAOs arise when cold and dry air flows over warmer and moist sea surfaces. CAOs occur frequently in the North Atlantic in boreal winter and spring, and the clouds associated with them are predominantly sLCCs (Fletcher et al., 2016; Papritz and Sodemann, 2018; Geerts et al., 2022; Mateling et al., 2023). While the influence of CAOs is also visible in the f_snow patterns across various seasons, with the exception of boreal summer, the f_snow values over land areas tend to be lower. Although it is possible that the higher minimum height bin of the radar used over land (∼1000 m) may lead to the 2C-SNOW-PROFILE product missing some of the precipitation over land, the variability in the f_snow values in regions where similar cloud types are expected suggests that this land–sea contrast is likely a representative feature. This discrepancy can be attributed to the limited moisture content and shallower boundary layer over land during winter, resulting in lower f_snow compared to ocean areas where CAOs provide deeper moisture-rich clouds. Once clouds embedded in CAOs reach a certain distance from the sea ice, they frequently produce precipitation (e.g., Abel et al., 2017). In contrast, during the winter months, northern Europe is often covered by non-precipitating supercooled stratus clouds (Cesana et al., 2012; McIlhattan et al., 2017). The models also show a lower f_snow over land, supporting the idea that the observed land–sea contrast is a genuine characteristic rather than an artifact of measurement techniques.

Figure 5 The figure presents the seasonal averages of f_snow in the NH mid-to-high latitudes. The layout and area-weighted averages are calculated the same as those shown in Fig. 1.

In contrast to the CloudSat–CALIPSO observations, both ERA5 and the CMIP6 model mean have much higher f_snow values (> 60 %), which are evenly distributed over the NH mid-to-high latitudes, without a significant seasonal variability. In regions where all three datasets have sLCCs, ERA5 and the CMIP6 model mean overestimate the f_snow by ∼50 % (Fig. C3). This indicates that ERA5 and the CMIP6 models produce snowfall much more frequently from sLCCs than observed. It is important to note that since f_snow accounts for the number of sLCCs, ERA5 and the CMIP6 models are too efficient at producing snowfall, regardless of the number of sLCCs that are simulated/assessed.

Figure 6 The figure presents the seasonal averages of the f_snow in the SH study region. The layout and area-weighted averages are calculated in the same way as for those shown in Fig. 1.

Regarding the SH mid-to-high latitudes, the CloudSat–CALIPSO observations show that there is no significant seasonal variation in f_snow, with spatially averaged values of ∼20 % for all seasons (Fig. 6a–d). However, when comparing the results from ERA5 and the CMIP6 model mean with observations from CloudSat–CALIPSO, we find that similar to the NH, the simulated f_snow values are much higher than those observed (Fig. 6). When only considering regions where CloudSat–CALIPSO observe f_snow (45–82° S), the area-weighted seasonal average of f_snow is approximately 60 %–70 % higher in ERA5 and the CMIP6 model mean (Fig. C4). The high simulated values of f_snow in both the NH and SH indicate that ERA5 and the CMIP6 models produce snowfall far too frequently from sLCCs compared to the satellite observations. In fact, when sLCCs are present, snow is produced practically all of the time (between 70 % and 90 %). This bears resemblance to a well-established bias in ESMs, namely the persistent “perpetual drizzle” problem for warm liquid clouds (Mülmenstädt et al., 2020; Lavers et al., 2021).

Figure 7 Annual cycle of f_snow for (a) the NH, which is defined here as the region between 45 and 82° N and (b) SH, defined here as the region between 45 and 82° S. Colors and boxplots are identical to Fig. 3.

To investigate the annual cycle of f_snow and the interannual variation between 2007 and 2010, we present the area-averaged data based on grid boxes in which CloudSat–CALIPSO is capable of making observations (|45°–82°|; Fig. 7). CloudSat–CALIPSO displays a relatively clear seasonal cycle in f_snow in the NH but not in the SH. The difference in the f_snow between boreal summer months and the rest of the year in the NH is 15 %, while the SH value is constant around 20 % throughout the year. In contrast, ERA5 and the CMIP6 model mean display a more significant seasonal cycle, with a difference of 25 % between boreal summer and winter months (Fig. 7a). Also, in the SH, there is more of a seasonal cycle evident for ERA5 and the CMIP6 model mean. Also worth noting is that ERA5 shows more pronounced interannual variability in both hemispheres of f_snow for a given month, meaning that the year-to-year fluctuations are larger than for CloudSat–CALIPSO and the CMIP6 model mean.

Figure 8 Magnitude of seasonal area-weighted averages (between |45–82°|) of f_snow for CloudSat–CALIPSO, ERA5, and the CMIP6 models (numbers). The heatmap colors correspond to the differences in the area-weighted averages of f_snow between CloudSat–CALIPSO and ERA5 and the CMIP6 models. Green (pink) values indicate underestimation (overestimation) of the individual model with respect to CloudSat–CALIPSO. Panel (a) shows the mid-to-high-latitude NH, and panel (b) shows the mid-to-high-latitude SH. Per season, the smallest (largest) seasonal and spatial differences are outlined with blue (red) lines.

The comparison of the magnitude and model differences for f_snow in Fig. 8 indicates how well the individual ESMs represent the frequency of occurrence for snow and if the same models that perform well for f_sLCC (Fig. 4) do so for surface snowfall as well. Among the models, AWI-ESM-1-1-LR and MPI-ESM1-2-LR come closest to matching CloudSat–CALIPSO observations, but their f_snow values are still much too high (NH: ∼30 %; SH: ∼45 %; Fig. 8). At the same time, IPSL-CM5A2-INCA deviates the most for f_snow in the NH and SH, having a negative difference of ∼65 % and ∼75 %, respectively. Interestingly, it appears that some models that perform well for f_sLCC perform poorly for f_snow. This could be an indication that models that are able to simulate the right frequency of occurrence for sLCCs can do so because they are converting cloud condensate to snow too readily. In other words, the models may be getting the right answer for the wrong reason.

4 Sensitivity tests

In the previous section, we examined the frequency of occurrence of sLCCs (f_sLCC) and the frequency of occurrence of surface snowfall (f_snow) from sLCCs in CloudSat–CALIPSO, ERA5, and CMIP6 data. We found that the reanalysis and ESMs overestimate sLCCs and the frequency of surface snowfall in comparison to CloudSat–CALIPSO. These biases in the CMIP6 mean values have potentially significant implications for our ability to predict how sLCCs and snowfall might change with future warming. At the same time, it is important to note that although the CMIP6 multimodel mean has this overestimation, some of the individual ESMs are much closer to the observations (green dots in Fig. 3), suggesting that these members may have more representative changes in clouds and snowfall in the future. It is therefore important to ensure that these findings are robust and not overly reliant on subjective decisions or limitations in the design of the comparison. In the following subsections, we investigate whether the identified discrepancies between the observations on the one hand and ERA5 and the CMIP6 models on the other hand can be explained by sampling biases or instrument sensitivity.

It is important to acknowledge that CloudSat–CALIPSO data come with inherent uncertainties, including those associated with snowfall retrievals (Stephens and Kummerow, 2007; Hiley et al., 2011; Schirmacher et al., 2023). CloudSat surface snowfall is subject to sampling biases and ground clutter issues and also relies on ECMWF temperature data to differentiate between snowfall and rainfall (Boisvert et al., 2020). Milani et al. (2018) found that applying adjustments and a temperature threshold to the CloudSat snowfall retrieval led to a decrease in the estimated occurrence of snowfall events by up to 30 %, primarily in the ocean regions surrounding Antarctica. Although these adjustments did not have the same effects everywhere, this highlights the sensitivity of the CloudSat retrievals to the assumptions made within them.

In addition to the uncertainties in the assumptions used in retrievals, comparing satellite observations with ESMs is challenging due to several factors linked to the sampling bias of the CloudSat–CALIPSO mission. Numerous considerations come into play when deriving climate statistics from satellite transects (Kotarba, 2022). These encompass the narrow swath coverage, infrequent revisits, and latitude-dependent ground track density (Kotarba, 2022; von Lerber et al., 2022). Varying the domain size and timescales of CloudSat observations showed that data on smaller scales or shorter time periods can introduce significant uncertainties associated with cloud, aerosol, or atmospheric properties (Henderson et al., 2013). However, instantaneous errors (using satellite data on smaller scales or shorter periods) tend to cancel out over longer time periods (Henderson et al., 2013; von Lerber et al., 2022). Another limitation of space-based remote-sensing observations is the uncertainty associated with the retrieval method employed to derive the snowfall rate from radar reflectivity measurements (Kulie and Bennartz, 2009; Milani et al., 2018). However, in our study, we examine the frequencies of occurrence and only use the surface snowfall amount to classify if an sLCC is snowing ($\mathrm{sf}\ge \mathrm{0.01}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) or not. Such a binary variable should be less sensitive to the exact snowfall rate retrieval. Edel et al. (2020) presented a CloudSat snowfall climatology, the frequency of snow, and snowfall rates over the Arctic. They found that the distribution of snowfall rate does not always match the distribution of the snowfall frequency in CloudSat. Nevertheless, future modeling studies should ideally use the Cloud Feedback Model Intercomparison Project Observation Simulator Package (COSP; Bodas-Salcedo et al., 2011) to make satellite observations and model datasets more comparable. At present, the COSP output is not available for a sufficient number of models and at a high enough temporal resolution to be useful for the present study.

Nevertheless, these satellite observations serve as valuable tools for identifying potential biases in reanalysis and ESMs, as demonstrated by previous studies (e.g., McIlhattan et al., 2017; Milani et al., 2018; Daloz et al., 2020; Heymsfield et al., 2020; Boisvert et al., 2020). To test the robustness of the overall findings of this study to our comparison approach, next we want to show how different temporal resolutions of ERA5, as well as different LWP thresholds applied to the CMIP6 data, impact the results.

4.1 Sensitivity to ERA5 temporal resolution

In the comparison by McIlhattan et al. (2017), the analysis was performed on a 6 h instantaneous model output to define the occurrence of LCCs in CESM-LE. As we compare ERA5 and CMIP6 models to CloudSat–CALIPSO, we only have daily means of LWP available. Thus, instead of four instantaneous values per day for a given model grid box, we use only one daily mean value of the LWP. In this way, we are comparing 30 daily values from ERA5 and CMIP6 models with thousands of CloudSat–CALIPSO profiles within each grid box to calculate f_sLCC. Due to the 16 d repeat cycle of CloudSat–CALIPSO, we cannot perform the satellite analysis on a daily temporal resolution at the given spatial resolution of the models. While CMIP6 output is not available at a higher temporal resolution, for ERA5, we can actually test the sensitivity using an even higher time resolution than McIlhattan et al. (2017). We take the hourly outputs from ERA5 from July and November 2007 to 2010 and create monthly mean values. July and November are chosen as these are the months with the lowest and highest f_sLCC in CloudSat–CALIPSO (Fig. 3a). We find that the f_sLCC values emerging from ERA5 data are not very sensitive to the output frequency, as shown in Fig. 9d–e and i–j. The area-weighted difference between CloudSat–CALIPSO and ERA5 hourly values is negligible in July (1 %) and ∼15 % in the inner Arctic and reduces from 12 % overestimation to 7 % in November. However, the conclusion that ERA5 overestimates the f_sLCC does not change by changing to a higher time resolution (Fig. 9). While the sensitivity analysis for ERA5 data indicates that changes in output frequency do not significantly affect the results, the same conclusions cannot be directly extrapolated to CMIP6 models. However, based on our analysis, it is reasonable to assume that similar results could be observed for individual ESMs within CMIP6, although this cannot be confirmed.

Figure 9 NH mid-to-high-latitude monthly averages of f_sLCC for July (a–e) and November (f–j) 2007 to 2010. The columns show the CloudSat–CALIPSO (CC) monthly means (a, f), ERA5 monthly means based on daily average values (b, g), and ERA5 monthly means based on hourly values (c, h). The fourth and fifth columns are the difference plots between CloudSat–CALIPSO and the ERA5 daily averages (d, i) and ERA5 hourly values (e, j), respectively. In the last two columns, green (pink) values indicate underestimation (overestimation) with respect to the reference used. Area-weighted averages for the study area (45–82° N) are located in the lower-left corner of each map and exclude the dotted area (b, c, g, h).

4.2 Sensitivity to LWP threshold for CMIP6 models

McIlhattan et al. (2017) conducted a similar analysis to that in this study by examining the frequency of Arctic LCCs and f_snow. They experimented with different LWP thresholds. Specifically, they used the 5 and 0.01 g m⁻² thresholds and based their LWP threshold of 5 g m⁻² on the approximate retrieval uncertainty in the LWP from ground-based microwave radiometer observations. We performed sensitivity tests in which we varied the threshold between 3, 5, 10, and 15 g m⁻² (Fig. 10). While some of these values are unreasonably high, it is nevertheless useful to test how adjusting the LWP threshold value affects the identified bias between the observations and simulations. As expected, Figs. 10 and E1 indicate that as the LWP threshold increases, the overestimation of f_sLCC decreases in the NH and SH. However, unless we increase the LWP threshold in ERA5 and in the CMIP6 models to an unrealistically high value, the general findings and conclusions hold – ERA5 and the CMIP6 model mean still overestimate the occurrence of sLCCs in the study area in comparison to CloudSat–CALIPSO.

Figure 10 Annual average difference plots (CloudSat–CALIPSO (CC) minus model) for f_sLCC in the NH mid-to-high latitudes. The first row represents ERA5 data (a–d), and the second row shows the CMIP6 model mean (e–h) difference. LWP thresholds of 3 g m⁻² (a, e), 5 g m⁻² (b, f), 10 g m⁻² (c, g), and 15 g m⁻² (d, h) were applied. The analysis in this study focuses on the LWP threshold of 5 g m⁻². The maps in each row are accompanied by area-weighted averages for the study area (located in the lower-left corner of each map). The averages are calculated for areas where CloudSat–CALIPSO observations have valid f_sLCC values.

5 Discussion and conclusion

We find that mid-to-high-latitude sLCCs and snowfall are produced more frequently in ERA5 and CMIP6 than in CloudSat–CALIPSO. While previous studies have focused on the underestimation of the supercooled liquid fraction (SLF) in mixed-phase clouds (Komurcu et al., 2014; Cesana et al., 2015; Tan and Storelvmo, 2016; Kay et al., 2016; Bruno et al., 2021; Shaw et al., 2022), this research focuses on the frequency of occurrence of sLCCs. This difference in the cloud-phase metric can lead to seemingly contradicting conclusions. We illustrate why our results do not necessarily contradict previous findings with the following example: if an ESM consistently predicts a low percentage of liquid in clouds, while satellite observations suggest that liquid occurs only 50 % of the time but with a high liquid percentage when it does occur, then it can create a discrepancy between the metric used here and other cloud-phase metrics like SLF. Even though the model would then have an sLCC frequency of occurrence of 100 %, the amount of supercooled liquid could still very well be underestimated. When combined with previous studies finding that ESMs generally underestimate SLF, our findings suggest that models too frequently simulate genuinely mixed-phase clouds with both phases present at the same time and point, while observations more frequently reveal conditionally mixed conditions in which ice and liquid are separated in time and/or space (Korolev and Milbrandt, 2022).

Previous studies not only use different metrics but also, in some cases, evaluate different generations of climate models, each with distinct cloud characteristics. For instance, many models included in CMIP5 were found to have too few sLCCs (Cesana et al., 2012, 2015; Tan and Storelvmo, 2016; Kay et al., 2016; Lenaerts et al., 2017; McIlhattan et al., 2017). This deficiency was addressed in several of the next-generation models used in CMIP6, resulting in an increase in the simulated SLF (Mülmenstädt et al., 2021). Although some CMIP6 models used in this study still over- or underestimate the traditional SLF metrics compared to CALIPSO observations, models such as CNRM-CM6-1 and CNRM-ESM2-1 are within 1 standard deviation of the satellite observations (Fig. F1). Studies focusing on high latitudes have shown that comparing individual CMIP5 models with CMIP6 models reveal improved cloud representation, more closely aligning with observations due to microphysical adjustments in the newer model versions (Lenaerts et al., 2020; McIlhattan et al., 2020). Indeed, McIlhattan et al. (2020) showed that when going from the CMIP5 to the CMIP6 version of CESM, the frequency of the LCC bias switched direction and overestimated in comparison to CloudSat–CALIPSO. Lenaerts et al. (2020) showed that cloud coverage increased and is slightly overestimated in the CMIP6 version of CESM but still underestimates the ice water path. Consequently, some of the differences between our findings and those from previous studies may be attributed to these improvements in the newer model generations.

5.1 Links between biases in f_sLCC and f_snow

We find that ERA5 and CMIP6 overestimate the magnitude of f_sLCC and f_snow in the mid-to-high latitudes of both hemispheres when compared to CloudSat–CALIPSO. In contrast, McIlhattan et al. (2017) showed that the CESM-LE underestimates the f_LCC by ∼17 % and overestimates the f_snow by ∼57 % in the Arctic. However, since we utilize a different metric (sLCC instead of LCC), there is no reason to expect the model biases to be identical. Nevertheless, in a further study by McIlhattan et al. (2020), it was shown that the LCC frequency in the newer CESM version is more aligned with the satellite observations, except for in the summer months, where it overestimates the LCC frequency.

Similarly, we generally observe that ERA5 and the CMIP6 mean overestimate the f_sLCC to various extents during all seasons. These overestimations are likely linked to the microphysical parameterizations of cloud processes that govern the cloud phase. This finding aligns with McIlhattan et al. (2020), indicating that while newer model versions have advanced in representing LCC frequencies more accurately, they can still overestimate cloud occurrences, particularly in specific seasons. Precipitation in the new CESM version is more frequent but lighter overall compared to the previous version, which is similar to our findings, indicating that the models's sLCCs produce continuous snowfall, analogous to the “perpetual drizzle” problem (Mülmenstädt et al., 2020; Lavers et al., 2021).

Furthermore, while McIlhattan et al. (2017, 2020) considered a single ESM, our study considers an ensemble of CMIP6 models. Nevertheless, the insights from McIlhattan et al. (2017, 2020) provide a relevant context for our finding that ERA5 and the CMIP6 model mean produce sLCCs more frequently than observed in the NH and SH mid-to-high latitudes, especially over the sea ice and land, depending on the season (Figs. C1 and C2). In these regions, not only is the frequency of occurrence of sLCCs too high, but the sLCCs are too efficient at producing snowfall in ERA5 and CMIP6 models (Figs. 5 and 6). The latter finding is consistent with the findings of McIlhattan et al. (2017), who showed that LCCs produce snow too frequently in the CESM-LE model.

McIlhattan et al. (2017) explained the overestimation of f_snow from LCC by exploring one potential microphysical pathway for removing supercooled cloud liquid – the WBF process. An overactive WBF process and the subsequent removal of cloud ice from the cloud through snow formation in ERA5 and CMIP6 models could also explain the large overestimation of f_sLCC and f_snow in ERA5 and the CMIP6 model mean. ERA5 and the CMIP6 models may thus have more frequent but lighter snowfall events that lead to a higher f_snow, while maintaining the clouds for longer compared to the CloudSat–CALIPSO observations. However, for some ESMs, there is no explicit simulation of the WBF process but rather a simple temperature-dependent partitioning of cloud condensate into liquid and ice. For such crude parameterization schemes, the identified biases are inevitable, as clouds at a given mixed-phase temperature would always have some liquid present (hence a very high f_sLCC) and at the same time always have ice available that could be converted to snow (hence a very high f_snow).

The difference in the simulated f_sLCC and f_snow compared to CloudSat–CALIPSO in the NH and SH could, in theory, result from differences in surface conditions and/or weather patterns and thereby differences in moisture availability for the formation of sLCCs (McIlhattan et al., 2017). However, the fact that very similar biases to those in the reanalysis are found in the CMIP6 models might suggest that surface conditions and weather patterns are not the main explanation for the identified biases. That being said, reanalysis products are known to be less reliable at high latitudes (Liu and Key, 2016). We therefore cannot rule out differences in weather patterns as a partial cause of the biases discussed above but maintain that deficiencies in the model representation of cloud microphysics is the more likely culprit. The inclusion of both ERA5 and CMIP6 models in this analysis is thus beneficial as it helps highlight potential root causes of the identified biases, although it should be noted that individual CMIP6 models do not all show biases that match ERA5 (Figs. 4 and 8). As stated above, we cannot entirely rule out a modest contribution to the biases from circulation differences, even for the reanalysis. For example, according to Boisvert et al. (2020), all reanalysis products show lower precipitation amounts in the southern Ross Sea and Weddell Sea embayments near the ice shelves. They explain that this is due to persistent cold and dry katabatic wind blowing from the continent across the ice shelves and out over the sea ice, which is not beneficial for precipitation formation and not captured well by the reanalyses.

Under the assumption that the overall weather patterns and surface conditions are accurate in ERA5, the most plausible root causes of the biases in f_sLCC and f_snow are linked to the microphysical parameterizations in the reanalysis and ESMs (Kiszler et al., 2024).

Boisvert et al. (2020) suggest that the difference in temperature and temperature threshold among various reanalysis microphysical products may explain the latitudinal inequality in the snowfall amount observed between CloudSat and the reanalysis products. But, here, no strong latitudinal dependence is observed in Fig. 6 as we investigate f_snow instead of the snowfall amount. The same principle should apply to ERA5 and the CMIP6 models as they rely on microphysical parameterizations. For the most part, the 0 °C isotherm shown in Figs. 1 and 2 supports our argument that the primary issue with the ERA5 and CMIP6 datasets lies not with the simulated temperature itself but with the representation of cloud properties and microphysics. This distinction highlights that the observed deviations in f_sLCC and f_snow are driven more by inaccuracies in cloud simulation than by temperature discrepancies. Exceptions occur over central Europe during DJF between ECMWF-AUX, ERA5, and the CMIP6 ensemble mean (Fig. 1a, e, i). However, more notable is the difference in the CMIP6 ensemble mean over the central Arctic during summer (Fig. 1k), where simulated temperatures appear to be too cold. In this specific case, the cloud bias could stem from a temperature bias, suggesting a potential link between temperature inaccuracies and cloud simulation in the CMIP6 ensemble mean for this region and season. Figures 4 and 8 show some models performing better than others for f_sLCC and f_snow in comparison to CloudSat–CALIPSO. As discussed above, many models have a simple temperature-dependent cloud phase that would almost certainly cause them to overestimate the f_sLCC and f_snow. A follow-up study will investigate the SLF by isotherms for mixed-phase clouds of several CMIP6 models. Figure F1 shows that some of the CMIP6 models have a more variable SLF in time and space (e.g., AWI-ESM-1-1-LR and MPI-ESM1-2-LR), which might be an indication that they could have more sophisticated microphysics and could, therefore, perform better in terms of sLCC and surface snowfall occurrence. However, it is not immediately clear which models use a purely temperature-dependent cloud-phase partitioning. From the interquartile range (IQR; Fig. F2) CNRM-CM6-1, CNRM-ESM2-1, and IPSL-CM6A-LR are assumed to have a simple temperature-dependent cloud phase due to the smaller IQR of the temperature at various SLF. It is interesting to note that these models with the simpler temperature-dependent schemes are the most poorly performing models for f_snow (Fig. 8).

Imura and Michibata (2022) studied the effect of changing the microphysical parameterization scheme from the traditional diagnostic scheme used in MIROC6, and participating in CMIP6, to a prognostic scheme. They found that the prognostic scheme improved cloud coverage and snowfall in the Arctic. In our study, MIROC6 is one of the models with the largest overestimation of f_sLCC in the SH and the model with the smallest difference in boreal spring (Fig. 4). Furthermore, MIROC6, with the diagnostic scheme used here, has a large negative difference in f_snow (45 %–60 %) that is independent of hemisphere and season (Fig. 8). However, the MIROC6 with the prognostic parameterization scheme in Imura and Michibata (2022) produced light snowfall too frequently, which is in agreement with the results presented here, despite the use of a different microphysical scheme.

ERA5 and the CMIP6 model mean are not able to reproduce f_snow in comparison to CloudSat–CALIPSO. The overestimation in ERA5 could be related to the global average wet bias of up to 0.27 mm d⁻¹ in the ECMWF product used in ERA5 (Lavers et al., 2021). If the modeled sLCCs are practically always snowing, then the sLCCs should have a shorter lifetime and subsequently have a lower f_sLCC than observations suggest. However, ERA5 and the CMIP6 model mean seem to have sLCCs that produce a little bit of snow all the time and, at the same time, maintain the sLCCs in the atmosphere. This bias may be the sLCC's counterpart to the “perpetual drizzle problem” that has been identified in ESMs for warm liquid clouds (Mülmenstädt et al., 2020).

5.2 Implications for modeling and future projections

An overestimation of the snowfall frequency can have significant effects on the Earth's radiative budget. The fact that the models seem to snow a little bit all the time instead of producing occasional heavy snowfall events, as well as periods with no snowfall, is therefore concerning. One implication is that fresh snow will be deposited too frequently on top of the snowpack, ice sheets, sea ice, and land surfaces, with consequences for the surface albedo. With time, the aging snow surfaces darken due to snow/ice metamorphism processes and the deposition of absorbing aerosols (e.g., Picard et al., 2012; Carlsen et al., 2017). Continuous light snowfall in the model could thus lead to a simulated overestimation of the surface brightness, which in turn could limit sea ice melt in ESMs. Indeed, only ESMs that simulate excessive Arctic warming are able to reproduce the observed sea ice loss in recent decades (Notz and Stroeve, 2016).

The overestimations of f_sLCC and f_snow could stem from how models handle cloud microphysical processes, particularly the WBF process and the partitioning of cloud condensates into the liquid and ice phases. Some identified biases in ESMs may be due to the specific parameterization schemes used in the models, which could rely on simple temperature-dependent cloud-phase partitioning. For example, models with a temperature-dependent cloud phase could simulate a fixed liquid fraction at a given temperature within the mixed-phase range, making it impossible for them to simulate an all-ice or all-liquid cloud in that range. Consequently, these models could tend to under- or overestimate the occurrence of sLCCs.

Interestingly, models with more sophisticated microphysics schemes (Figs. F1 and F2) do not necessarily perform better (Figs. 4 and 8). Larger biases in cloud phase can still occur, for example, if these models do not accurately represent INPs. This underscores how generational advancements and implementations of observational constraints can significantly influence the comparison and interpretation of cloud and precipitation patterns.

Continuing work should be performed on the comparison between ESMs, reanalysis, and observations at specific locations to better understand microphysical processes and their representation in ESMs and reanalysis data. Cloud microphysical scheme development should, therefore, focus on improving these processes in order to reduce cloud-phase and surface snowfall biases. More sophisticated schemes should be explored, and ensemble studies with varied microphysical model parameterizations can provide insights into how different parameterizations affect the representation of sLCCs and snowfall in climate simulations.

Ground and aircraft observations in connection with field campaigns at specific locations should, therefore, be increased, especially to improve our understanding of cloud phase, liquid and ice water content, and cloud-top temperatures in regions where sLCCs form. These field campaigns should also focus on the processes responsible for snowfall from these clouds.

Finally, as shown in the present study, satellite observations are an essential validation tool for assessing the representations of clouds and precipitation in ESMs. To ensure that future model generations are as accurate as possible, continuous improvements to satellite retrievals are required to further refine our understanding of cloud and snowfall properties, especially in high-latitude regions where limitations in the Sun angle and ground track are prevalent. The accurate representation of these cloud and snowfall processes in models, validated through enhanced satellite observations, is essential for reliable simulations of weather and climate. This accuracy is particularly important for making future predictions of the hydrological cycle, which are essential for understanding and mitigating the impacts of climate change.

Appendix A: CMIP6 models and LWP calculation

The daily averages of the total mass of liquid water (LWP) in a grid box are calculated using the hydrostatic equation in the CMIP6 models, apart from the models from the modeling institute, MOHC (the UK Meteorological Office). First, isobaric pressure levels are retrieved with the individual CMIP6 hybrid sigma model formula found in Table A1. The terms hyam and hybm represent the vertical coordinate formula term along the dimension Nlev. p_sfc(i,j) is the surface air pressure for each latitude (i) and longitude (j) grid box, and p₀(Nlev) is the vertical coordinate formula term for the reference pressure (both in Pa) for the individual model. For CMIP6 models with orographic vertical coordinates (UKESM1-0-LL and HadGEM3-GC31-LL; Table A1) representing bm(Nlev), the vertical coordinate formula term and orog(i,j) the surface altitude are given in units of meters.

After retrieving the pressure coordinate on half-isobaric-pressure levels ($p(i,j,\mathrm{Nlev}-(k{\pm }^{\mathrm{1}}{/}_{\mathrm{2}}\left)\right)$) and $\mathrm{clw}(i,j,\mathrm{Nlev})$ on the full-isobaric-pressure levels (Nlev) in CMIP6 the $\mathrm{lwp}(i,j,\mathrm{Nlev})$ in each vertical grid box can be calculated with the hydrostatic equation.

$${\displaystyle \frac{\mathrm{\Delta }p}{\mathrm{\Delta }Z}}=-{\mathit{\rho }}_{\mathrm{air}}(i,j,\mathrm{Nlev})\cdot g$$

${\mathit{\rho }}_{\mathrm{air}}(i,j,\mathrm{Nlev})$ is the density of the air mass at full-isobaric-pressure level. ΔZ depicts the height difference between the half-isobaric- pressure levels, and g=9.81 m s⁻² is the gravity acceleration. After the American Meteorological Society (AMS) Glossary of Meteorology, the $\mathrm{lwp}(i,j,\mathrm{Nlev})$ is defined as

$$\begin{array}{ll}{\displaystyle }\mathrm{lwp}(i,j,\mathrm{Nlev})& {\displaystyle }=\underset{z=\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathit{\rho }}_{\mathrm{air}}(i,j,\mathrm{Nlev})\cdot \mathrm{clw}(i,j,\mathrm{Nlev})\cdot \mathrm{d}z\\ {\displaystyle }& {\displaystyle }=\underset{\mathrm{0}}{\overset{p=p_{\mathrm{0}}}{\int }}{\mathit{\rho }}_{\mathrm{air}}(i,j,\mathrm{Nlev})\cdot \mathrm{clw}(i,j,\mathrm{Nlev})\\ {\displaystyle }& {\displaystyle }\cdot \left(-{\displaystyle \frac{\mathrm{d}p}{{\mathit{\rho }}_{\mathrm{air}}\cdot g}}\right)\\ {\displaystyle }\mathrm{lwp}(i,j,\mathrm{Nlev})& {\displaystyle }=\underset{p=p\left(i,j,\mathrm{Nlev}-\left(k{-}^{\mathrm{1}}{/}_{\mathrm{2}}\right)\right)}{\overset{p=p\left(i,j,\mathrm{Nlev}-(k+{}^{\mathrm{1}}{/}_{\mathrm{2}})\right)}{\int }}{\displaystyle \frac{\mathrm{clw}(i,j,\mathrm{Nlev})}{g}}\cdot \mathrm{d}p.\end{array}$$

The LWP(i,j) can then be calculated by summing the individual $\mathrm{lwp}(i,j,\mathrm{Nlev})$ per vertical grid box from the surface to the top of the atmosphere. And the following applies for the daily mean total liquid water path per column air in each pixel:

$$\begin{array}{}\text{(A1)}& {\displaystyle }\mathrm{LWP}(i,j)=\sum _{k=\mathrm{0}}^{\mathrm{NLEV}+\mathrm{1}}\mathrm{lwp}(i,j,\mathrm{Nlev}),\text{(A2)}& {\displaystyle }\begin{array}{rl}& =-{\displaystyle \frac{\mathrm{1}}{g}}\sum _{k=\mathrm{0}}^{\mathrm{NLEV}+\mathrm{1}}\mathrm{clw}(i,j,\mathrm{Nlev})\\ & \cdot \left[p\left(i,j,\mathrm{Nlev}-(k-{}^{\mathrm{1}}{/}_{\mathrm{2}})\right)\right.\\ & \left.-p\left(i,j,\mathrm{Nlev}-(k+{}^{\mathrm{1}}{/}_{\mathrm{2}})\right)\right].\end{array}\end{array}$$

Tatebe et al. (2019)Swart et al. (2019a)Ackermann et al. (2020)Mauritsen et al. (2019)Sellar et al. (2019)Roberts et al. (2019)Voldoire et al. (2019)Séférian et al. (2019)Boucher et al. (2020a)Sepulchre et al. (2020)

Table A1 This study uses the Coupled Model Intercomparison Project Phase 6 (CMIP6) models. Models marked with an asterisk (*) have some “historical” simulations only until 2009. Isobaric levels are calculated with a formula given by the individual model. The three equations used are (⊗) $p(i,j,\mathrm{Nlev})=\mathrm{hyam}\left(\mathrm{Nlev}\right)\cdot p_{\mathrm{0}}+\mathrm{hybm}\left(\mathrm{Nlev}\right)\cdot p_{\mathrm{sfc}}(i,j)$, (◦) $p(i,j,\mathrm{Nlev})=\mathrm{hyam}\left(\mathrm{Nlev}\right)+\mathrm{hybm}\left(\mathrm{Nlev}\right)\cdot p_{\mathrm{sfc}}(i,j)$, and (△) $z(i,j,\mathrm{Nlev})=\mathrm{am}\left(\mathrm{Nlev}\right)+\mathrm{bm}\left(\mathrm{Nlev}\right)\cdot \mathrm{orog}(i,j)$. Further description of the ESMs is available on the ESDoc interface (https://explore.es-doc.org/cmip6/models/, last accessed 25 October 2023).

Appendix B: Variable unit transformation

To compare ERA5 and CMIP6 snowfall ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$; Table 1) to CloudSat–CALIPSO (sf_CC), we apply the following multiplication to achieve the satellite surface snowfall rate (in mm h⁻¹; Table 1). The density of water (ρ_water) is about 1000 kg m⁻³ and 1000 mm are 1 m.

$$\begin{array}{}\text{(B1)}& {\mathrm{sf}}_{\mathrm{CC}}={\displaystyle \frac{\mathrm{mm}}{\mathrm{h}}}\cdot {\displaystyle \frac{\mathrm{1}\mathrm{m}}{\mathrm{1000}\mathrm{mm}}}\cdot {\mathit{\rho }}_{\mathrm{water}}={\displaystyle \frac{\mathrm{m}}{\mathrm{h}}}\cdot {\displaystyle \frac{\mathrm{1000}\mathrm{kg}}{\mathrm{1000}{\mathrm{m}}^{\mathrm{3}}}}={\displaystyle \frac{\mathrm{kg}}{{\mathrm{m}}^{\mathrm{2}}\mathrm{h}}}\end{array}$$

In the ERA5 reanalysis data, the daily mean of msr has units of $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (Table 1). To make ERA5 snowfall daily means comparable to CloudSat and CMIP6 snowfall rate, we apply the following multiplication to achieve the mean snowfall rate per hour (msr), where we know that 1 h has 3600 s.

$$\begin{array}{}\text{(B2)}& {\displaystyle }{\mathrm{sf}}_{\mathrm{ERA}}{\displaystyle }={\displaystyle \frac{\mathrm{kg}}{{\mathrm{m}}^{\mathrm{2}}\mathrm{s}}}\cdot {\displaystyle \frac{\mathrm{3600}\mathrm{s}}{\mathrm{h}}}={\displaystyle \frac{\mathrm{kg}}{{\mathrm{m}}^{\mathrm{2}}\mathrm{h}}}\cdot \mathrm{3600}\end{array}$$

The same is done for the CMIP6 snowfall flux parameter (prsn; Table 1).

$$\begin{array}{}\text{(B3)}& {\displaystyle }{\mathrm{sf}}_{\mathrm{CMIP}\mathrm{6}}{\displaystyle }={\displaystyle \frac{\mathrm{kg}}{{\mathrm{m}}^{\mathrm{2}}\mathrm{s}}}\cdot {\displaystyle \frac{\mathrm{3600}\mathrm{s}}{}}\mathrm{h}={\displaystyle \frac{\mathrm{kg}}{{\mathrm{m}}^{\mathrm{2}}\mathrm{h}}}\cdot \mathrm{3600}\end{array}$$

Appendix C: Spatial distribution of f_sLCC and f_snow

Figure C1 Seasonal averages of f_sLCC in the NH mid-to-high latitudes between 2007 and 2010. Combined CloudSat and CALIPSO observations are shown in panels (a)–(d). The last two rows are the difference plot. They are CloudSat–CALIPSO (CC) observations minus ERA5 (e–h) or the CMIP6 model mean (i–l), where valid data occur, with green (pink) values showing underestimation (overestimation) in ERA5 and the CMIP6 model mean concerning the satellite observations. Areas where the difference between CloudSat–CALIPSO and CMIP6 model mean is not significant (< 95 %) are marked with hatches. The area-weighted averages for the study area where CloudSat–CALIPSO have observations are displayed in the lower-left corner of each map. The red line(a–d) shows the average sea ice edge of the 20 % sea ice concentration (SIC) between 2007 and 2010 for the given season.

Figure C2 Seasonal averages of f_sLCC in the SH mid-to-high latitudes. Layout and differences are identical to Fig. C1.

Figure C3 The figure presents the seasonal average of f_snow in the NH. The layout and presentation resemble Fig. C1.

Figure C4 Seasonal averages of f_snow in the SH. The layout and presentation resemble Fig. C1.

Appendix D: Zonal mean of 2 m temperature

Figure D1 Zonal mean 2 m temperature of ECMWF-AUX (black) and ERA5 (dashed orange) for mid-to-high latitudes.

Appendix E: LWP threshold sensitivity – SH

Figure E1 Annual average difference plots (CloudSat–CALIPSO (CC) minus model) for f_sLCC in the SH mid-to-high latitudes. The first row represents ERA5 data (a–d), and the second row shows the CMIP6 model mean (e–h) difference. LWP thresholds of 3 g m⁻² (a, e), 5 g m⁻² (b, f), 10 g m⁻² (c, g), and 15 g m⁻² (d, h) are applied. The analysis in this study focuses on the LWP threshold of 5 g m⁻². The maps in each row are accompanied by area-weighted averages for the study area that are located in the lower-left corner of each map. The layout and area-weighted averages are calculated in the same way as those in Fig. 10.

Appendix F: Supercooled liquid fraction in CMIP6 models

Figure F1 Supercooled liquid fraction (SLF) as a function of temperature for CALIPSO (black) and the CMIP6 models used in this study (color). Error bars on the CALIPSO SLF values correspond to 1 standard deviation. All values represent an area-weighted average for lat≥45° in the Northern Hemisphere (NH; solid lines) and Southern Hemisphere (SH; dashed lines).

Figure F2 Interquartile range of the temperature at given supercooled liquid fractions (e.g., SLF 50 %; representing 50 % ice and 50 % liquid) for the different models used in this study for lat≥45° over the Northern Hemisphere (NH; a) and Southern Hemisphere (SH; b).