Analysis of the cloud fraction adjustment to aerosols and its dependence on meteorological controls using explainable machine learning

Abstract

Aerosol–cloud interactions (ACI) have a pronounced influence on the Earth's radiation budget but continue to pose one of the most substantial uncertainties in the climate system. Marine boundary-layer clouds (MBLCs) are particularly important since they cover a large portion of the Earth's surface. One of the biggest challenges in quantifying ACI from observations lies in isolating adjustments of cloud fraction (CLF) to aerosol perturbations from the covariability and influence of the local meteorological conditions. In this study, this isolation is attempted using 9 years (2011–2019) of near-global daily satellite cloud products in combination with reanalysis data of meteorological parameters. With cloud-droplet number concentration ($N_{\mathrm{d}}$) as a proxy for aerosol, MBLC CLF is predicted by region-specific gradient boosting machine learning (ML) models. By means of SHapley Additive exPlanation (SHAP) regression values, CLF sensitivity to $N_{\mathrm{d}}$ and meteorological factors as well as meteorological influences on the $N_{\mathrm{d}}$–CLF sensitivity are quantified. The regional ML models are able to capture, on average, 45 % of the CLF variability. Based on our statistical approach, global patterns of CLF sensitivity suggest that CLF is positively associated with $N_{\mathrm{d}}$, particularly in the stratocumulus-to-cumulus transition regions and the Southern Hemispheric midlatitudes. However, $N_{\mathrm{d}}$ retrieval bias may contribute to non-causality in these positive sensitivities, and hence they should be considered upper-bound estimates. CLF sensitivity to estimated inversion strength (EIS) is ubiquitously positive and strongest in tropical and subtropical regions topped by stratocumulus and within the midlatitudes. Globally, increased sea-surface temperature (SST) reduces CLF, particularly in stratocumulus regions. The spatial patterns of CLF sensitivity to horizontal wind components in the free troposphere may point to the impact of synoptic-scale weather systems and vertical wind shear on MBLCs. The $N_{\mathrm{d}}$–CLF relationship is found to depend more on the selected thermodynamical variables than dynamical variables and in particular on EIS and SST. In the midlatitudes, a stronger inversion is found to amplify the $N_{\mathrm{d}}$–CLF relationship, while this is not observed in the stratocumulus regions. In the stratocumulus-to-cumulus transition regions, the $N_{\mathrm{d}}$–CLF sensitivity is found to be amplified by higher SSTs, potentially pointing to $N_{\mathrm{d}}$ more frequently delaying this transition in these conditions. The expected climatic changes in EIS and SST may thus influence future forcings from the CLF adjustment. The novel data-driven framework, whose limitations are also discussed, produces a quantification of the response of MBLC CLF to aerosols, taking into account the covariations with meteorology.

1 Introduction

The emission of aerosols into the atmosphere affects the Earth's climate in particular by masking part of the warming effect from greenhouse gases by reflecting solar radiation and changing cloud properties. Aerosol–cloud interactions (ACIs) can strongly influence the Earth's energy distribution and thus also contribute a substantial uncertainty to past and future climate projections. The effective radiative forcing due to ACI (ERF_aci) is assessed to be $-$1.0 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, with an uncertainty range of $-$1.7 to $-$0.3 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ albeit decades of effort and headway have been made in understanding the complex system of aerosols, clouds, and their environmental controls. The correct representation of ACI in Earth system models (ESMs) remains a tremendous challenge because of the lack of accurate global quantification of the cloud-related fine-scale processes and the lack of larger-scale constraints from the existing measurement systems at the ESM spatiotemporal resolution .

Marine boundary-layer clouds (MBLCs) cover over 23 % of the global ocean surface . Due to relatively small temperature differences between MBLC top and the sea surface, they only weakly impact outgoing longwave radiation but greatly reflect incoming shortwave radiation, leading to a strong net cooling effect . MBLCs play a critical role in the Earth's radiative balance and, in this regard, are the most important cloud type . Furthermore, MBLCs are especially susceptible to aerosol perturbations due to their relatively low optical depths and their formation in environments typically characterized by lower anthropogenic aerosol loading than continental clouds . Therefore, a deeper understanding of the aerosol–MBLC interactions is crucial to reduce the uncertainties in climate predictions. Atmospheric aerosols are critical for the formation of clouds as cloud condensation nuclei (CCN). Increases in aerosols are associated with increases in cloud-droplet number concentration ($N_{\mathrm{d}}$). As the cloud water is distributed among more droplets, cloud-droplet effective radius ($r_{\mathrm{e}}$) shrinks at constant liquid water content, resulting in an enhancement of cloud brightness and a negative instantaneous radiative forcing . The likelihood of collision and coalescence subsequently decreases due to smaller drop sizes, hampering rainfall formation, which can prolong cloud lifetime and thus increase cloud fraction (CLF) . However, the aerosol–CLF relationship is complex, and the sign of the CLF adjustment can also be the opposite. This has been found in particular for non-precipitating clouds, stemming from enhanced entrainment mixing with ambient air over the clouds owing to shorter evaporation timescales or reduced sedimentation because of smaller droplet sizes.

From the perspective of observations at satellite scales, though there are studies suggesting a negative relationship between aerosols and CLF , it has been documented by multiple studies that the overall CLF increases in response to increasing aerosols e.g.. Likewise, studies based on ESMs reported substantial negative ERF_aci due to liquid water path (LWP) and CLF adjustments e.g.. In spite of the attribution of such adjustments in ESMs primarily to LWP adjustments , a global satellite-based study by suggested that LWP adjustments are overestimated in ESMs and that aerosol impact on CLF dominates the negative aerosol forcing. This is supported by observational evidence presented by , who also reported an overestimation of LWP adjustment in climate models, and by , who recently highlighted the role of CLF increases due to aerosols from a large volcano eruption as the main cause of the associated forcing. Some large-eddy simulations have, however, suggested a negative response of CLF of trade wind cumulus to aerosol perturbations . While most studies, from both observational and model points of view, are in agreement that generally CLF increases with increasing aerosols due to a prolonged lifetime , the magnitude of the response of CLF to aerosols and its corresponding adjustments are still highly uncertain. For satellite-based analyses, one of the most challenging aspects in the quantification of CLF adjustment is isolating the influence of the aerosol loading on cloud properties from confounding covariations with meteorological parameters paired with aerosol retrieval issues related to aerosol swelling and 3D radiative effects in the vicinity of clouds . Recent observational studies have utilized different methods to tackle this issue. A first approach is to stratify the data by meteorological factors, therefore accounting for local meteorology in the relationships e.g.. Secondly, using $N_{\mathrm{d}}$ as a mediating variable was proposed by to analyse the causal pathway between aerosol optical depth and CLF. Another approach is to use a sampling strategy that applies a cloud–aerosol pairing algorithm . However, these methods do not account for aerosol retrieval issues, meteorological influencing factors, and confounders at once, which is essential to constrain the CLF adjustment. Recently, several studies have successfully used machine learning (ML) to account for non-linearities and meteorological factors to quantify ACI . ML regression algorithms allow for the prediction of CLF (predictand) on the basis of aerosol and meteorological factors at the same time and treat the aerosol–cloud–meteorology system as a whole. In addition, ML models can represent non-linear interactive systems, which can be analysed in sensitivity analyses with explainable ML techniques. Explainable ML refers to the techniques explaining the predictions of a trained ML model by explicitly quantifying the relationships, which helps improve the understandability, transparency, and trustworthiness of the ML models .

In this study, we set up region-specific ML models at a global scale using satellite and reanalysis data sets to predict CLF to analyse $N_{\mathrm{d}}$-induced changes in MBLCs. The goal of the explainable ML framework is to quantify the global sensitivity patterns of CLF to $N_{\mathrm{d}}$ and meteorological factors. In addition, we aim to estimate the magnitude of the dependence of $N_{\mathrm{d}}$–CLF sensitivity on the meteorological factors using SHapley Additive exPlanation (SHAP) interaction values, providing a new and insightful pathway to more profound knowledge of the physical processes relevant to the CLF adjustment and, hence, to a global constraint on aerosol-induced CLF changes accounting for meteorological covariations. The hypothesis of this study is that the response of cloud fraction of MBLCs to aerosol perturbations is positive but buffered, i.e. reduced or amplified, by ambient meteorology and that both the sensitivities and the interactions with meteorological factors have distinct regional patterns.

2 Data and methods

2.1 Data sets

This work combines 9 years (2011–2019) of satellite retrievals from Moderate Resolution Imaging Spectroradiometer (MODIS) and reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF) from 60° N to 60° S. In this study, MBLCs are defined as single-layer warm cloud fields with cloud top temperatures higher than 268 $\mathrm{K}$. To achieve this, the information on CLF (Cloud_Retrieval_Fraction_1L_Liquid product), $r_{\mathrm{e}}$ (Cloud_Effective_Radius_1L_Liquid_Mean product), cloud optical depth (${\mathit{\tau }}_{\mathrm{c}}$; Cloud_Optical_Thickness_1L_Liquid_Mean product), cloud top temperature (CTT; Cloud_Top_Temperature_Mean product), and satellite viewing geometry are obtained from MODIS level-3 collection-6.1 atmosphere daily products on the Terra platform (MOD08_D3), which are gridded into 1° $\times$ 1° globally from level-2 atmospheric products. CLF serves as the predictand in this study. The computation of $N_{\mathrm{d}}$ relies on ${\mathit{\tau }}_{\mathrm{c}}$ and $r_{\mathrm{e}}$, with filtering criteria based on CTT, solar zenith viewing angle, and satellite zenith angle, as elaborated in the following.

The equation used to calculate the MODIS $N_{\mathrm{d}}$ is from depends on the retrievals of $r_{\mathrm{e}}$ and ${\mathit{\tau }}_{\mathrm{c}}$ and so do the uncertainties in the errors propagated from $r_{\mathrm{e}}$ and ${\mathit{\tau }}_{\mathrm{c}}$:

1$$N_{\mathrm{d}}=\mathit{\alpha }{\mathit{\tau }}_{\mathrm{c}}^{\mathrm{0.5}}{r}_{\mathrm{e}}^{-\mathrm{2.5}},$$ where $\mathit{\alpha }=\mathrm{1.37}\times {\mathrm{10}}^{-\mathrm{5}}$ ${\mathrm{m}}^{-\mathrm{0.5}}$ is a constant related to adiabatic growth rate. The uncertainties in $N_{\mathrm{d}}$ retrievals are exhaustively evaluated by , which suggests that the uncertainties in averaged $N_{\mathrm{d}}$ over a 1° $\times$ 1° grid box (spatial resolution of the MODIS products used in this study) decrease by over 50 % compared to pixel-level uncertainties. This derivation approach relies on the assumed adiabaticity in global marine warm clouds where liquid water content and $r_{\mathrm{e}}$ increase monotonically and $N_{\mathrm{d}}$ is distributed as constant vertically. Departure from the adiabatic assumption (e.g. due to entrainment) would result in $N_{\mathrm{d}}$ retrieval biases . The uncertainty related to the estimation of $N_{\mathrm{d}}$ from MODIS also depends on liquid CLF. $N_{\mathrm{d}}$ is less biased in the regions of larger CLF, where clouds are more homogeneous, while in the regions with lower CLF $N_{\mathrm{d}}$ retrievals are sparser and less reliable . In such heterogeneous cloud fields, subpixel effects in the retrieval of $r_{\mathrm{e}}$ can negatively bias the retrieved $N_{\mathrm{d}}$ values . Such retrieval biases could cause a bias in the $N_{\mathrm{d}}$–CLF relationship as well. Furthermore, the interpretation of the causal effect of $N_{\mathrm{d}}$ on CLF can also be obscured by small-scale sampling issues. In particular, apart from the retrieval errors in $r_{\mathrm{e}}$ and ${\mathit{\tau }}_{\mathrm{c}}$, the natural spatial variability in cloud fields can also propagate to the $N_{\mathrm{d}}$ estimate and distort the $N_{\mathrm{d}}$–CLF relationship .

Following the screening criteria for a more reliable $N_{\mathrm{d}}$ demarcated by , only clouds restricted to a single layer in the liquid phase with a CTT higher than 268 $\mathrm{K}$ are considered. As suggested by , samples with $r_{\mathrm{e}}<\mathrm{4}$ $\mathrm{\mu }\mathrm{m}$ and ${\mathit{\tau }}_{\mathrm{c}}<\mathrm{4}$ are excluded to cope with the high $r_{\mathrm{e}}$ retrieval uncertainties at low ${\mathit{\tau }}_{\mathrm{c}}$. In addition, solar and sensor viewing zenith angles respectively greater than 65° and 55° are removed to avoid the large biases in $r_{\mathrm{e}}$ and ${\mathit{\tau }}_{\mathrm{c}}$ retrievals as in. The pixels selected according to the above sampling strategies generate more reliable $N_{\mathrm{d}}$ estimates.

Atmospheric and oceanic variables are taken from the fifth-generation ECMWF atmospheric reanalysis of the global climate (ERA5) at an hourly frequency (Table ) . The ERA5 data sets are harmonized to fit the level-3 MODIS data by first being resampled to 1° $\times$ 1° from their default 0.25° $\times$ 0.25° spatial resolution using bilinear interpolation; they are subsequently collocated with Terra MODIS by extracting hourly data to align with the UTC overpass times of the Terra satellite for each grid cell, yielding a spatiotemporally matched MODIS-ERA5 combined data set for training the ML models. For $N_{\mathrm{d}}$ retrievals, only samples within 1st–99th percentiles are retained to exclude potential unrealistic outliers from $r_{\mathrm{e}}$ and ${\mathit{\tau }}_{\mathrm{c}}$ retrievals . Furthermore, the explanation of ML models in this study relies on using linear regressions to capture the distribution of individual prediction instances, and the extreme values may excessively magnify or reduce the sensitivity or interactive effects quantified by SHAP (shown in Fig.  and discussed in Sect. ). The threshold of 1st–99th percentiles for each predictor is thus adopted to remove the values at the very tails of the specific distribution and to improve the robustness of the estimated sensitivities. To define the sensitivities of CLF and the interactive effects of meteorological factors, the natural logarithm of $N_{\mathrm{d}}$ is taken (see Sect.  in detail). Estimated inversion strength (EIS) is calculated based on the formulation from , and in this study, it is dependent only on atmospheric temperatures at 700 hPa and at the level of 1000 hPa.

All input predictors for each Extreme Gradient Boosting (XGB) model (i.e. for each 5° $\times$ 5° window aggregated from 1° $\times$ 1° grid boxes, as detailed in Sect. ) are standardized by centring around the mean and scaling to have unit variance as in . suggested that the standardization process is a standard practice when aiming for comparability of sensitivity estimates across predictors. This process eliminates the influence of units and aligns data on the same scale instead of the original natural ones, thereby ensuring the comparability of the quantified sensitivities and interactive effects with meteorology among different variables. This standardization procedure has been applied in other studies investigating different cloud sensitivities to various cloud-controlling factors e.g.. This procedure, however, may result in reduced spatial comparability due to variations in mean and standard deviation values across different 5° $\times$ 5° windows. To assess the trade-off between comparability among different predictors and comparability in space, we provide results without standardization in the Supplement (Figs. S2 to S7 therein) as done by . In terms of spatial patterns, the results are nearly identical to their corresponding ones presented in the following sections of the main text, suggesting that standardizing the data based on the local mean and standard deviation for each window has only a small impact on comparability across each window. Therefore, we primarily benefit from achieving comparability among different predictors while making only a minor compromise in spatial comparability.

Table 1

Summary of the predictors from ERA5 reanalysis.

Predictor name	Abbreviation	Unit
Instantaneous pressure-level parameters (at 700 $\mathrm{hPa}$, 850 $\mathrm{hPa}$)
Relative humidity	RH700, RH850	%
Specific humidity	SH700, SH850	$\mathrm{kg}{\mathrm{kg}}^{-\mathrm{1}}$
Temperature	$t_{\mathrm{700}}$, $t_{\mathrm{850}}$	$\mathrm{K}$
Vertical velocity	${\mathit{\omega }}_{\mathrm{700}}$, ${\mathit{\omega }}_{\mathrm{850}}$	$\mathrm{Pa}{\mathrm{s}}^{-\mathrm{1}}$
Eastward wind component	$u_{\mathrm{700}}$, $u_{\mathrm{850}}$	$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$
Northward wind component	$v_{\mathrm{700}}$, $v_{\mathrm{850}}$	$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$
Surface and single-level parameters (instantaneous or mean rates/fluxes)
Eastward and northward wind component at 10 m	$u_{\mathrm{10}}$, $v_{\mathrm{10}}$	$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$
Boundary-layer height	BLH	$\mathrm{m}$
Convective available potential energy	CAPE	$\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}$
Sea surface temperature	SST	$\mathrm{K}$
Total column water vapour	TCWV	$\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$
Mean large-scale precipitation fraction	PF	proportion
Mean surface sensible/latent heat flux	SHF/LHF	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
Calculated
Estimated inversion strength	EIS	$\mathrm{K}$

2.2 Machine learning model setup

Extreme Gradient Boosting (XGB) is a distributed tree boosting algorithm aiming to provide a scalable, portable, and flexible library under the gradient boosting framework . The state-of-the-art XGB algorithm can be implemented efficiently in Python and has been recently used to study clouds and ACI . As an extension of previous gradient boosting methods, XGB has incorporated regularization techniques which help prevent overfitting and improve model generalization. Besides, the subsampling on training subsets and column (feature) subsampling techniques can shorten the running time and also avert overfitting and hence elevate model performance . Relevant regularization and subsampling hyperparameters are tuned using Bayesian optimization to determine the best combination; see Table  for the search space.

Table 2

Overview of the hyperparameters tuned for regional Extreme Gradient Boosting models using Bayesian optimization.

Hyperparameter name	Search space
learning_rate	0.01–0.5
max_depth	3–10
min_child_weight	1–10
subsample	0.5–1
colsample_bytree	0.5–1
gamma	0–10
alpha	0–10
lambda	0–10

Data from 2011 to 2016 are used for training and data from 2017 to 2019 for testing (ratio of independent train to test split of about 67 % / 33 %). By chronologically splitting the training and test sets without random shuffling, we ensure that the training data does not see future information and the autocorrelation in data does not lead to overoptimistic evaluation of the model's performance . As suggested by , a single ML model may not perform well across all regions due to the heterogeneity of relevant processes. Therefore, data at a 1° $\times$ 1° spatial resolution are aggregated into 5° $\times$ 5° geographical windows, where an individual independent XGB model is trained and tested for each “window”. Hereby, a region-specific ML framework is established to potentially capture regional relationships and characteristics and thus the regional patterns of CLF adjustment. The coarser 5° $\times$ 5° spatial resolution of the modelling grid increases the sample size by a factor of $\approx$ 25, which is helpful to establish robust sensitivity estimates. In addition, at the spatial resolution of 1° $\times$ 1° summarized in 5° $\times$ 5° windows, the spatial scale is adequate for ACI sensitivity estimation . To ensure a sufficient data amount for training and testing the XGB models, only the geographical windows with over 6000 available data points are retained. Consequently, 34 out of 1190 oceanic windows have been excluded. These windows located between 47.5° W and 122.5° E and 52.5 and 57.5° S in the Southern Ocean (Fig. ) contain fewer than 6000 valid samples due to the screening for $N_{\mathrm{d}}$ retrievals. For each model, the hyperparameters are tuned by implementing Bayesian optimization, which uses a Gaussian process prior distribution over hyperparameters to initialize a probabilistic model for the objective function to be optimized. After the initialization, the probabilistic model is updated iteratively, and Bayesian optimization suggests the optimal combination of hyperparameters to try for the next iteration according to the previous one and samples gathered from the search space (Table ) . Each iteration is evaluated by five-fold cross-validation using the root mean square error (RMSE) as score. The number of boosting rounds (the number of trees) for each XGB model is then determined by the early stopping technique to further avoid overfitting; i.e. the training of the model stops early once it is monitored, so the score of cross-validation does not improve within 20 iteration rounds.

2.3 Explaining the machine learning models

2.3.1 SHapley Additive exPlanation (SHAP) values

SHAP values were proposed by on the basis of cooperative game theory to explain the outputs of ML models. The SHAP approach has been implemented with XGB in Python, and it has been reported that outputs from XGB models with various number of trees can be well explained by the SHAP framework in different subject areas e.g.. The contribution of a predictor value to a specific model prediction is calculated as the difference between the predictions of the model in the presence and absence of this particular predictor for all possible combinations of predictor values. Since this is performed at a “local” level (i.e. for this specific instance's prediction), it allows for insights into how a certain model outcome is achieved, thereby complementing more traditional “global” (considering all instances) feature importance measures (e.g. partial dependence plot).

The base value in the context of SHAP values is what would be predicted in the absence of any feature information , and it is typically computed as the average of all predictions by ML models over the entire training data set. Positive (negative) SHAP values indicate that the specific feature value increases (decreases) the prediction compared to this base value. In other words, the base value serves as the reference point against which the contributions of individual features are measured. SHAP values for all features always sum up to the difference between the base value and the final model prediction so that SHAP values are additive and internally consistent. The base value could be analogous to the climatological CLF for a given geographical window, assuming no information about the input parameters is known. In this context, the SHAP values of input features indicate the extent to which knowing information about each feature value would deviate the prediction from the climatological CLF (base value).

Furthermore, the quantification of the influence of meteorology on the $N_{\mathrm{d}}$–CLF relationship can be analysed using SHAP interaction values, which are an extension of SHAP values. They measure the difference between the SHAP values for a feature when another (secondary) feature is included versus when it is not included, offering a potential tool for insights into feature interactions captured by the tree ensembles. SHAP values have been applied to study atmospheric aerosols in the context of air pollution and have been used by to explore satellite-observed $N_{\mathrm{d}}$–LWP relationship in MBLCs in the southeast Atlantic, finding that meteorological variables have considerable influences on the $N_{\mathrm{d}}$–LWP relationship using SHAP interactive values. Moreover, the use of SHAP interaction values in these studies allows for a more profound and in-depth comprehension of the underlying processes with respect to local meteorology. SHAP values provide insights into the behaviour of the XGB models, and as all statistical/ML models, they may not necessarily reflect real-world physical causality. Nevertheless, this state-of-the-art technique allows us to account for meteorological covariations when deriving sensitivities and to appraise to what extent the meteorological predictors interact with and influence the $N_{\mathrm{d}}$–CLF relationship beyond traditional global-level feature attributions.

2.3.2 Quantification of sensitivities and interactive effects

Figure  is an exemplary graph for a regional XGB model at a specific 5° $\times$ 5° window (27.5–32.5° S, 122.5–127.5° W). SHAP values and SHAP interaction values are used to explain this XGB model and to quantify and isolate the CLF sensitivity to $N_{\mathrm{d}}$ and the interactive effects of meteorological factors (here sea-surface temperature, SST). Each dot in Fig.  represents an individual data instance (i.e. a single observation at a specific grid cell and time step) and shows how individual $N_{\mathrm{d}}$ or ln$N_{\mathrm{d}}$ values impact the CLF prediction.

Plotting SHAP values of $N_{\mathrm{d}}$ against $N_{\mathrm{d}}$ values without the standardization process (Fig. a) for each data sample illustrates that increased $N_{\mathrm{d}}$ values lead to an increase in the predicted CLF, while the rate of the increase (dSHAP $/$ $\mathrm{d}{N}_{\mathrm{d}}$) drops with $N_{\mathrm{d}}$ as shown by the orange line. For each 20 ${\mathrm{cm}}^{-\mathrm{3}}$ wide bin of $N_{\mathrm{d}}$, dSHAP $/$ $\mathrm{d}{N}_{\mathrm{d}}$ is calculated as the slope of the linear regression between $N_{\mathrm{d}}$ and $N_{\mathrm{d}}$ SHAP values. The non-linear positive association between $N_{\mathrm{d}}$ and predicted CLF aligns well with findings of prior studies e.g. that the aerosol impact on CLF saturates at relatively high aerosol loading. This relationship also resembles the one reported by , which is attributed to the precipitation suppression effect due to a relatively high $N_{\mathrm{d}}$.

Expressing the sensitivity logarithmically in $N_{\mathrm{d}}$ is ideal because cloud processes are prone to respond to a relative change in $N_{\mathrm{d}}$ rather than an absolute one . Furthermore, the log-transformed $N_{\mathrm{d}}$ facilitates the application of simple linear regressions to capture the relationship between the contribution of $N_{\mathrm{d}}$ and the predicted CLF ($N_{\mathrm{d}}$ SHAP values) and its feature values. As depicted in Fig. b, the contribution of ln$N_{\mathrm{d}}$ to the predicted CLF increases almost linearly with a rising ln$N_{\mathrm{d}}$. Thus, the CLF sensitivity to $N_{\mathrm{d}}$ is estimated as the slope of the linear regression between ln$N_{\mathrm{d}}$ SHAP values and ln$N_{\mathrm{d}}$ values (0.098 CLF ${\mathit{\sigma }}^{-\mathrm{1}}$). A similar method to estimate sensitivity has also been used by , where it is also suggested that this method can enhance the robustness of the sensitivity estimation. Because it can leverage the benefits of an XGB model, including bagging techniques and no need for distribution assumptions, along with the advantages of SHAP, which provides global interpretations consistent with local explanations . It should be noted that the notably linear relationship in Fig. b does not hold across all geographical windows. Figure S1 displays additional exemplary windows where the relationships exhibit less linearity. Our approach also captures non-linearity in the system; in these cases, the linear regression helps decrease the convolved relationships as in . Note that unlike $N_{\mathrm{d}}$ (${\mathrm{cm}}^{-\mathrm{3}}$) in panel (a), ln$N_{\mathrm{d}}$ and SST in (b) and (c) have been standardized, and thus sensitivities and interaction indices (IAIs) are expressed with the unit of cloud fraction change per standard deviation (CLF ${\mathit{\sigma }}^{-\mathrm{1}}$). Standardizing all predictors ensures that the results become comparable across all of them. We also present the SHAP dependence plots for the same example window in Fig. S2 where non-standardized ln$N_{\mathrm{d}}$ and SST are used to plot panels (b) and (c). The patterns are alike and only the magnitudes of the example sensitivity and IAI are different because they are no longer expressed on a physical scale.

The vertical dispersion around the ln$N_{\mathrm{d}}$–CLF relationship captured by the SHAP dependence plot is due to the dependence of the ln$N_{\mathrm{d}}$ contribution to the predicted CLF on meteorological factors (e.g. SST) in the model, which is captured by SHAP interaction values, as displayed in Fig. c. The colouring of the data points by SST illustrates how interactions with SST split up the ln$N_{\mathrm{d}}$–CLF relationship, with low SST values amplifying the ln$N_{\mathrm{d}}$ contribution and vice versa. To quantify this interaction effect, the meteorological data are then divided into a group of above-average feature values and a group of below-average feature values. A linear regression is fit to the ln$N_{\mathrm{d}}$ values and the SHAP interaction values in each group. An interaction index (IAI) is derived from these regression fits and defined as the slope for the high-value group ($>$ mean) with the slope for the low-value group ($<$ mean) subtracted:

2$$\mathrm{IAI}={\mathit{\beta }}_{x,\mathrm{high}}-{\mathit{\beta }}_{x,\mathrm{low}},$$ where $\mathit{\beta }$ is the slope of the linear regression between SHAP interaction values and ln$N_{\mathrm{d}}$ values and the subscripts denote the high-value group and the low-value group for a specific meteorological variable $x$ (SST in the example) respectively. At the exemplary geographical window, the influence of SST on the $N_{\mathrm{d}}$–CLF sensitivity is quantified by IAI $=-\mathrm{0.029}$ CLF ${\mathit{\sigma }}^{-\mathrm{1}}$ (Fig. c). Similar to sensitivities, the unit of IAIs is also CLF ${\mathit{\sigma }}^{-\mathrm{1}}$. Therefore, for a positive sensitivity such as the $N_{\mathrm{d}}$–CLF sensitivity shown in Fig. b, a negative IAI value means that the $N_{\mathrm{d}}$–CLF sensitivity is larger with low feature values, as shown in Fig. c (the positive relationship is weakened by high SST values). On the contrary, a positive IAI value corresponds to a larger positive sensitivity with high feature values.

Figure 1

SHAP dependence plots for the cloud-droplet number concentration ($N_{\mathrm{d}}$) in the region from 27.5 to 32.5° S and from 122.5 to 127.5° W. (a) Dots show $N_{\mathrm{d}}$ SHAP values versus $N_{\mathrm{d}}$ values. The orange line shows the change rate of $N_{\mathrm{d}}$ SHAP values with respect to $N_{\mathrm{d}}$ (dSHAP $/$ $\mathrm{d}{N}_{\mathrm{d}}$) versus $N_{\mathrm{d}}$ values for each $N_{\mathrm{d}}$ bin of 20 ${\mathrm{cm}}^{-\mathrm{3}}$ wide. Panel (b) is similar to panel (a) but shows the relationship between ln$N_{\mathrm{d}}$ SHAP values and ln$N_{\mathrm{d}}$ with the corresponding sensitivity defined as the slope of the linear regression. Panel (c) shows SHAP interaction values coloured by sea-surface temperature (SST) showing the dependence of ln$N_{\mathrm{d}}$–CLF relationship on the interactive effects of SST. The interaction values are further divided into two groups by the mean feature value of SST. Linear regressions are performed respectively for the high-value group and low-value group and the interaction index (IAI) is defined as the slope for the high-value group by subtracting the slope for the low-value group. The horizontal dashed lines are a demarcation between negative and positive SHAP (interaction) values. Note that $N_{\mathrm{d}}$ in (a) is not standardized, while ln$N_{\mathrm{d}}$ and SST in (b) and (c) are standardized.

[Figure omitted. See PDF]

2.3.3 Limitations of observation-based machine learning of aerosol-cloud processes

In this section, limitations of this study are discussed. A fundamental limitation of our study is that the assertion of causality from the statistical relationships of aerosols/$N_{\mathrm{d}}$ and cloud fraction/properties is not easily done. While causal inference approaches exist and have been applied in the field of aerosol–cloud interactions , we employ a more traditional approach of analysing statistical relationships of instantaneous observations (i.e. correlations). Unless nonetheless explicitly incorporating such causal inference approaches, studies utilizing statistical or ML models to explore observational aerosol–cloud processes contend with this common limitation. For instance, some studies assessed satellite-based statistical relationships between CLF and $N_{\mathrm{d}}$ , between LWP and $N_{\mathrm{d}}$ , and between $N_{\mathrm{d}}$ and other aerosol proxies , all resting on statistically inferring sensitivities of cloud quantities to aerosol proxies . While we interpret the derived relationships with respect to the known physical relationships, uncertainties regarding the physical interpretation are mainly driven by two sources: uncertainties in the data and uncertainties from the methods.

• 1.

  Data. Uncertainties exist for each satellite/reanalysis quantity, but may be particularly large in $N_{\mathrm{d}}$. For example, the subpixel effect can introduce more bias in the $N_{\mathrm{d}}$ retrieval process within broken-cloud regimes due to increased heterogeneity. The $N_{\mathrm{d}}$ retrieval biases are discussed in Sect. . Also, $N_{\mathrm{d}}$ and CLF observations are not fully independent, which may introduce a spurious positive correlation between the two variables. As such, we expect the physical relationship of $N_{\mathrm{d}}$ and CLF to be weaker than our estimate so that the derived sensitivities present an upper bound of the physical relationship.

  Another caveat in our data is that $N_{\mathrm{d}}$ values in our study are computed using MODIS level-3 large-scale mean $r_{\mathrm{e}}$ and ${\mathit{\tau }}_{\mathrm{c}}$ values instead of joint histograms as in . This may introduce additional biases considering the non-linearity of the $N_{\mathrm{d}}$ calculation. In future work, $N_{\mathrm{d}}$ data calculated from underlying joint histograms or pre-filtered data by could be applied to be compared with the results in this study.

• 2.

  Methods.

  • a.

    The exact quantification of sensitivities is dependent on the choice of the statistical/machine learning model. While for (more linearly related) monthly data, have shown that XGB, artificial neural networks, and linear models tend to lead to very similar results, this is not expected for more instantaneous data. Here, non-linear relationships are expected, and a more complex non-linear model is a more appropriate choice. XGB and other tree ensemble methods are a particularly popular choice because of their interpretability, high accuracy considering computational efficiency , and ability to model the interactions between predictors . They have been frequently used to study aerosols and clouds in the past . Besides, the Tree SHAP algorithm, specifically tailored for tree-based models to compute exact Shapley values, can even further enhance their interpretability and has been applied in this field as well .

  • b.

    The quantification of sensitivities with SHAP values depends on details: the choice of the algorithm to effectively estimate Shapley values is application-specific and comes to the trade-off between being true to the data and true to the model, which relies on an observational and interventional conditional expectation respectively . The true to the model approach is preferable when trying to understand how an ML model makes a prediction, which requires assuming feature independence. In this study, we focus on potential mechanisms behind CLF sensitivities, and thus we tend to respect the correlations spread among input features (true to the data) . Consequently, we suffer from the disadvantage of being true to the data: entangled importance attributions of correlated features, e.g. a feature not explicitly used by the model for the prediction task, might be assigned a non-zero contribution. Yet we refrain from the drawback of being true to the model – unrealistic input instances . Despite the inherent trade-off, SHAP approach has been employed in the context of being true to the data e.g..

The derived estimates of sensitivities and interactive effects in this paper should thus be interpreted with these limitations and uncertainties in mind.

3 Results and discussion

3.1 Model performance

The skills of the region-specific XGB models in predicting CLF are evaluated by the coefficient of determination ($R^{\mathrm{2}}$) on the unseen hold-out test data. The global weighted mean $R^{\mathrm{2}}$ is 0.45 (about 45 % on weighted average and up to 73.57 % of the variability in CLF prediction is explained) and the standard deviation 0.10. While this means that, on average, about half of the variability in CLF cannot be explained by the machine learning models, this is expected as previous studies have shown that the performance of statistical models decreases when going from monthly to daily data , and the performance is on par with that reported by , who used machine learning models to predict $N_{\mathrm{d}}$ with daily reanalysis data. The models in tropical regions in the Indian Ocean and the western Pacific relatively poorly explain the variability in CLF, while XGB models perform well in the stratocumulus regions in the subtropics near the continents and in the midlatitudes, particularly the Southern Hemispheric midlatitudes. The high skill of predicting CLF in the Southern Hemispheric midlatitudes is in contrast to a recent study where this region has been found to be particularly difficult to model statistically with monthly data . In this region, the day-to-day CLF variability is high due to the large influence of synoptic-scale weather systems, and hence data at the daily resolution are more adequate to represent the CLF variability in these regions.

Figure 2

$R^{\mathrm{2}}$ score of regional Extreme Gradient Boosting models predicting the cloud fraction of marine boundary-layer clouds in the independent test data set (2017–2019).

[Figure omitted. See PDF]

3.2 CLF sensitivity: global perspectives and regional characteristics

3.2.1 Global overview of CLF sensitivities

Figure  summarizes the means and distributions of the near-global sensitivities of CLF to all predictors. The sensitivities are estimated as described in Sect. . The sequence is sorted by descending mean values of the absolute sensitivities (i.e. by feature importance) of the predictor variables. A strong and consistently positive $N_{\mathrm{d}}$–CLF sensitivity is found. The fact that CLF is the most sensitive to $N_{\mathrm{d}}$ is to be expected, as cloud observations from the same sensor are more directly related than a reanalysis product, so their overall magnitude should not be compared . The entrainment of relatively dry air from the free troposphere into the MBL is impeded by a stronger inversion (i.e. higher EIS), resulting in a shallower, better-mixed, and more humid MBL conducive to stratocumulus clouds . The salient positive sensitivity to EIS is in accordance with the links found in previous studies e.g., suggesting that EIS is a crucial controlling factor for low marine cloud cover. Note that in some studies, the strength of the inversion over the boundary layer is measured by lower tropospheric stability, which can be regarded as a similar metric outperformed by EIS . Precipitation fraction is the fraction of the original ERA5 grid box covered by large-scale precipitation. The strong positive CLF sensitivity to precipitation fraction is likely caused by the ML model learning that precipitation can be viewed as a proxy for cloudiness rather than being an indicator of the physical processes via which precipitation exerts controls on the macrophysics of MBLCs. Humidity shows positive CLF sensitivities greater at 850 hPa, where cloud tops are often located , than at 700 hPa, which is typically in the free troposphere above the MBLCs . Likewise, the atmospheric temperature at 850 hPa ($t_{\mathrm{850}}$) presents stronger CLF sensitivity than the temperature at 700 hPa ($t_{\mathrm{700}}$). Nonetheless, in the case of winds the 700 hPa pressure level is more relevant than that at 850 hPa. A relatively pronounced negative sensitivity to the eastward wind component at 700 hPa ($u_{\mathrm{700}}$) seems to indicate that clouds are depleted due to more westerlies at this level. CLF exhibits negative sensitivities to vertical pressure velocities at both 850 and 700 hPa, showing that large-scale ascending motion is connected to increases in MBLCs . In general, the global averages of CLF sensitivity in terms of dynamical predictors (i.e. 3D winds at surface and pressure levels) vary in sign and are less strong. A marked negative sensitivity of CLF to SST is found, which is in agreement with many prior studies e.g., where increases in SST have been found to lead to low cloud breakup and dissipation due to a number of processes as described in, for example, . One of these is that the associated enhancement of mean surface latent heat flux (LHF) deepens MBL and facilitates buoyancy and thus the entrainment of dry free-tropospheric air . However, CLF is much less sensitive to LHF than to SST, which may indicate that this mechanism is less important at the spatial scale and timescale considered in this study. CLF exhibits a considerable negative sensitivity to mean surface sensible heat flux (SHF), which quantifies an increase in CLF with increasing SHF (upward SHF is negative). While increased SHF can promote the transition from decks of stratus or stratocumulus clouds (high CLF) to more convective clouds (low CLF) due to the deepening of the boundary layer , potentially leading to a positive SHF–CLF relationship, increased SHF is associated with situations of cold air advection where turbulent surface fluxes are enhanced, which could lead to marked increases in CLF .

Figure 3

The distribution of the sensitivities of the cloud fraction to all predictors as depicted in Table . Boxes represent the interquartile range, which is extended by whiskers to up to 1.5 interquartile ranges, with outliers shown as points outside the range. The solid line and white dot in each box show the median and mean values of the sensitivities respectively. Predictors are sorted by the mean values of absolute sensitivity values. The dashed line across the figure separates positive and negative sensitivity values.

[Figure omitted. See PDF]

3.2.2 Spatial patterns of the CLF sensitivity to $N_{\mathrm{d}}$

The sensitivity of the MBLC fraction associated with the aerosol proxy, $N_{\mathrm{d}}$, is ubiquitously positive in accordance with the global correlations or sensitivities found in, for example, and . This is presumably due to the lifetime effect but could also partially result from $N_{\mathrm{d}}$ retrieval biases discussed in Sect. . The global weighted mean value of the $N_{\mathrm{d}}$–CLF sensitivity is 0.074 CLF ${\mathit{\sigma }}^{-\mathrm{1}}$, with a standard deviation of 0.036 CLF ${\mathit{\sigma }}^{-\mathrm{1}}$. The relationship between CLF and $N_{\mathrm{d}}$ is found to be particularly strong in the regions of frequent stratocumulus-to-cumulus transition off the western continental coasts. These marked positive $N_{\mathrm{d}}$–CLF sensitivities may be caused by high $N_{\mathrm{d}}$, delaying the transition from stratocumulus to cumulus clouds . However, as this cloud regime transition involves clouds shifting from more overcast to more broken, the strong relationships in these regions may be more affected by $N_{\mathrm{d}}$ retrieval errors. The $N_{\mathrm{d}}$–CLF sensitivity is also pronounced in the Southern Hemispheric midlatitudes, where stratiform clouds dominate. The $N_{\mathrm{d}}$–CLF sensitivity is weak and close to zero in the tropics, in particular in the deep convective warm-pool region. These spatial patterns of $N_{\mathrm{d}}$–CLF sensitivity resemble those found by , in particular the ones where they mediated the aerosol optical depth–CLF relationship by $N_{\mathrm{d}}$ but are more pronounced in the Southern Hemispheric midlatitudes. This difference in estimated sensitivity seems noteworthy and should thus be investigated in future work. As $N_{\mathrm{d}}$ retrievals tend to negatively bias at lower CLF and positively bias at higher CLF, the $N_{\mathrm{d}}$–CLF sensitivity may be overestimated and, at the scales considered here, should be interpreted as an upper bound to the physical $N_{\mathrm{d}}$–CLF sensitivity. The global weighted average of the CLF–ln$N_{\mathrm{d}}$ sensitivity without standardization is 0.112 (unitless), and its spatial pattern is shown in Fig. S4. This value is higher than the upper bound of 0.1 reported by , which is based on global climate models and large-eddy simulations. This may be partly due to the aforementioned bias. However, it is important to note that our non-standardized CLF–$N_{\mathrm{d}}$ sensitivity, shown in Fig. a, closely mirrors that from , with a similar range. In addition, the high lnCLF–ln$N_{\mathrm{d}}$ values estimated in and suggest that values exceeding the upper bound of 0.1 might be plausible. These recent observational studies, including quantifying cloud fraction adjustment based on ship tracks , volcano aerosol perturbations , and our SHAP approach using global satellite observations, indicate that the 0.1 upper bound may be extended. In future work, estimating a radiative forcing using the SHAP-based sensitivities will make our study more comparable with other research on cloud fraction adjustment.

Figure 4

Sensitivity of the marine boundary-layer cloud fraction to ln$N_{\mathrm{d}}$.

[Figure omitted. See PDF]

3.2.3 Spatial patterns of the CLF sensitivity to thermodynamical drivers

There has been a strong consensus that EIS and SST are the two important determinants of cloud fraction of marine boundary clouds and their corresponding radiative effects across different geographical regions and on varying timescales e.g. . Stronger inversions capping MBL (i.e. higher EIS) will hamper the entrainment of aloft dry air from the troposphere and thus lead to a shallower MBL and more moisture trapped within MBL, promoting the development and maintenance of low-level clouds . The regional EIS–CLF sensitivity patterns (Fig. a) show that low marine cloud fraction increases ubiquitously in response to stronger EIS, in particular in the tropical and subtropical stratocumulus-capped regions and within the midlatitudes. The sensitivity pattern is in good agreement with that found by and , related studies at different timescales .

MBLC cover reduces globally in response to increased SST, particularly pronounced in the stratocumulus regions over eastern oceanic basins (Fig. b), consistent well with . SST can favour MBLC dissipation through increasing surface latent heat fluxes and deepening MBL, facilitating dry entrainment and eventually desiccating the MBL and clouds . Yet as stated in Sect. , the weak CLF sensitivity to LHF in relation to the strong sensitivity to SST may imply that the other process makes more substantial contributions – namely, that the higher moisture gradient between the troposphere and MBL arising from the increased SST makes the entrained air more efficient in evaporating cloud water . This process has been shown to be the driving mechanism for the observed reduction in marine low cloud cover near the coast of Baja California .

Figure c shows that low marine cloud fraction increases with negative (upward) SHF most markedly in the stratocumulus regions. CLF can increase in response to increased surface fluxes in situations of cold advection . Over the south Indian Ocean, a marked SHF–CLF sensitivity is also found. Here, enhancements of SHF due to the subtropical anticyclone and midlatitude storm-track activity have been found to increase CLF . The results may be a hint that the increase in CLF presumably due to increased SHF (e.g. due to cold advection) outweighs the influence of SHF on CLF by controlling the transition from marine stratocumulus to open-cellular marine clouds in the core stratocumulus regions. Consequently, the SHF–CLF sensitivity is less pronounced in regions of frequent closed- to open-cell and cumulus transitions. Relative humidity at 850 hPa (RH₈₅₀) is positively related to marine low liquid cloud fraction across the globe. The positive sensitivity is particularly strong in the trade cumulus regions, where the 850 hPa level is representative of the boundary layer. In the coastal stratocumulus regions, clouds are frequently below this level , so that clouds are not as sensitive to variability in RH at that level.

Figure 5

Sensitivity of the marine boundary-layer cloud fraction to the estimated inversion strength (EIS), sea-surface temperature (SST), sensible heat flux (SHF), and relative humidity at 850 hPa (RH₈₅₀). Note that the range of colour bars of SHF and RH₈₅₀ ($-$0.075 to 0.075) is narrower than EIS and SST ($-$0.15 to 0.15).

[Figure omitted. See PDF]

3.2.4 Spatial patterns of the CLF sensitivity to dynamical drivers

Large-scale circulations and dynamical conditions play an essential role in controlling cloud fraction and the indirect effects of aerosols . The large-scale dynamics are represented by the horizontal and vertical winds at 700 and 850 hPa, which display clear and distinct regional patterns of CLF sensitivity (Fig. ). It can also be seen that at the considered scales and pressure levels, horizontal wind vectors have stronger CLF sensitivities than large-scale vertical motion. There is a coherent pattern of negative CLF sensitivity to the zonal wind at 700 hPa in the stratocumulus-dominated regions (also apparent at 850 hPa), and the Southern Hemispheric midlatitudes, indicating a decrease in MBLCs with westerly anomalies at this pressure level. Recently, a study using monthly data has also found a similar sensitivity pattern of stratocumulus clouds to zonal wind at 700 hPa, finding that the reduced CLF is related to increased vertical wind shear (as the boundary-layer flow is easterly), leading to increased turbulence and dry-air entrainment . However, using monthly data, did not find a similar CLF sensitivity to zonal winds in the Southern Hemispheric midlatitudes. As the CLF sensitivity to $u_{\mathrm{700}}$ in the Southern Hemispheric midlatitudes is only apparent using daily data and only at 700 hPa, it seems likely that it is related to synoptic variability that drives day-to-day variability in MBLCs in this region . Positive CLF sensitivities to $u_{\mathrm{700}}$ (higher CLF with westerly anomalies) and, to a lesser degree, $u_{\mathrm{850}}$ are found off the eastern Asian and North American continents. CLF increases due to cold-air outbreaks in NW Atlantic and NW Pacific may be the reason for these positive sensitivities. Cold-air outbreaks occur during winter as cold continental air moves over warmer SSTs, increasing moisture and heat fluxes into the MBL so that the formation of MBLCs is favoured . This leads to wintertime maxima in CLF in these regions .

The sensitivity of CLF to the meridional winds at 700 hPa exhibits two bands straddling the subtropical regions between about 15 and 35° in both hemispheres but opposite in sign (positive in the Northern Hemisphere and negative in the Southern Hemisphere), illustrating that in these regions, the poleward winds are associated with an increase in low cloud fraction. The bands are still apparent at 850 hPa, while the negative band in the Southern Hemisphere extends northward to tropical areas. These hemispheric sensitivity bands to the $v$ wind component at 700 hPa closely resemble those found in , with their analysis suggesting that the poleward winds on the eastern side of midlatitude cyclones may be related to warm and moist advection, increasing CLF. However, they also find a strong correlation of these free-tropospheric poleward winds with large-scale ascending air motion making the assertion of causality difficult. Poleward winds are also found to decrease CLF over the Southern Hemispheric midlatitudes.

CLF is negatively connected to the vertical pressure velocity at both 700 and 850 hPa (${\mathit{\omega }}_{\mathrm{700}}$ and ${\mathit{\omega }}_{\mathrm{850}}$) over the entire Earth, indicating that ascending large-scale air motion enhances the cover of MBLCs globally. It is shown in the bottom of Fig.  column (a) that the CLF sensitivity to ${\mathit{\omega }}_{\mathrm{700}}$ is larger in the midlatitude ocean basins, whereas the CLF sensitivity to ${\mathit{\omega }}_{\mathrm{850}}$ is larger in the subtropical oceans, where subsidence is climatologically prevalent . This seems indicative of CLF being the most sensitive to large-scale ascending motion at the typical altitude of the clouds. It is interesting to note that between 30° N and 30° S, no marked CLF sensitivity to ${\mathit{\omega }}_{\mathrm{700}}$ is found, contrasting the finding of enhanced subsidence at this level reducing MBLCs by . This effect is likely better described in the ${\mathit{\omega }}_{\mathrm{850}}$ data, which is more related to the altitude of the cloud top.

Figure 6

Sensitivity of cloud fraction to wind component vectors $\mathit{u}$ and $\mathit{v}$ and vertical velocities at 700 hPa (column a) and 850 hPa (column b). Note that the range of the colour bars is in general smaller ($-$0.04–0.04) than in Fig. .

[Figure omitted. See PDF]

3.3 Dependence of $N_{\mathrm{d}}$–CLF relationship on meteorology

3.3.1 Global overview of the interaction indices

In this section, we use the IAI as defined in Sect.  to quantitatively show how the response of MBLC fraction attributed to the aerosol proxy $N_{\mathrm{d}}$ varies with the meteorological factors. As discussed in Sect. , since the sensitivity related to $N_{\mathrm{d}}$ is positive across the globe (Fig. d), a positive IAI can be interpreted as an amplification of the $N_{\mathrm{d}}$–CLF sensitivity with high (above-average) feature values of a meteorological variable, whereas a negative IAI signifies an amplification of the sensitivity at low feature values.

In Fig. , analogous to Fig. , the features along the $x$ axis are arranged in descending order based on their averaged absolute IAIs, that is, by the strength of the impact of each meteorological feature on the $N_{\mathrm{d}}$–CLF sensitivity. Similar to the feature importance summarized by Fig. , EIS, SST, RH₈₅₀, and SHF have relatively large strength of interaction effect and can thus be regarded as critical controlling factors for not only marine low cloud cover but also their response to changes in $N_{\mathrm{d}}$ (and in extension aerosols). Compared to the CLF sensitivities, the IAIs associated with atmospheric temperatures at 700 and 850 hPa have greater strengths. Furthermore, it can also be seen that the vertical and horizontal winds at the surface and different pressure levels are generally ranked lower. In general, the thermodynamical factors seem to have a stronger influence on the $N_{\mathrm{d}}$–CLF sensitivity than the dynamical factors.

Figure 7

Similar to Fig.  but for the interaction effect of $N_{\mathrm{d}}$ with all environmental parameters, quantified by the interaction index (CLF ${\mathit{\sigma }}^{-\mathrm{1}}$).

[Figure omitted. See PDF]

3.3.2 Spatial patterns of the interaction indices

Coherent and distinct spatial distributions of the impact of selected meteorological parameters on the $N_{\mathrm{d}}$–CLF relationship can be observed. Hereafter, we show the regional characteristics of the interaction effects of EIS and SST, which are the two most important meteorological factors for CLF in MBLCs and have the greatest absolute strengths of IAI. EIS exerts the most noticeable positive IAIs over the midlatitude oceanic areas (Fig. a), reflecting that stronger temperature inversions capping the MBL over these regions may amplify the positive $N_{\mathrm{d}}$–CLF relationship. The interpretation of possible underlying physical mechanisms of these interaction effects is difficult and remains speculative. The results seem to suggest that in these regions, potentially through hampering the entrainment of drier air from the free troposphere, the stronger inversion and more stable conditions are capable of trapping more moisture within a shallower MBL and could thus weaken the evaporation–entrainment feedback. As a result, it may ultimately favour a more positive $N_{\mathrm{d}}$–CLF relationship . It is interesting to note that these interactions are not apparent in the stratocumulus regions, where EIS is a strong control of CLF, and in the stratocumulus-to-cumulus transition regions, where found the aerosol effect on this transition to be confined to stable atmospheric conditions. This may imply that the suggested entrainment effect is dependent on the EIS and stronger at slightly lower EIS values typically found in the midlatitudes . The observed impact of EIS on the $N_{\mathrm{d}}$–CLF relationship found in the midlatitudes may also have implications within the context of climate change. While in the subtropics global climate models predict an increase in EIS with a warming climate, in the midlatitudes EIS is predicted to decrease , potentially decreasing the sensitivity of CLF to $N_{\mathrm{d}}$ there.

Figure b shows that higher SSTs are found to amplify the positive $N_{\mathrm{d}}$–CLF relationship (positive IAI) in the regions of frequent stratocumulus-to-cumulus transition . The physical interpretation could be the following: here, higher SSTs tend to lead to the transition from stratocumulus clouds to shallow convective clouds ; however, this transition has been found to be delayed when aerosol is increased . Tentatively, the positive IAIs in these transition regions may thus point to increased control of $N_{\mathrm{d}}$ on CLF at higher SST values as these are the situations where transitions typically occur and when increased $N_{\mathrm{d}}$ can act to delay this transition. In these regions, higher SSTs in the future might thus increase the sensitivity of MBLC CLF to aerosols. It should be noted that the quantification of the dependence of the $N_{\mathrm{d}}$–CLF relationship on meteorological factors (EIS, SST discussed in this section) is also likely subject to the biases in the $N_{\mathrm{d}}$–CLF sensitivity caused by the $N_{\mathrm{d}}$ retrieval biases as a function of CLF. This would potentially contribute to the non-causal facets of the relationships and interactive effects quantified by SHAP values.

Figure 8

Patterns of the interaction index showing the dependence of the $N_{\mathrm{d}}$–CLF relationship on estimated inversion strength (EIS) (a) and sea-surface temperature (SST) (b).

[Figure omitted. See PDF]

4 Conclusions

In this study, 9 years (2011–2019) of daily satellite and reanalysis data have been analysed to better understand the effect of $N_{\mathrm{d}}$ on CLF in MBLC and its dependence on meteorological factors. We have established a near-global machine learning framework to predict the cloud fraction of marine boundary clouds using regionally specific XGB regression models. Including many confounding and influencing factors as a whole, the explainable machine learning technique of SHAP regression values has been used to explain the regional XGB models; to quantify the CLF sensitivity to all cloud controlling factors with a specific focus on $N_{\mathrm{d}}$; and, moreover, to quantify the meteorological influence on the $N_{\mathrm{d}}$–CLF relationship at a global scale. The statistical sensitivities and interactive effects are interpreted with the guidance of hypothesized causal pathways and the state-of-the-art physical understanding of the system. The main findings of this study, which should be interpreted in light of the data and methodology limitations discussed in Sect. ), are summarized as follows:

• The marine boundary-layer cloud fraction shows a notable positive sensitivity to $N_{\mathrm{d}}$ (a surrogate for aerosols) in the regions of stratocumulus-to-cumulus transition, which may arise from the high $N_{\mathrm{d}}$ delaying this transition. The $N_{\mathrm{d}}$–CLF sensitivity in the Southern Hemispheric midlatitudes is observed to be higher than in previous studies, which should be investigated in future work. The estimated $N_{\mathrm{d}}$–CLF sensitivity and its magnitude suggest that aerosols likely have a considerable impact on MBL cloudiness although this may partially result from an overestimation caused by the effect of a positive retrieval bias of $N_{\mathrm{d}}$ at high CLF.

• Consistent with the literature, our statistical method shows that EIS and SST are two important determinants for low marine clouds by regulating surface fluxes and dry-air entrainment processes. In addition, strong negative CLF sensitivity and spatial patterns for SHF are also found, suggesting that the effect of cold air advection might surpass the SHF enhancement of closed-to-open-cell and cumulus transitions. Dynamic drivers (meridional and zonal winds) indicate that midlatitude synoptic-scale disturbances and vertical wind shear seemingly make considerable contributions to marine low cloud amounts.

• In general, thermodynamical parameters exert a more important influence on the $N_{\mathrm{d}}$–CLF relationship than dynamical parameters. EIS, RH₈₅₀, SST, and temperatures at 700 and 850 hPa have the strongest effect on the $N_{\mathrm{d}}$–CLF sensitivity. In the midlatitudes, higher EIS is found to amplify the positive $N_{\mathrm{d}}$–CLF sensitivity, which may be related to a reduced entrainment feedback in these conditions, whereas higher SST is found to amplify the $N_{\mathrm{d}}$–CLF sensitivity in stratocumulus-to-cumulus transition regions, which is potentially because the transition induced by higher SSTs may be delayed by increased $N_{\mathrm{d}}$. These findings have potential implications for possible future changes in the sensitivity of CLF to aerosols.

• For the dynamical and thermodynamical factors shown here, both CLF sensitivities and interactive effects (dependence of $N_{\mathrm{d}}$–CLF relationship on meteorology) exhibit distinct regional patterns. These coherent spatial patterns indicate that the proposed explainable machine learning framework not only is capable of skilfully predicting CLF for marine low clouds but also has the potential to capture regional characteristics of the relation between CLF and $N_{\mathrm{d}}$ as well as meteorological influences.

In the future, the observation-based sensitivities and interactive effects quantified by the ML framework here will be compared to those in ESMs, which have the potential to evaluate ESM parameterizations related to ACI and even help gain insights into how the models could be tuned in this respect. In addition, incorporating causal approaches for SHAP, such as those proposed by and , would help to test to which extent the observed statistical relationships and interaction effects represent physical processes.

Code availability

Code is available from the corresponding author upon reasonable request.

Data availability

All data sets used in this study are publicly available. The MODIS data set (10.5067/MODIS/MOD08_D3.061, ) was acquired from the Level-1 and Atmosphere Archive and Distribution System (LAADS) Distributed Active Archive Center (DAAC) (NASA: MODIS Data Collection, https://ladsweb.modaps.eosdis.nasa.gov/search/, last access: 17 November 2024); the hourly reanalysis data at single levels (10.24381/cds.adbb2d47, ) and pressure levels (10.24381/cds.bd0915c6, ) are obtained from the Copernicus Climate Change Service (C3S) Climate Date Store.

The supplement related to this article is available online at: https://doi.org/10.5194/acp-24-13025-2024-supplement.

Author contributions

HA and JC designed the initial research idea. YJ, HA, and JC developed the study concept and methodology. YJ and HA obtained and analysed the data sets. YJ implemented the explainable machine learning framework, performed the visualization, and wrote the original draft. All authors contributed to interpreting the results and reviewing and improving the paper.

Competing interests

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

Disclaimer

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

Acknowledgements

The (co-)authors have received funding from European Union’s Horizon 2020 research and innovation programme under grant agreement no. 821205 (FORCeS) and the Deutsche Forschungsgemeinschaft (DFG) as part of the project Constraining Aerosol-Low cloud InteractionS with multi-target MAchine learning (CALISMA; project no. 440521482). We thank three anonymous reviewers whose helpful comments contributed to improving the manuscript.

Financial support

This research has been supported by Horizon 2020 (grant no. 821205) and the Deutsche Forschungsgemeinschaft (grant no. 440521482).The article processing charges for this open-access publication were covered by the Karlsruhe Institute of Technology (KIT).

Review statement

This paper was edited by Yuan Wang and reviewed by three anonymous referees.