2.1 Data Set

The ε and δ data set used in this study is derived from global observations from June–August 2017 by the Visible Infrared Imaging Radiometer Suite (VIIRS) sensor onboard the Suomi National Polar-orbiting Partnership (SNPP) satellite. Additionally, the 6-hourly NCEP FNL reanalysis data (National Centers for Environmental Prediction, 2000) and National Oceanic and Atmospheric Administration (NOAA) sea surface temperature (SST) data (Huang et al., 2021), both on a 1° × 1° grid, were used in the retrieval process. The cloud product retrievals and the derivation of the ε and δ data sets are detailed in Zhu et al. (2025). The uncertainty analysis indicates that ε and δ exhibit relatively small uncertainty under the current retrieval framework (see Text S1 in Supporting Information S1 for details). Cloud droplet number concentration (N_d) was retrieved using cloud optical thickness (COT) and cloud droplet effective radius (r_e) (Szczodrak et al., 2001; Wang et al., 2021, 2023). The RH_e was obtained from the reanalysis data, while cloud base vertical velocity (w) was retrieved using the algorithm proposed by Zheng and Rosenfeld (2015). A total of 83,240 cloud samples were collected on a global scale, exhibiting broad spatial coverage across both land and ocean areas with no significant seasonal and regional biases (Figures S1 and S2 in Supporting Information S1). To ensure the reliability and generalizability of the parameterizations, 80% of the data was randomly selected for parameterization development, while the remaining 20% was reserved for independent evaluation. In the ML approach, the 80% portion used for model construction was further split into training and validation sets in an 8:2 ratio. The sampling process was repeated ten times with different random seeds to ensure robustness and quantify uncertainty.

2.2 Machine Learning Approach

Machine learning has found widespread applications in data analysis across multiple disciplines due to its ability to automatically learn data patterns and optimize models. Among various ML algorithms, gradient boosting, an ensemble learning method, combines multiple weak learners into a robust predictive model. Studies have demonstrated that it outperforms random forests, artificial neural networks, and convolutional neural networks (Abdollahi-Arpanahi et al., 2020; Ma et al., 2018; Yan et al., 2021). Light Gradient Boosting Machine (LightGBM), developed by Microsoft, is a gradient-boosting algorithm optimized for efficiently handling large-scale structured data sets (Ke et al., 2017). It offers advantages such as high training speed, low memory usage, and the ability to perform parallel computations. Gao et al. (2024) compared four ML algorithms for predicting cloud droplet number concentration and cloud droplet relative dispersion and found that LightGBM achieves the best performance in both computational efficiency and prediction accuracy, emphasizing its potential for advancing research on entrainment-mixing processes.

To achieve the study's objective of predicting ε and δ during entrainment-mixing processes with ML, LightGBM was selected as the preferred algorithm. Model input features were chosen based on their relevance to ε and δ and their availability in numerical models. For predicting ε, the selected features included environmental variables (RH_e, specific humidity (q_ve), environmental temperature (T_e), and environmental pressure (P_e)), in-cloud variables (w, liquid water content (LWC), N_d, r_e, liquid water path, and COT), and meteorological conditions (cloud fraction, lower tropospheric stability (LTS), SST, surface relative humidity (RH_surf), and terrain height (H_terrain)). For predicting δ, the observed ε was included as an additional input feature. Note that all variables used are temporally matched to the satellite observation time and spatially matched to a 1° × 1° grid, as detailed in Zhu et al. (2025).

One of the most widely used methods for interpreting ML is the SHapley Additive exPlanations (SHAP) (Lundberg & Lee, 2017), which has been widely applied for ML interpretation in geoscientific applications (García & Aznarte, 2020). SHAP values provide a fair allocation of each feature's contribution to the prediction, facilitating a deeper understanding of how ML makes predictions. For any given feature, a positive SHAP value indicates a positive contribution relative to the baseline prediction, while a negative value indicates a negative contribution. The mean absolute SHAP value is a standard metric for ranking feature importance. Given that entrainment-mixing is a small-scale process characterized by complex nonlinear interactions among variables, SHAP offers a robust framework for elucidating the complex physical mechanisms underlying ML prediction of ε and δ (Gao et al., 2024).

3.1 Parameterization of Entrainment and Detrainment Rates Based on the Physical Approach

Previous studies have proposed parameterizations through curve fitting and other physical approaches (Li et al., 2025; Xu et al., 2021; Zhang et al., 2016). In this section, parameterizations for ε and δ are developed using physical approaches. Following similar methods in earlier research (Lu et al., 2016; Xu et al., 2021), ε is parameterized as a function of RH_e and w. Figure 1 presents the fitted functions derived from a randomly selected 80% of the data, as well as the correlation coefficient (R) and root mean squared error (RMSE) between the fitted and calculated values. Additionally, it shows the R and RMSE between the fitted and calculated ε for the remaining 20% of the data, with calculated ε serving as the ground truth for evaluation. The parameterization based on RH_e (Figure 1a) outperforms that based on w (Figure 1b), achieving higher R and lower RMSE across both the fitting and test data sets. Due to the complexity of cumulus processes (Li et al., 2019; Wu et al., 2023; Yeom et al., 2025), ε is influenced by multiple factors (Lu et al., 2016; Xu et al., 2021). Therefore, a multi-variable parameterization using both RH_e and w, which are independent of each other, is presented in Figure 1c. The results indicate that the multi-variable parameterization (Figure 1c) yields higher R and lower RMSE than the single-variable cases (Figures 1a and 1b) in both fitting and test data sets. Overall, for the parameterization of ε, the multi-variable parameterization based on RH_e and w performs better than single-variable parameterizations.

For the parameterization of δ, similar to Gregory and Rowntree (1990), a linear fit based on ε demonstrates strong performance, with high R and low RMSE in both the fitting and test data sets (Figure 1d). Moreover, except for a few outliers, the fitted and calculated values of δ align closely with the 1:1 line. Compared to the linear fitting, a power-law fitting based on ε achieves comparable R yet lower RMSE, and the agreement between fitted and calculated values is closer to the 1:1 line (Figures 1d and 1e). The simple parameterization of δ is fine. However, given the significant influence of RH_e on δ and its frequent inclusion in δ parameterizations, this study also explores two alternative formulations (δ = ε(a − RH_e) and δ = aε(1 − RH_e), where a is a fitting parameter), inspired by Bechtold et al. (2014) and Stratton and Stirling (2012). Figures 1f and 1g illustrate the performance of these two parameterizations. Although the two forms are similar, the parameterization based on δ = ε(a − RH_e) outperforms the alternative, showing higher R and reduced RMSE in both the fitting and test data sets. Considering that the single-variable power-law parameterization exhibits enhanced performance, this study further optimizes the two-variable parameterizations by introducing a power-law exponent, resulting in two new parameterizations: δ = ε^b (a − RH_e) and δ = aε^b (1 − RH_e). Compared to the original schemes from Bechtold et al. (2014) and Stratton and Stirling (2012) (Figures 1f and 1g), the new parameterizations show comparable or improved R and reduced RMSE in both the fitting and test data sets. Moreover, the comparison between fitted and calculated values under the new parameterizations aligns more closely with the 1:1 line (Figures 1h and 1i). Overall, the δ parameterization based on δ = ε^b (a − RH_e) demonstrates the strongest performance.

3.2 Machine Learning-Based Prediction of Entrainment and Detrainment Rates

Due to the inherent complexity of the cumulus entrainment-mixing process, ε and δ are influenced by multiple factors, including environmental meteorological conditions, cloud microphysical properties, and macrophysical characteristics. As a result, in addition to developing traditional parameterizations, ML offers an alternative method for predicting ε and δ. Previous studies have shown that LightGBM outperforms other algorithms in investigating entrainment processes (Gao et al., 2024). Accordingly, this section utilizes the LightGBM to predict ε and δ.

Using LightGBM trained on 80% of the data, Figures 2a and 2b presents the validation results obtained using the remaining 20% of the data. The ε values predicted by the ML approach align well with the calculated values, with a high R and a low RMSE. Similarly, the predicted δ values exhibit strong correlations with the calculated values and align closely along the 1:1 line, also showing high R and low RMSE. Overall, the ML approach highlights its effectiveness in accurately predicting both ε and δ.

Although the ML approach shows strong performance in predicting ε and δ, interpreting how the ML generates these predictions is equally important. This not only enhances the physical interpretability of the ML approach but also offers useful insights for future parameterizations. Through SHAP analysis, the contributions of each input feature to the predictions can be explored. Figure 2 also presents the scatter plots of SHAP values for the features contributing to the predictions of ε and δ, as well as their rankings. From the SHAP scatter plots for predicting ε (Figure 2c), higher q_ve and RH_e contribute positively to the prediction, while higher T_e, w, P_e, and LWC have negative contributions. The feature ranking (Figure 2e) indicates that the top six features influencing ε are q_ve, RH_e, T_e, w, P_e, and LWC, followed by H_terrain, LTS, and SST. The remaining features have negligible impact. Among these, q_ve and RH_e are the most influential factors, with SHAP value contributions accounting for 26.6% and 20.6%, respectively. These results are consistent with physical expectations: higher RH_e reduces buoyancy and vertical velocity within the cloud core, allowing more time for cloud-environment mixing and thus increasing ε (Lu et al., 2016, 2018; Zhu et al., 2021; Zhu, Lu, Chen, et al., 2024). The global satellite observations used for training include cloud samples with sufficient variability (relevant to the natural selection mechanism) and are not confined to convection-favorable environments (associated with the cloud-size mechanism). Therefore, the ε-RH_e relationship learned by the ML approach is more likely to reflect the natural selection mechanism rather than the cloud-size mechanism, as discussed in Section 1 (Derbyshire et al., 2011). The ML approach effectively identifies the key factors influencing ε. Beyond RH_e and w, the interpretation of the ML also highlights the significant influence of T_e and P_e.

For the prediction of δ, Figures 2d and 2f also present the SHAP-based interpretation. The results reveal that an increase in ε contributes positively to the prediction of δ, and ε is the most significant one among the input features, with a SHAP value contribution as high as 55.9%. Other features that play relatively important roles in predicting δ include P_e, w, LWC, RH_e, T_e, SST, q_ve, and RH_surf. Among these, P_e, w, LWC, and T_e contribute positively to δ prediction, while RH_e, SST, q_ve and RH_surf have negative contributions. The remaining features exhibit negligible correlation with δ and thus can be neglected. The dominant contribution of ε to δ prediction is consistent with physical expectations, as entrainment and detrainment processes are closely related. Beyond ε, the contributions of environmental variables and w to the prediction of δ are also notable.

Note that the ML approach identifies several physically meaningful variables not directly involved in the retrieval of ε and δ as important predictors (such as w and RH_e), and SHAP analysis confirms that their relationships with ε and δ align with those in physical parameterizations, indicating that the ML approach likely captures underlying physical processes rather than retrieval-related patterns. Moreover, results show that the ML approach consistently yields accurate predictions of both ε and δ across diverse spatial and environmental conditions (Figures S3, S4, and S5 in Supporting Information S1), supporting the generalizability of the developed ML-based parameterizations. Nevertheless, we still acknowledge the possibility that the ML approach may learn the systematic and structural biases embedded in the retrieval algorithm.

3.3 Discussion: Comparison Between Physical and Machine Learning Approaches

The results from the previous two sections indicate that both physical and ML approaches achieve strong performance in predicting ε and δ. Table 1 provides a comprehensive comparison of the ε and δ parameterizations developed using the two approaches, based on statistical evaluation metrics (see Text S2 for details) from ten independent runs with different random seeds. Compared to the physical approach, the LightGBM trained on all features (LGB1) yields predictions with lower mean absolute error (MAE), mean squared error (MSE), and RMSE, indicating that ML outperforms the physical approach in accurately predicting ε. Notably, the physical approach relies on only one or two parameters for parameterization (e.g., the ε parameterization based on RH_e is referred to as P1_a; based on w as P1_b; and based on both RH_e and w as P1_c). In contrast, LGB1 uses 14 features for training and predicting ε. To ensure a fair comparison, this study further trained ML using the same input parameters as those employed in the physical approach. Specifically, the LightGBM trained on RH_e is labeled as LGB1_a, the LightGBM trained on w as LGB1_b, and the LightGBM trained on both RH_e and w as LGB1_c. As shown in Table 1, LGB1_a, LGB1_b, and LGB1_c yield lower MAE, MSE, and RMSE compared to P1_a, P1_b, and P1_c, respectively. This demonstrates that, even when using the same parameters for parameterization, ML consistently outperforms the physical approach in accurately predicting ε. These results demonstrate the superiority of ML over the physical approach in developing ε parameterizations.

Table 1 also summarizes the performance evaluation of the δ parameterizations. Compared to the physical approach, the LightGBM trained on all features (LGB2) significantly reduces the MAE, MSE, and RMSE between the predicted and calculated δ, demonstrating its superior predictive accuracy. Similarly, ML for δ prediction was developed using the same parameters as those used in the physical parameterizations. Specifically, the LightGBM trained on ε is labeled as LGB2_a, and the LightGBM trained on ε and RH_e is labeled as LGB2_b. The results show that LGB2_a and LGB2_b yield lower MAE, MSE, and RMSE compared to their respective physical parameterizations. This indicates that ML exhibits superior accuracy in predicting δ compared to the physical approach.

Compared to the physical approach, the ML approach inherently possesses more flexible functional forms, making the comparison structurally unbalanced. To ensure a more balanced comparison, a generalized additive model (GAM) with comparable complexity was applied to predict ε and δ using all 14 input variables. As shown in Table 1, the ML approach yields lower MAE, MSE, and RMSE than the GAM, confirming its superior predictive performance. Based on the evaluation of the prediction performance for ε and δ (Table 1), the parameterizations developed using ML exhibit significantly higher accuracy than those developed using physical approaches. This advantage likely arises from the fact that physical approaches are limited in their ability to account for the combined effects of multiple influencing factors, whereas ML can capture complex nonlinear relationships between variables through data-driven training. Furthermore, Figure S6 in Supporting Information S1 shows that the training and validation loss curves of the abovementioned ML models exhibit convergence without signs of overfitting (Gao et al., 2024; Seifert & Rasp, 2020). These findings offer valuable guidance for the future development of ε and δ parameterizations and support their applicability in broader climate and weather modeling contexts.