Evaluation of Global Climate and Storm-Resolving Model Representations of Mixed-Phase Clouds and Their Hemispheric Contrasts

1. Introduction

Mixed-phase clouds, which contain both liquid droplets and ice crystals, exert a strong influence on Earth’s radiation budget. Their radiative effect depends critically on the partitioning between liquid and ice: for the same water content, liquid clouds are optically thicker and reflect more solar radiation than ice clouds. However, the physical processes governing liquid–ice partitioning remain incompletely understood, and their representation in weather and climate models is still highly uncertain [1].

Uncertainty arises from the interplay of several factors, including temperature, vertical motion, the availability of aerosols acting as cloud condensation nuclei or ice-nucleating particles, and microphysical processes such as nucleation, riming, and secondary ice production [2,3]. These processes determine the supercooled liquid fraction (SLF), defined as the ratio of supercooled liquid water to the total condensate. Because SLF strongly affects cloud radiative forcing, evaluating how well models reproduce SLF is essential for reducing uncertainties in climate projections.

Satellite observations have provided valuable insights into mixed-phase clouds. Active sensors such as CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation) have been widely used to assess cloud phase occurrence and phase partitioning [4,5]. Passive sensors, despite their limitations, offer broader spatial and temporal coverage. Using observations from both active and passive satellite sensors, Bruno et al. (2021) [6] reported systematic hemispheric contrasts in SLF: values were generally larger in the Southern than in the Northern Hemisphere, except for continental low-level clouds where the opposite pattern was found. Whether global and storm-resolving models can reproduce these contrasts has not yet been fully assessed.

To evaluate model performance against observations, we adopt the framework introduced by Cesana et al. (2015) [7], who used joint histograms of SLF as a function of cloud top temperature (CTT). This approach provides a consistent way to compare satellite retrievals and model outputs despite differences in native variables and resolution. In this work, we extend this method by introducing joint histograms of SLF and averaged liquid effective radius ($\overline{r_{liq}}$). These histograms offer additional insight into microphysical processes by linking phase partitioning with droplet size.

Building on Bruno et al. (2021) [6], the present study expands the analysis in two ways. First, we evaluate not only a global climate model—the Community Atmosphere Model version 6, Oslo configuration (CAM6-Oslo)—but also three storm-resolving models: the ICOsahedral Non-hydrostatic model (ICON), the Simple Cloud-Resolving Energy Exascale Earth System Model Atmosphere Model (SCREAM), and the Goddard Earth Observing System model (GEOS). Second, we propose the SLF–$\overline{r_{liq}}$ diagnostic as a new benchmark for assessing cloud microphysics.

The objectives of this study are therefore:

• To provide updated reference values for mixed-phase cloud occurrence based on satellite observations and models;
• To assess the ability of global and storm-resolving models to reproduce hemispheric contrasts of mixed-phase clouds;
• To explore the potential of SLF–$\overline{r_{liq}}$ relationships as a benchmark for evaluating model microphysics schemes.

2. Datasets

In this paper, datasets derived from observations, a global climate model, and three global storm-resolving models were analysed.

The observations included two passive satellite sensor-based datasets: Cloud_cci AVHRR-PMv3 (Cloud Climate Change Initiative Advanced Very High Resolution Radiometer, [8,9]), with "PM" indicating the equatorial crossing time in the afternoon; and

MODIS_D3_MODIS (Moderate Resolution Imaging Spectroradiometer (Aqua/Terra) Cloud Properties Level 3 daily [10,11]).

In Bruno et al. (2021) [6], cloud property datasets derived from the passive sensor AVHRR onboard polar-orbiting NOAA satellites were compared to the active sensor CALIOP (Cloud-Aerosol Lidar with Orthogonal Polarization) onboard the polar-orbiting CALIPSO satellite (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation). The AVHRR-based datasets showed a warm bias in cloud top temperature compared to CALIOP; furthermore, a thermodynamic phase mismatch was found with $\sim 9\%$ of liquid cloud tops in CALIOP retrieved as ice in AVHRR-based datasets. Nevertheless, AVHRR-based datasets were shown to be able to reproduce the hemispheric contrasts from the active satellite sensor CALIOP. Because of their agreement and because active sensors provide less spatial coverage than passive sensors for the same time frame, Cloud_cci v3 — the most recent AVHRR-based dataset analysed in Bruno et al. (2021) [6] — has been included in this paper (hereafter AVHRR).

MCD06COSP_D3_MODIS [10,11] is a product specifically designed to facilitate the comparison of MODIS observations and climate model output produced with the Cloud Feedback Model Intercomparison Project (CFMIP) Observation Simulator Package 2 (COSP2) simulator [12], which is included in the global climate model (GCM) dataset used in this study. The MCD06COSP_D3_MODIS dataset provides MODIS measurements from the polar-orbiting Aqua and Terra satellites, aggregated to a coarser spatial resolution ($1^{\circ }$ x $1^{\circ }$) during dataset processing. The cloud optical properties are retrieved by using only one channel (3.7µm) instead of six. For brevity and readability, we hereafter refer to this dataset as MODIS.

The GCM dataset we analysed here is the atmospheric component of the Norwegian Earth System Model version 2 (NorESM2), which is based on the Community Atmosphere Model Version 6 (CAM6, [13]) and is hereafter referred to as CAM6-Oslo. CAM6-Oslo uses the OsloAero5.3 [14] aerosol scheme and aerosol-radiation-cloud interactions from [15], and is the only model in this work to treat aerosols interactively as possible cloud condensation nuclei and ice nucleating particles.We specifically used MODIS-like cloud fields from CAM6-Oslo produced by the COSP2 satellite simulator to provide the best comparison with passive satellite retrievals. This output is still referred to as CAM6-Oslo in this paper.

The three global storm-resolving model datasets we analysed are part of the DYnamics of the Atmospheric general circulation Modeled On Non-hydrostatic Domains (DYAMOND) Winter project [16]: ICON (ICOsahedral Nonhydrostatic Model, [17,18]), SCREAM (Simple Cloud-Resolving E3SM (Energy Exascale Earth System Model) Atmosphere Model, [19]), and GEOS (Goddard Earth Observing System, [20]). For consistency, the atmosphere-only experiments of all models were used in this study. As the study includes an analysis of mixed-phase clouds and the passive satellite sensors provide only cloud-top information, only the cloud top information of model outputs was used, therefore reproducing the satellite view of clouds. The following cloud parameters were considered in this study: cloud-top phase (CPH) considered either as a liquid or ice pixel (for AVHRR), or as liquid (CLW) and ice (CLI) cloud fraction (for MODIS and CAM6-Oslo), or as mass mixing ratio of ice ($q_i$) and liquid water ($q_l$, for the DYAMOND datasets); cloud-top pressure (CTP); cloud-top temperature (CTT); cloud optical thickness (COT); and liquid cloud effective radius ($r_{liq}$) at the cloud top. Not all datasets provide all variables. These are listed in Table ?? for each dataset together with the most relevant features for our study.

3. Methods

All datasets provide global data with daily or higher time resolution, here analysed near-globally (${60}^{\circ }$N to ${60}^{\circ }$S), thus excluding satellite data with low confidence due to the presence of sea ice [21] and high solar zenith angle [22]. In addition to the near-global analysis, continental and marine regions were considered separately, as were different cloud types and latitudes: Northern Hemisphere mid-latitudes (NH) from ${60}^{\circ }$ N to ${30}^{\circ }$ N; Southern Hemisphere mid-latitudes (SH) from ${30}^{\circ }$ S to ${60}^{\circ }$ S. We analysed 4 years (from June 1, 2009 to May 31, 2013, see [6]) of AVHRR, MODIS, and CAM6-Oslo. The DYAMOND datasets, however, which are computationally very costly, are only available for a 30-day analysis period (January 30 - February 28, 2020) following a 10-day spin-up period. In this study, only the February portion is used.

3.1. Cloud Identification

AVHRR, MODIS and the COSP outputs from CAM6-Oslo provided the cloud top information in 2-dimensional composites as a function of latitude and longitude, while the DYAMOND datasets provided the 3-dimensional view of the atmosphere. For this reason, a top-down approach was used to locate the cloud top in the DYAMOND models. In ICON and GEOS, the cloud top was defined as the highest atmospheric layer for which $q_w+q_i>{10}^{-8}$ kg ${\mathrm{kg}}^{-1}$, following [23]; to make data from SCREAM consistent with ICON and GEOS, the cloud top was defined as the highest atmospheric layer for which COT $>0$.

We also wish to account for the diurnal and COT-dependent sampling characteristics of AVHRR and MODIS when producing comparable model fields with CAM6-Oslo. Specifically, the analysis of AVHRR and MODIS was limited to daytime measurements, as the COT involved in the cloud-type classification could only be detected by satellite sensors in the visible range. To remove very thin clouds for which retrievals are more uncertain, only satellite measurements with COT > 0.3 and effective radius > 4 µm were used for this study [24]. For consistency, the same filters were applied to CAM6-Oslo, with the cloud top being defined as the highest cloud layer below an overlying COT of 0.3 for the COSP outputs.

3.2. SLF Derivation

The SLF was computed at the top of mixed-phase clouds. As in AVHRR the thermodynamic phase of a single cloudy pixel can be either liquid or ice, a reduction in spatial resolution was required to compute the SLF within gridboxes: the original latitude-longitude daily composites consisting of 3000x7200 pixels ($0.{05}^{\circ }$ spatial resolution) were reorganised into 500x1200 gridboxes of 6x6 pixels ($0.3^{\circ }$ spatial resolution), with each gridbox considered cloudy only if it contained at least two cloudy pixels (i.e., more than $\sim 5.56\%$ of the gridbox contained clouds). Furthermore, the gridboxes were only included in the analysis if the largest difference between the coldest and the warmest cloud top temperature of all valid pixels was less than ${20}^{\circ }$C, filtering out gridboxes including clouds of very different heights (such as low- and high-level clouds, [6]). Finally, gridboxes were termed "mixed-phase" if they included coexisting liquid and ice pixels with a cloud top temperature between $-{40}^{\circ }$C and $0^{\circ }$C. For each gridbox the SLF was then calculated using the number of liquid and ice pixels in the gridbox, following the Equation (1):

$$SL{F}_{n^{\circ }pixel}={\displaystyle \frac{{\mathrm{n}}^{\circ }\mathrm{liquid}\mathrm{pixels}}{{\mathrm{n}}^{\circ }\mathrm{liquid}\mathrm{pixels}+{\mathrm{n}}^{\circ }\mathrm{ice}\mathrm{pixels}}}$$

(1)

The model datasets and MODIS provide information on both liquid and ice content or cloud fraction in each pixel. For this reason, no gridboxes were needed to compute the SLF, as the individual pixels could already contain mixed phase. Therefore, the SLF was calculated as the ratio of the liquid cloud fraction and the sum of the liquid and ice cloud fractions (for MODIS and CAM6-Oslo, Equation (2)) or as the ratio of the mass mixing ratio of liquid water and the sum of the mass mixing ratios of ice and liquid water (for ICON, SCREAM, and GEOS, Equation (3)).

$$SL{F}_{cloudfraction}={\displaystyle \frac{\mathrm{liquid}\mathrm{cloud}\mathrm{fraction}}{\mathrm{liquid}\mathrm{cloud}\mathrm{fraction}+\mathrm{ice}\mathrm{cloud}\mathrm{fraction}}}$$

(2)

$$SL{F}_{massfraction}={\displaystyle \frac{{\mathrm{q}}_{\mathrm{w}}^{\phantom{0}}}{{\mathrm{q}}_{\mathrm{w}}^{\phantom{0}}+{\mathrm{q}}_{\mathrm{i}}^{\phantom{0}}}}$$

(3)

Hoose et al. (2018) [25], used a similar approach to generate cloud top liquid fraction from models and found that $SL{F}_{n^{\circ }pixel}$ was $1^{\circ }$C to $2^{\circ }$C warmer than $SL{F}_{massfraction}$ at SLF ∼ 0.5. This result suggests that these metrics are generally comparable, though $SL{F}_{massfraction}$ may be slightly higher than $SL{F}_{n^{\circ }pixel}$.

ICON, SCREAM, and GEOS were analysed after coarse graining to $0.5^{\circ }$ x $0.5^{\circ }$ by averaging for better comparison with AVHRR.

3.3. Classification of Cloud Regimes and k-Means Clustering

The clustering of clouds into groups called regimes is a widely used method over the last few decades. In some works the classification is based on vertical velocity [26,27], on sea level pressure [28], or on sea surface temperature and large-scale vertical motion [29,30,31], in some others on CTP, COT, and total cloud cover [32,33,34,35,36,37,38]. Our cloud-type classification was based on standard ISCCP (International Satellite Cloud Climatology Project) cloud top pressure (CTP$=[0,440,680,1000]$ hPa) and cloud optical thickness (COT$=[0,3.6,23,379]$) thresholds, as in Bruno et al. (2021) [6]. Each COT-CTP combination was identified by the cloud name (e.g, cumulus, altocumulus, etc., [39]). As a classification of cloud types based solely on COT and CTP has limitations [40], the cloud names are mainly used as a reference to compare clouds within the same COT-CTP ranges. ICON and GEOS do not provide the COT, therefore their clouds were sorted by using CTP thresholds of $[0,440,680,1000]$ hPa to classify them into high-, mid-, and low-level clouds. SCREAM provides the COT for each vertical level, therefore, the column integral of the COT for each location was considered.

As the spatial resolution of AVHRR and DYAMOND datasets has been coarsened (Section 3.2), pixel-based COT and CTP information in the original resolution has been collected in COT-CTP joint histograms. While the clouds from the MODIS and CAM6-Oslo could be directly classified on a pixel basis by using COT-CTP joint histograms, a different approach was required for AVHRR and DYAMOND datasets as COT and CTP were already collected in joint histograms. Therefore, k-means clustering was applied to the joint histograms in order to classify cloud regimes in the AVHRR and DYAMOND datasets.

The k-means algorithm [41] is a method used to sort data (points or vectors) into clusters resulting from the minimum sum of squared Euclidean distances between each cluster and its nearest cluster mean called centroid. This algorithm has limitations occurring, for example, when the distance between clusters is too large and therefore the centroids cannot move between clusters, or when stable clusters prevent the movement of centroids (e.g., high-density data regions may accumulate more centroids than necessary, with the risk of losing less dense clusters). Another difficulty with cluster analysis is finding the best number of clusters in which the data can be considered well represented. These limitations can be overcome by a better initialisation method and the iteration of the k-means with different initial solutions [42].

Initialising 42 centroids for AVHRR covered all possible COT-CTP combinations (excluding those for COT $<0.3$ as one of the applied filters for satellite observations and simulator, the possible COT-CTP combinations were 6 times 7, from COT$=[0.3,1.3,3.6,9.4,23,60,379]$ and CTP$=[0,180,310,440,560,680,800,1000]$ hPa bins, respectively) and made the k-means algorithm deterministic. In this way, the risk of neglecting some cloud types was prevented, clusters with low density were taken into account unless they were empty, and there was no need to repeat the k-means clustering since the initial centroids would be the same, and hence the final clusters too. Similarly to AVHRR, 7 prescribed initial centroids were used to cluster CTP histograms for ICON and GEOS (7 as the number of CTP bins), while 49 initial clusters (7 COT times 7 CTP combinations) were used for SCREAM, as no filters on COT has been applied to this dataset (see Section 3.1).

The final clusters needed to be merged to obtain the final 9 cloud types [39] or to organise them into 3 cloud heights. For this, the most frequent COT-CTP combination in each of the resulting clusters was used to determine to which cloud type the clusters belonged to. After this step, each of the 9 final COT-CTP joint histograms and 3 CTP histograms was individually normalised and represented a cloud type.

3.4. Geographical Distributions and Joint Histograms

Once the cloud classification was complete and the SLF computed, we derived the geographical distribution, and both the SLF-CTT and SLF-$\overline{r_{liq}}$ joint histograms for each cloud type.

We defined the relative frequency of occurrence of all clouds ($H_{all}$, Equation (4)) and mixed-phase clouds ($H_{MP}$, Equation (5)) as percentage:

$$H_{all}={\displaystyle \frac{{\mathrm{n}}^{\circ }\mathrm{cloudy}\mathrm{pixels}}{{\mathrm{n}}^{\circ }[\mathrm{cloudy}+\mathrm{invalid}\mathrm{cloudy}+\mathrm{clear}\mathrm{sky}]\mathrm{pixels}}}\times 100$$

(4)

$$H_{MP}={\displaystyle \frac{{\mathrm{n}}^{\circ }\mathrm{mixed}\text{-}\mathrm{phase}\mathrm{cloudy}\mathrm{pixels}}{{\mathrm{n}}^{\circ }[\mathrm{cloudy}+\mathrm{invalid}\mathrm{cloudy}+\mathrm{clear}\mathrm{sky}]\mathrm{pixels}}}\times 100$$

(5)

where the pixels are replaced by gridboxes in AVHRR.

We also defined the geographical relative frequency of occurrence of a cloud category $H_{ge{o}_{i,j}}$, normalised by the number of total available observational or model output time steps (including cloud-free pixels) $N_{to{t}_{i,j}}$, for each latitude-longitude position and expressed it in percentage:

$$H_{ge{o}_{i,j}}={\displaystyle \frac{N_{i,j}\left(\mathrm{cloud}\mathrm{category}\right)}{N_{to{t}_{i,j}}}}\times 100,$$

(6)

where i and j represented the position in latitude and longitude of the geographical distribution, while the cloud category is the combination of height (high, middle, or low clouds) and thermodynamic phase (all or mixed-phase clouds).

In the SLF-CTT joint histograms, the SLF changing with CTT is shown for all cloud types, with an SLF bin width of 0.05 and a CTT bin width of $1^{\circ }$C. For MODIS, CTT was replaced by CTP ranging from 0 hPa to 1000 hPa with a bin width of 25 hPa. We normalised these histograms, for each cloud type, with their own total frequency of occurrence ($H_{k,l}$, Equation 7, with $N_{k,l}$ being the absolute frequency of occurrence at bin combination (k,l), k referring to CTT or CTP, and l referring to SLF). The SLF-$\overline{r_{liq}}$ joint histograms show how $\overline{r_{liq}}$ varies with SLF in all cloud types, which we could use to interpret cloud dynamics and microphysics. We normalised these histograms individually for each cloud type (Equation 7, this time with k referring to $\overline{r_{liq}}$).

$$H_{k,l}={\displaystyle \frac{N_{k,l}}{{\sum }_{k,l}{N}_{k,l}}}$$

(7)

For consistency to Bruno et al. (2021) [6], we additionally normalised the histograms using an approach similar to CFADs (Contoured Frequency by Altitude Diagrams, $H_{CFADk,l}$, Equation (8)). CFAD histograms were introduced by Yuter and Houze (1995) [43] and later used in other forms for other studies (e.g., [44,45,46,47]). To highlight hemispheric contrasts in SLF over land and ocean, we computed the difference of SLF-CTT histograms using the CFAD-like histograms.

$$H_{CFA{D}_{k,l}}={\displaystyle \frac{N_{k,l}}{{\sum }_k{N}_{k,l}}}$$

(8)

A significance test based on the p-value had been applied to these histograms, showing dots where the bins showed a statistically significant difference between NH and SH. The p-value had been computed in the highest resolution possible, i.e., for each bin and its eight surrounding bins. For this reason, no significance test could be applied to the bins at the edges of the histograms. The threshold value for significance can be arbitrary; in this work, a threshold value of 0.05 was considered [48].

4. Results

4.1. Near-Global Cloud Cover and Mixed-Phase Occurrence

The relative frequencies of occurrence of all ($H_{all}$) and mixed-phase ($H_{MP}$) clouds retrieved from the datasets between ${60}^{\circ }$ North and ${60}^{\circ }$ South are summarised in Table 2.

The use of different sensors and models with individual characteristics, the application of different filters, the use of only daytime data (providing COT) in AVHRR, MODIS, and CAM6-Oslo, and the use of customised equations can lead to differences in the frequency of occurrence of clouds, so that: $H_{all}$ is more consistent between the datasets than $H_{MP}$, partly due to resolution and the definition of "mixed-phase"; $H_{all}$ of AVHRR decreases from 59.2% to 52.9% when the dataset is coarse-grained from ${0.05}^{\circ }$ to $0.3^{\circ }$, as a cloudy gridbox must contain at least two cloudy pixels; ICON shows the closest $H_{all}$ (65.5%) and $H_{MP}$ (6.9%) to the observations of AVHRR (52.9% and 6.2%, respectively). Figure 4-42 in the Product Validation and Intercomparison Report - ESA Cloud_cci (2020)[49] shows the time series of the monthly cloud fraction (here $H_{all}$) for afternoon satellite retrievals, where it can be seen that the largest monthly differences within one year are mainly confined to a 5% range (and this is always the case considering e.g., NOAA-19 only), so the large differences in $H_{all}$ between the DYAMOND datasets (including only February) and the other datasets (including full years) are not expected to depend on the different time periods analysed.

The observational datasets (AVHRR and MODIS) show large differences in both $H_{all}$ and $H_{MP}$. This may be due to the use of only the 3.7 µm channel in MODIS to retrieve cloud optical properties, as well as the combination of daily time frames of Terra and Aqua observations to provide cloud top related parameters at the final $1^{\circ }$x$1^{\circ }$ spatial resolution. Previous validation studies including AVHRR (e.g.,[6,8]) have shown the strong performance of this dataset compared to others derived from passive and active satellite sensors. In addition, $H_{all}$ from AVHRR is closer to previous satellite studies (ranging from $61.9\%$ to $67.2\%$ for afternoon satellites; [49]) than MODIS. For these reasons, we subsequently focus on the observational data from AVHRR and using the MODIS dataset for its intended comparisons with COSP outputs from CAM6-Oslo.

4.2. Geographical Distribution and Occurrence of Clouds

Figure 1 shows the relative frequency of occurrence of all clouds derived from AVHRR and MODIS, alongside four model datasets.

Cloud types at the same height are combined to obtain high-, mid-, and low-level clouds such that ICON and GEOS, which do not provide COT, can be included in the comparison. The relative frequency of occurrence $H_{all}$ (calculated with Equation 4) refers to CTP > 680 hPa for low-level clouds, 440 hPa < CTP ≤ 680 hPa for mid-level clouds, and CTP < 440 hPa for high-level clouds, and it varies significantly between the datasets. While AVHRR and MODIS retrieve more low-level clouds (respectively, $H_{all}=37.83\%$ and $50.59\%$) than mid- ($H_{all}=7.9\%$ and $19.3\%$) and high-level clouds ($H_{all}=7.15\%$ and $17.41\%$), most of the models show more high-level ($H_{all}$ ranging from $27.17\%$ to $52.06\%$) than low-level clouds ($H_{all}$ ranging from $13.35\%$ to $27.10\%$). Despite all these differences, which may also be related to the filter applied to COT for the cloud top definition used in AVHRR, MODIS, and CAM6-Oslo, the analysed datasets mostly agree on the regions where clouds at certain heights occur. The geographical distributions in Figure 1 show that low-level clouds are mostly marine clouds while mid- and high-level clouds are found over both continents and oceans, consistent with other studies [33,35,36,50]. In particular, high- and mid-level clouds are mainly found in the extratropics and around the Intertropical Convergence Zone (ITCZ), which are areas often subject to synoptic- and mesoscale events, and low-level clouds are found mainly over the oceans, and especially over their eastern edges, where they form below a capping inversion created by the atmospheric subsidence above a low sea surface temperature.

However, the geographical distributions also show differences between the datasets: for example, low-level clouds are more widely distributed over the western sides of the oceans in AVHRR and MODIS than in the models (e.g., in the Pacific Ocean near South America), while high-level clouds are present on the eastern side of the Pacific Ocean near North America in SCREAM and GEOS, with lower frequency in AVHRR and MODIS, but these are absent in CAM6-Oslo and ICON. The differences in cloud locations between datasets, together with the differences found in $H_{all}$, suggest a possibly wrong representation of the vertical structure of clouds in models.

Similarly to Figure 1 and Figure 2 shows the geographical distributions of the mixed-phase clouds at the three height-levels for the different datasets. The datasets show consistent patterns, with mixed-phase clouds mainly located in the mid-latitudes, in agreement with Hu et al. (2010) [51], as supercooled liquid water has been observed in both convective and stratiform clouds in these regions.

The main contributors to mixed-phase clouds are low-level clouds, except for AVHRR and CAM6-Oslo, for which mixed-phase clouds are mostly in the mid-level height, with $H_{MP}$ of AVHRR in the mid-level being only 0.02% larger than in the low-level height. However, the MODIS product includes more mixed-phase clouds than the other datasets, also present in the tropics in areas where strong updrafts lead to the formation of high-level clouds with ice tops (e.g., the ITCZ). This may be due to the MODIS cloud retrieval, which averages MODIS Aqua and Terra measurements and provide cloud properties daily at $1^{\circ }$ spatial resolution. Therefore, some of the additional mixed-phase clouds in MODIS may be the result of averaging. The differences between the models in the geographical distributions in Figure 2 may imply different dynamics and therefore the possibility of regionally dependent cloud micro- and macrophysics within the same cloud height. Similarly to Figure 1 and Figure 2 shows different cloud locations depending on the datasets. Specific patterns in the DYAMOND models may be due to seasonal effects (e.g., boreal Winter storm track) or could suggest a possible misrepresentation of the cloud vertical structure in the models.

4.3. Temperature Dependence of Mixed-Phase Cloud Types

Figure 3 shows the SLF-CTT joint histograms of the different cloud types for AVHRR normalised as in Equation 7. Each cloud type has a unique SLF-CTT joint histogram, indicating individual cloud processes that could be driven by cloud dynamics or regional effects, e.g., different aerosol composition. Furthermore, CTT decreases from low- to mid- and high-level clouds, as expected from the direct dependence of CTT on the cloud top height, and SLF decreases with decreasing CTT. All datasets display these trends (Figures S1 to S5) although the histograms appear very different at first sight. In the temperature range between $-{40}^{\circ }$C and $0^{\circ }$C, the decrease of SLF at lower temperatures could be driven by heterogeneous ice nucleation as well as secondary ice formation and the Wegener-Bergeron-Findeisen (WBF) process. For example, about 30% of all stratocumulus clouds have an SLF $<0.5$, while only about 10% of stratus clouds have an SLF $<0.5$, indicating very inefficient glaciation processes in stratus clouds compared to stratocumulus clouds or the possibility for stratus clouds to completely glaciate (very efficient glaciation processes) and quickly dissipate due to ice crystal sedimentation or snow. The presence of cirrus clouds in the mixed-phase regime may arise from the cloud-type classification based on COT–CTP thresholds [40]. Previous studies (e.g., [52]) have shown the presence of supercooled liquid droplets in cirrus clouds. However, the presence of cirrus clouds at temperatures characteristic of heterogeneous glaciation processes, as in our results, suggests the presence of lower clouds in the cirrus regime and implies that this regime represents a mixture of clouds with CTP$<440$ hPa and COT $<3.6$.

4.4. Hemispheric Differences of Continental and Maritime Mixed-Phase Clouds

Figure 4 shows, for all datasets, the hemispheric contrast for continental and maritime areas for only cirrostratus and stratocumulus clouds, which we found to generally be the most frequent cloud types in the mixed-phase for high- and low-level clouds, respectively (Table S1). An example of derivation of these difference plots between CFAD-like histograms is shown in Figure S6 for AVHRR.

For a given CTT $<-{15}^{\circ }$C in AVHRR (Figure 4, top row), cirrostratus clouds in the SH show higher SLF than in the NH over both land and the ocean, as displayed by the turquoise colors at high SLF in the difference plot. The same holds only over marine regions for stratocumulus clouds, while the opposite is true over land and at CTT $>-{20}^{\circ }$C, in agreement with Bruno et al. (2021) [6]. CAM6-Oslo shows agreement with AVHRR for stratocumulus clouds over ocean but not over land, while cirrostratus clouds show inconsistent patterns in hemispheric contrasts, with higher SLF in the SH than in the NH only for SLF $<0.35$ when considering statistically significant differences. CAM6-Oslo shows similar results to MODIS only over land for cirrostratus clouds and over both land and ocean for stratocumulus clouds (Figure 4, SLF-CTP histograms, last two rows). The DYAMOND datasets fail in representing the hemispheric differences over continental and maritime regions, with few exceptions at specific temperature ranges.

The presence of an aerosol scheme and aerosol-cloud interactions in CAM6-Oslo (see Section 2) may explain the agreement of this model with observations. However, the simulation of continental stratocumulus clouds failed to capture the hemispheric differences found in AVHRR. This suggests that the micro- and macrophysics of low-level clouds may be more influenced by the Earth’s surface than high-level clouds, making the simulation of low-level clouds more complex, especially over land where orography and aerosol type heterogeneity, distribution and transport can vary significantly in space and time. Nevertheless, CAM6-Oslo shows good agreement with MODIS in representing the hemispheric contrast of continental stratocumulus clouds. However, this is not the case for CTP $>825$ hPa, where warmer temperatures can occur, reflecting the stronger disagreement with AVHRR for CTT$>-{20}^{\circ }$C in the hemispheric contrast for continental stratocumulus clouds. Differences between NH and SH in cloud dynamics and cloud microphysics of the simulated cloud types, as well as the use of both day- and nighttime data and possible seasonal effects influencing results from the DYAMOND datasets (only available for February 2020) may explain the disagreements between the datasets as well as the contrast with the observations. Unfortunately, this study would not benefit from constraining the analysis to observations for February only, as the data would be insufficient due to the temporal (daily) and spatial resolution of AVHRR (coarse-grained by using gridboxes) and MODIS.

4.5. Thermodynamic Dependence of the Liquid Effective Radius in Mixed-Phase Clouds

Figure 5 shows SLF-$\overline{r_{liq}}$ joint histograms of the different cloud types for AVHRR (a), MODIS (b), CAM6-Oslo (c), and GEOS (d). ICON and SCREAM do not provide information on the effective radius, therefore these model datasets could not be included in the comparison. Focusing in particular on AVHRR (Figure 5a) the main results can be summarised as follows:

• Low-level and altostratus clouds are more likely to be liquid.
• In altostratus, altocumulus and high-level clouds, $\overline{r_{liq}}$ increases with low SLF, in agreement with Coopman et al. (2021) [53]. The presence of this trend in cirrus and cirrostratus clouds, which are expected to consist of ice crystals, may be due to a cloud misclassification [40]. In fact, it should be stressed again that the cloud type classification is based on the COT-CTP joint histograms, therefore the same cloud type may undergo different processes depending on where it forms or the atmospheric or dynamical conditions, possibly classifying some clouds in the nearby COT-CTP histogram, i.e., as a different cloud type.
• High frequency of $\overline{r_{liq}}$ between 12 and 13 µm is found in stratocumulus and altostratus clouds, and less but still present in cirrostratus clouds, most likely due to the a priori value assigned to $\overline{r_{liq}}$, which is 12 µm [9].

Some limitations have to be considered for AVHRR. In Bruno et al. (2021) [6], a cloud-phase mismatch between AVHRR-based datasets and CALIOP (Cloud-Aerosol Lidar with Orthogonal Polarization) dataset was shown: AVHRR retrieved some clouds as ice which CALIOP retrieved as liquid. This concerned $9\%$ of the total collocated cloudy pixels. This implies that the liquid effective radius may be affected by errors due to incorrect phase assignment of cloudy pixels. In Stengel et al. (2017) [54] a validation study between several AVHRR datasets is presented. It shows that several datasets derived from passive satellite sensors overestimate the occurrence of liquid clouds at the expense of ice clouds. The sensitivity of passive retrievals to optically very thin ice cloud layers above liquid cloud layers is very low and the filter on COT applied in this work, which removes clouds for COT $<0.3$, may not be sufficient to overcome this problem. In addition, there could always be some ice in the primarily liquid pixels, and these large ice crystals may bias $\overline{r_{liq}}$.

Although the histograms of MODIS (Figure 5(b)) and CAM6-Oslo (Figure 5(c)) look very different from those of AVHRR, these datasets also show an increase in $\overline{r_{liq}}$ with a decrease in SLF for high and mid-level clouds (CAM6-Oslo also for cumulus and stratocumulus clouds and, in general, reaching larger $\overline{r_{liq}}$ values than AVHRR and MODIS). This outcome is likely attributable to the propensity of cold regions, which frequently contain pristine air masses with minimal CCN, enabling the available particles to grow larger than in regions with higher CCN under equivalent humidity. GEOS (Figure 5 (d)) does not show this trend, nor the large variability in radius size shown by the other datasets.

I think the main explanation for this is that cold regions which are more likely to have low SLF are also pristine and have few CCN.

Coopman et al. (2020) [55] has shown that the glaciation of clouds with larger droplets is more likely to occur at warmer temperatures. As SLF decreases from 1 to 0 with decreasing CTT, warmer temperatures can be associated with SLF approaching 1, and this may explain why smaller values of $\overline{r_{liq}}$ are found as SLF approaches 1. This result can be explained by several microphysical processes: in clouds with strong updrafts, ice crystals and liquid droplets experience a decrease in temperature with height and the vapour pressure over ice and over liquid are larger than the in-cloud vapour pressure, leading to the simultaneous growth of liquid droplets and ice crystals at the expense of available water vapour [56]; during ascent in the cloud (where a decrease in SLF is expected due to the decreasing CTT), droplets can collide and aggregate, resulting in larger droplets [3]; collisions in the mixed phase temperature range can cause droplets freezing [57,58] and the formation of secondary ice [59], which can explain the decrease in SLF with the increase in $\overline{r_{liq}}$. The satellite retrievals limit the size of the liquid effective radius to 35µm and 30µm in AVHRR and MODIS, respectively. This is particularly evident in Figure 5(b) for MODIS, where the $\overline{r_{liq}}$ size reaches values below 31µm, although the $\overline{r_{liq}}$ distributions in AVHRR and MODIS are centred on much lower values. In CAM6-Oslo, the $\overline{r_{liq}}$ distribution reaches values up to 39µm. This difference between the observational datasets and CAM6-Oslo may be due to the possible absence of mechanisms in the cloud microphysics scheme in CAM6-Oslo that should limit droplet growth, although the presence of a fixed maximum $\overline{r_{liq}}$ value in the observational retrievals may also lead to a bias. The decrease in $\overline{r_{liq}}$ with decreasing SLF in GEOS may reflect features of the model’s cloud parameterisation. In addition, the GEOS configuration for the DYAMOND Winter project simulations uses a single moment cloud microphysics scheme and the Goddard Chemistry Aerosol Radiation and Transport (GOCART) model, in which the aerosol mass of some of the key aerosol types found in the troposphere with their sources, sinks, and transport, are estimated, but the aerosol-cloud interaction is absent, with possible consequences in the simulated $\overline{r_{liq}}$. In addition, CAM6-Oslo clearly shows lower $\overline{r_{liq}}$ in optically thicker clouds. This result may be due to the Twomey effect: clouds with a higher number concentration of liquid droplets are optically thicker and, for the same amount of water in the cloud, the droplets are smaller as they all compete for water.

4.6. Geographical Distribution of Liquid Droplets in Mixed-Phase Low-Level Clouds

Figure 6 and Figure 7 show the geographical distributions of $\overline{r_{liq}}$ of mixed-phase high- and low-level clouds, respectively, for AVHRR, MODIS, CAM6-Oslo, and GEOS. In both figures, continental clouds in AVHRR and MODIS generally have smaller $\overline{r_{liq}}$ than marine clouds, while this contrast is less pronounced in CAM6-Oslo and absent in GEOS.

Compared to the observations, CAM6-Oslo simulate larger $\overline{r_{liq}}$ in shallow cumulus (Figure 7) and all high-level clouds (Figure 6) and shows lower $\overline{r_{liq}}$ in optically thicker clouds. This result is in line with Figure 5c and possible reasons for that have been discussed in Section 4.5. GEOS underestimates $\overline{r_{liq}}$ and simulates very uniform distributions. This is probably related, again, to the use of a single moment microphysics scheme and the absence of aerosol-cloud interaction.

In Figure 7, the observations show larger values of $\overline{r_{liq}}$ also over some mountainous regions (e.g., the Andes Mountains in South America, the North American Cordillera, and the Himalayas). CAM6-Oslo shows larger values of $\overline{r_{liq}}$ over the Andes Mountains in South America, while GEOS shows much less variability in $\overline{r_{liq}}$ (Figure 7). However, the regions where larger $\overline{r_{liq}}$ is found can vary not only with the dataset but also with the cloud type. For example, the continental region in the SH with larger droplets in AVHRR is found in shallow cumulus clouds, while in MODIS larger droplets are found in stratus clouds. Also, the continental region in the NH with larger droplets in AVHRR includes northeastern China, while in MODIS larger droplets are found in southern Russia and part of the east coast of China. These differences may be related, again, to the averages in the MODIS retrieval, which also causes the presence of mixed-phase clouds in the tropics (as already discussed in Section 4.2). Smaller $\overline{r_{liq}}$ in continental than marine regions may be an indirect effect of the number concentration of aerosols acting as cloud condensation nuclei or ice-nucleating particles, being larger over land than over the ocean, so that more droplets compete for water vapour over land than over the ocean. Larger values of $\overline{r_{liq}}$ found over mountainous regions support the hypothesis that orographic updrafts may cause droplets to grow to larger sizes, although not all the mountainous chains show large $\overline{r_{liq}}$. The continental area from ${30}^{\circ }$ to ${60}^{\circ }$ N is larger than the continental area from ${30}^{\circ }$ to ${60}^{\circ }$ S, which is mainly mountainous. If the presence of larger $\overline{r_{liq}}$ over land in the SH than in the NH for low-level clouds found in AVHRR and partially in MODIS is taken as true, the results in Bruno et al. (2021) [6] with AVHRR showing a lower SLF over land in the SH than in the NH for the same cloud types could be explained by the higher probability of larger droplets (likely grown by orographic updrafts) colliding and freezing at negative temperatures, as well as colliding with ice particles with the possibility to rime and produce secondary ice.

5. Conclusions

In this paper, we have used near-global (${60}^{\circ }$ N to ${60}^{\circ }$ S) data from passive satellite sensors (AVHRR and MODIS), three storm-resolving models from the DYAMOND project (ICON, SCREAM, and GEOS), and one global climate model (CAM6-Oslo) to analyse the cloud top thermodynamic phase, expressed in terms of supercooled liquid fraction (SLF), in relation to cloud top temperature (CTT) and mean liquid effective radius $\overline{r_{liq}}$. We compared the datasets for different cloud types, classified using cloud top pressure and cloud optical thickness thresholds. The comparison included the geographical cloud cover and focused on hemispheric midlatitude (${30}^{\circ }$ to ${60}^{\circ }$) contrasts of continental and maritime regions based on 4 years analysis (from June 1, 2009 to May 31, 2013), except for the storm-resolving models, for which only February 2020 was available. To deal with the different parameters provided by the datasets for the thermodynamic phase of the cloud, we used different metrics for SLF. The main findings are summarized here:

• We derived new reference values for the fraction of mixed-phase clouds up to 60° latitude. While models exhibited a wide spread in occurrence, they generally agreed with observations on their geographical distribution, with mixed-phase clouds mainly located in the midlatitudes.
• Models showed an increase of SLF with CTT, in agreement with the observations, regardless of the cloud types and regions.
• Observations confirmed larger SLF values in the Southern Hemisphere than in the Northern Hemisphere, except for continental stratocumulus, consistent with Bruno et al. (2021) [6]. CAM6-Oslo reproduced these contrasts most closely, while the DYAMOND storm-resolving models showed only limited agreement.
• We introduced SLF–$\overline{r_{liq}}$ joint histograms to link phase partitioning to droplet size. Observations showed an increase in $\overline{r_{liq}}$ as SLF decreased. CAM6-Oslo reproduced this trend, although with a tendency to overestimate droplet size, while GEOS showed a flatter relationship and less variability.
• Geographical distributions of $\overline{r_{liq}}$ from observations show smaller $\overline{r_{liq}}$ over land than over the ocean, whereas models exhibit a less pronounced land/ocean contrast.

Several reasons may explain the disagreements found between models and observations. These include an unrealistic representation of the aerosol influence (especially for low-level clouds) in models, an insufficient representation of the cloud macro- and microphysics in models, but also uncertainties in the observational datasets (e.g., the presence of snow in mountainous regions and the use of only one channel in MODIS to retrieve the cloud optical properties may affect the quality of observations), the differences between datasets in SLF computation and the substantial differences in the satellite retrievals. We have shown that like-for-like comparisons between storm-resolving models and observations are challenging. Future high-resolution modeling experiments could employ satellite simulators or similar tools to enable fairer comparisons between models and observations. Despite these caveats, our results demonstrate that high spatial resolution alone does not guarantee a realistic representation of mixed-phase clouds (specifically, the inability of storm-resolving models to represent SLF hemispheric and land/ocean $\overline{r_{liq}}$ contrasts). Instead, the inclusion of interactive aerosol–cloud coupling and more advanced cloud microphysics schemes, as implemented in CAM6-Oslo, appears to be decisive. The systematic differences among models in reproducing hemispheric contrasts and SLF–$\overline{r_{liq}}$ relationships also have implications for simulated cloud radiative effects, since biases in phase partitioning can alter the cloud radiative effect and thus influence modeled climate feedbacks. By linking supercooled liquid fraction and droplet size, the diagnostic framework developed here provides quantitative benchmarks for evaluating and tuning cloud schemes, offering a practical path toward reducing uncertainties in the representation of mixed-phase clouds and their feedbacks in global and storm-resolving climate models.