The day-to-day co-variability between mineral dust and cloud glaciation: A proxy for heterogeneous freezing

To estimate the global co-variability between mineral dust aerosol and cloud glaciation, an aerosol model reanalysis was combined with satellite retrievals of cloud thermodynamic phase. We used the CALIPSO-GOCCP and DARDAR products from the A-Train satellite constellation to assess whether clouds are composed of liquid or ice and the MACC reanalysis to estimate the dust mixing-ratio in the atmosphere. Night-time retrievals within a temperature range from 10 +3°C to −42°C for the period 2007-2010 were included. The results confirm that the cloud thermodynamic phase is highly dependent on temperature and latitude. However, at midand high latitudes, at equal temperature and within narrow constraints for humidity and static stability the average frequency of fully glaciated clouds increases by +5 to +10% for higher mineral dust mixing-ratios. The differentiation between humidity-stability regimes reduced the confounding influence of meteorology on the observed relationship between dust and cloud ice. Furthermore, for similar mixing-ratios of mineral 15 dust, the cloud ice occurrence-frequency in the Northern Hemisphere was found to be higher than in the Southern Hemisphere at −30°C but lower at −15°C. This may suggest a difference in the susceptibility of cloud glaciation to the presence of dust. Based on previous studies, the differences at −15°C could be explained by higher feldspar fractions in the Southern Hemisphere, while the differences at −30°C may be explained by the higher freezing efficiency of clay minerals in the Northern Hemisphere. 20


DARDAR-MASK
The DARDAR-MASK v1.1.4 product Hogan, 2008, 2010) available at the ICARE data center combines the 125 attenuated backscatter from CALIOP (at 532 nm; sensible to small droplets), the reflectivity from the CPR (at 94 GHz; sensible to larger particles) and the temperature from the ECMWF-AUX product to assess cloud thermodynamic phase. The radar pixels have a horizontal resolution of 1.4 km (cross-track) × 3.5 km (along-track) and a vertical resolution of 500 m, with a nadir angle of 0.16° of the radar beam. A decision is made for each pixel with a 60 m vertical resolution to take advantage of the lidar resolution. These pixels are collocated with the CloudSat footprints (1.1 km horizontal resolution). If 130 the backscatter lidar signal is high (>2·10 -5 m −1 sr −1 ), strongly attenuated (down to at least 10% in the next 480 m) and penetrates less than 300 m into the cloud, it is assumed that supercooled droplets are present. In this case, the pixel is categorized as supercooled or mixed-phase depending on the radar signal, which is assumed a priori to indicate the presence of ice particles. Otherwise, the pixel is categorized as ice (Delanoë et al., 2013;Mioche et al., 2014). For reasons that will become clear later, we will coerce the mixed-phase category into the liquid category. 135

MACC and ERA-Interim
The Monitoring Atmospheric Composition and Climate reanalysis (MACC, Eskes et al., 2015) is based on ECMWF's Integrated Forecasting System (IFS), and simulates the emission, transport and deposition of various aerosol species and trace gases with an output resolution of 1.125° × 1.125° and 60 vertical levels. In this study, we use the dust mixing-ratio and large-scale vertical velocity from the daily MACC reanalysis product on model levels provided by the ECMWF. 140 Additionally, the Relative Humidity (RH) from the ERA-Interim reanalysis daily product (Dee et al., 2011) will be used in Sect. 5. The cloud properties in the MACC reanalysis are derived from the ECMWF Integrated forecast system (IFS Cycle 36r1 4D-Var). This atmospheric model is analogous to the one used in the ERA-Interim reanalysis (IFS Cycle 31r2 4D-Var).
At the time of this study, the new generation of reanalysis based on IFS Cycle 41r was not yet publicly available. However, it is expected that future studies will use the new CAMS (Copernicus Atmosphere Monitoring Service) and ERA5 reanalysis 145 instead of the MACC and Era-Interim reanalysis.
Dust emission in the MACC model is parameterized as a function of the 10 m wind, vegetation, soil moisture and surface albedo. The dust loadings are corrected by the assimilation of the total column AOT at 550 nm retrieved from the MODIS instrument on board NASA's Aqua and Terra satellites. Dust sinks are simulated including dry and wet deposition, as well as in-cloud and below-cloud removal. 150 The freezing efficiency of INPs depends mainly on their surface area concentration (Atkinson et al., 2013;Hartmann et al., 2016;Murray et al., 2011;Niedermeier et al., 2011Niedermeier et al., , 2015Price et al., 2018). While the number concentration of dust aerosol is generally dominated by fine-mode dust (particle diameter < 0.5 µm), the surface area concentration is often determined by both fine and coarse (particle diameter > 1 µm) dust particles (Mahowald et al., 2014). Moreover, the atmospheric lifetime of fine-mode dust is longer than that of coarse-mode dust due to the lower dry deposition rates of finer particles (Mahowald et 155 al., 2014;Seinfeld and Pandis, 1998). Additionally, it has been shown that the MACC model underestimates the coarsemode dust fraction in relation to fine-mode dust (Ansmann et al., 2017;Kok, 2011). In the MACC reanalysis, dust aerosols are represented by three size bins, with size limits of 0.03, 0.55, 0.9 and 20 µm diameter. In this work, we will define the size bin between 0.03 and 0.55 µm as fine-mode dust. Because the fine mode contributes to both the number and surface area concentration it will be used as a proxy for the concentration of dust INP. Although mostly focused on the Northern 160 Hemisphere, several studies have evaluated the simulated dust mixing-ratios from the MACC reanalysis with observations. A mean bias of 25% was found between MACC and LIVAS, a dust product based on CALIPSO observations over Europe, northern Africa and Middle East (Georgoulias et al., 2018). Additionally, the correlation between MACC and AERONET was found to range from 0.6 over the Sahara and Sahel to 0.8 over typical regions of dust transport (Cuevas et al., 2015).
Finally, using shipborne measurements of long-range dust transport, it was found that the MACC model significantly 165 overestimates the fine-dust fraction compared to observations (Ansmann et al., 2017).

Methods
In this section, the different steps in our processing of the datasets presented in Sect. 2 will be described. Fig. 1 presents a flow chart of the data processing and a roadmap for the following subsections.

Selection of cloud profiles 170
To exclude the effects of the scattering of sunlight on the detection of the CALIOP lidar signal, only night-time retrievals were used. Additionally, to avoid biases in the radar retrievals at pixels where the lidar is fully attenuated only pixels where the CALIOP retrieval was classified as cloudy (SR > 5) were used. This warranted a dataset free of lidar-attenuated pixels.
To avoid biases on the radar reflectivity due to rain droplets, only non-precipitating clouds were included in Sect. 4.1-4.3. the products algorithms, the temperature at each profile pixel is interpolated using the information from the reanalyses and between -42°C and +3°C. Then, using this information, we agreggate the pixels into 3 K intervals.
For each of these 3 K intervals, the data is then averaged horizontally using nearest-neighbor regridding. We regridded the 185 dataset first into a Gaussian T63 grid, then aggregated every 16 gridboxes along the longitude (1.875°×30°; lat×lon gridboxes) to better fill the horizontal gaps between the satellite orbits. The Gaussian T63 grid is commonly used in Global Climate Models (Randall et al., 2007) and facilitates comparisons with global simulations of the cloud thermodynamic phase. In section 4.4 and onwards, latitude bands of 30°×360° are used to allow a direct comparison with previous studies (Zhang et al., 2018). 190

3.3
Frequency Phase Ratio and cloud volume weighting Huang et al. (2015) compared the cloud thermodynamic phase retrieved by the DARDAR-MASK and the CALIOP level 2 cloud layer product (Hu et al., 2009) -195 MASK. This discrepancy is attributed partly to the instrumental differences and partly to differences in the classification algorithms.
To assess the differences between the cloud phase from the DARDAR-MASK and CALIPSO-GOCCP products, we defined a new phase ratio based on the DARDAR-MASK classification. In this alternative definition, which we will call ALT-DARDAR, only gridboxes (1.875°×30°×3 K) filled with ice pixels are considered as ice (fully glaciated), so that just a single 200 liquid pixel is enough to define a gridbox as liquid (not fully glaciated). One advantage of this marginal definition is that it ignores cloud ice in mixed-phase clouds, which is mostly only detected as such by the DARDAR-MASK product and neglected by the CALIPSO-GOCCP product. However, this neglection of ice in mixed-phase clouds is only carried out to clarify the differences between the products.
For FPRGOCCP and FPRDARDAR, the FPR is calculated as the ratio of ice pixels to the total number of pixels within each 205 gridbox. The FPR_ALTDARDAR uses gridboxes instead: The gridbox phase is set to one (fully glaciated cloud) if all cloud pixels in the gridbox are classified as ice and zero (not fully glaciated cloud) otherwise. The differences in FPRGOCCP, FPRDARDAR and FPR_ALTDARDAR will be studied in detail in Sect. 4.1-4.2. Here it will be shown that the FPR_ALTDARDAR definition does well in mimicking the limitations of the CALIPSO-GOCCP product.
We expect cloudy overcast gridboxes to be statistically more robust than partly cloudy gridboxes. Cloud cover is generally 210 defined for the whole atmosphere or certain levels (low, mid or high). Therefore, we use instead the cloud volume fraction (also known as three-dimensional or 3-D cloud fraction; Chepfer et al., 2010Chepfer et al., , 2013Li et al., 2017a;Yin et al., 2015) retrieved by the CALIPSO-GOCCP product as a weight for the averages used in Sect. 4.1-4.3. The cloud volume fraction is defined as the number of cloudy pixels divided by the total number of pixels within a given length-height domain along the satellite swath. The length is the segment of the satellite swath crossing a given gridbox, and the height interval corresponds 215 to each temperature bin (3 K) in this study. The main benefit from using the cloud volume fraction instead of the cloud cover is that the former is defined for each temperature bin. This allows to differentiate between vertically thick and shallow clouds. Using the cloud volume fraction as a weight results in a higher representation of clouds with larger spatial extension (vertical as well as horizontal). It also introduces a bias towards the cloud tops for thick clouds because the lidar signal is attenuated at higher cloud depths. More details on the spatiotemporal variability of the cloud volume fraction can be found 220 on the supplement (S8) to this article.
In Sect. 4.1, the adjusted ice volume fraction * = (2 · − 1) · (3.1) is used instead of the traditional FPR, with cvf the cloud volume fraction obtained from the GOCCP product. The adjusted FPR* helps to visualize the cloud thermodynamic phase of the significant (high cvf) clouds in the retrieval. In Sect 4.

Meteorological regimes 230
Dust aerosol can produce or be accompanied by changes in atmospheric stability and relative humidity. To disentangle such effects, we constrain the cloud environment in Sect. 4.4 using the air relative humidity, the lower troposphere static stability With Tx and Px the temperature and pressure at the surface or at x hPa using the pressure levels of the ERA-Interim reanalysis. R is the gas constant and Cp the specific heat capacity of air (Klein and Hartmann, 1993). The relative humidity is obtained directly from the ECMWF-AUX reanalysis.
To increase the sample size prior to the regime classification, we included back gridboxes containing precipitating or convective clouds back into the dataset. However, most of such clouds are expected to fall into high RH and low LTSS 240 regimes and, therefore, could still be excluded later on.

Classification of dust loads and day-to-day correlation
In contrast to previous studies, in this work we want to isolate the day-to-day correlation between dust aerosol and cloud phase. To exclude the spatial component of the correlation, the complete time span 2007-2010 is used to determine the timedeciles of dust mixing-ratio using the MACC reanalysis for each volume gridbox (latitude, longitude, temperature). These 245 deciles are used to sort the daily data depending on the daily dust mixing-ratio into 10 different decile ranks. These ranks can be also understood as dust mixing-ratio bins (from now on simply deciles).
Next, we also need to exclude the seasonal component of the temporal correlation. With this purpose, for each 3 K temperature bin and each gridbox the daily data is averaged within each dust decile and each month of the year. This is done as a multiyear average (e.g., January containing Jan'07, Jan'08, Jan'09 and Jan'10). The resulting field contains one extra 250 dimension for each gridbox (month, dust decile, temperature, latitude, longitude). Fig. 2 present a visualization of this process.

3.6
Data availability and averaging order latitudes and about 500 to 1500 datapoints in the high-latitudes, with the lowest sample size found for the high southern latitudes. Potential reasons for missing data are: • The satellite swaths (orbits) produce a different density of retrieved profiles at different latitudes.
• Using only night-time data, the sample size in the meteorological summer time (shorter nights) is lower.
• The cloud phase retrievals are less frequent for seasons, regions and heights with low cloud cover. 260 • At high latitudes, relatively warm temperatures (e.g., -15°C) exceeding the surface temperature can be found, and therefore no information is available for such temperatures (e.g., over Antarctica in winter).
To avoid artefacts arising from the averaging of dimensions containing missing values, the averaging order of dimensions was defined (going from the first to the last dimension to be averaged) as: longitude, month, decile, latitude, temperature.
Latitude and temperature are averaged last because of the higher associated correlations with cloud phase (Sect. 4.2 -4.3 of 265 this study; Choi et al., 2010;Tan et al., 2014). Each 1.875°×30° gridbox contains on average 100 to 200 datapoints at −15°C (within a 12 K range) in the mid-latitudes. Meanwhile, in the subtropics and in the high latitudes, the sample size is much more heterogeneously distributed and can drop below 50 datapoints near the poles and in subsidence regions. A detailed view of the spatiotemporal distribution of the sample size for stratiform clouds can be found in the supplement (S14) to this article. 270

Case study
To better understand the differences between the ice-to-liquid ratio retrieved in the DARDAR-MASK and CALIPSO-GOCCP product, this section provides a detailed case study of a stratiform cloud scenario, in which four stratiform cloud types from the CloudSat classification are included -stratocumulus(low-level clouds), altostratus and altocumulus (mid-275 level clouds), and cirrus (high-level clouds). Although not present in the case study, Nimbostratus are included in the analysis of cloud phase as well and are particularly important in the high latitudes. Stratus clouds are defined for temperatures above 0°C; therefore, they are not relevant for this study. Finally, the horizontal extension of Cumulus and Deep Convection clouds is very low compared to the stratiform clouds and can be therefore ignored in our study, especially outside the tropics (Sassen and Wang, 2008). 280 Fig. 4 shows a case study at 9:50 UTC on Dec 14, 2010 over the Southern Ocean for temperatures between −42°C and +3°C.
This A-train segment has been already chosen for a previous case study by Huang et al. (2015) due to the variety of cloud types it contains. Fig. 4a-b show for the same segment the cloud volume cover (CALIPSO-GOCCP) of clouds classified (2B-CLDCLASS) as cirrus or altocumulus ( Fig. 4a) and as altostratus or stratocumulus (Fig. 4b). These cloud types are frequently thin enough to be penetrated by lidar and radar systems and are therefore a good target to study cloud glaciation 285 processes (Bühl et al., 2016;D.Zhang et al., 2010b). Moreover, stratiform clouds have also simpler microphysics compared to convective clouds, where the dynamical forcing is usually stronger. Fig. 4c shows the mixing-ratio of fine (0.03µm-0.55µm) dust aerosol (MACC reanalysis) for the same vertical plane. 40°S and +3°C to −6°C, the ice virgae falling from the cloud (FPRDARDAR) are missed in the FPRGOCCP. Because this study aims at assessing the occurrence-frequency of fully glaciated clouds, such mixed-phase clouds are then reclassified in FPR_ALTDARDAR as liquid clouds. A similar case is observed for the stratocumulus clouds at 50-55°S and +3°C to −6°C, and for the altostratus at 35-45°S below the −20°C isotherm (at higher temperatures). Finally, the cirrus clouds above −33°C remain nearly unaffected by the reclassification in FPR_ALTDARDAR as it is classified as fully glaciated. Clouds between -295 38°N and -44°N, ranging from -6°C to -33°C in temperature, are classified mostly as altostratus by the 2B-CLDCLASS product. These altostratus clouds offer a good opportunity to compare the three FPR variables in more detail.

FPRGOCCP:
The detected ice virgae below the liquid cloud top suggest that the cloud top did not fully attenuate the lidar signal (not optically thick enough). The number and/or size of the ice particles near the cloud top probably was not enough to increase the depolarization ratio above the threshold value for the GOCCP algorithm and was therefore classified as liquid. 300 FPRDARDAR: In the decision tree of the DARDAR algorithm there are multiple alternatives for a mixture of cloud droplets and ice particles (e.g., at cloud top) to be classified as ice only (Mioche et al., 2014): a) If the lidar backscatter signal (β) is lower than 2.10 −5 m −1 sr −1 b) If not a): If it is weakly attenuated (less than 10 times) or not rapidly attenuated (at a depth larger than 480 m).

c) If not b):
If the layer thickness of the cloud is larger than 300 m. This is equivalent to 5 pixels with a lidar vertical 305 resolution of 60 m. Therefore, there are many cases where a mixed-phase cloud (and especially an optically thin stratiform cloud) can be missclassified as ice only in the DARDAR product and consequently in the FPR_DARDAR variable. In this specific case, we speculate that c) is the most probable cause because of the large vertical extent of the clouds around 1 to 5 km using a moist adiabatic lapse rate of −6 K/km for the estimation). 310 FPR_ALTDARDAR: In the case of droplets and ice particles coexisting at cloud top, we expect that at some location the cloud droplets will be enough in number for the pixel to be classified as liquid (strong attenuation) in the DARDAR-MASK algorithm. If this is the case, the entire gridbox value of FPR_ALT_DARDAR will be LIQUID. This should be interpreted as a non-completely glaciated cloud.
In summary, the GOCCP algorithm is unable to detect ice in mixed-phase clouds and the DARDAR algorithm tends to 315 classify mixed-phase clouds as ice. Therefore, we avoid using the frequency of cloud ice (FPR) to compare the GOCCP and DARDAR products. Instead, we introduced FPR_ALTDARDAR, which has been defined to address the limitations of both products. In FPR_ALTDARDAR, a significant portion of mixed-phase clouds that would otherwise be classified as ICE are now classified as LIQUID. This however partly reintroduces the inability of the GOCCP algorithm to detect ice in mixed-phase clouds. Therefore, the frequency of completely glaciated clouds, which is represented by FPR_ALTDARDAR and FPRGOCCP, 320 allows a better comparison of both algorithms, mostly by ignoring ice virgae in FPR_ALTDARDAR when cloud droplets are also present in the same gridbox. This idea is summarized in Table 1.

Temperature dependence
Mixed-phase clouds between 0°C and −25°C are usually topped by a liquid layer (Ansmann et al., 2008;De Boer et al., 2011;Westbrook and Heymsfield, 2011). Below this layer, there is often a thicker layer containing ice particles. Because the 325 CPR is more sensitive to larger particles, this results in a large fraction of the cloud classified as ICE in the DARDAR-MASK. In contrast, the CALIOP backscatter signal is usually already strongly attenuated at such depths and often cannot detect large ice particles. Therefore, the CALIPSO-GOCCP algorithm usually classifies the whole cloud layer as liquid (Huang et al., 2012;. As a result, FPRDARDAR tends to be higher than FPRGOCCP. about 20 % at −1.5°C and down to 0 % at +1.5°C. This temperature dependence between −42°C and 0°C is also observed for a wide range of parameterizations in global climate models (Cesana et al., 2015), in ground-based measurements (Kanitz et al., 2011), in spaceborne lidar measurements (Tan et al., 2014) and in aircraft measurements (McCoy et al., 2016). However, for the same temperature range, the FPRDARDAR only decreases down to 60 % at 1.5°C. This is partly due to the higher sensitivity of the radar to ice particles, especially falling ice. Additionally, in the DARDAR algorithm water can be still 335 classified as ice at +1.5°C due to the melting layer being set to a wet-bulb temperature (Tw) of 0°C. This allows the detection of ice at temperatures slightly above 0°C dry-bulb temperatures (named simply temperature in this work). For instance, at a relative humidity of 50%, a temperature of about +2.5°C would correspond to a Tw of −2.5°C. Nevertheless, this last effect is not relevant for temperatures below freezing. In contrast, FPR_ALTDARDAR follows very closely the pattern of the FPRGOCCP down to −1.5°C. The absolute differences of the global averaged FPR_ALTDARDAR and FPRGOCCP are less 340 than 10 % between −42°C and 0°C. This shows that the temperature dependence of the alternative phase ratio FPR_ALTDARDAR and FPRGOCCP agree better than for FPRDARDAR. Therefore, for the rest of the study, only FPR_ALTDARDAR and FPRGOCCP will be considered.
Additionally, the average fine-mode dust mixing-ratio is also shown in Fig. 5. At the height of the 0°C isotherm, the mixingratio is on average higher than at the −42°C isotherm (note the logarithmic right y-axis). This reflects the fact that, on 345 average, dust mixing-ratios tend to be higher near the dust sources at the surface. However, this does not imply any general relationship between dust and temperature. Moreover, instant vertical profiles of dust loading and temperature may differ greatly from this average, especially in the long-range transport of dust plumes.

4.3
Latitude dependence Fig. 6 shows the latitudinal dependence of dust and cloud thermodynamic phase at −30°C (averaged from −36°C to −24°C; 350 Fig5a) and at −15°C (averaged from −21°C to -9°C; Fig 5b). For both temperature ranges shown in Fig. 6 the absolute maximum of FPR is located near the Equator. At −30°C, the maximum is 85% for FPRGOCCP and 78% for FPR_ALTDARDAR and at −15°C the FPR from both products peak at 44%. These maxima are probably associated with the enhanced homogeneous freezing in the tropics at temperatures below −40°C and the resulting downward transport of cloud ice.
Similarly, minima are observed towards the high latitudes. At −30°C, the FPRGOCCP has two local maxima with values of 76 355 % and 84 % near 39°S and 39°N, respectively. Similar local maxima are observed for the FPR_ALTDARDAR but at higher latitudes, at 61°S and 61°N with values 69 % and 74 %. At −30°C, both products show a higher FPR in the Northern Hemisphere than in the Southern Hemisphere, in particular for the high latitudes. This higher FPR coincides with the higher average dust mixing-ratio in the Northern Hemisphere. Such positive spatial correlations between FPR and dust aerosol have been already pointed out using the dust occurrence-frequency derived from CALIOP (Choi et al., 2010;Tan et al., 2014;360 Zhang et al., 2012).
In comparison, the differences between FPRGOCCP and FPR_ALTDARDAR at −15°C are much lower than at −30°C as shown in Fig. 6b. Moreover, the FPRGOCCP at −15°C is lower than the FPR_ALTDARDAR at the southern mid-latitudes and northern high-latitudes. In the southern high latitudes, for both variables, a local minimum near 73°S is followed by a steep increase at 84°S. The larger standard deviation in these latitudes is possibly a result of the low sample size in the region, as mentioned 365 in Sect. 2. However, the higher FPR in the southern than in the northern polar region is consistent with the fraction of ice clouds reported previously in the literature at −20°C (Li et al., 2017). On the other hand, it has been shown that the orographic forcing in Antarctica can lead to high ice water contents for maritime air intrusions (Scott and Lubin, 2016). In other words, maritime air intrusions associated with higher temperatures, higher concentrations of INP and stronger vertical motions could explain the observed pattern in the southern polar regions. However, the low sample size near the South Pole 370 ( Fig. 3 and supplement material S.14.b) and the low altitude of the -15°C isotherm (S.12.b) result in a lower confidence in the results for this region. For example, at −15°C, the zonal standard deviation of the FPR significantly increases from 60°S towards the South Pole -from about ±0.08 to ±0.16 in Fig.6a -at the same time that the sample size decreases from 2200 to 300 (Fig.3).
For the clouds studied, the time-averaged large-scale vertical velocity (from the MACC reanalysis) is somewhat regionally 375 correlated with the FPR at −15°C -with a pearson correlation coefficient of 0.47 using zonal averages and of 0.31 using the 30°×1.875° gridbox averages. Moreover, in another study, the spatial correlation between large-scale updraft velocity at 500 hPa was also found to be positively correlated (spatially) to the occurrence-frequency of ice clouds at −20°C (Li et al., 2017a). In other words, both the dust mixing-ratio and the large-scale vertical velocity appear to be to some extent correlated (spatially) to the FPR. There are some plausible explanations for this: 380 • The spatial correlation can be a result of an enhanced transport of water vapor to higher levels at temperatures below −40°C and the subsequent sedimentation of ice crystals from the homogeneous regime (cloud seeding).
• The updrafts are associated with higher availability of INP at the cloud level (from below the cloud) and the effect is large enough to mask the enhanced droplet growth typically associated with updrafts.
• The updrafts enhance a certain type of heterogeneous nucleation requiring saturation over liquid water (e.g., 385 immersion freezing). Updrafts generate a local adiabatic cooling, possibly activating INPs that may have not been active before at higher temperatures.
However, to understand which (if any) of these explanations influence the freezing processes inside the cloud remains a complex challenge of ongoing debate (Sullivan et al., 2016).

4.4
Constraining the influence of static stability and humidity on the dust-cloud-phase relationship 390 To study the temporal correlation between mineral dust mixing-ratio and cloud ice occurrence-frequency (from now on referred to as dust-cloud-phase relationship) it is crucial to systematically classify weather regimes to constrain the meteorological influence. By doing so, the resulting dust-cloud-phase relationship may offer a good insight into how heterogeneous nucleation of mineral dust may affect the day-to-day average cloud thermodynamic phase. In other words, to extract the specific influence of mineral dust on cloud glaciation, it is necessary to identify and constrain relevant 395 meteorological confounding factors (Gryspeerdt et al., 2016). The atmospheric relative humidity and static stability are good candidates for such a confounding factor (Zamora et al., 2018). Both are correlated with the transport of mineral dust and vary between different cloud regimes. Additionally, relative humidity is, next to the temperature, one of the main factors in the initiation of ice nucleation in laboratory studies (Hoose and Möhler, 2012;Welti et al., 2009). The static stability of the atmosphere (see equations 3.4 and 3.5) represents the gravitational resistance of an atmospheric column to vertical motions 400 and is defined as the difference of potential temperature between two pressure levels (Klein and Hartmann, 1993). Because such vertical motions are traduced in a temperature change rate within the air parcel, the static stability can have an important impact on the heterogeneous freezing rates, especially on immersion freezing. We note that the dynamic component of the atmospheric stability is not included in the static stability. Especially in the upper troposphere, atmospheric gravity waves occurring during stable thermal conditions may also result in vertical motions affecting ice production. 405 In general, moist and unstable conditions are associated with enhanced lifting of air that likely causes nucleation of hydrometeors. The effect of humidity and static stability on ice production is not straightforward. In moist convective conditions (high humidity and low static stability) between 0°C and -40°C, the supersaturation of water vapor over liquid increases the liquid formation because ice growing (deposition nucleation) is rather inefficient at this conditions. However, at temperatures below −40°C, ice production by deposition and homogeneous nucleation are favoured, which can result in a 410 higher occurrence of cloud ice in the mixed-phase regime below due to ice sedimentation. To constrain both the atmospheric stability and humidity, a subset of the data must be found within a narrow range of these variables. At the same time, enough data points must still be available to assess the dust-cloud-phase relationship. For this purpose, we use a probability histogram to define the regime bounds such that at least 10% of the data is included in each regime. For example, if the probability density between 4−6 K and 70−80% is 0.01, then 20% of the data is contained between these bounds. The magenta boxes represent the different stability-humidity regimes used for the lower and upper troposphere. To 420 maximize the sample size, for this classification and for the following results precipitating and convective clouds are also included. Fig. 8 shows the dust-cloud-phase relationship for the mid-and high-latitudes separated by humidity and LTSS at −15°C using the FPRGOCCP product and MACC reanalysis. For dust mixing-ratios between 0.1 and 2.0 µg kg −1 , the cloud-phase curve in the mid-latitudes follows a similar logarithmic increase of cloud ice occurrence-frequency of about +6% for low-425 LTSS and +4% for high-LTSS conditions. After analysing 11 years of ground-based lidar measurements in Leipzig, Seifert et al. (2010) reported a slightly higher increase by about +10% between −10°C and −20°C for dust concentrations between 0.001 to 2 µg m −3 (note the different units). In our results at −15°C, the cloud ice occurrence-frequency tends to be higher for higher relative humidity and the LTSS seems to have a major effect on the dust-cloud-phase relationship. For high-LTSS conditions (Fig. 8a-b), a positive dust-cloud-phase correlation can be observed at all four latitude bands. The logarithmic, this means that for the high-latitudes small increases in dust mixing-ratio are associated with a high increase in cloud ice occurrence-frequency. However, the range of ice occurrence-frequency values is higher for the high latitudes.
Particularly for the low-RH regime, the ice occurrence-frequency in the southern high latitudes increases by +8%, and at the mid-latitudes, the increase is only about +4%. In both mid-and high latitudes, the cloud ice occurrence-frequency for the 435 same dust mixing-ratio is about +2% to +8% higher in the Southern than in the Northern Hemisphere. This could point to a factor other than dust aerosol causing an increased ice occurrence-frequency in the Southern Hemisphere. It could also suggest a potential difference in the sensitivity of cloud glaciation to mineral dust between hemispheres. The difference between the Northern and Southern Hemisphere is reduced in the high-RH regime, as well as the standard deviation of the ice occurrence-frequency, possibly due to the higher sample size density in the high-RH regime. For the low-LTSS regime 440 ( Fig. 8c-d), the cloud thermodynamic phase in the high-latitudes remains mostly constant for increasing dust mixing-ratios, and the dust-cloud-phase curves for the mid-latitudes coincide so that the maximum cloud ice occurrence-frequency in the southern mid-latitudes is similar to the minimum in the northern mid-latitudes. Fig. 9 show a similar constraint on humidity and UTSS at −30°C. For all regimes, the cloud ice occurrence-frequency in the southern high latitudes remains almost constant for increasing dust mixing-ratios. For the high-RH regime, the cloud ice 445 occurrence-frequency tends to be higher, especially for the southern high latitudes for which the cloud ice occurrencefrequency is about +4% higher at the high-RH regime. For dust mixing ratios between 0.1 and 1.5 µg kg −1 , the cloud ice occurrence-frequency at −30°C increase by about +5%. The highest increase is found for the northern latitudes. However, the results from the southern mid-latitudes contradict the notion that the INP activity of mineral dust is of secondary importance in the Southern Hemisphere due to low dust aerosol concentrations (Burrows et al., 2013;Kanitz et al., 2011). 450 Nevertheless, recent studies have acknowledged that the importance of mineral dust in the southern latitudes still cannot be ruled out (Vergara-Temprado et al., 2017). Fig. 10 shows the dust-cloud-phase relationship for different UTSS and relative humidity at −22°C. Similar to the results at −15°C and −30°C, the cloud ice occurrence-frequency is higher in the high-RH regime. For high-UTSS conditions, the dust-cloud-phase curves are in closer agreement between the Northern and Southern Hemisphere. Overall, at −22°C the four latitude bands show the best agreement between Northern/Southern Hemisphere and 455 mid-/high-latitudes. Combining the results from all mid-and high latitudes, the ice occurrence-frequency increases by about 25% at high-UTSS conditions and by about 20% at low-UTSS conditions for mixing ratios between 0.01 and 1.0 µg kg −1 at −22°C. This suggests that the dust mixing-ratio may explain both the day-to-day differences in cloud ice occurrencefrequency and the differences between latitudes.
At all temperatures studied, higher humidity values were associated with a higher cloud ice occurrence-frequency. For 460 similar dust loadings, the cloud ice occurrence-frequency was found to be higher at the mid-latitudes than at the highlatitudes. However, against our expectations, for similar dust loadings the cloud ice occurrence-frequency at −15°C was higher in the Southern than in the Northern Hemisphere.

Discussion
Some studies have already suggested that the lower occurrence-frequency of cloud ice in the higher latitudes may be 465 associated with lower INP concentrations (Li et al., 2017a;Tan et al., 2014;Zhang et al., 2012). This hypothesis has been supported mainly by the spatial correlation between the dust relative aerosol frequency and the occurrence-frequency of ice clouds retrieved from satellite observations. However, evidence of the global temporal co-variability between INP and ice occurrence-frequency on a day-to-day basis was lacking up to now. Furthermore, by studying the temporal correlation between mineral dust and cloud ice occurrence-frequency it is possible to extract new information about the differences in 470 cloud glaciation at different latitudes and to connect these differences to previous studies of heterogeneous freezing.
Particularly, our results may be used to evaluate our current knowledge of the global differences in the mineralogy of dust aerosol and its freezing efficiency.

North-South contrast
We have found that the ice occurrence-frequency can vary at different latitudes even for similar mixing-ratios of mineral 475 dust. This could be explained by differences in the mineralogical composition of the mineral dust aerosol at the Southern and Northern Hemisphere. It has been suggested that the freezing efficiency of clay minerals from the Northern Hemispheres (composed mostly from Illite and Smectite minerals) can be well represented by the mineral Montmorillonite while the Southern clay minerals are better represented by the mineral Kaolinite (Claquin et al., 1999;Hoose et al., 2008). The freezing efficiency of Kaolinite and Montmorillonite are known for the immersion and contact freezing mode (Diehl et al., 480 2006;Diehl & Wurzler, 2004). Following this simplification, the immersion freezing rates at −30°C would be 300 times higher in the Northern than in the Southern Hemisphere. This could explain the higher ice occurrence-frequency in the Northern Hemisphere relative to the Southern Hemisphere for similar dust mixing-ratios at −30°C. Below −25°C, the contact freezing is expected to dominate over immersion freezing. However, for contact freezing between −25°C and −16°C the freezing rate is similar for Kaolinite and Montmorillonite. This again may explain why the ice occurrence-frequency in the 485 Northern Hemisphere is only slightly higher for similar dust mixing-ratios at −22°C. Finally, between −15°C and −4°C, the contact freezing efficiency of Montmorillonite is slightly higher than for Kaolinite. However, this fails to explain the higher ice occurrence-frequency found in the Southern Hemispheres at −15°C. Nevertheless, at such high temperatures, other dust minerals like feldspar mineral are much more efficient as ice nucleating particles than clay minerals (Atkinson et al., 2013).
Moreover, it could be that the effect of such feldspar minerals dominates over the effect of clay minerals at high 490 temperatures. Indeed, such efficient minerals are believed to be quickly depleted trough heterogeneous freezing, so that only few would reach lower temperatures. Therefore, they are likely more relevant at temperatures above −20°C, where the immersion efficiency of clay minerals quickly decay (Boose et al., 2016;Broadley et al., 2012;Murray et al., 2011). If feldspar minerals do dominate the heterogeneous freezing due to mineral dust above −20°C, then the higher cloud ice occurrence-frequency in the Southern Hemisphere may be due to a higher fraction (or higher efficiency) of feldspar minerals in the southern dust particles. Some evidence for this has been already found by comparing the immersion freezing efficiency of dust particles from different deserts worldwide (Boose et al., 2016). In these results, the immersion efficiency of dust particles lays mostly between Kaolinite and K-feldspar. The dust samples from sources in the Southern Hemisphere (Australia, Etosha and Atacama milled) have a higher freezing efficiency than most of the samples from the Northern Hemisphere sources including Saharan sources for temperatures below −24°C. Although only four of these samples were 500 studied for higher temperatures, between −23°C and −11°C it was again a sample from the Southern Hemisphere (Atacama milled) which exhibited the highest freezing efficiency. Assuming that the higher freezing efficiency of the southern dust sources may be extrapolated to temperatures above −20°C, the higher immersion efficiency of southern mineral dust, possibly due to higher feldspar fractions, may explain the higher ice occurrence-frequency in the Southern Hemisphere at −15°C. The highly efficient particles, most likely feldspar minerals, would be quickly depleted at temperatures around 505 −15°C and would therefore not interfere with the Kaolinite-Illite(Montmorillonite) differences at −30°C. Furthermore, such a depletion of highly efficient INP during the transport of dust aerosol may also explain the higher ice occurrence-frequency at the mid-latitudes compared to the high-latitudes for similar mixing-ratios of mineral dust, especially at higher temperatures.
The ageing (e.g., internal mixing with sulfate or "coating") of dust particles may also reduce the freezing efficiency of dust aerosol during the transport from low to high latitudes. The hypotheses explaining the differences in the freezing behaviour 510 of dust between the Northern and Southern Hemisphere are summarized in Table 2.

Assumptions and uncertainties
In the analysis presented above, certain assumptions were made to assess the potential effect of mineral dust on cloud thermodynamic phase. In this section, these assumptions and the uncertainties that arise from them, as well as the subsequent 515 limitations of the resulting interpretation will be discussed.
Concerning the vertical resolutions of the different products, the choice of 3 K bins is based on the resolution of the CPR (480 m) -used in the DARDAR-MASK product -and the original 3 K bins of the CALIPSO-GOCCP product. Using a coarser vertical resolution (e.g., 6 K bins) would hinder the assessment of the role of dust as INP. For example, a decrease of 3 K in temperature is roughly equivalent to a fivefold increase in INP concentrations (e.g., Niemand et al., 2012). Because at 520 the mid-and high-latitudes the typical standard deviation of the day-to-day dust mixing-ratio corresponds to roughly a fourfold increase from the mean (See supplement figure S.5), we expect that the variability of dust loading should dominate over temperature variations, given a temperature constraint of 3 K or less. The statistical distribution of the phase ratio also limits the resolution options. The cloud phase values for single pixels in the DARDAR and GOCCP products are binary (1 or 0). Therefore, a minimal sample size is required for the averaged cloud phase ratio -within a certain temperature range, 525 gridbox and percentile of dust -to achieve a normal distribution along time and space, which allows interpreting the correlation with dust loading directly. For this reason, temperature bins smaller than 3 K result in a less normally-distributed cloud phase ratio.
As mentioned in section 3.5, we excluded the seasonal component of the dust-cloud-phase correlation by calculating the deciles independently for each month of the year. However, shorter cycles (e.g., weather variability) may still have an 530 influence in the variabiliry of dust and cloud phase. Although the cloud phase may be affected by such cycles (e.g., more liquid clouds at convective fronts and more cirrus clouds at the detrainment regions), it is still possible to distinguish between dusty and non-dusty conditions at each point of the weather cycle. Therefore, once we average over the weather cycle -using monthly means inside each dust percentile -we expect the dust-cloud-phase relationship to be dominated by the microphysical effect of dust on cloud phase. 535 Despite the long period (2007)(2008)(2009)(2010) used in the study, a significant fraction of the 5-dimensional space used for our analysis (10 dust deciles, 12 months, 15 temperature bins, 96 latitudes, and 12 longitudes) is sparsely sampled or even contains missing values. In the high-latitudes, a sampling bias exists towards the respective winter seasons because very few nighttime retrievals are available in summer. However, the seasonal variability was not found to be a dominating factor in the day-to-day impact of dust mixing-ratio on the FPR (See S.19 in the Supplement to this article). Furthermore, many factors 540 may contribute to higher standard deviations for the ice occurrence, including: • Changes in dynamical forcing (e.g., updrafts) and cloud regimes • Temperature changes after cloud glaciation (e.g., latent heat release) • Ice sedimentation from above (cloud seeding), and INPs other than dust • Cloud vertical distribution within the studied temperature ranges 545 • Turbulence favouring aerosol mixing and sub-grid temperature fluctuations • Differences in dust mineral composition, electric charge and/or size • Coatings (e.g. Sulfate) affecting aerosol solubility and freezing efficiency • Subsetting of the data (e.g., only night-time retrievals) Additionally, some issues arise from the coarse spatial resolution used in our study. A high dust mixing-ratio simulated in a 550 gridbox indicated as cloudy by the satellite observations does not ensure that the dust is actually mixed with the cloud. The subgrid-distribution of dust relative to the exact cloud position remains unresolved. Higher dust mixing-ratios should be interpreted as an indicator or a higher probability that a significant amount of dust was mixed with a collocated cloud. This mixing may have happened during or before the observation by the satellite. The latter is only true if both the cloud and the dust followed the same trajectory up to the moment of the observation. Overall, at coarse resolutions, the combination of 555 modelled dust concentrations with satellite-retrieved cloud properties cannot guarantee the mixture of aerosol and clouds (R. Li et al., 2017b). Similarly, the atmospheric parameters obtained from the reanalysis may not match the conditions for the exact position of the clouds in the satellite retrievals. However, the atmospheric parameters are expected to match on average the large-scale conditions influencing the aerosol-cloud interactions.
In general, we expect that the assimilation of the total AOD from MODIS in the MACC reanalysis produce a fair estimation 560 of the large-scale aerosol conditions on a day-to-day basis. At least for the Northern Hemisphere, this has been already validated with in situ measurements (Cuevas et al., 2015). As for the consistency among the different reanalyses, both the ERA-Interim and the MACC reanalysis are based on the IFS model and use a similar assimilation algorithm. Among the different satellite products, both CALIPSO-GOCCP and DARDAR-MASK rely on CALIOP to determine the presence of clouds. Nevertheless, the reader should be aware that several uncertainties remain, for example, between the meteorology in 565 the reanalysis and in the real atmosphere, particularly on the sub-grid scale. In the worst-case that the reanalyses are entirely inconsistent with the retrievals of cloud phase, we expect the result would be the lack of correlation between dust and the ice occurrence (Fig 8-10). In other words, given the large dataset included in the study, we expect that mismatches between reanalysis and cloud retrievals would cause an underestimation -and not an overestimation -of the dust-cloud-phase correlation. 570 Concerning the interpretation of our results, it cannot be ruled out that the increase in ice cloud occurrence in the Southern Hemisphere for higher dust loading arises from other types of INP such as biogenic aerosol (Burrows et al., 2013;O'Sullivan et al., 2018;Petters and Wright, 2015) or background free-tropospheric aerosol (Lacher et al., 2018) , which could be misclassified as mineral dust in the reanalysis. Similarly, a possible correlation between ice cloud occurrence and the atmospheric conditions leading to the emission and transport of mineral dust should be further investigated (e.g., dusty 575 air masses from land are usually warmer and drier). Another interesting explanation of the results presented in this study could involve the mixing of mineral dust particles with ice nucleation active macromolecules . Such particles are in the size of few 10 nm (Fröhlich-Nowoisky et al., 2015) and would therefore not be detected if mixed with dust aerosol. Furthermore, biases such as the overestimation of the fine-mode dust aerosol in the MACC reanalysis (Ansmann et al., 2017;Kok, 2011) may shift the mixing-ratios shown in Sect. 4.4. However, as long as such 580 biases are not limited to certain meteorological conditions, the cloud phase averaged inside each dust decile should remain unaffected.
In general, meteorological parameters have a larger impact on cloud properties than aerosols do (Gryspeerdt et al., 2016).
For example, different updraft regimes can change the aerosol-cloud interactions in warm clouds by an order of magnitude.
Therefore, it is important to study how such meteorological parameters relate to the dust aerosol loading. With this purpose, 585 Fig. 11 shows the mean relative humidity, cloud height and large-scale updraft at −15°C for the different fine-mode dust mixing-ratio deciles and for the four latitude bands studied in Sect. 4.4. Firstly, the correlation between fine-mode dust mixing-ratio from the MACC reanalysis and the RH from the ERA-Interim reanalysis -weighted by cloud volume fraction -was found to be negative (Fig. 11a). We note that the RH from the ERA-Interim reanalysis represents the conditions at a large-scale and not the conditions at a specific location and the moment of the interaction between dust aerosol and 590 supercooled cloud droplets. Still, this relationship is consistent with the intuition that dust is mostly associated with drier air masses. Second, The significant positive correlation found between dust aerosol mixing-ratio and the height of the isotherms (weighted by cloud volume fraction) points to an important source of uncertainty (Fig. 11b). This could be due to clouds being detected in a higher temperature bin after being glaciated at lower temperatures, thus erroneously suggesting an enhanced glaciation occurrence frequency at higher temperatures. Therefore, it is crucial for future studies to take into 595 account this possibility when studying the occurrence of ice clouds at a certain isotherm. More details on the spatiotemporal variability of the cloud height can be found in the supplement (S12) to this article. Lastly, Fig. 11c shows a positive correlation between the fine-mode dust and the large-scale vertical velocity from the MACC reanalysis at −15°C. Updrafts favour saturation over liquid water and therefore CCN activation, droplet growth and inhibition of the WBF (Wegener-Bergeron-Findeisen) process. Therefore, a positive dust-updraft correlation could lead to an underestimation of the dust-600 cloud-phase relationship.
In summary, much of the co-variability between dust, humidity, updrafts, temperature and cloud ice occurrence-frequency is still poorly understood. However, we expect that the constrains on humidity and static stability minimized most of the biases discussed in this section.

Conclusions 605
For the first time, the MACC reanalysis was combined with satellite-retrieved cloud thermodynamic phase to investigate the potential effect of mineral dust as INP on cloud glaciation on a day-to-day basis at a global scale. Satellite products of cloud thermodynamic phase for the period 2007-2010 were included. We focused on clouds observed at night-time in the mid-and high latitudes. Our main findings can be summarized as follows: 1. Between −36°C and −9°C, day-to-day increases in fine-mode dust mixing-ratio (from lowest to highest decile) were 610 mostly associated with increases in the day-to-day cloud ice occurrence-frequency (FPR) of about 5% to 10% in the mid-and high-latitudes.
2. The response of cloud ice occurrence-frequency to variations in the fine-mode dust mixing-ratio was similar between the mid-and high-latitudes and between Southern and Northern Hemispheres. Moreover, increases in FPR from first to last dust decile were also present in the northern and southern high-latitudes, even though dust aerosol 615 is believed to play a minor role in cloud glaciation in the Antarctic region.
3. Using constraints on atmospheric humidity and static stability we could partly remove the confounding effects due to meteorological changes associated with dust aerosol.
4. The results also suggest the existence of different sensitivities to mineral dust for different latitude bands. The north-south differences in ice occurrence-frequency for similar mineral dust mixing-ratios agree with previous 620 studies on the mineralogical differences between Southern and Northern Hemisphere. A larger fraction of feldspar in the Southern Hemisphere could explain the differences at −15°C, and the higher freezing efficiency of Illite and Smectite (more abundant in the Northern Hemisphere) over Kaolinite (more abundant in the Southern Hemisphere) could explain the differences at −30°C.
We believe these new findings may have an important influence on improving the understanding of heterogeneous freezing 625 and the indirect radiative impact of aerosol-cloud interactions. The authors hope that the results of this work will also motivate further research, including field campaigns in remote regions to study the day-to-day variability of cloud thermodynamic phase and the role of mineral dust in ice formation, satellite-based studies of associated changes in the radiative fluxes, and modelling studies to test the representation and relevance of specific processes involved in ice formation and mineral dust transport. Such studies could help to further improve our understanding of the influence of 630 mineral dust on cloud glaciation and the climate system.

Author contribution
DV, IT, BH and PS contributed to the design of the study. DV processed the datasets, performed the analysis, designed the figures and drafted the manuscript. All authors contributed valuable feedback throughout the process. All authors helped with the discussion of the results and contributed to the final manuscript. 635 Calculate 11 day-to-day (2007)(2008)(2009)(2010) percentiles (0 th , 10 th , … , 100 th ) of dust mixing-ratio at each gridbox (x, y, T°) and for each month Create new dimension "Dust level" by selecting the days for which the dust mixing-ratio falls between each pair of deciles (10) Aggregate into multiyear-monthly values (12) The result is a 5-dimensional dataset of dimensions: time(12)="Months" level(10)="Dust percentiles" y (96)="Latitude" x (12)="Longitude" z (15)="Temperature range"