Quantifying Parametric Uncertainty Effects on Tropical Cloud Fraction in an AGCM

The underestimation of cloud fraction, especially the low stratus cloud (LSC) fraction over the eastern oceans, remains a problem in most atmospheric global climate models. This study investigated potential improvements through perturbing nine moist physical parameters, using uniform sampling and Latin hypercube sampling methods, and quantified the parametric uncertainty and effects of non‐linear interaction between parameters on the tropical cloud fraction in the Grid‐Point Atmospheric Model of the Institute of Atmospheric Physics/State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics, version 2. Results showed that the uncertainty ranges of the tropical total cloud fraction and LSC associated with multiple‐parameter perturbation were larger than those from any single‐parameter perturbation. The total cloud fraction was significantly improved with multiple‐parameter perturbation and the LSC also increased notably when using parameter values optimized for the total cloud fraction because of the indirect parametric effect on lower‐tropospheric stability. Non‐linear effects between parameters on simulating LSC were much stronger than those on simulating the total cloud fraction. Furthermore, the non‐linear interaction reduced large values of total cloud fraction but had a strong incremental effect on large values of LSC over the southeast Pacific. The findings demonstrate the feasibility of improving the simulation of LSC over the eastern oceans by tuning moist parameters, and provide further insight into the non‐linear effects between parameters on the simulation of cloud fraction.

In addition to the parameterization schemes themselves, the uncertain parameters within cloud and convective parameterizations are another important factor in determining model performance, and they introduce large uncertainties into climate simulations and projections (Murphy et al., 2004;Stainforth et al., 2005). Quaas (2012) revealed that the cloud fraction is sensitive to the choice of the relative humidity (RH) threshold. Similarly, M.  and Y.  reported the large sensitivity of the RH threshold for low clouds (RHMINL) and demonstrated that RHMINL tuning can improve the simulation of clouds in the National Center for Atmospheric Research (NCAR) Community Atmosphere Model. The importance of other parameters in the convection, cloud microphysics and macrophysics, and boundary layer cloud schemes in the simulation of cloud fraction has also been highlighted (Mauritsen et al., 2012;Pathak et al., 2020;Qian et al., 2018). Aside from the cloud fraction, much effort has also been made to quantify the parametric uncertainty or improve the performance of models in simulating phenomena and variability closely related to tropical cloud systems-for example, the Pacific Walker circulation (Xie et al., 2017, El Niño-Southern Oscillation (Watanabe et al., 2011), the Madden-Julian Oscillation (Boyle et al., 2015;X. W. Liu et al., 2019), the East Asian summer monsoon (Yang, Zhang, Qian, Huang, & Yan, 2015;Yang, Zhang, Qian, Wu, et al., 2015), and the diurnal cycle of tropical precipitation (Del Genio & Wu, 2010;Slingo et al., 2004;Stratton & Stirling, 2012;Y. Q. Wang et al., 2007). Moreover, for quantifying uncertainty and calibrating parameters in an effort to reduce model biases, numerous different approaches have been applied (C. S. Jackson et al., 2008;Mauritsen et al., 2012)-for instance, Monte Carlo or quasi-Monte Carlo method (C. Zhao et al., 2013), Latin hypercube selection (Stein, 1987), the downhill simplex algorithm (Severijns & Hazeleger, 2005; T. Zhang et al., 2015), multiple very fast simulated annealing algorithm (Ingber, 1989;Mu et al., 2004;Yang et al., 2012;Yang, Zhang, Qian, Huang, & Yan, 2015;Yang, Zhang, Qian, Wu, et al., 2015), and simulated stochastic approximation annealing algorithm (Yang et al., 2013).
Despite the considerable body of work outlined above, influences of parameters from different moist physical processes on the simulation of tropical cloud fraction, particularly on clouds at different altitudes and their geographical distributions, have not been sufficiently studied or quantified. This is an issue because being able to quantify the influences of parameters may help guide the direction and extent of changes in a model's ability to simulate cloud-related quantities or processes. Therefore, it is of great importance to quantify the parametric uncertainty effects on cloud fractions, not only for providing useful evidence and information for model optimization, but also for a better understanding of the physical processes within the climate system. In particular, it is well known that most coupled GCMs (CGCMs) suffer from warm sea surface temperature (SST) biases and too little low cloud in the southeast Pacific (e.g., Gordon et al., 2000;Kiehl & Gent, 2004;Large & Danabasoglu, 2006;Lin, 2007;Ma et al., 1996;McAvaney et al., 2001;Mechoso et al., 1995;Richter, 2015;Wittenberg et al., 2006). A model's inability to reproduce the abundance of low stratus cloud (LSC) can lead to excessive absorption of shortwave radiation in the upper ocean and an overly warm SST, which in turn induces positive feedback to the initial error (Richter, 2015;Zuidema et al., 2016). Previous studies have demonstrated that these simulation biases over the southeast Pacific are mainly caused by the uncorrected air-sea fluxes originating from the errors in the atmospheric model rather than the ocean model (G. Li & Xie, 2014;Xiang et al., 2017). Considering the poorly simulated low cloud over the eastern oceans in most CGCMs, we were also inspired to determine whether the simulation of low cloud or even marine stratus cloud could be improved through tuning multiple moist parameters in an atmospheric GCM (AGCM) as a first step toward tuning or calibrating a CGCM, even if they do not control the marine stratus cloud directly. Furthermore, given that the impacts of various parameters may interact (C. Zhao et al., 2013), another important aspect of parametric uncertainty is how to extract the non-linear interactions between parameters and quantify their impacts on the simulation of cloud fraction as well as its optimization. Such issues have not yet been investigated in depth, but they can be addressed effectively by using multiple-parameter perturbation combined with single-parameter perturbation, as developed by Haerter et al. (2009). That is, taking the results from single-parameter perturbation as a baseline, the effects of non-linear interaction between parameters can be derived and separated from the results of multiple-parameter perturbation (see Section 2.3.2 for further details).
Hence, the objectives of this study were (a) to quantify the sensitivity of tropical cloud-related quantities and the impacts of different moist parameters on the simulation of cloud fraction by applying single-parameter perturbation, (b) to explore the potential for improving the simulation of tropical cloud fraction by tuning multiple parameters, and (c) to investigate the effects of non-linear interaction between parameters. The remainder of the paper is organized as follows. The model, the selected moist physical parameters, the experimental design, and the observational data are introduced in Section 2. Section 3 analyzes the results, including the isolated effects of individual parameters, the impacts of optimization, and the effects of non-linear interaction among multiple parameters on the simulation of tropical cloud fraction. Section 4 provides a summary and some further discussion.

GAMIL2 Model
This study used the Grid-Point Atmospheric Model of the Institute of Atmospheric Physics, National Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics, version 2 (GAMIL2, L. J. Li, Wang, et al., 2013), which is the atmospheric component of the Flexible Global Ocean-Atmosphere-Land System model, Grid-point version 2 (L. J. Li, Lin, et al., 2013). GAMIL2 uses a finite difference dynamical core that conserves the total mass and effective energy under the standard stratification approximation (B. Wang et al., 2004), with a 128 × 60 horizontal grid and 26 sigma-levels in the vertical resolution.
GAMIL2 includes the following moist physical processes: deep convection, shallow convection, cloud fraction, and cloud macro-/microphysics. The G. J. Zhang and McFarlane (1995) and G. J. Zhang and Mu (2005) parameterizations are used for deep convection, which is a mass flux cumulus convection scheme. The Hack (1994) parameterization is used to characterize the shallow convective activities. There are three types of cloud in GAMIL2: convective cloud, low-level marine stratus, and layered cloud (Rasch & Kristjánsson, 1998). The layered cloud, which dominates the cloud fraction output, is diagnosed by a quadratic function of RH (i.e., it forms in the layer where the RH exceeds its threshold); the convective cloud is calculated by detrainment in the convective updraft mass flux; and the low-level marine stratus cloud is estimated by the lower-tropospheric stability defined by the potential temperature difference between 700 hPa and sea level. The total cloud fraction is estimated by the "maximum-random" overlap assumption, whereby vertically continuous clouds are assumed to be maximally overlapped while clouds at different heights that are separated by an entirely cloud-free model level are randomly overlapped (Hogan & Illingworth, 2000). The stratiform cloud process parameterization is formulated by the cloud macrophysics (M. H. Zhang et al., 2003) and cloud microphysics (Morrison & Gettelman, 2008). The macrophysical processes control the mass exchanges between the water vapor and cloud condensates, and the temperature changes associated with the phase changes, which together with the diagnostic cloud relationship that represent the subgrid-scale variability of RH, form a closure to calculate the fractional condensation rate (M. H. Zhang et al., 2003). The microphysical processes describe the conversion from condensate to precipitate and the evaporation process (Morrison & Gettelman, 2008). The boundary layer scheme is a nonlocal K-profile scheme, which calculates diffusion coefficients based on the diagnosed boundary layer height and the related characteristic velocity scale (Holtslag & Boville, 1993;Vogelezang & Holtslag, 1996).

Selected Moist Physical Parameters
The focus in this study was on nine tunable parameters from the moist physical processes. A description of the parameters along with their default values and perturbed ranges is given in Table 1. The parameter ranges were based on the experience of the model developers and referred to previous studies (C. S. Jackson et al., 2008;Yang et al., 2013;Yang, Zhang, Qian, Huang, & Yan, 2015;. Of these nine parameters, five were from the deep convection scheme: the precipitation efficiency for deep convection (C0_deep), the evaporation efficiency for deep convection (KE), the threshold value for convective available potential energy (CAPE) for deep convection (CAPELMT), the initial cloud downdraft mass flux (ALFA), and the threshold value for RH for deep convection (RHCRIT). The timescale for the consumption rate of shallow CAPE (CMFTAU) and the precipitation efficiency for shallow convection (C0_shc) were shallow convective parameters. The threshold RH for low clouds (RHMINL) and high clouds (RHMINH) were layered cloud fraction parameters closely related to the stratiform cloud processes.

Sampling Methods
To explore the parametric uncertainty effects and the non-linear interactions between parameters, two sampling methods were utilized. One was uniform sampling (US) in which we modified only one parameter at a time and maintained all others at their default values (Table 1). Sampled values were equally spaced between the lower and upper bounds sampled for all parameters. The number of samples was not the same for all parameters. For example, RHCRIT had the least (61) and CMFTAU the most (127). These different sample numbers were used to ensure the sampled points contained the default parameter values, to facilitate calculation of the perturbed deviation of each individual parameter. A total of 807 integrations were performed for the US-experiments.
The other method was Latin hypercube sampling (LHS), a technique that is effective for a large number of dimensions, in which all parameters are varied simultaneously within the nine-dimensional parameter space, making it possible to improve the model skill by optimizing multiple parameters and explore the joint effects of non-linear interactions between parameters. It uses a stratification strategy to produce good uniformity in each dimensional projection. This method divides the range of each parameter into several equal bins and makes each bin appear once at an optimal distance, thus ensuring an even distribution of the sampled points in the multi-dimensional space. Loeppky et al. (2009) suggested that the number of sampled points should be at least 10 times the number of parameters. In this study, 130 perturbed points were sampled for the nine parameters-a ratio exceeding 14.  (Yang et al., 2013). First, it reduced the computational costs. More importantly, however, 2003 was not an El Niño or La Niña year, and was thus less likely to have been affected by interannual variabilities. All of the 807 runs in the US-experiments were valid, as were 128 out of the 130 runs in the LHS-experiments. Instability occurred in two runs in the LHS-experiments, so were not included in the analyses.

Analysis Method
Isolated effects of individual parametric uncertainty were represented by the model responses to a single one of the nine parameters in the US-experiments, which were measured as regression coefficients between the output variables and parameters (Section 3.1). Note that, in order to quantify and compare sensitivities of different parameters, each parameter value was normalized to [0, 1] by subtracting the minimum within its range and then scaling by the difference between the maximum and minimum.
To describe model variable changes due to the perturbed parameters in the LHS-experiments, we express the output variable as follows: where LHS ( ) is the output variable simulated by the LHS-experiments; is the parameter vector {C0_deep, KE, CAPELMT, ALFA, CMFTAU, RHCRIT, C0_shc, RHMINL, RHMINH}; 0 is the default parameter vector (normalized to [0, 1]); and x is the perturbation from 0 .
In order to derive and separate the joint effects of non-linear interaction between parameters (Section 3.3), four calculation steps were used following Haerter et al. (2009) with some modifications. First, third-order polynomial regression-fitting equations for each parameter were obtained by regressing the output variables in the US-experiments to a single parameter to generate the isolated effect of the individual parameter. Second, each parameter perturbed in the LHS-experiments was substituted into its corresponding regression-fitting equation to gain the fitted values as f it ( 0 + ) in Equation 2, where i = 1, …, 9 refers to the nine parameters. Third, the default output value was subtracted from the fitted values for each parameter, and then the values for each parameter were summed and superimposed on the default output value to obtain the reconstructed values as the right-hand side of Equation 2, which represent the linear superposition of the isolated effect of each of the parameters: Here, Recons ( ) is the reconstructed variable value, (3)

Observational Data
The observational data used in this paper were the monthly mean cloud fraction from the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) GCM-oriented CALIPSO Cloud Product; resolution (GOCCP): 2.0° × 2.0° data set (Chepfer et al., 2010) for the period January 2007-December 2020. The multi-annual mean cloud fraction from CALIPSO-GOCCP during 2007-2020 was calculated as the climatological mean state for evaluating the model.

Results
The annual mean total cloud, high cloud, mid-level cloud, and low cloud fractions, from observations and simulated by the standard GAMIL2 with default parameter sets, are presented in Figure 1. The high, mid-level, and low clouds in GAMIL2 are defined as clouds within P ≤ 440, 440 < P ≤ 700, and P > 700 hPa, respectively, which are close to the CALIPSO-GOCCP definitions of P ≤ 440 hPa for high clouds, 440 < P ≤ 680 hPa for mid-level clouds, and P > 680 hPa for low clouds. The model captures general characteristics of the observed spatial pattern of total cloud fraction reasonably well, but with some biases; in particular, the total and high cloud fractions over 30°N-30°S in GAMIL2 are systematically smaller than those in CALIPSO-GOCCP. GAMIL2 also evidently underestimates the low cloud maxima near the eastern coastlines of tropical and subtropical oceans. This underestimation is a widespread problem for most GCMs, especially over the marine stratus regions near the eastern oceans (e.g., Karlsson    . Note that the mean biases in high cloud, mid-level cloud, and low cloud fractions over 30°N-30°S are no more than 0.039, whereas the bias of the total cloud fraction reaches 0.144. These discrepancies between the model and observations may be associated with the cloud vertical overlap scheme and the fact that satellite cloud simulators were not used. However, the focus of this paper is the effects of moist parameters on these biases. In this section, we first quantify and analyze the sensitivity of cloud-related quantities over the tropics to the nine selected moist physical parameters and how these parameters affect the model simulation, including the moisture, cloud fraction, and grid-box cloud liquid water path (GLWP, a vertical integral of liquid water content between two levels in the atmosphere). Next, the potential for improving the simulation of tropical cloud fraction by tuning multiple parameters is explored. And finally, we quantify the joint effects of non-linear interaction between parameters on the total cloud fraction by comparing the results reconstructed by the US-experiments and those simulated by the LHS-experiments.

Isolated Effects of Individual Parameters in the US-Experiments
The isolated effects of individual parametric uncertainties were quantified and analyzed using the 807 US experiments. Table 2 lists the cloud-related physical variables over the tropics that are usually considered as important in measuring AGCM performance. The responses of these variables (see Table 2) to the nine moist physical parameters are summarized in Figure 2.
In terms of the variables, the sensitivity ranking and the sign of the response of total precipitation to most parameters (except C0_shc) are consistent with those of surface latent heat flux (LH), indicating that the greater the total precipitation, the stronger the surface evaporation. The response of convective precipitation is anticorrelated with that of stratiform precipitation, except for CAPELMT and RHMINL. The sign of the response of net shortwave radiation at the top of the atmosphere (FSNT_TOA) is opposite to that of the total cloud fraction, as an increase (decrease) in total cloud fraction could induce stronger (weaker) shortwave cloud radiative forcing (SWCF). Such linkages between the responses of total cloud fraction, outgoing longwave radiation at the top of the atmosphere (OLR_TOA), and longwave cloud radiative forcing (LWCF) are also found for most parameters, except KE and CAPELMT. The reason is that larger KE induces more absorption of longwave radiation by moisture and less clear-sky OLR_TOA (OLRC_TOA) (not shown), which exceeds the extent of the decrease in OLR_TOA, thus leading to weaker LWCF. As for CAPELMT, the responses based on the regression coefficients are not statistically significant, due to its extremely weak impact, which gives inconsistent responses of total cloud fraction, OLR_TOA, and LWCF. In terms of the parameters, RHCRIT is the one with the greatest impact, with strong influences on most of the variables, except FSNT_TOA and SWCF. CMFTAU can also strongly affect most variables but has a smaller impact on specific humidity at 850 hPa. Both RHCRIT and CMFTAU affect most of the variables throughout the troposphere. In contrast, C0_deep has a more significant influence on mid-and high-level model variables, including mid-level and high cloud fractions, specific humidity at 400 hPa, and temperature at 200 hPa, and KE has a more significant influence on the lower-level variables, including the low cloud fraction, specific humidity and temperature at 850 hPa. One possible reason for this is that the maximum height that deep convective activities can reach exceeds 300 hPa ; the role of C0_deep is to affect the cloud liquid water conversion to rain and water vapor   Table 2). Sensitivity ranking was calculated based on the regression coefficients between output variables and input parameters from the US-experiments. Positive and negative correlations are denoted by plus and minus symbols, respectively.
condensation throughout most of the troposphere and the role of KE-a coefficient multiplied by the precipitation flux at each model layer-is to evaporate convective precipitation into the grid-scale environment more in the lower layer owing to the larger convective precipitation flux near the surface. RHMINL is very important for cloud-related variables in the lower layer, including the low cloud fraction, LH, LWP, and SWCF, as well as FSNT_TOA. In contrast, RHMINH mainly affects the mid-level and high cloud fraction, as well as LWCF. The influences of CAPELMT, ALFA, and C0_shc are generally weak, and these three parameters will not be included in the subsequent analyses.
The cloud fraction output in GAMIL2 is dominated by the layered cloud that is directly diagnosed by RH. The vertical profiles of RH responses and geographical distributions of cloud responses are shown in Figures 3 and 4, respectively. RHCRIT is the parameter with the strongest impact on the simulated RH (Figure 3d), but has a relatively weaker impact than RHMINL on the total cloud fraction (Figures 4m and 4q). Larger RHCRIT greatly suppresses deep convection by making the trigger conditions more difficult to attain, resulting in less water vapor being transported from the lower layer to the mid and high layers and more water vapor accumulating in the lower layer, which correspondingly means that the RH increases significantly below 700 hPa and decreases significantly above 700 hPa. Furthermore, owing to the diagnostic relationship between RH and the layered cloud fraction, low cloud increases while mid-level and high cloud decrease (Figures 4n-4p). However, larger RHCRIT decreases the low cloud fraction over the marine stratus regime. This is possibly a response to weakened convection in remote regions as indicated by S. Liu et al. (2021), who pointed out that convective parameter perturbations in remote regions, especially those in the Intertropical Convergence Zone (ITCZ), can induce considerable impacts on the simulated cloud fraction over the southeast Pacific by modulating subsidence and moisture there. In contrast, the effect of RHMINL on the low cloud fraction results primarily from the role as a threshold in modulating low layered cloud, rather than from modulating moisture. An increase in RHMINL makes it difficult for low layered cloud to generate, resulting in a significant reduction in low cloud. The areas where RHMINL has its strongest effect on the low cloud fraction are in regions with strong subsidence over the eastern tropical and subtropical oceans, suggesting that parameter optimization involving RHMINL could potentially help reduce the model bias for low cloud there. Similarly, RHMINH regulates the generation of mid-level and high cloud through its role as a threshold (Figures 4w and 4x). Larger CMFTAU can induce more total cloud by increasing the low, mid-level and high clouds (Figures 4j-4l). The underlying mechanism is that enhanced CMFTAU decreases the intensity of shallow convection by raising the adjustment time and strengthens deep convection through the interactions among moist processes , meaning more water is transported upwards by the strengthened deep convection, making the layer above (below) 850 hPa wetter (drier) (specific humidity not shown). Although the specific humidity in the lower layer decreases, the RH increases because the lower level is colder owing to larger CMFTAU (not shown). Thus, the RH increases throughout the troposphere, which contributes to the increase in the total cloud fraction. Larger C0_deep converts more liquid cloud water to precipitation and leads to less mid-level and high cloud (Figures 4c and 4d), corresponding to reduced moisture at 850-400 hPa over the tropics (Figure 3a). In contrast to C0_deep, the effect of larger KE is to enhance the evaporation efficiency of convective precipitation in the lower layer, which increases the RH below 700 hPa (Figure 3b), leading to more low cloud (Figure 4f).
Marine stratus is an important part of low cloud and highly sensitive to the low-level thermodynamic structure-in particular, the strength of the inversion at the cloud top (Gettelman & Sherwood, 2016;Klein & Hartmann, 1993;Wood & Bretherton, 2006). As seen in Figure 4n, larger RHCRIT causes a decrease in the low cloud fraction over the eastern coastlines of tropical and subtropical oceans where marine stratus persistently prevails, opposite to its impact in other regions where it causes an increase in the low cloud fraction. The responses of LSC and the potential temperature difference between 700 hPa and sea level to some of the parameters are illustrated in  Larger RHCRIT leads to less LSC over the eastern tropical and subtropical oceans by reducing the potential temperature difference between 700 hPa and sea level. RHMINL has a similar but weaker impact on LSC. The areas affected by CMFTAU are located mainly over the eastern oceans in the Southern Hemisphere, whereas the responses of LSC in the Northern Hemisphere are weak. In summary, moist parameters not only directly affect the simulated layered cloud fraction by modulating moisture, but also indirectly affect the LSC through modulating the lower-tropospheric stability. Therefore, parameters in moist processes could potentially be tuned or calibrated to enhance the model's ability in simulating the LSC, and even the low cloud fraction, over marine stratus regions.
In addition to cloud fraction, vertical profiles of the GLWP response are also affected by the parameters (Figure 6). In general, when the response of cloud fraction to a certain parameter is strong (weak), the response of GLWP is also strong (weak) (Figures 4 and 6). RHMINL and CMFTAU are the most important parameters,  having strong but opposite influences on GLWP. With a maximum at 900 hPa, GLWP increases below 500 hPa in response to larger CMFTAU but changes little above 500 hPa (Figure 6c), whereas the mid-level and high cloud increases (Figures 4k and 4l). Larger RHMINL leads to less GLWP below 850 hPa (Figure 6e). In contrast, larger RHMINH could cause less GLWP between 800 and 300 hPa (Figure 6f). Like the moisture and cloud fraction responses, the GLWP response to increased RHCRIT shows apparent differences between lower and upper levels, with more GLWP below 650 hPa but less above (Figure 6d). The impacts of C0_deep and KE are relatively weak (Figures 6a and 6b). Larger KE induces more precipitation evaporation and more GLWP below 500 hPa (Figure 6b), but the height of the strongest GLWP response is a little higher than that of the strongest cloud fraction response (Figure 4f).

Impacts of Optimization on the Cloud Fraction in the LHS-Experiments
To investigate the potential for improving the simulation of cloud over the tropics by parameter tuning, 128 simulations from multiple-parameter perturbation were carried out, each with a different parameter set generated by the LHS method. Figure 7 shows the frequency distribution of "good" simulations, defined as simulations that produced a smaller root-mean-square error (RMSE) of tropical total cloud fraction based on the CALIPSO-GOCCP data than the standard simulation that used the default values for individual parameters. RHMINL and CMFTAU have the most significant effects on model RMSEs. When RHMINL is smaller than 0.89, the probability of obtaining a good experiment is higher than 80%, but when RHMINL is larger than 0.94, the frequency of good simulations is always below 20%. When CMFTAU becomes larger than 8,000 s, the probability of a good simulation reaches 80%, except when the value of CMFTAU is around 10,500 s. Some parameters have peaks in their frequency distribution. For example, when KE is around 0.7 × 10 −5 (kg m −2 s −2 ) −1/2 s −1 , there is a very high probability of a good simulation. CAPELMT has double peaks in the frequency of good simulations, at values of around 50 and 140 m −2 s −2 . For other parameters, the frequency distributions are either more uniform or irregular, implying large uncertainties in those parameters based on skill measures defined by RMSE. In particular, although the sensitivity to RHCRIT in the US-experiments is evidently strong, the frequency distributions of RHCRIT are uniform, which may result from the skill measure chosen here or complex non-linear interaction between parameters.
The main bias of the standard GAMIL2 simulation is the underestimated total cloud fraction shown in Figure 1.
To measure whether the bias can be reduced through multiple-parameter perturbation, metrics including the average tropical cloud fraction, RMSE, and spatial correlation coefficient of cloud fraction between simulations in the LHS experiments and the CALIPSO-GOCCP data were calculated. These metrics are arranged in ascending order of the average tropical total cloud fraction in Figure 8. The bias of the underestimated total cloud fraction is apparently reduced in the LHS-experiments. The largest simulated total cloud fraction can reach the observed value of 0.617 (Figure 8a), which is contributed by the increase in the low cloud fraction (Figure 8d). With the reduced mean bias, the corresponding RMSE of the total cloud fraction is also smaller (Figure 8e), whereas the spatial correlation coefficient changes little (Figure 8i). For the low cloud fraction, despite the RMSE increasing to some degree, the spatial correlation coefficient unexpectedly increases. This implies that when the simulation of the total cloud fraction is improved, that of the spatial distribution of low cloud could also potentially be improved at the same time via parameter optimization. To investigate how the simulated spatial distributions of cloud fraction can be improved, a sub-ensemble of the top 15 simulations (SET15) with the least mean bias of total cloud fraction in the LHS-experiments was selected (the last 15 simulations ranked in Figure 8a). Distributions of the ensemble mean of SET15 and its difference from the standard GAMIL2 simulation are shown in Figure 9. Generally, SET15 outperforms the standard GAMIL2 simulation. The tropical mean total cloud fraction substantially increases from 0.473 in the standard GAMIL2 simulation to 0.591 in SET15, which is much closer to the observed value of 0.617 (Figures 1a and 9a). Also, the high cloud fraction increases from 0.310 in the standard GAMIL2 simulation to 0.347 in SET15, very close to the observed value of 0.348 (Figures 1b and 9b). The increase in total cloud fraction is mainly from the increase in low cloud fraction. More importantly, the observed low cloud maxima near the eastern coastlines of the tropical and subtropical oceans, which remain a challenge for the vast majority of current AGCMs, are successfully reproduced by SET15. Specifically, the LSC over the eastern oceans is enhanced significantly, contributing to the reduced bias in low cloud fraction ( Figure 10). That is, a remarkable improvement in the simulated low cloud maxima near the eastern coastlines of the tropical and subtropical oceans can be simultaneously obtained through improving the model's performance in simulating the total cloud fraction via parameter optimization.
To illustrate which parameters play important roles in the model's improvement, the values of the nine parameters were also arranged in ascending order of the tropical average total cloud fraction ( Figure 11). It is apparent that RHMINL has an obvious descending trend, indicating that it is the most important parameter for reducing the biases of total cloud (Figure 11h). The much smaller RHMINL in SET15 than in the default simulation can lead to an increase in total cloud and low cloud fraction over the tropics, especially an increase in LSC over the eastern oceans (Figures 4q, 4r, and 5c). Meanwhile, with a generally larger value than in the default simulation (Figure 11e), CMFTAU also contributes to the improvement by increasing the total cloud and low cloud fraction ( Figures 4i and 4j). The changes in these two parameters can increase the high cloud fraction slightly (Figures 4l and 4t). Another important cause of the increased high cloud fraction can be attributed to the smaller than default value of RHCRIT (Figure 11f). Despite leading to less low cloud (Figure 4n), a smaller RHCRIT induces a remarkable increase in LSC over the marine stratus regions, as expected (Figure 5b). Other parameters either vary randomly around their default values (i.e., KE, CAPELMT and RHMINH) or have less of an impact, as shown in Figures 4 and 5 (i.e., C0_deep, ALFA, and C0_shc).

Effects of Non-Linear Interaction Between Parameters
The non-linear interaction between parameters could be identified by comparing the results from the US-experiments with those from the LHS-experiments using Equations 1-3. Figure 12 shows the responses of the tropical total cloud fraction to the nine parameters from the US-experiments. The total cloud fraction is strongly correlated with every individual parameter perturbation, and is even linearly related to some parameters (i.e., KE, RHMINL, and RHMINH). All nine parameters have significant effects on the total cloud fraction. Considering that the responses to some parameters (i.e., ALFA and C0_shc) are not linear, we applied third-order polynomial regression-fitting to the perturbed parameters. All regression-fitted lines are statistically significant at the 95% confidence level, and the coefficients of the third-order polynomial regression-fitting are given in Table 3. It is clear that more total cloud can be obtained with larger (smaller) values for KE, CAPELMT, ALFA, and CMFTAU (C0_deep, RHCRIT, RHMINL, and RHMINH).  (Figure 13a). This dependence means that the effects of the individual parameter perturbations manifest in a joint manner in the LHS-experiments. As expected, the uncertainty range of F CLDTOT (0.337-0.608) is larger than that found for any of the individual variations (0.40-0.55) shown in Figure 12, indicating that the uncertainty of F CLDTOT due to parameter combinations is amplified by multiple-parameter perturbation. Furthermore, to explore whether this large uncertainty range is caused by non-linear interaction between changes due to the nine parameters or represents a linear superposition of the isolated effects of each individual parameter, the reconstructed total cloud fraction was compared with that obtained from the LHS-experiments by adopting the method introduced in Section 2.3.2 (Figure 13b). The larger the difference between the two, the larger the non-linear effect. The points of F CLDTOT in Figure 13b are generally scattered around the diagonal. However, there are overall biases for values of F CLDTOT between 0.4 and 0.53, and for large values of F CLDTOT . When F CLDTOT is between 0.4 and 0.53, the reconstructed values of total cloud fraction are often smaller than the values from the LHS-experiments. In contrast, as F CLDTOT becomes larger, the reconstructed values are larger than those obtained from the LHS-experiments. That is, to some degree, the non-linear effect among multiple parameters tends to increase the simulated values of total cloud fraction between 0.4 and 0.53 but reduce the large simulated values of total cloud fraction compared with the linear superposition of the isolated effects of individual parameters. The reduction effect of non-linear interaction between parameters on large values of F CLDTOT is also demonstrated by the smaller slope of the total cloud fraction from the LHS-experiments depending on ( ) than the slope of the reconstructed values (Figure 13a). In addition, differences between multiple-parameter perturbation and the linear reconstruction of large F CLDTOT (blue box in Figure 13b), corresponding to the improved total cloud fraction in SET15, imply that non-linear interaction between parameters plays an important role in the improvement.  perturbation ( Figure 14e). The points of F LSC in the LHS-experiments deviate more from the diagonal than those of F CLDTOT , indicating a much stronger non-linear effect between parameters on the simulated LSC ( Figure 14f). Opposite to F CLDTOT , when F LSC is smaller than 0.25, most of the reconstructed values of LSC are larger than the values from the LHS-experiments; but when F LSC is larger than 0.3, the reconstructed values are smaller than those obtained from the LHS-experiments (Figure 14f). The slope of the LSC from the LHS-experiments depending on ( ) is larger than that of the reconstructed values ( Figure 14e). That is, unlike the reduction effect on the simulated large total cloud fraction, the non-linear interaction among multiple parameters evidently increases the simulated large values of LSC over the southeast Pacific. Moreover, the incremental effects of non-linear interaction between parameters on the simulated large values of LSC over the southern Atlantic, northern Atlantic, and northeast Pacific are clear to see ( Figure S1).
What remains unclear is why the effects of non-linear interaction between parameters on the LSC over the southeast Pacific are stronger than those on the total cloud fraction. In this regard, we speculate that the probable cause is that LSC simulation incorporates a more indirect effect of moist parameters than total cloud fraction. First, the strong link between the formation of LSC and lower-tropospheric stability implies that LSC is tightly coupled with the general circulation (Bony et al., 2015;Bretherton & Hartmann, 2009;Wood, 2012). Prevalent over the cold parts of the eastern tropical and subtropical oceans, LSC is favored in the descending branches of    large-scale atmospheric circulations such as the Hadley and Walker circulations (Klein & Hartmann, 1993;Lipat et al., 2017;Ma et al., 1996;Wood, 2012), which are not directly linked with moist parameters. Xie et al. (2017) revealed that parameter changes indirectly affect the strength of the descending branches of the Pacific Walker circulation through changes in vertical velocity linked to diabatic heating and emphasized the importance of lower-tropospheric stability responses in modulating vertical velocity. Such indirect parametric effects associated with complex dynamics are thought to introduce more non-linear interactions between parameters into LSC simulations. In contrast, the main regions where the total cloud fraction is sensitive to parameter changes are where convection activities are strong, such as the equatorial western Pacific and the ITCZ, which are directly affected by moist parameters, leading to weaker non-linear interaction between parameters. Second, the LSC regions are located to the west of adjacent continents with substantial orography. Such high and steep orography has a strong impact on the behavior of LSC through modulating local circulation (Richter & Mechoso, 2006). For example, removal of the Andes Mountains in a regional climate model significantly reduced the incidence of stratocumulus off the Peruvian coast owing to the intrusion of dry, warm air from the South American continent within the planetary boundary layer (Xu et al., 2004). A similar orographic impact on stratocumulus but from local circulation changes above the planetary boundary layer was revealed by Richter and Mechoso (2004). Therefore, orographically induced processes, not directly affected by moist parameters, will probably involve stronger non-linear interaction between parameters in LSC simulations. By contrast, smaller-scale orography within the region of large total cloud fraction would result in weaker non-linear interaction between parameters.

Summary and Discussion
It remains a challenge for current AGCMs to successfully reproduce the observed tropical cloud fraction, especially over the marine stratus regions near eastern oceans. In this study, using the GAMIL2 model, we investigated the influences of nine moist physical parameters on the simulated tropical cloud fraction and cloud-related quantities, and further examined the potential of improving the simulation of cloud fraction by tuning multiple parameters and exploring the effects of non-linear interaction between parameters. Two sets of experiments were carried out. In the US-experiments, we varied only one parameter, keeping all others at their default values. In the LHS-experiments, the nine parameters were varied simultaneously using a Latin hypercube sample approach.
Results showed that the standard GAMIL2 had large biases in simulating the total cloud fraction over 30°N-30°S, in particular underestimating the low cloud fraction maxima over the eastern coastlines of the tropical and subtropical oceans. The tropical cloud-related variables were sensitive to six moist physical parameters. Changes in the convective parameters of C0_deep, KE, CMFTAU, and RHCRIT affected the layered cloud fraction by modulating the vertical transport of moisture and condensation/evaporation associated with moist physical processes, while the effects of parameters from the layered cloud fraction scheme (i.e., RHMINL and RHMINH) resulted mainly from their roles as thresholds in generating layered cloud. The cloud fraction was very sensitive to CMFTAU, RHCRIT, and RHMINL. Besides their direct impacts on the layered cloud fraction, the indirect effects of these parameters on marine stratus over the eastern oceans, linked to lower-tropospheric stability, contributed to the responses of the low cloud fraction and even the total cloud fraction. Using appropriate parameter values in the LHS-experiments, the ensemble mean of SET15 exhibited much better performance in simulating the total cloud fraction and succeeded in reproducing the low cloud fraction maxima over the eastern coastlines of the tropical and subtropical oceans. The significantly enhanced LSC over the eastern oceans that contributed to the increased low cloud fraction indirectly resulted from smaller values of RHMINL and RHCRIT, and larger values of CMFTAU, than their defaults. This indicated that a remarkable improvement in the simulation of the low cloud maxima near the eastern coastlines of the tropical and subtropical oceans can be simultaneously achieved through improving the model's performance in simulating the total cloud fraction via parameter optimization.
The multiple-parameter perturbation in the LHS-experiments led to larger uncertainties in the tropical total cloud fraction by about 0.27, and in the LSC over the southeast Pacific by about 0.59, than those found for any of the individual variations in the US-experiments, indicating that the perturbation for multiple parameter combinations could amplify the uncertainty. The results reconstructed from the US-experiments and LHS-experiments suggested that, even though changes due to isolated effects of individual parameter perturbation could be linearly superimposed, substantial non-linear interaction between parameters did indeed occur, which tended to increase the total cloud fraction where it had values within 0.4-0.53, but reduce it to some degree for large values of total cloud fraction. Compared with total cloud, the effects of non-linear interaction between parameters on marine stratus cloud were much stronger and tended to reduce the LSC over the southeast Pacific within values of 0.04-0.25 but increase it for large values of LSC. These findings provide further insight into the non-linear interaction between parameters within cloud fraction simulations and motivate further studies on this non-linear interaction within other important physical processes or phenomena in climate simulations.
Several limitations of this study should be taken into consideration and deserve further investigation. One of the important purposes of uncertainty quantification is to improve model performance by parameter tuning or calibration. Selecting the parameter for tuning is crucial to the optimization result. Taking the parameter sensitivities in Figure 2 as an example, RHCRIT had the most effect, whereas the impacts of CAPELMT and C0_shc were extraordinarily weak. If only these three parameters were selected for tuning, it is likely that the perturbation of RHCRIT would dominate the variation of the simulated climate system, weakening or even eliminating the impact of the other two. This is not what we desire. Therefore, more attention should be paid to parameter selection before model optimization. Besides, model biases induced by parameters are only a portion of the total bias. Thus, studies on other limitations or deficiencies in parameterizations, such as closure assumptions and systematic structure errors in models, are also important for identifying the source of model biases (M. . In addition, the intercomparison of clouds from direct model output and from CALIPSO-GOCCP is not completely accurate. In future work, using a satellite simulator could help to achieve a more rigorous evaluation of the clouds simulated by GAMIL2. Another issue that deserves further investigation is the interaction among deep convective, shallow convective, and stratiform cloud processes in AGCMs. Since precipitation is a complex product of these interacting processes, it cannot be guaranteed that the improved cloud simulation is accompanied by better model performance in precipitation. It is therefore useful to evaluate the specific behaviors of these moist processes in generating precipitation. Moreover, interaction among moist processes may be a temporary problem for the modeler and climate researcher, which results from the artificial separation of moist physical processes in current AGCMs. To address this issue in the future, Arakawa (2004) emphasized the importance of developing a conceptual framework for a "unified cloud parameterization" or "unified model physics" for AGCMs. Yang et al. (2021) coupled the Zhang-McFarlane deep convection scheme with a unified parameterization for boundary-layer turbulence, shallow convection, and cloud macrophysics, which improves the simulation of precipitation in terms of the mean state and variability at various timescales.
It is not quantitatively explicit whether the non-linear interactions between parameters in simulating LSC will differ when orography with different altitudes is applied in an AGCM. Studies on parametric sensitivity using orography with different altitudes are potentially needed to better understand the orographic effect on the non-linear interaction between parameters in simulating LSC. Moreover, the model biases over the southeast Pacific are more prominent in CGCMs than in AGCMs, because the biases are amplified via coupled feedbacks (Lin, 2007;Zheng et al., 2011). Since biases in CGCMs could originate partially from the misrepresented atmosphere due to air-sea interactions, future efforts with CGCMs should involve studying the LSC uncertainty effect induced by atmospheric moist parameters on the simulated SST and its potential role in reducing the long-standing warm SST biases.

Data Availability Statement
The simulated data used in this study is available from zenodo at https://doi.org/10.5281/zenodo.7608625. The NCAR Command Language is used to generate the figures in this study and can be available following the procedures described in http://www.ncl.ucar.edu. For details please contact ljli@mail.iap.ac.cn referencing this paper.