An effective parameter optimization with radiation balance constraint in CAM5 (version 5.3)

Uncertain parameters in physical parameterizations of general circulation models (GCMs) greatly impact model performance. In recent years, automatic parameter optimization has been introduced for tuning model performance of GCMs, but most of the optimization methods are unconstrained optimization methods under a given performance indicator. Therefore, the calibrated model may break through essential constraints that models have to keep, such as the radiation balance at the top of the model. The radiation balance is known for its importance in the conservation of model energy. In this study, an automated and efficient parameter optimization with the radiation balance constraint is presented and applied in the Community Atmospheric Model (CAM5) in terms of a synthesized performance metric using normalized mean square error of radiation, precipitation, relative humidity, and temperature. The tuned parameters are from the parameterization schemes of convection and cloud. The radiation constraint is defined as the absolute difference of the net longwave flux at the top of the model (FLNT) and the net solar flux at the top of the model (FSNT) of less than 1 W m⁻². Results show that the synthesized performance under the optimal parameters is 6.3 % better than the control run (CNTL) and the radiation imbalance is as low as 0.1 W m⁻². The proposed method provides an insight for physics-guided optimization, and it can be easily applied to optimization problems with other prerequisite constraints in GCMs.

1 Introduction

The subgrid-scale physical processes in general circulation models (GCMs) are represented by parameterization schemes, which may exist with several uncertain parameters. Inappropriate parameters can seriously affect the overall performance of the model. The Intergovernmental Panel on Climate Change Fifth Assessment Report (IPCC AR5) pointed out that studies on parameter uncertainty are critical to improve climate simulation capabilities (Mastrandrea et al., 2011). Bauer et al. (2015) also indicated that small errors in the physical parameterization schemes could lead to large-scale systematic errors. Traditionally, to achieve better performance, the uncertain parameters are tuned based on the experience of model experts and statistical analysis. This is a labor-intensive job, and the tuning results make it difficult to achieve local or global optimality in complex climate models.

To efficiently reduce parameter-introduced uncertainty, quite a few automated parameter calibration methods have been proposed. These calibration methods can be categorized into two types. One attempts to obtain the probability distributions of the parameters by likelihood and Bayesian estimation methods. Cameron et al. (1999) exploited the generalized likelihood uncertain estimation to obtain parameter ranges with a specific confidence level. An adaptive Markov Chain Monte Carlo (MCMC) was used to calibrate the uncertain parameters in the fifth-generation atmospheric general circulation model (ECHAM5) (Järvinen et al., 2010). Edwards et al. (2011) proposed a simplified procedure of Bayesian calibration to make a quantification of uncertainty in climate forecasting. This type of method has also been successfully applied to the CAM3.1 model and the third Hadley Centre Climate Model (HadCM3) (Jackson et al., 2008; Williamson et al., 2013).

The other method is to adjust parameters using optimization methods to minimize the errors between model simulations and observations, which are formulated with a given performance indicator. Many intelligent evolutionary optimization algorithms were applied to model tuning. For example, both simulated stochastic approximation annealing (SSAA) (Yang et al., 2013) and multiple very fast simulated annealing (MVFSA) (Zou et al., 2014) were used for uncertainty quantification and parameter calibration.

Both methods can consider the interaction of parameters, achieve automatic optimization, and avoid the subjectivity and experientiality of manual calibration. However, they also share high computation cost challenges due to the hundreds and thousands of required simulations. This is usually unacceptable, especially for high-resolution climate models. To overcome the computational issues, the surrogate model, which is a way to replace the real climate model with a cheaper statistical model for faster optimization, has been recently introduced. Applications of these methods in climate models include the work presented by Neelin et al. (2010) and Wang et al. (2014). However, training a precise surrogate model for a complicated climate model such as the Community Earth System Model (CESM) is very challenging. Moreover, capturing the climatic characteristics of extreme events is difficult for the cheap statistical model. To make it possible to optimize parameters efficiently and quickly in the complex and highly nonlinear earth system models, an improved simplex algorithm was presented by Zhang et al. (2015). This method can overcome the shortcomings of the traditional simplex downhill method, and the computing efficiency of the algorithm is improved compared with evolutionary optimization algorithms.

The application of various automatic parameter optimization methods in climate models has gradually received more attention; however, the optimization algorithms mentioned above are mostly unconstrained, and they lack emphasis on the physical mechanisms of the model itself. This paper takes radiation balance as an example. According to the Earth's energy conservation theory, the absorbed solar radiation is approximately equal to outgoing longwave radiation at the top of the model. Forster et al. (2007) proposed that radiative balance is critical to the Earth's system, and the bias of radiation has a big impact on the change in surface temperature. Hourdin et al. (2017) pointed out that a 1 W m⁻² change in global energy balance may result in a global mean surface temperature change of 0.5 to 1.5 K in coupled simulations. Additionally, Wild (2008) indicated that radiation biases in the GCMs may influence climate sensitivity, thus possibly distorting the prediction of future climate conditions. Lin et al. (2010) showed the importance of climate energy imbalance and stressed that long-term high-precision measurements of top of the atmosphere (TOA) radiation are necessary.

Radiation balance is critical for GCMs, but its deviation can still exceed 1 W m⁻² in some CMIP5 models (Smith et al., 2015). To better understand this problem, many studies have tried to determine the cause of radiation deviation by analyzing the influence of uncertain parameters and making corresponding adjustments. Zhao et al. (2013) concluded that cloud microphysics and emission-related parameters have statistically important impacts on the global mean net radiation flux. Qian et al. (2015) indicated that net radiation flux is very sensitive to some parameters in cloud microphysics and convection. The improvement of the simulation performance of the climatology and variability based on the radiation balance is very meaningful. However, the constrained optimization methods used to calibrate parameters with physical constraints in climate models remain to be further studied. Cheng et al. (2018) showed that penalty functions and the separation of objective and constraint methods are popular for solving constrained problems. Penalty methods encourage the search toward feasible regions by increasing the objective function value with a penalty value for the points that violate the constraints. The exterior penalty method is relatively easy to implement, and it can be widely used in various algorithms. The separation of objective and constraint is commonly used by transforming constraints into objectives, but it is limited by the convergence of the multi-objective algorithms when the optimization problem has a high computational cost.

The purpose of this paper is to propose an effective constrained optimization method and demonstrate its feasibility in the calibration of uncertain parameters under the premise of ensuring the balance of radiation. This paper is organized as follows. Section 2 describes the details of the model and experimental design. Section 3 introduces the new constrained parameter calibration method. Evaluations and analysis of the optimization results are presented in Sect. 4. The last section contains the conclusion and discussion.

Table 1 Parameters description of CAM5. The default, final optimal values by constrained and unconstrained calibrations, as well as the ranges of parameters. CAPE means the convective available potential energy.

2 Model and experiment

2.1 Model description

The model used in this study is CAM5 (release v5.3), which is the atmospheric component of the CESM 1.2.1. The dynamical core uses the finite-volume method developed by Lin and Rood (1996) and Lin (2004).

More details on CAM5 can be found in the work of Neale et al. (2010). Deep convection is handled by a parameterization scheme developed by Zhang and McFarlane (1995) with the further modifications of Richter and Rasch (2008), as well as Neale et al. (2008). The original parameterization of stratiform cloud microphysics is handled by Morrison and Gettelman (2008). Modifications of ice nucleation and ice supersaturation can be found in Gettelman et al. (2010). The parameterization of fractional stratiform condensation is described by Zhang et al. (2003) as well as Park et al. (2014). Radiation scheme uses the rapid radiative transfer method for GCMs (RRTMG) (Mlawer et al., 1997; Iacono et al., 2008).

2.2 Experiment design

Table 1 shows the parameters to be adjusted, the ranges, and the default values. These parameters were identified as sensitive to cloud and convection processes in previous studies. Qian et al. (2018) showed that deep-convection precipitation efficiency zmconv_c0_lnd and zmconv_ c0_ocn have significant impact on the variance of shortwave cloud forcing (SWCF) over the land and ocean, respectively. Thresholds of relative humidity for high and low stable clouds (cldfrc_rhminh and cldfrc_rhminl) are regarded as the important parameters for cloud and radiation (Zhang et al., 2018). In addition, the relative humidity threshold for low clouds is also one of the strongest parameters impacting the global mean precipitation, and it makes a huge contribution to the TOA net radiative fluxes in CAM5 (Qian et al., 2015). The timescale for the consumption rate of deep convective available potential energy (CAPE) (zmconv_tau) is identified as the most sensitive parameter to the convective precipitation in the Zhang–McFarlance scheme by Yang et al. (2013). The cloud ice sedimentation velocity (cldsed_ai) has a significant effect on cloud radiative forcing (Mitchell et al., 2008), and it has been identified as a high-impact parameter in sensitivity experiments related to temperature, radiation, and precipitation, etc. (Sanderson et al., 2008). The ranges of these parameters are based on previous studies (Qian et al., 2015; Zhang et al., 2018).

Table 2 The output variables used to evaluate performance metric index and the source of the corresponding observations.

The output variables used to synthesize a performance indicator are longwave cloud forcing (LWCF), SWCF, precipitation (PRECT), humidity at 850 hPa (Q850), and temperature at 850 hPa (T850), shown in Table 2. The observations of LWCF and SWCF are from CERES-EBAF (Clouds and the Earth's Radiant Energy System-Energy Balanced and Filled; Loeb et al., 2014). PRECT is from GPCP (Global Precipitation Climatology Project; Adler et al., 2003), and Q850 and T850 are from ERA-Interim, which was produced by the ECMWF (Dee et al., 2011).

In this study, we use 1.9^∘ latitude × 2.5^∘ longitude resolution CAM5 with 30 vertical layers. Each simulation is a 5-year atmospheric model intercomparison project (AMIP) from 2000 to 2004 with the observed climatological sea surface temperature (SST) and sea ice (Rayner et al., 2003). The simulations in the last 3 years are used to evaluate the synthesized performance metric and constraint.

3 Method

A constrained automatic optimization method is proposed based on Zhang et al. (2015). The synthesized metrics used to evaluate the performance of overall simulation skills are shown in Eq. (3):

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{\left({\mathit{\sigma }}_{\text{m}}^F\right)}^{\mathrm{2}}=\sum _{i=\mathrm{1}}^{I}w\left(i\right){\left(x_{\text{m}}^F\left(i\right)-x_{\text{o}}^F\left(i\right)\right)}^{\mathrm{2}},\text{(2)}& {\displaystyle }& {\displaystyle }{\left({\mathit{\sigma }}_r^F\right)}^{\mathrm{2}}=\sum _{i=\mathrm{1}}^{I}w\left(i\right){\left(x_r^F\left(i\right)-x_{\text{o}}^F\left(i\right)\right)}^{\mathrm{2}},\text{(3)}& {\displaystyle }& {\displaystyle }{\mathit{\chi }}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{N^F}}\sum _{F=\mathrm{1}}^{N^F}{\left({\displaystyle \frac{{\mathit{\sigma }}_{\text{m}}^F}{{\mathit{\sigma }}_r^F}}\right)}^{\mathrm{2}}.\end{array}$$

$({\mathit{\sigma }}_{\text{m}}^F{)}^{\mathrm{2}}$ represents a criterion for the simulation skill of the models with modified parameters; $({\mathit{\sigma }}_r^F{)}^{\mathrm{2}}$ is an evaluation of the default experiment simulation skill. If the indicator χ² is less than 1, this means that the simulation with tuned parameters is better than the control run (CNTL). The smaller the index, the better performance of model. The model outputs are represented by $x_{\text{m}}^F\left(i\right)$, and $x_{\text{o}}^F\left(i\right)$ denotes the corresponding reanalysis or observation data. The expression $x_r^F\left(i\right)$ represents model outputs from the CNTL. The weight of the different grids on the sphere is represented by ω. I denotes the total number of grids in the model. The number of the output variables in the performance index is represented by N^F.

The radiation balance is defined as the absolute value of the difference between net solar flux (FSNT) and net longwave flux (FLNT) in climatology at the top of the model of less than 1 W m⁻², which is the maximum deviation of radiation observations in the decade before 2014 (Trenberth et al., 2014).

Coupled with the radiation balance constraint, the optimization problem of this study can be expressed as Eqs. (4) and (5):

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }\min {\mathit{\chi }}^{\mathrm{2}},\text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{subject}\mathrm{to}\mathrm{ABS}\left({\text{FSNT}}_{\text{m}}-{\text{FLNT}}_{\text{m}}\right)<\mathrm{1}.\end{array}$$

Converting the unconstrained problem into a constrained problem using the penalty function method, it can be transformed into the augmentation function as Eq. (6):

$$\begin{array}{}\text{(6)}& F\left(x\right)={\mathit{\chi }}^{\mathrm{2}}+\mathit{\beta }\times \text{ABS}\left({\text{FSNT}}_{\text{m}}-{\text{FLNT}}_{\text{m}}\right).\end{array}$$

The penalty factor β for the constraint in Eq. (6) is chosen to be 10 000 if the constraint in Eq. (5) is not satisfied; otherwise it is equal to 0. The purpose of this choice is to optimize the search space to avoid the possibility of radiation imbalance. This penalty function method is easy and effective when dealing with this tightly constrained optimization.

We use the improved simplex downhill method, proposed by Zhang et al. (2015), to optimize the augment function. Firstly, the single-parameter perturbation sample method (SPP) is used to obtain several better initial values, while ensuring that the initial geometry of simplex downhill is well-conditioned. The initial value preprocessing mechanism ensures that the algorithm starts from a good basis. This is important for the simplex algorithm, which may easily fall into local optimum. Next, the simplex downhill algorithm is applied to search for better performance.

Figure 1 The change in augmentation function F(x) across the optimization iterations. The x axis is the number of iterations. The y axis is the value of F(x) in Eq. (6). The black line shows the value of F(x) in a given iteration step, while the red line shows the best F(x) value up to the current iteration step.

Figure 2 Comparison of results between the constrained optimization algorithm and the unconstrained optimization algorithm. The 15 red squares and 15 black triangles are optimized solutions found by the unconstrained optimization algorithm and constrained algorithms, respectively. The blue diamond is the result of the CNTL experiment. The x axis is the synthesized metric index in Eq. (3). The y axis is the radiation bias at the top of the model.

Figure 3 Taylor diagram of the climate mean state of each output variable from 2002 to 2004 between the model run with optimal parameters and the CNTL run. The number (1) in the diagram stands for EXP, and (2) stands for CNTL.

F(x) gradually converges as shown in Fig. 1. There are some cases in which the radiation balance is not satisfied at the beginning of optimization. However, as the iteration step increases, the search space of the algorithm is constrained within the feasible range. The goal is then to make the synthesized performance metric smaller. In addition, a comparative experiment with unconstrained algorithms is done to verify our doubts about unconstrained methods. Figure 2 shows the performance indices and radiation deviations corresponding to the first 15 solutions after the two algorithms converge. The constrained optimization algorithm can find solutions that are more radiation-balanced; however, the final solution metrics are not as good as the unconstrained optimization algorithm. Compared to the CNTL experiment, we can find quite a few solutions with better metrics and smaller radiation biases.

4 Result

4.1 The optimal model

The best uncertainty parameters obtained by the unconstrained optimization method optimize the overall performance of the simulation by 10.1 %, but they have a radiation deviation of up to 3.8 W m⁻². When considering the converged constrained optimization algorithm, the optimal parameters can improve the model performance by 6.3 %, and the radiation imbalance is as low as 0.1 W m⁻². The corresponding results of the optimal solutions with the two methods are shown in Table 3. Both unconstrained optimization and constrained optimization can further improve the simulation performance, but unconstrained optimization may encounter an optimal solution that does not satisfy the radiation balance, thus leading to meaningless optimization. The optimization results discussed below are based on the proposed constrained optimization method.

Table 3 Synthesized performance metric index and radiation bias in the CNTL run and the optimal model run with unconstrained and constrained methods.

Figure 4 Meridional distribution of the difference between EXP/CNTL and observed data of (a) LWCF, (b) SWCF, (c) PRECT, (d) Q850, and (e) T850. The position of the dark blue line is 0; the red and black solid lines represent the difference between EXP ∕ CNTL and the observations.

Table 4 Performance metric index of each variable in the optimal model run with unconstrained and constrained methods.

The optimization of each output variable is shown in Table 4. In addition, a Taylor diagram is used to estimate the model performance through the standard deviation and correlation (Fig. 3). By combining the results of Table 4 and Fig. 3, it can be concluded that SWCF and Q850 receive most optimization, as they achieve a better performance index. Also, compared to the default experiment, their standard deviations have improved. Table 5 shows the standard deviations of the variables, which are important for the model but not used as evaluation criteria. It is noteworthy that they are also close to the default experiment.

Table 5 The percentage of standard deviation of the eight fields between the CNTL run and the optimal model run with constrained optimization according to the corresponding observations.

Figure 5 The spatial distribution of TOA SW cloud forcing of (a) observation, (b) CNTL–observation, (c) EXP, (d) EXP–observation, (e) CNTL, and (f) EXP–CNTL.

For a more comprehensive analysis of the spatial variation of the output variables, the zonal distribution of the difference between the control (labeled as CNTL)/the optimized (labeled as EXP) simulations and observations of all metric variables are shown in Fig. 4. SWCF and Q850 have been obviously improved over low and middle latitudes, but the changes in PRECT and T850 are not particularly notable. Further, LWCF only showed significant improvement near the Equator, and it slightly deteriorated over the middle and high latitudes.

4.2 Interpretation of the results

The optimized parameters values are provided in the “Constrained tune” column of Table 1. The deep-convection precipitation efficiency over land and ocean is reduced relative to the default values. The timescale for the consumption rate of CAPE for deep convection is smaller than the default value, and both relative humidity thresholds for high and low clouds are increased. Additionally, the sedimentation velocity of cloud ice is larger. Next, we will explain how the changes in these parameters are related to the results of the simulations.

Figure 6 The spatial distribution of specific humidity at 850 hPa of (a) observation, (b) CNTL–observation, (c) EXP, (d) EXP–observation, (e) CNTL, and (f) EXP–CNTL.

The relative humidity threshold for low clouds is larger in optimization experiments than the default value, which will obviously lead to the decrease in low-cloud fraction. The decreased low-cloud fraction is consistent with the increase in SWCF. The CNTL experiment has excelled at simulating the spatial distribution of SWCF (Fig. 5c), but it has a negative bias over the ocean in the low latitudes, where the improvement is significant in the optimal experiment.

The zonal mean specific humidity at 850 hPa is significantly improved, and its spatial distribution is presented in Fig. 6. In the optimal experiment, the atmosphere is drier in the tropics and middle latitudes, which is closer to the observation than the CNTL experiment. Meanwhile, the middle to low troposphere is also slightly drier in these areas (Fig. 7), which may be related to the increased convective precipitation. A quasi-equilibrium closure is used in the deep-convection scheme in CAM5, which is based on CAPE. The adjustment timescale represents the denominator of the cloud bottom convective mass flux. When the timescale is shorter with less changed CAPE, the increased cloud-base mass flux would help to enhance the convective precipitation. Additionally, compared to the CNTL experiment, the lower troposphere gets warmer and the middle troposphere is colder, which exacerbates the instability of the temperature structure (Fig. 8) and leads to more convective precipitation. The spatial distribution of convective precipitation over the tropics where convection occurs most frequently can be seen in Fig. 9. The increase in convective precipitation may be related to the decrease in specific humidity at 850 hPa. However, the increase in total precipitation is not particularly significant and is dominated by the changes in convective precipitation. The main reason is likely associated with the decreased precipitation efficiency parameters, which could reduce the convective precipitation as compensation. Therefore, the decreases in precipitation efficiency partially offset the precipitation change caused by tau and temperature structure.

Figure 7 Pressure–latitude distributions of specific humidity of (a) EXP–OBS, (b) CNTL–OBS, and (c) EXP–CNTL.

Figure 8 Pressure–latitude distributions of temperature of (a) EXP–OBS, (b) CNTL–OBS, and (c) EXP–CNTL.

Figure 9 The spatial distribution of convective precipitation over the tropics of EXP (a), CNTL (b), and EXP–CNTL (c).

Figure 10 Pressure–latitude distributions of cloud fraction of EXP (a), CNTL (b), and EXP–CNTL (c).

It is difficult for all variables to be optimized, due to the strong interaction among parameters and the complex relationship among output variables. The simulations of T850 between optimal and CNTL experiments are very similar. It is likely the result of the combined effects of all relevant parameterizations. In the optimal experiment, LWCF is closer to the observation in the tropics, but it becomes slightly smaller at middle to high latitudes compared to the CNTL experiment, which implies the larger bias. The relative humidity threshold for high clouds and the sedimentation velocity of ice crystals are correspondingly increased, and both of them would lead to the reduction in high clouds. High-cloud fraction changes compared to the CNTL experiment can be seen in Fig. 10c. The reduced high cloud is consistent with the reduction in LWCF. Cloud changes also inevitably affect SWCF. It can be seen that the middle cloud has increased relative to the default experiment (Fig. 10c), and the increase in the middle cloud may be related to the decrease in precipitation efficiency over the ocean.

Note that three of six parameters hit their lowest allowable limit with the TOA balance constraint. We found that the incoming shortwave radiation flux is more sensitive to tuning parameters than the outgoing longwave radiation flux. Thus, to reduce the TOA imbalance and keep the reasonable model performance, the shortwave radiation flux should be reduced largely via increasing low-cloud fraction and liquid water content. These three variables can help achieve this by setting to the lowest bounds. This suggests that getting both the TOA balance and reasonable model performance is a relatively complex and difficult problem due to model structure problems, as pointed out by Qian et al. (2018) and Yang et al. (2019). Meanwhile, finding out how to pick parameters with a similar sensitivity to both longwave and shortwave radiation flux might be a potential approach to overcome the bound limit and it warrants further studies.

In conclusion, the increase in SWCF is consistent with the decrease in cloud fraction for the sake of a larger relative humidity threshold of low clouds. Changes in the Q850 are related to increased convective precipitation. Precipitation only slightly increases in the tropics, and the global total precipitation has changed very little, which is related to the comprehensive effect of the changes in the convection adjustment timescale, the precipitation efficiency parameter, and the vertical temperature structure. T850 simulated by the optimization experiment is similar to the default experiment. The reduced LWCF is related to the decreased high clouds caused by the increased relative humidity threshold for high clouds and the increased sedimentation velocity of ice crystals.

5 Conclusion and discussion

Radiation balance is a crucial factor for the long-term energy balance of GCMs, but it has not received enough attention in automatic parameter optimization. First of all, this paper points out that the previous parameter optimization algorithms do not consider radiation balance a necessary condition, and the obtained optimization parameters are likely to break this important physical constraint, which may lead to unacceptable calibrated parameters. Thus we propose an efficient constrained automatic optimization algorithm to calibrate the uncertainty parameters in CAM5 with the constraint of the absolute value of the difference of net solar flux and net longwave flux at the top of the model (less than 1 W m⁻²). In the parameter calibration, we use the comprehensive performance with the five fields of LWCF, SWCF, PRECT, Q850, and T850 as the performance indicators. We choose the uncertain parameters in cloud and convection parameterizations, including the deep-convection precipitation efficiency over land and ocean, thresholds of relative humidity for stable high and low clouds, the timescale for consumption rate deep CAPE, and the ice falling speed. Each simulation in our optimization experiments is a 5-year AMIP experiment forced with prescribed seasonal climatology of SST and sea ice.

The optimal parameters found by our method can increase the overall performance of the model by 6.3 %, and the radiation imbalance is as low as 0.1 W m⁻². The most optimized variables are SWCF and Q850. The increase in SWCF is consistent with the decrease in cloud water due to a larger relative humidity threshold value for low clouds. The reduction in the Q850 in the troposphere may be related to the increase in convective precipitation. The change in global total precipitation is not particularly obvious, which is likely the comprehensive effect of the changes in the convection adjustment timescale, the precipitation efficiency parameter, and the structure of temperature over the troposphere. The change in T850 is very small, and the result is slightly better than that of the default experiment. Meanwhile, under the constraint of energy balance, LWCF has deteriorated in the middle and high latitudes. This also reflects some issues that may exist in the structure of the model.

The unconstrained optimization methods calibrate the uncertain parameters in climate models without a consideration of the principles that model have to hold, this creates challenges in maintaining the physics constraints and improving the structure of models. Perhaps a more physics-guided optimization is a better way to understand the principles of climate systems and to best use these principles in tuning processes. In the future, we will apply this method to coupled models, where the radiation imbalance has a more significant impact on long-term simulation stability. In addition, we will also try to introduce more constraints, such as the surface energy balance, into automatic parameter calibration.