Assessing the performance of climate change simulation results from BESM-OA2.5 compared with a CMIP5 model ensemble

The main features of climate change patterns, as simulated by the coupled ocean–atmosphere version 2.5 of the Brazilian Earth System Model (BESM), are compared with those of 25 other CMIP5 models, focusing on temperature, precipitation, atmospheric circulation, and radiative feedbacks. The climate sensitivity to quadrupling the atmospheric CO₂ concentration was investigated via two methods: linear regression (Gregory et al., 2004) and radiative kernels (Soden and Held, 2006; Soden et al., 2008). Radiative kernels from both the National Center for Atmospheric Research (NCAR) and the Geophysical Fluid Dynamics Laboratory (GFDL) were used to decompose the climate feedback responses of the CMIP5 models and BESM into different processes. By applying the linear regression method for equilibrium climate sensitivity (ECS) estimation, we obtained a BESM value close to the ensemble mean value. This study reveals that the BESM simulations yield zonally average feedbacks, as estimated from radiative kernels, that lie within the ensemble standard deviation. Exceptions were found in the high latitudes of the Northern Hemisphere and over the ocean near Antarctica, where BESM showed values for lapse rate, humidity feedback, and albedo that were marginally outside the standard deviation of the values from the CMIP5 multi-model ensemble. For those areas, BESM also featured a strong positive cloud feedback that appeared as an outlier compared with all analyzed models. However, BESM showed physically consistent changes in the temperature, precipitation, and atmospheric circulation patterns relative to the CMIP5 ensemble mean.

1 Introduction

The effects of increased atmospheric CO₂ concentrations on the climate system have been studied over the last 120 years (Arrhenius, 1896; Callendar, 1938; Plass, 1956; Kaplan, 1960; Manabe and Wetherald, 1967, 1975; Manabe and Stouffer, 1980; IPCC, 2007, 2013; Pincus et al., 2016; Good et al., 2016,  and many others). The human-induced increase in atmospheric greenhouse gas (GHG) concentrations, sometimes given as the CO₂-equivalent concentration, contributes to a radiation imbalance at the top of the atmosphere (TOA) that causes less outgoing radiation to leave the Earth system. The trapping of infrared radiation results in a temperature rise at the lower levels of the troposphere and an increase in ocean heat content. In addition, the increased GHG concentration can trigger climate feedback processes that either amplify or damp the initial radiative perturbation (Cubasch and Cess, 1990). Earth system models (ESMs) are the most advanced tools available for analyzing the coupled climate system (atmosphere, ocean, land, and ice) physical processes and their interactions, although even these models still exhibit important uncertainties in their projections of climate change (IPCC, 2013).

The equilibrium global mean surface temperature change induced by doubling the CO₂ concentration in the atmosphere, referred to as the equilibrium climate sensitivity (ECS), remains a centrally important measure of a model's climate response to CO₂ forcing. In the fifth Intergovernmental Panel on Climate Change (IPCC) assessment report (AR5), climate model ECS estimates range from 2.0 to 4.5 K. For more than 40 years, this inter-model spread has been considered one of the most critical uncertainties for the evaluation of future climate change (IPCC, 2013). This inter-model dispersion arises principally from differences in how the climate models simulate climate feedback processes. Among them, cloud feedback constitutes the largest source of variation in the climate sensitivity estimates (Cess et al., 1989, 1990; Dufresne and Bony, 2008; Vial et al., 2013; Caldwell et al., 2016).

Beyond the ECS, the response of precipitation patterns to anthropogenic GHG emissions is a topic of great interest in climate science given their potential effects on both societies and ecosystems. Changes in precipitation can generally be decomposed into two processes: a thermodynamic component due to increased moisture with no circulation change and a dynamic component due to circulation change with no moisture change (Bony et al., 2006; Seager et al., 2010). The thermodynamic component gives rise to the well-known “wet-gets-wetter” and “dry-gets-drier” patterns of precipitation change first described by Held and Soden (2006), which are associated with the Clausius–Clapeyron relation (i.e., a temperature-dependent exponential increase in the saturation specific humidity) (Marvel and Bonfils, 2013). The dynamic component, which is associated with circulation change, sometimes yields strong deviations from the thermodynamic pattern of precipitation, and this component is known to dominate the inter-model deviation in estimates of total precipitation due to uncertainties in the regional circulation change (Xie et al., 2015).

In this study, we assess the main features of climate change patterns as simulated by the Brazilian Earth System Model ocean–atmosphere coupled version 2.5 (BESM-OA2.5), with a focus on temperature (climate sensitivity and feedbacks), precipitation, and atmospheric circulation. The recent development of the BESM-OA2.5 has been a coordinated effort at the National Institute for Space Research (INPE) in Brazil intended to advance the understanding of the causes of global and regional climate changes and their effects on the socioeconomic sector. We evaluate how the BESM-simulated climate change prediction compares with those from Coupled Model Intercomparison Project Phase 5 (CMIP5) models, also discussing peculiarities in the BESM-OA2.5 climate response. This paper is structured as follows: Sect. 2 presents the description of the new features of BESM-OA2.5, Sect. 3 presents the methodology, Sect. 4 presents the results, and Sect. 5 presents the summary and conclusions.

2 Model description

2.1 BESM-OA2.5

BESM-OA2.5 is the result of coupling the Center for Weather Forecast and Climate Studies (CPTEC/INPE) Brazilian Atmospheric Model (BAM; Figueroa et al., 2016) and the Geophysical Fluid Dynamics Laboratory (GFDL) Modular Ocean Model version 4p1 (Griffies et al., 2004) via the Flexible Modular System (FMS) (also from GFDL). The dynamical core and physical parameterizations of the atmospheric component of BESM-OA2.5 are the same as those discussed in Veiga et al. (2019). BAM is a hydrostatic model, with its dynamical core based on the spectral transform method, which employs global spherical harmonic basis functions. The Eulerian advection scheme option is used in this study, with a two-time-level semi-Lagrangian scheme for the transport of moisture and microphysics prognostic variables, which are carried out completely within the model grid space. Simplified fast physical parameterizations are used here to increase the computational efficiency of long integrations, thus resulting in a decreased computational demand compared with that required by the operational Numerical Weather Prediction (NWP) model. A summary of the main differences in the physical parameterizations between BAM (as used in this paper) and the BAM NWP operational model is provided in Table 1. The dynamical equations in BAM are discretized following a spectral transform with horizontal resolution truncated at triangular wavenumber 62 (an equivalent grid size of approximately 1.875^∘) and 28 layers unevenly spaced in the vertical sigma coordinate with the top level at approximately 2.73 hPa. The oceanic component uses a tripolar grid at a horizontal resolution of 1^∘ in longitude, and in the latitudinal direction the grid spacing is 1∕4^∘ between 10^∘ S and 10^∘ N, decreasing uniformly to 1^∘ at 45^∘ and to 2^∘ at 90^∘ in both hemispheres. The ocean grid has 50 vertical levels with a 10 m resolution in the upper 220 m, decreasing gradually to approximately 370 m at deeper levels.

(Iacono et al., 2008)(Tarasova and Fomin, 2007)(Iacono et al., 2008)Harshvardhan et al. (1987)(Morrison et al., 2005)Ferrier et al. (2002)Foley et al., 1996Kubota, 2012(Xue et al., 1991)Mellor and Yamada (1982)Holtslag and Boville (1993)(Park and Bretherton, 2009)Tiedtke (1984)Grell and Dévényi (2002)Grell and Dévényi (2002)Webster et al. (2003)Alpert et al. (1988)Slingo (1987)

Table 1 Atmospheric physical parameterizations used in BAM (Figueroa et al., 2016) BESM-OA2.5.

Veiga et al. (2019) showed that BESM-OA2.5 can capture the general mean climate state; however, substantial biases appeared in the simulation associated with a double Intertropical Convergence Zone (ITCZ) over the Pacific and Atlantic oceans and regional biases in the precipitation over the Amazon and Indian regions. BESM-OA2.5 can also reproduce the most important large-scale interannual and decadal climate variability patterns. The Atlantic Meridional Mode (AMM) (Nobre and Shukla, 1996) is well simulated by the model in terms of its spatial pattern and temporal variability, whereas this mode is poorly represented in most CMIP5 models (IPCC, 2013; Liu et al., 2013; Richter et al., 2014; Amaya et al., 2017). The maximum strength of the Atlantic Meridional Overturning Circulation (AMOC) represented by BESM-OA2.5 is 14 Sv, which is lower than the value determined within the RAPID project (17 Sv; McCarthy et al., 2015) but in the range of the observed root mean square variability, and this value is comparable to the ensemble AMOC simulated by the CMIP5 models. Moreover, the spatial structures of both the North Atlantic Oscillation (NAO) and the Pacific Decadal Oscillation (PDO) variability are well captured (Veiga et al., 2019).

2.2 Comparison to a previous model version

The recently developed BESM-OA2.5 is an advancement of BESM-OA2.3, which was presented by Nobre et al. (2013). The main differences between BESM-OA2.5 and the previous version (BESM-OA2.3) lie in the atmospheric model and how some surface layer variables are estimated. The total energy balance at the TOA is better represented in BESM-OA2.5 than in BESM-OA2.3, with a reduced global mean bias of approximately −4 W m⁻² compared with −20 W m⁻² for the latter. It should be noted that BESM-OA2.5 has a new set of parameterizations, mainly related to an improved representation of microphysical processes. For instance, the precipitation in the previous model was parameterized only in terms of the large-scale condensation. Moreover, BESM-OA2.5 underwent improvements in the representation of wind, humidity, and temperature in the surface layer with the use of the similarity functions method presented by Jiménez et al. (2012). Based on Monin–Obukhov theory, the wind (u_10 m), humidity (q_2 m), and temperature (θ_2 m) are estimated from the values of the first atmospheric model level and the surface, as described in Eqs. (24), (25), and (26) of Jiménez et al. (2012). Furthermore, the similarity functions ψ_m and ψ_h depend on the stability regimes (Businger et al., 1971). For BESM-OA2.5, those regimes are associated with stable ($\mathit{\zeta }/L>\mathrm{0}$) and unstable ($\mathit{\zeta }/L\le \mathrm{0}$) conditions (Arya, 1988). These diagnostic variables are important for BESM because they are used in the ocean–atmosphere coupling strategy.

Both versions reproduce the main climate variability, particularly over the Atlantic, such as the AMOC and the AMM, but simulate a weak El Niño–Southern Oscillation (ENSO) interannual variability over the equatorial Pacific (Nobre et al., 2013; Veiga et al., 2019). Concerning the general mean present-day climate state, BESM-OA2.5 shows improvements in reproducing the Intertropical Convergence Zone (ITCZ), and it reduces both the global precipitation root mean square error (RMSE) and the sea surface temperature (SST) RMSE compared with those modeled by BESM-OA2.3.

Global simulations that were 1 year long and 6-hourly outputs were performed with BAM configured with surface layer schemes based on Arya (1988) and Jiménez et al. (2012), here called BAM-Arya (the original scheme) and BAM-Jimenez (the new scheme), respectively. The normalized RMSE was computed with respect to the reanalysis NCEP-DOE (National Centers for Environmental Prediction – Department of Energy) version 2 (Kanamitsu et al., 2002). The normalized RMSE of the wind at 10 m and the temperature and humidity at 2 m for the two surface layer schemes was investigated. Consistent improvements of BAM-Jimenez relative to BAM-Arya were noted in all three variables over the oceanic regions. The normalized RMSE analysis over the continents yielded less consistent results, with an improved BAM-Jimenez representation of both winds and temperature but with an inferior representation of the humidity field (figures not shown).

3 Methodology

3.1 Experimental design

For this study, climate simulations were performed using BESM-OA2.5 (hereinafter BESM) for the piControl (preindustrial control scenario, run for 300 years with atmospheric CO₂ concentration invariant at 274 ppmv) and abrupt4xCO2 (run for 150 years after the abrupt quadrupling of atmospheric CO₂ at year 151 of the piControl simulation) scenarios; therefore, both experiments were run in parallel for 150 years. These two scenarios are commonly employed in CMIP5 studies for climate sensitivity assessment (Taylor et al., 2012; Eyring et al., 2016). Climate change is evaluated as the difference between the abrupt4xCO2 and piControl experiments. In addition, BESM's results were compared with a selection of 25 CMIP5 models listed in Table 2. All models, including BESM, were interpolated at 2.5^∘ × 2.5^∘ longitude–latitude horizontal resolution. All CMIP5 model data are available from the Earth System Grid Federation (ESGF).

Table 2 Models belonging to the CMIP5 ensemble used in this study.

3.2 Climate change sensitivity estimates

Here we estimate the climate feedback with two different methods using either a regression according to Gregory et al. (2004) or radiative kernels (Soden et al., 2004, 2008). The Gregory method is more straightforward computationally; however, it returns only a global mean value. Moreover, the ECS can be estimated with this method. On the other hand, it is possible to obtain the seasonal feedback for every latitude–longitude point with the radiative kernel method; furthermore, the feedback can be decomposed into different processes.

3.2.1 Linear forcing–feedback regression analysis

The regression method for computing the thermal response to radiative forcing was applied for 26 CMIP5 models, including BESM. The method consists of linear regression between the annual change (considering abrupt4xCO2 minus piControl) in the global mean near-surface temperature (ΔT_as) and the net radiation flux change (ΔR) at the TOA.

If G is the radiative forcing imposed on the climate system (here associated with an abrupt increase in atmospheric CO₂ concentration) and ΔR is the resulting radiative imbalance in the global mean net radiative budget at the TOA, then at any time, the response of the climate system to this radiative imbalance corresponds to the radiative forcing according to the following equation:

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }R=\mathit{\lambda }\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}+G,\end{array}$$

where λ (<0) is the climate feedback parameter and $\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}$ is the global mean near-surface temperature change. In a sufficiently long simulation (coupled atmosphere–ocean models take millennia), the climate system reaches a new equilibrium when ΔR=0. As G can be approximated via backward regression towards $\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}=\mathrm{0}$, ECS can be estimated as $\mathrm{ECS}=-G/\mathit{\lambda }$. As the ECS is the theoretical equilibrium temperature for doubling CO₂, in a quadrupling of CO₂ it is common to divide the result derived from 4xCO₂ simulations by 2 (Andrews et al., 2012).

By using this linear forcing response framework, we can estimate climate sensitivity, radiative forcing, and feedback parameters following the method proposed by Gregory et al. (2004). The values of λ (slope) and G (y intercept) are estimated via the ordinary least-squares regression of the global annual mean of ΔR against $\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}$ under all-sky conditions. Using the same linear technique, we decompose the feedback parameter into shortwave (SW) and longwave (LW) radiation components, and we extract the clear-sky radiative flux components from the BESM and CMIP databases to estimate the cloud radiative forcing or cloud radiative effect (ΔCRE) defined as the difference between the all-sky and clear-sky feedback parameters (Andrews et al., 2012). Estimates of G, λ, ΔCRE, and ECS for all models are presented in the next section.

3.2.2 Separating individual climate feedbacks using radiative kernels

The radiative kernel technique (as in Soden and Held, 2006; Soden et al., 2008; Vial et al., 2013) is used to separate the feedback parameter λ into contributions from the temperature response (λ_T), water vapor (λ_lnq), surface albedo (λ_a), and cloud (λ_c) feedbacks plus a residual term (Re) (Vial et al., 2013) as expressed in Eq. (2).

$$\begin{array}{}\text{(2)}& \mathit{\lambda }={\mathit{\lambda }}_T+{\mathit{\lambda }}_{\ln q}+{\mathit{\lambda }}_a+{\mathit{\lambda }}_{\mathrm{c}}+\mathrm{Re}\end{array}$$

We used both GFDL (Soden et al., 2008) and National Center for Atmospheric Research (NCAR) (Shell et al., 2008) radiative kernels to estimate climate feedbacks. Such radiative kernels represent the impact on the radiative balance at the TOA via arbitrary increases in the atmospheric temperature, albedo, and specific humidity. For calculating the temperature kernel, an increment of 1 K is added for all model levels (including the surface). For the albedo kernel, the albedo value is increased by 0.01 (1 %). Finally, for the water vapor kernel, the specific humidity is increased by a value equivalent to a 1 K atmosphere temperature increase, with the relative humidity remaining constant. Furthermore, we used the ln q instead of q, considering the quasi-proportionality of the absorption of radiation by water vapor to ln q.

Following Soden and Held (2006), Jonko et al. (2013), and Vial et al. (2013), we decompose the total feedback parameter (λ) into contributions from λ_T, λ_lnq, λ_a, and λ_c as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathit{\lambda }=& \sum _x{\mathit{\lambda }}_x+\mathrm{Re}=\sum _x{\displaystyle \frac{\partial R}{\partial x}}{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}}}+\mathrm{Re}\\ =& \sum _x{K}_x{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}}}+\mathrm{Re}\\ \mathit{\lambda }=& \left(K_{T_{\mathrm{s}}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{s}}}{\mathrm{d}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}}}+K_T{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}}}\right)+\left(K_{\ln q}{\displaystyle \frac{\mathrm{dln}q}{\mathrm{d}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}}}\right)\\ +& \left(K_a{\displaystyle \frac{\mathrm{d}a}{\mathrm{d}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}}}\right)+{\mathit{\lambda }}_{\mathrm{c}}+\mathrm{Re}.\end{array}\end{array}$$

The temperature feedback has been separated into the Planck feedback (the vertically uniform tropospheric warming equal to the surface warming) and the lapse rate feedback (the deviation from the tropospheric uniform warming):

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathit{\lambda }}_T={\mathit{\lambda }}_p+{\mathit{\lambda }}_{\mathrm{lr}}& =\left(K_{T_{\mathrm{s}}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{s}}}{\mathrm{d}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}}}+K_T{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{s}}}{\mathrm{d}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}}}\right)\\ & +\left(K_T{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}}}-K_T{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{s}}}{\mathrm{d}{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}}}\right).\end{array}\end{array}$$

In Eq. (3), K_x (the radiative kernel for a variable x) and x (temperature (T_s and T; K), the natural logarithm of humidity (lnq; kg kg⁻¹), and the albedo, which is represented by “a” and is dimensionless) are functions of the longitude, latitude, and pressure vertical coordinates in the monthly climatology. To obtain tropospheric averages, the water vapor and temperature feedbacks are vertically integrated from the surface up to the tropopause, defined as 100 hPa at the Equator, and varying linearly to 300 hPa at the poles. The stratospheric temperature and water changes are not accounted for in calculating the feedbacks, and they are shifted to the residuum.

It is worth noting that in the regression method, the radiative feedback is consistent with the actual radiative transfer scheme used in the climate model, while in the radiative kernel, the feedback is not necessarily consistent. In fact, the kernel is obtained from another climate model that is not among the CMIP5 models analyzed. Model intercomparison is easily achieved via this method, as the same kernel can be applied to all models (Soden and Held, 2006; Soden et al., 2008). However, the resulting kernel-derived feedbacks can only be assumed to reflect the actual feedback in the considered models under the premise of small differences between the radiative transfer codes.

Due to the nonlinearities involving clouds and net radiation at the TOA (Soden et al., 2008), the cloud feedback is not calculated directly from these radiative kernels, which represents one of the key limitations of the kernel method. Instead, the cloud feedback is approximated using the cloud radiative forcing (ΔCRE) corrected by removing the non-cloud feedback effects as in Soden et al. (2004, 2008). After the calculation of non-cloud feedbacks for both all-sky and clear-sky (superscript cs) conditions, we thus estimate the cloud feedback (λ_c) as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\mathrm{\Delta }\mathrm{CRE}& =\mathrm{\Delta }R-\mathrm{\Delta }{R}^{\mathrm{cs}}\\ \mathrm{\Delta }{\mathrm{CRE}}_k& =(G-G^{\mathrm{cs}}{)}_{{\mathrm{CO}}_{\mathrm{2}}}-\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}\sum _x({\mathit{\lambda }}_x-{\mathit{\lambda }}_x^{\mathrm{cs}})\\ \mathrm{\Delta }{\mathrm{CRE}}_{\mathrm{a}}& =\mathrm{\Delta }\mathrm{CRE}-\mathrm{\Delta }{\mathrm{CRE}}_k\\ {\mathit{\lambda }}_{\mathrm{c}}& ={\displaystyle \frac{\mathrm{\Delta }{\mathrm{CRE}}_{\mathrm{a}}}{{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}}},\end{array}\end{array}$$

where ΔR^cl is the clear-sky net radiation flux at the TOA. Following Soden et al. (2008), $(G-G_{\mathrm{cl}}{)}_{{\mathrm{CO}}_{\mathrm{2}}}$ was defined as 2×0.69 W m⁻². The index k represents the change in the ΔCRE due to the non-cloud feedbacks, while the index “a” means the adjusted ΔCRE. Finally, a 30-year mean relative to the period from the 120th to 150th years of each scenario was used for all feedback estimations. It is worth mentioning that the term ${\sum }_x({\mathit{\lambda }}_x-{\mathit{\lambda }}_x^{\mathrm{cs}})$, introduced in Eq. (5), represents the cloud masking on the non-cloud radiative feedbacks. It can be physically interpreted as the differences in the distribution of the temperature, vapor water, and albedo between an all-sky and a clear-sky atmosphere (Soden et al., 2004).

3.3 Changes in the atmospheric circulation and precipitation

Monthly mean climatologies were computed for the last 30 years of the piControl and abrupt4xCO2 runs, and the projected climate response to CO₂ increase was evaluated from the difference between these abrupt4xCO2 and piControl monthly mean climatologies. The statistical significance of this difference was calculated based on the Student's t test with a significance level of 90 %. Furthermore, to evaluate how similar two spatial patterns are, we used the spatial inner product calculated as $\sum (A_i\cdot B_i)/\left(\right|A|\cdot |B\left|\right)$, where A and B are the 2-D variables and i is the spatial index related to their latitude–longitude coordinates.

Figure 1 Annual global mean linear regression between $\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}$ and (a) net radiation, (b) $\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}$ and ΔLW (clear sky), and (c) $\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}$ and ΔSW (clear sky) for BESM-OA2.5.

4 Results

4.1 G, λ, ΔCRE, and ECS estimated via the Gregory method

Figure 1 shows the linear regressions of ΔR , ΔLW (clear sky), and ΔSW (clear sky) against $\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}$ from BESM. The linear regressions based on all-sky radiative flux are used to estimate ECS, G, and λ, while the regressions based on clear-sky data are used to obtain ΔCRE (as mentioned in the previous section). BESM features G=8.62 W m⁻², $\mathit{\lambda }=-\mathrm{1.45}$ W m⁻² K⁻¹, $\mathrm{\Delta }\mathrm{CRE}=-\mathrm{0.13}$ W m⁻² K⁻¹, and ECS=2.96 K.

The parameters G, λ, ΔCRE, and ECS as computed for all models are shown in Table 3. The climate sensitivities of the 26 CMIP5 coupled models (including BESM-OA2.5) were compiled in an extension of the previous work by Andrews et al. (2012), who evaluated 15 CMIP5 coupled models. In the present work, we added the following models: ACCESS1-0, ACCESS1-3, bcc-csm1-1, BESM-OA2.5, BNU-ESM, CCSM4, FGOALS-g2, FGOALS-s2, GISS-E2-H, GISS-E2-R, and inmcm4. In Andrews et al. (2012) the ECS ranges from 2.07 to 4.74 K for the 15 models analyzed there, which is largely confirmed by our analysis. The small differences can possibly be attributed to the interpolation of the data. G and λ vary from 5.01 to 8.95 W m⁻² and from −1.66 to −0.60 W m⁻² K⁻¹, respectively. The inter-model spread in G among the models is due to differences in the radiative codes used and the rapid adjustment processes in the troposphere and at the surface (Collins et al., 2006; Gregory and Webb, 2008; Andrews and Forster, 2008). The spread in the ECS is more robustly influenced by λ than by G (Fig. 2), as was also suggested by Andrews et al. (2012). The correlation coefficient between the ECS and λ is −0.82, which is significant at a 1 % significance level (Fig. 2b). On the other hand, the correlation between the ECS and G is −0.01, which is not statistically significant (Fig. 2a). Thus, inter-model variation in the balance of feedbacks explains the dispersion in the ECS better than the initial radiative imbalance triggered by the CO₂ increase (related to G). Although BESM yielded one of the highest G values among all of the CMIP5 models, it showed a warming response to CO₂ doubling that is well within the range of 3.30±0.76 K as presented by the ensemble models.

Figure 2 (a) Equilibrium climate sensitivity (ECS, in red) and radiative forcing (G, in blue) values with ECS values increasing from left to right; (b) ECS (red) and climate sensitivity (λ, in green) with ECS values increasing from left to right.

Table 3 CO₂ forcing (W m⁻²) (G), net feedback (W m⁻² K⁻¹) (λ), climate response (W m⁻²) (ΔCRE), and equilibrium climate sensitivity (K) (ECS) values.

The ΔCRE from BESM is −0.13 W m⁻² K⁻¹, while the CMIP5 models yield ΔCRE values ranging from −0.50 to 0.70 W m⁻² K⁻¹. Unlike the ΔCRE_a, this term does not consider the masking effects of clouds as estimated via the radiative kernel method (Eq. 5). Therefore, the ΔCRE cannot be interpreted to reflect a change in the cloud properties alone.

4.2 Climate feedbacks estimated via the radiative kernel method

Figure 3 shows the global mean feedbacks for lapse rate, water vapor, lapse rate plus water vapor, albedo, and cloud (SW, LW, and total) for each CMIP5 model. Both radiative kernels are used to test whether the results are sensitive to the choice of radiative kernel and whether the inter-model deviation is greater than the distribution of the radiatively active constituents (temperature, water vapor, albedo, and cloud) of the base model. It is worth clarifying that positive (negative) values of feedbacks contribute to the amplification (damping) of global warming. The strongest positive feedback (Fig. 3) is due to water vapor (mean value: 1.39 W m⁻² K⁻¹), followed by clouds (mean value: 0.96 W m⁻² K⁻¹) and then surface albedo (mean value: 0.32 W m⁻² K⁻¹). The global mean Planck feedback is negative, with an average value of −3.60 W m⁻² K⁻¹ (not shown in Fig. 3), followed by a lapse rate feedback of −0.77 W m⁻² K⁻¹. BESM yields values near the ensemble mean for the albedo and cloud feedbacks, i.e., 0.27 and 0.95 W m⁻² K⁻¹, respectively. For the lapse rate feedback BESM yields a value of −0.71 W m⁻² K⁻¹, slightly underestimating the ensemble mean value in magnitude. In turn, BESM is among the models with the lowest global water vapor feedback average, with a value of around 1.24 W m⁻² K⁻¹.

Figure 3 Global mean feedbacks for 25 CMIP5 models and BESM-OA2.5 (circle). Changes in abrupt4xCO2 relative to the piControl were averaged over years 120–150. The triangles represent the mean estimated feedback values calculated using the NCAR radiative kernel, whereas the upside-down triangles represent the estimated feedback values calculated using the GFDL radiative kernel.

Figure 4 shows the latitudinal profiles for the global mean feedback values in Fig. 3, allowing us to identify the regions that induce deviations of BESM results from the CMIP ensemble. In Fig. 4a–b, there is a nearly constant Planck feedback of approximately −3.4 W m⁻² K⁻¹ from 90 ^∘S to 60^∘N, with a notable increase in the ensemble standard deviation in the sub-Antarctic latitudes (around 60 ^∘S), which is in accordance with Rieger et al. (2017). The exception is in the Arctic region, where the mean value reaches −10 W m⁻² K⁻¹ with a similarly increased standard deviation value. In the sub-Antarctic and Arctic latitudes, BESM yields one of the most negative values for the Planck feedback, revealing that BESM has a stronger vertically homogeneous warming (corresponding to large surface warming) among the CMIP5 models. Furthermore, for those same regions, BESM showed more positive lapse rate feedback than the ensemble. Therefore, BESM yields a warming relatively larger at the surface and relatively weaker at the upper troposphere, resulting in a stronger vertical temperature gradient in comparison to the other models.

Figure 4 Feedbacks for the CMIP5 multi-model ensemble mean (solid line) and BESM-OA2.5 (solid line with dots) for the Planck feedback (a, b), lapse rate feedback (c, d), water vapor feedback (e, f), albedo feedback (g, h), and cloud feedback (i, j). Inter-model standard deviations for each latitude are in yellow. In blue are the feedback limits based on the maximum and minimum values for each latitude among the models, excluding BESM-OA2.5. All feedbacks are based on the averaged over years 120–150.

In the tropics, where there is an intense moist convection, atmospheric warming almost follows a moist adiabat (temperature increase is larger in the upper troposphere compared to that at the surface), implying a negative lapse rate feedback (Fig. 4c–d) (Manabe and Stouffer, 1980). In accordance with this upper tropospheric warming in the tropics, an increase in the specific humidity occurs (Manabe and Wetherald, 1975), which causes a reinforcement of the greenhouse effect, reflected by a positive water vapor feedback as shown in Fig. 4d–e. Because of this close link between the lapse rate and water vapor feedbacks, it is common to consider their effects as a sum, as displayed in Fig. 3. BESM shows a lapse rate feedback near the ensemble mean for the tropics. The greatest BESM deviations are observed near the Antarctic and over the Arctic (Fig. 4c–d), where this feedback became positive for all models. For the water vapor feedback, greater dispersion of the models was observed in the tropics, with BESM systematically yielding values below the ensemble mean for the same latitude band (Fig. 4e–f). These lower values extend throughout the Northern Hemisphere, consistent with the low global mean water vapor feedback value relative to the ensemble (shown in Fig. 3).

The albedo feedback profiles from BESM and the CMIP5 models are compared in Fig. 4g–h. Nonzero results mostly occur over the polar regions, where there is a reduction in sea ice and snow cover (Chung and Soden, 2015; Block et al., 2020). The positive albedo feedback signals yielded by all of the models in such regions imply that the reduction in the albedo corresponds to an increase in the radiation budget at the TOA due to reduced upward shortwave radiation. As expected, the polar regions present a large dispersion among the models, which is related to how fast the sea ice melts in the different climate models. The regions over the Arctic and the ocean near the Antarctic show the largest surface warming, and this positive albedo, together with the positive lapse rate feedbacks, is the main factor responsible for a phenomenon known as polar amplification (Pithan and Mauritsen, 2014; Block et al., 2020). BESM yielded an albedo feedback greater than the ensemble standard deviation over the Southern Ocean at around 60^∘ S. This same latitude is where BESM shows negative Planck and positive lapse rate feedbacks outside the models limits, as previously discussed.

Finally, regarding cloud feedback, most of the inter-model spread arises from the SW component (Figs. 3 and 4i–j). This dispersion is also reflected in the standard deviation and in the limit between the minimum and maximum of the zonally averaged cloud feedback shown in Fig. 4i–j. The SW cloud feedback ranges from −0.28 to 1.40 W m⁻² K⁻¹, while the LW cloud effect ranges from 0.10 to 0.96 W m⁻² K⁻¹. The combined SW and LW cloud effects result in positive cloud feedback ranging from 0.35 to 1.69 W m⁻² K⁻¹. This result is similar to that found by Soden et al. (2008) for CMIP3 (IPCC AR4, IPCC, 2007) models, who also reported a near-zero to positive cloud feedback. BESM presents positive values of around 0.5 W m⁻² K⁻¹ for both SW and LW cloud feedback, resulting in a total cloud feedback of about 1.0 W m⁻² K⁻¹ (as shown in Fig. 3). The highest positive values are in regions with strong albedo feedback (Fig. 4i–j).

Although BESM yielded global area-averaged feedbacks near the model ensemble mean values, differences are found mainly at high latitudes. In fact, for cloud feedback, BESM is an outlier due to a strong shortwave component response over both the Arctic and the Southern Ocean near the Antarctic. This effect is evident via the decomposition of the cloud feedback into the SW and LW components for the ensemble and BESM values, as shown in Fig. 5. To assess the analytical causes of this strong BESM shortwave cloud feedback departure from the ensemble, we separately computed the contributors to the SW cloud feedback, i.e., the SW CRE (as described by Cess et al., 1989) and the feedback cloud masks, as in Eq. (5). For the shortwave component, the feedback cloud masks are obtained for the albedo and SW water vapor feedbacks by calculating the all-sky minus clear-sky radiation flux for each feedback. We find that the higher BESM cloud feedback values (Fig. 6g–h) are mainly consequences of the sum of the SW CRE (Fig. 6a–b) and the effect of cloud masking for the albedo feedback (Fig. 6c–d) in the sub-Antarctic and Arctic regions. In turn, the cloud masking for the SW water vapor (Fig. 6e–f) does not contribute to the higher positive BESM values. As shown in Fig. 6, it is possible to attribute BESM's status as an outlier over the Arctic region to the SW CRE, while in the Southern Ocean (around 60^∘ S), the major contribution comes from the albedo feedback cloud mask.

Figure 5 Cloud feedbacks calculated using the NCAR radiative kernel for the CMIP5 ensemble (a, c, e) and BESM-OA2.5 (b, d, f). These results are based on the average from years 120–150.

Figure 6 Shortwave cloud radiative effect (a, b), the albedo (c, d) and shortwave water vapor (e, f) feedback cloud maskings, and shortwave cloud feedback (g, h) for the CMIP5 multi-model ensemble mean (solid line) and BESM-OA2.5 (solid line with dots). Inter-model standard deviations for each latitude are in yellow. In blue are the feedback limits based on the maximum and minimum values for each latitude among the models, excluding BESM-OA2.5 (physical units of these feedbacks: W m⁻² K⁻¹).

A deepened physical analysis to understand the BESM cloud feedback behavior in high latitudes is obtained by examining the zonal mean of the change in the cloud vertical profile for BESM, as shown in Fig. 7. Over the Arctic and near the Antarctic, BESM showed an increase in the cloud fraction above 850 hPa and a decrease below that level, indicating an upward shift of low-level clouds (Fig. 7a). However, the increase in cloud cover above 850 hPa is stronger than the reduction below. Because of this increase in the total cloud fraction, a negative SW CRE (not adjusted) appears in those regions (Fig. 6a–b), consistent with an increase in sun shading (Fig. 7b). Moreover, the SW cooling is smaller than the heating provided by LW radiation due to the upward shift of the low-level clouds, as evident in the net effect (Fig. 7d). Furthermore, the positive albedo and lapse rate feedbacks (Fig. 4c–f) are consistent with this vertical cloud shifting. In this manner, a loss of SW energy at the surface associated with an increase in the total cloud fraction explains the negative SW CRE, and the gain of LW energy is responsible for the sea ice melting. Consequently, this gain of LW energy is indirectly linked to the albedo feedback cloud mask for BESM, since the mask ($\mathrm{\Delta }a/\mathrm{\Delta }{\stackrel{\mathrm{‾}}{T}}_{\mathrm{as}}(K_a-K_a^{\mathrm{cs}})$) is proportional to the albedo change (Δa). As discussed before, both the SW CRE and albedo feedback cloud mask contribute to the large positive cloud feedback over the Arctic and sub-Antarctic areas observed in BESM.

Figure 7 Vertical profiles of the zonal mean of the abrupt4xCO2–piControl mean difference for the following variables: (a) cloud fraction; radiative heating or cooling rate (dT∕dt; K d⁻¹) of (b) shortwave, (c) longwave, and the (d) sum of the shortwave and longwave for BESM-OA2.5.

4.3 Changes in temperature, atmospheric circulation, and precipitation

Figure 8 shows the annual mean surface temperature differences between the abrupt4xCO2 and piControl scenarios for the ensemble of 25 CMIP5 models and BESM. As clearly shown in Fig. 8, despite the generalized increase in the surface temperature over most of the globe, BESM shows a generally lower temperature increase, principally over the continental areas. The CMIP5 ensemble yields a mean continental temperature increase of 6.78 K, while BESM yields a value of 5.57 K. Nevertheless, the spatial patterns of the temperature increases are similar, as measured by the spatial inner product (as described in the previous section) between Fig. 8a and b, which results in a value of 0.96 (values near 1 indicate that both variables have a similar spatial pattern, whereas values near 0 mean that there is hardly any pattern correlation between variables). One point of interest within the scientific community is the relatively low temperature increase over the subpolar North Atlantic, also referred to as the warming hole (Drijfhout et al., 2012). In the CMIP5 ensemble mean, the North Atlantic does not show a temperature decrease, although it is the region with the smallest temperature increase globally, while BESM shows an area of temperature decrease in this region. Such a decrease is also present in six other analyzed models (CSIRO-Mk3-6-0, FGOALS-s2, GFDL-ESM2G, GFDL-ESM2M, GISS-E2-R, and inmcm4). These results are consistent with those of Drijfhout et al. (2012), who showed that both observations and CMIP5 models present maximum cooling in the center of the subpolar gyre. The authors argue that there is evidence that both the subpolar gyre and the AMOC adjust in concert with different time lags.

Figure 8 Differences (averaged over years 120–150) in surface temperature between the abrupt4xCO2 and piControl simulations in (a) the CMIP5 ensemble and (b) BESM-OA2.5. (c) A scatter plot of the global average surface temperature for the CMIP5 models used in the ensemble and BESM-OA2.5 (black dot). The shaded areas in (a) and (b) have confidence levels greater than 90 %; the black line represents the contour line of zero temperature difference.

The regions with the largest temperature increase in the abrupt4xCO2 scenario are the polar regions, particularly over the Arctic. The equatorial Pacific shows a relative maximum in warming when the abrupt4xCO2 scenario is compared with the piControl both in the CMIP5 ensemble and BESM. Such changes in the Pacific mean state are consistent with the IPCC AR5, which reports that the Pacific Ocean becomes warmer near the Equator compared to the subtropics in the CMIP5 projections (Liu et al., 2005; Gastineau and Soden, 2009; Cai et al., 2015). The scatter plot of the global average under abrupt4xCO2 conditions versus the piControl conditions presented in Fig. 8 provides additional information that helps to understand the dispersion around the mean value among the different models. Even though the outputs of most of the models lie in either quadrant 1 or 3 (top right and bottom left, respectively), it is not possible to claim any robust linear relationship. This result indicates that models with warmer or cooler mean climates in the piControl runs apparently do not show a corresponding warmer or cooler climate in the abrupt4xCO2 experiments. BESM yields a temperature near that of the ensemble in both the piControl and abrupt4xCO2 runs; consequently, it also showed a temperature increase near the ensemble mean, consistent with its Planck feedback (Figs. 3 and 4a–b).

Figure 9 shows the precipitation changes between the abrupt4xCO2 and piControl scenarios for the multi-model ensemble and the BESM. The results are, in general, similar to those of Held and Soden (2006), with wet regions becoming wetter (near-equatorial and subpolar regions) and dry regions becoming drier (centered around 30^∘ in both hemispheres). The precipitation pattern in the CMIP5 ensemble shows increased precipitation over the equatorial Pacific, which could be related to the equatorial Pacific warming pattern shown in the temperature change (Fig. 8). Furthermore, the CMIP5 ensemble shows a decrease in precipitation in northern South America. The BESM precipitation pattern is similar to the spatial patterns in the CMIP5 ensemble but with some notable exceptions. For example, the decreased precipitation over the South Pacific shown in the CMIP5 ensemble plot is extended into the Indonesian region in BESM. It is also worth noting that in the BESM simulation, the South Pacific convergence zone (SPCZ) shifts southward in the abrupt4xCO2 scenario compared with its position in the piControl scenario. In both the multi-model ensemble and BESM, the precipitation change pattern over South America is similar to that which occurs during El Niño years (Kayano et al., 1988; Marengo and Hastenrath, 1993; Grimm and Tedeschi, 2009), with increased precipitation over southeastern South America and decreased precipitation over northern–northeastern South America. The scatter plot in Fig. 9 emphasizes a linear relationship between the experiments, indicating that models with higher (lower) global average precipitation in the piControl scenario show higher (lower) precipitation in the abrupt4xCO2 scenario. As is obvious from Fig. 9, the BESM performance perfectly matches the ensemble mean behavior in global mean.

Figure 9 Differences (averaged over years 120–150) in precipitation (mm month⁻¹) between the abrupt4xCO2 and piControl simulations in (a) the CMIP5 ensemble and (b) BESM-OA2.5. (c) A scatter plot of the precipitation global averages for the CMIP5 models used in the ensemble and BESM-OA2.5 (black dot). The shaded areas in (a) and (b) have confidence levels greater than 90 %; the black line represents the contour line of zero precipitation difference.

Figure 10 Scatter plot of the ECS and ΔPr (%) values for all of the ensemble models. The solid black line shows the linear fit between the ECS and the perceptual precipitation change. As in Fig. 2, the models are sorted according their ECS value. The dashed lines represent the error limits considering the residual standard error.

Figure 10 shows a scatter plot of the ECS versus the precipitation difference between the abrupt4xCO2 and piControl scenarios (ΔPr) for all of the considered models. It is worth noting that all the models show increased global mean precipitation upon a quadrupling of atmospheric CO₂ with piControl preindustrial CO₂ concentrations (positive values on the y axis in Fig. 10). An apparent linear relationship between these differences (abrupt4xCO2 minus piControl) in the global mean precipitation and ECS is also evident in Fig. 10, in which the warmest models tend to have the largest precipitation changes. The slope of the linear regression reflects a 2.5 % precipitation change per Kelvin, which is close to that found by Held and Soden (2006). This slope is much lower than that predicted by the Clausius–Clapeyron relation, i.e., an approximately 6.5 % change in precipitation per Kelvin. In fact, such precipitation increases are not governed by moisture availability but rather by the surface and tropospheric energy balance, which incorporates the surface radiative heating, surface latent heat flux, and radiative cooling of the troposphere (Allen and Ingram, 2002).

MRI-CGCM3, ACCESS1-0, and HadGEM2-ES show greater deviations from the linear fit shown in Fig. 10. Furthermore, BESM is marginally out of the residual standard error interval, with a 9.5 % increase in precipitation (the error limit is 9.2 %). ACCESS1-0 and HadGEM2-ES use the same atmospheric model (Bi et al., 2013; Dix et al., 2013), which could explain the lower increase in precipitation in both coupled models.

In addition to temperature and precipitation changes, we are also interested in understanding the changes in the BESM atmospheric circulation (compared to other models) following a quadrupling of the CO₂ concentration. The sea level pressure (SLP) response patterns shown in Fig. 11 depict a poleward shift in the subtropical high-pressure cells in both the CMIP5 ensemble and BESM. Furthermore, when the models are subjected to the increased atmospheric CO₂ concentration, a decrease in the SLP over the polar regions is evident. This, connected with the increase in the midlatitudes, indicates a positive trend in Arctic oscillation (AO) and Antarctic oscillation (AAO) episodes, which have already been reported by Fyfe et al. (1999), Cai et al. (2003), and Miller et al. (2006). It is also interesting to note the statistically significant SLP decrease (increase) over the eastern (western) Pacific, a pattern that might indicate an ENSO-like pattern in scenarios with an increased CO₂ concentration. This pattern is consistent with those depicted in Fig. 8 for the SST changes in a 4xCO₂ scenario.

Figure 11 Difference (averaged over years 120–150) in sea level pressure (SLP; hPa) between the two scenarios (abrupt4xCO2 minus piControl; shaded) and SLP under piControl conditions (contours) in the CMIP5 model ensemble (a) and BESM-OA2.5 (b). The white areas have confidence levels less than 90 %.

The results from the piControl scenario (the contours in Fig. 12) show that the Southern Hemisphere subtropical jet, reflected by the core of the maximum eastward zonal wind, is localized around 35^∘ S at 200–150 hPa in both the CMIP5 ensemble and BESM. In Fig. 12a and b (BESM and the CMIP5 ensemble), we note that the regions with the strongest increases in westerly winds at all levels show a southward jet displacement. This observation is consistent with the poleward displacement of the high SLP center shown in Fig. 11. Furthermore, as the high-pressure centers experienced a poleward shift, the pressure gradients intensified in the subpolar areas; consequently, the near-surface wind velocity increased, following the geostrophic approximation ($u\approx -(\mathrm{1}/f\mathit{\rho })(\partial p/\partial y)$), where f is the Coriolis parameter and ρ is the air density.

Figure 12 Vertical profile of the difference (averaged over years 120–150) in the zonal wind (m s⁻¹) between the two scenarios (abrupt4xCO2 minus piControl; shaded) and the piControl conditions (contours) for (a) the ensemble of CMIP5 models and (b) BESM-OA2.5. The white areas have confidence levels less than 90 %.

Figure 13 shows the average differences from 5^∘ N to 5^∘ S (Walker circulation) between the abrupt4xCO2 and piControl scenarios for vertical motion (omega; shaded) and zonal wind and vertical motion (vectors). According to the vertical motion pattern in the piControl (contours), the multi-model ensemble and BESM show subsidence over an extensive area in the Pacific (150^∘ E–90^∘ W) whose intensity is lower in the abrupt4xCO2 simulation, as shown in Fig. 13 (blue). This finding is coherent with near-surface temperature patterns (Fig. 8), which show an equatorial warming pattern in the mean state (e.g., during El Niño years, a weakening of the Walker circulation occurs). Furthermore, there is enhanced subsidence for the difference between the two scenarios over South America (around 75^∘ W), consistent with the decrease in precipitation in tropical South America, in both BESM and the CMIP5 ensemble (Fig. 9).

Figure 13 Difference (averaged over years 120–150) between the abrupt4xCO2 and piControl conditions for generalized vertical motion (omega (Pa s⁻¹); shaded) and the mean zonal and vertical motion (vectors) averaged between 5^∘ S and 5^∘ N for (a) the CMIP5 ensemble and (b) BESM-OA2.5. The contours represents the averaged piControl vertical motion (omega) in the same region. The white regions have confidence levels less than 90 %.

5 Conclusions

The piControl and abrupt4xCO2 scenarios for 25 CMIP5 models have been compared with those generated by BESM based on their key sensitivity parameters, such as the equilibrium climate sensitivity (ECS) and climate feedbacks. Furthermore, changes in the temperature, atmospheric circulation, and precipitation patterns were investigated.

Applying the linear regression method (Gregory et al., 2004), we obtained ECS values for the 25 CMIP5 models analyzed that ranged from 2.07 to 4.74 K, with BESM showing 2.96 K, close to the ensemble mean value (3.30±0.76). BESM has one of the biggest radiative forcing (G) values, i.e., 8.62 W m⁻² K⁻¹, which is related to the radiative transfer model and the rapid adjustment process (Collins et al., 2006; Gregory and Webb, 2008; Andrews and Forster, 2008). Both G and the climate sensitivity (λ) define the ECS values calculated with this method; however, only λ shows a significant correlation with the ECS, corroborating the results of Andrews et al. (2012).

To go further in the analysis, the radiative kernel method was used to separate the climate feedback into Planck, lapse rate, water vapor, albedo, and cloud feedbacks. Two regions presented considerable inter-model variability for the Planck, lapse rate, and albedo values, i.e., the Arctic region and over the ocean near the Antarctic. Over these regions, the BESM zonal mean cloud feedback ranges outside the standard deviation for the analyzed models, reaching approximately 3 W m⁻² K⁻¹, while the zonal mean was close to zero. BESM showed an upward shift of low-cloud cover and an increase in cloud cover between 850 and 700 hPa, and these features were responsible for sun shading at the surface, which increased the reflected solar radiation at the TOA. Moreover, BESM presented a greater albedo change compared with those of the other models, especially in the sub-Antarctic area. Despite the loss of SW energy at the surface, which results in a negative SW cloud radiative effect, this effect was overcome by the albedo feedback cloud mask, which contributes to positive cloud feedback over those regions.

The atmospheric circulation patterns in BESM were similar to the patterns of the multi-model ensemble and those of other studies regarding the near-surface temperature (IPCC, 2007, 2013). For precipitation, the thermodynamic component reflects the well-known “wet-gets-wetter” and “dry-gets-drier” patterns of precipitation change (Held and Soden, 2006). BESM and the CMIP5 ensemble show consistent weakening of the Walker circulation, principally in the Pacific and over northern South America, which has been reported in previous studies (Collins et al., 2010; DiNezio et al., 2012; Huang and Xie, 2015; Cai et al., 2015). Regarding SLP, both BESM and the CMIP5 ensemble indicate a poleward displacement of the subtropical high-pressure systems, as shown in other studies (Fyfe et al., 1999; Cai et al., 2003; Miller et al., 2006). In line with such displacement, the subtropical jet also shifted polewards, and this effect was more distinct in the Southern Hemisphere.

Summarizing, we conclude that BESM-OA2.5 is a climate model that can reproduce proven physical processes that determine and modify changes in the global climate system. In this sense, the analysis methods used here have the potential to explain remaining process uncertainties causing inter-model spread in the cloud feedback in future work. This notwithstanding, the BESM team continues is effort to improve the cloud parameterization of the model as well as its land surface model in subsequent versions. Furthermore, it is important to mention that the radiative energy imbalance of −4 W m⁻² at the TOA, arising from our ocean–atmosphere coupling, is seen as an issue that is to be tackled in ongoing model development work. We hope that the next version will include improved energy flow diagnostics and that it will include physical parameterizations of atmosphere–ocean interactions that lead to better agreement with other models and with observations.