The response of the global climate system to atmospheric CO2 concentration increase in time is scrutinized employing the Brazilian Earth System Model Ocean–Atmosphere version 2.3 (BESM-OA2.3). Through the achievement of over 2000 yr of coupled model integrations in ensemble mode, it is shown that the model simulates the signal of recent changes of global climate trends, depicting a steady atmospheric and oceanic temperature increase and corresponding marine ice retreat. The model simulations encompass the time period from 1960 to 2105, following the phase 5 of the Coupled Model Intercomparison Project (CMIP5) protocol. Notwithstanding the accurate reproduction of large-scale ocean–atmosphere coupled phenomena, like the ENSO phenomena over the equatorial Pacific and the interhemispheric gradient mode over the tropical Atlantic, the BESM-OA2.3 coupled model shows systematic errors on sea surface temperature and precipitation that resemble those of other global coupled climate models. Yet, the simulations demonstrate the model’s potential to contribute to the international efforts on global climate change research, sparking interest in global climate change research within the Brazilian climate modeling community, constituting a building block of the Brazilian Framework for Global Climate Change Research.
The conclusions of the Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report (Solomon et al. 2007) are unequivocal regarding the human interference on the global climate via the buildup of greenhouse gases (GHG) into the atmosphere. The steady increase of GHG concentration in the atmosphere, with CO2 being the most notable of them, has led the scientific community and governments alike to debate the issue to exhaustion. Therefore, a detailed scrutiny of the influence of atmospheric GHG concentration on the climate system is a must and has been delved in by several research groups worldwide. Such scrutiny is presently applied to matters ranging from the development of state-of-the-art coupled climate models up to complex earth system models (ESM) that incorporate the complexity of the many components of the earth system. Examples of such complex system models are those developed by the largest climate research centers in the world, such as the National Center for Atmospheric Research (NCAR) Community Climate System Model, version 4 (CCSM4) (Bates et al. 2012; Grodsky et al. 2012); Geophysical Fluid Dynamics Laboratory Climate Model version 2.0 (GFDL CM2.3) (Delworth et al. 2012); and the Hadley Centre Global Environmental Model, version 2 (Earth System) (HadGEM-ES). Yet, other climate centers around the world are also developing their own climate models, either by adapting a particular state-of-the-art global ESM or by combining ocean, atmosphere, and surface component models from various sources, creating their “own” global climate model. This is the path chosen by Brazil, which has decided to develop its ESM, whose initial steps are underway and are detailed in this article.
The effort to develop Brazil’s ESM is part of a national coordinated effort to build a multidisciplinary research framework geared at understanding the causes of global climate change, its effects, and its impacts on society. In this coordinated effort, the Brazilian National Institute for Space Research (INPE) is leading the construction of a global climate model capable of predicting the consequences of the buildup of GHG in the atmosphere. This article describes the first version of this model, here represented by the global ocean–atmosphere coupled model, named the Brazilian Earth System Model Ocean–Atmosphere version 2.3 (BESM-OA2.3). For this article, we have followed the criteria for participation in phase 5 of the Coupled Model Intercomparison Project (CMIP5) protocol (Taylor et al. 2009), aiming at simulating the behavior of the coupled ocean–atmosphere system on decadal time scales under varying GHG concentrations in the atmosphere.
In so doing, we have developed a model that is capable of performing seasonal to decadal predictions, moving in the direction of a seamless prediction system (Palmer et al. 2008; Shukla 2009), as other centers have already embraced worldwide. With this goal in mind, this article presents the ability of BESM-OA2.3 to simulate global trends and global climate variability related to historical variations of GHG concentrations in the atmosphere. Such exercise constitutes a benchmark of the model’s capability to be used for studies involving more complex phenomena in the future.
BESM-OA2.3 is an evolution of previous versions of the Brazilian Center for Weather Forecasting and Climate Studies (CPTEC) coupled ocean–atmosphere model, which have been used to study the impacts of Amazon deforestation on climate (Nobre et al. 2009), the dynamics of the South Atlantic convergence zone (SACZ) over cold waters (Chaves and Nobre 2004; Nobre et al. 2012), the importance of ocean–atmosphere coupling to represent the Atlantic thermocline depth (Siqueira and Nobre 2006), and regional seasonal predictions with the Eta regional model over South America (Pilotto et al. 2012). All these studies highlighted the importance of ocean–atmosphere coupling to better predict the phenomena of study.
The challenge in this article is to demonstrate that the coupled ocean–atmosphere model developed at INPE is suitable to simulate climate phenomena along the last few decades and into the near future, for the period from 1960 to 2035. This article shows that the model is capable of reproducing global trends of ocean heat content and top-of-atmosphere radiation balance, as well as important modes of climate variability. The effects of observed variations of atmospheric CO2 during the period of study and its projection into the near future are also investigated. Regardless of the simplicity of BESM-OA2.3 compared to full ESMs (e.g., Muñoz et al. 2012; Delworth et al. 2012), its results are comparable to other global climate models.
2. Model description
The BESM-OA2.3 coupled climate model used in this research is constituted by the CPTEC/INPE atmospheric general circulation model (AGCM) coupled to GFDL’s Modular Ocean Model version 4p1 (MOM4p1) oceanic general circulation model (OGCM) via GFDL’s Flexible Modular System. The component models and coupling are described as follows.
a. Atmosphere model
The CPTEC/INPE AGCM has been constantly modified during the last few years. The implementation of two new dynamic cores (Eulerian and semi-Lagrangian) and the algorithm for the use of a reduced grid are now employed operationally at CPTEC. For the present work, the Eulerian core and reduced grid have been used. The options for physical parameterizations include the Simplified Simple Biosphere Model (SSiB) module surface (Xue et al. 1991), turbulence level 2 module (Mellor and Yamada 1982), gravity wave module (Anthes 1977), deep convection module (Arakawa and Shubert 1974; Grell and Devenyi 2002), shallow convection module (Tiedtke 1984), shortwave radiation codes (Tarasova et al. 2006; Chou and Suarez 1999), longwave module (Harshvardhan and Corsetti 1984; Harshvardhan et al. 1987), and the model for interaction of radiation with clouds (Slingo 1987; Hou 1990; Kinter et al. 1997). In the shortwave and longwave radiation schemes the carbon dioxide is globally specified, assuming either a 370-ppm fixed value (for the control integrations, labeled “physics2” for the CMIP5 protocol) or a time-dependent function (for transient integrations, labeled “physics1”).
The new deep convection scheme implemented in the CPTEC AGCM is based on the convective parameterization developed by Grell and Devenyi (2002) and referred to here as the Grell ensemble (GRE) scheme. This scheme includes several alternative closure assumptions to estimate the cloud base mass flux. Both the first type of closure (Grell closure) and the second type (Arakawa–Schubert closure) assume a quasi-equilibrium state between large-scale forcing and convection. Other closures used are the Kain–Fritsch closure, in which the atmospheric stability is removed by the convection within the specified convective time, moisture convergence closure, and low-level vertical velocity mass flux closure. The set of convection scheme ensembles used in the CPTEC AGCM includes 144 different versions of the cumulus parameterization, 16 from dynamic control, and 9 from static control.
The exchanges of heat, moisture, and momentum between the surface and atmosphere in the CPTEC AGCM over the ocean and continents are computed differently by various physical processes that define the surface fluxes. The model uses an empirical relationship that approximates the theory of Monin–Obukhov similarity to determine the transfer and the drag coefficient. The relationships over the ocean differ from those over land. Over the oceans, the momentum, heat, and moisture fluxes are parameterized considering the global aerodynamic scheme (bulk aerodynamic schemes), where the fluxes are assumed to be proportional to wind speed and the surface potential (difference in temperature or moisture between the ocean surface and the air near the surface). The proportionality or drag coefficient is calculated by an empirically determined functional fit. On the surface of the ocean the functional fit is similar to the one used by Sato et al. (1989). With regard to momentum, heat, and moisture fluxes over the continent, the parameterization used is more sophisticated. The surface scheme on the continent considers a variety of different types of vegetation, soil, and geographical formation. The SSiB (Xue et al. 1991) scheme used in the CPTEC AGCM, despite being a simplified model, has many physical processes that represent the interactions between land surface and atmosphere.
The AGCM has spectral horizontal resolution truncated at triangular wavenumber 62. This gives an equivalent grid size of 1.875° and 28 sigma levels unevenly spaced in the vertical (i.e., T062L28).
b. Ocean model
The ocean model used in this work is the MOM4p1 (Griffies 2009) from GFDL, which includes the Sea Ice Simulator (SIS) built-in ice model (Winton 2000). The horizontal grid resolution is set to 1° in the longitudinal direction, and in the latitudinal direction the grid spacing is ¼° in the tropical region (10°S–10°N), decreasing uniformly to 1° at 45° and to 2° at 90° in both hemispheres. For the vertical axis, 50 levels are adopted with a 10-m resolution in the upper 220 m, increasing gradually to about 370 m of grid spacing at deeper regions. The variant of the ocean model not coupled to an atmospheric model, which is referred to here as an OGCM or forced model, is forced by winds, temperature, specific humidity at 10 m, sea level pressure, radiative fluxes, sea ice cover and thickness, and river discharges. The OGCM spinup was done in a manner similar to that used by Nobre et al. (2012)—starting from zero currents velocities and Levitus (1982) ,T–S structure of the World Ocean, running the forced model for 13 yr with climatological atmospheric forcing fields, and followed by an additional 58 yr (1950–2007) forced by interannual fields from Large and Yeager (2009)—except for river discharges and the sea ice variables, which were kept climatological in both cases. During the forced model run, the dynamic and thermodynamic structure of the ocean was saved at regular periods of time in order to be used as initial states of the ocean for the coupled model experiments.
The model settings are similar to the Coordinated Ocean-Ice Reference Experiments (COREs) experiment settings (Griffies et al. 2009) apart from modifications made to obtain better results for our forced and coupled experiments at the chosen ocean grid, such as changing the smoothing of the sea surface height (“Laplacian” instead of “biharmonic”), disabling the “mixing downslope” module (which treats the mixing of tracer between dense shallow parcels and deeper parcels downslope), disabling the tidal forcing module, enabling the “convect” module (which vertically adjusts gravitationally unstable columns of ocean fluid), changing parameters of the frazil module (which computes the heating of seawater due to frazil ice formation), disabling the restoring masks, and other restoring options, with the exception of the salinity restoring, which is kept on with a time scale of 60 days. The turbulent flux calculation described by Large and Yeager (2009) is used in the forced model in order to compute the surface drag coefficients.
The coupled version of the ocean model has additional changes in relation to the forced model. The background vertical diffusivity was defined by the Bryan–Lewis scheme (Bryan and Lewis 1979), zonally constant but varying both meridionally and vertically. Between 10°N and 10°S, the Bryan–Lewis diffusion coefficient span from 2 × 10−5 m2 s−1 near the surface to 7.8 × 10−5 m2 s−1 at the bottom. Such values increase linearly poleward to compensate for the coarser grid, measuring 2.3 × 10−5 m2 s−1 at the surface and 8.7 × 10−5 m2 s−1 at the bottom near the poles. The consequent near surface eddy kinetic energy is coherent with the observations (not shown), considering the limitations of the simulations setup. The BESM-OA2.3 simulations develop up to 1.2 cm2 s−2 average eddy kinetic energy over a few degrees around the equator. This value rapidly decays poleward because of the gradually coarser grid resolution near the poles. Moreover, in the coupled model all restoring options including the salinity restoring are disabled.
c. Coupling strategy
The coupled model uses GFDL’s Flexible Modeling System (FMS) coupler to couple the CPTEC AGCM to GFDL MOM4p1 (Griffies 2009) and SIS (Winton 2000). The main characteristics of the CPTEC AGCM are described in the previous subsection. The original configuration of MOM4p1 is similar to that used in the COREs release of the model for the forced experiments using Large and Yeager (2009) fields. For the coupled model, all restore flags were turned off. Wind stress fields are computed, using Monin–Obukhov scheme within MOM4p1, from the wind field 10 m above the ocean surface. Adjustments were done to the Monin–Obukhov boundary layer scheme, whose parameters were tuned according to the wind fields output by the CPTEC AGCM. The hydrological cycle is balanced by an indirect method, as suggested in Griffies (2009). The AGCM receives the following two fields from the coupler: sea surface temperature (SST) and ocean albedo from ocean and sea ice models at an hourly rate (coupling time step). Adjustments were also made to ocean shortwave penetration parameters due to the CTPEC AGCM supply of visible and infrared short wave radiation. The coupling variables supplied by the AGCM are as follows: freshwater (liquid and solid precipitation), specific humidity, heat, vertical diffusion of velocity components, momentum fluxes, and surface pressure.
d. Experiment design
The numerical experiments design followed the CMIP5 protocols (Taylor et al. 2009) for decadal and short-term simulations and near-future predictions using two types of physics for the coupled model differing on the concentration of atmospheric CO2. Integrations labeled “physics1” used atmospheric CO2 concentrations derived from in situ air samples collected at Mauna Loa Observatory, Hawaii. Historical data from Bacastow et al. (1985) were regressed onto a parabolic equation defined by
where the CO2 concentration is in ppm, t is the number of fractional years since 2000, a = 0.011 669 6, b = 1.799 84, and c = 369. Integrations labeled “physics2” used atmospheric CO2 concentrations fixed at 370 ppm and are referred to as control run throughout the text.
Ensemble members were integrated for 10 yr each with initial conditions (IC) on 1–10 December of years 1960, 1965, 1970, 1975, 1980, 1985, 1990, 1995, 2000, and 2005, for each of physics1 and physics2. Three of these ensembles (viz., 1960, 1980, and 2005) were extended for an extra 20 yr for each of the 10 members, completing 30-yr-long integrations each, referred to in the text as decadal runs, and an extra 70 yr for only the IC conditions of 1 December, completing 100-yr-long single simulations, referred to in the text as centennial runs. In addition, one member CMIP5 “1pct-co2” run was computed, assuming a 1% yearly increase in the atmospheric CO2 concentration, starting at the 20th year of the 1960 physics2 run and completing 80 yr.
The CPTEC AGCM initial conditions for each ensemble member used the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) reanalysis fields (Kalnay et al. 1996) for the 0000 UTC of each day from 1 to 10 December of the chosen years. The ocean initial states were chosen for the same dates from a spinup run of MOM4p1 that used prescribed atmospheric fields of momentum, solar radiation, air temperature, and freshwater described earlier in this section.
The outputs from the ocean model were generated directly at the proper frequencies requested according to the CMIP5 experiment specification. The atmospheric data were output at a 3-hourly frequency and later processed using the Climate Model Output Rewriter version 2 (CMOR2) software (Doutriaux and Tayler 2011) to satisfy all CMIP5 output requirements.
a. Atmospheric CO2 concentration, ocean heat content, and TOA fluxes
The residual SST between the physics1 experiment and the physics2 control run show a uniform global SST warming for the 100-yr-long experiments (Fig. 1a) showing a trend of 0.054 K decade−1 (p < 0.01). This corresponds to approximately 50% of the observed global SST warming in the 1980–2011 period, as estimated from averaging the second Met Office Hadley Centre Sea Surface Temperature Dataset (HadSST2) anomalies globally (figure not shown).
The propagation of the CO2-induced extra heat into the surface layers of the ocean can be illustrated in Fig. 1b, which shows the global average physics1 minus physics2 mean temperature as a function of depth and time for the centennial2005 run. In Fig. 1b the surface warming penetrating characteristic with warmer physics1 ocean during the 2005–2105 period is noteworthy. The residual temperature propagation due to atmospheric CO2 concentration changes reaches the upper 400 m of the oceans in approximately 50 yr. Such ocean warming is consistent with the global ice reduction shown in Fig. 1c in terms of physics1 − physics2, acting to reduce the surface albedo.
Yet, the analyses of decadal trends in SST are only weakly indicative of changes in the top-of-atmosphere radiation balance (TOA). The importance of ocean heat content (OHC) as a diagnostic for changes in the climate system has been emphasized in many recent papers (Domingues et al. 2008; Levitus et al. 2012; Palmer et al. 2011). Trends in total ocean heat content strongly constrain TOA since the ocean is the primary heat store in the earth system. The essential balance in the earth system is between total OHC (the primary heat storage term) and (the total radiation flux entering or leaving the system). Therefore, we use conservation of energy to write the following equation:
where is the net inward radiative flux through the top of the atmosphere, is the area of the earth, and is the ocean heat content.
Numerous assumptions underlie this equation, particularly the neglect of the thermal energy of the other components of the climate system (atmosphere, land, cryosphere, and geothermal fluxes). Therefore, changes in OHC should approximately balance the time-integrated TOA. Climate modeling and observations indicate that, to understand the global energy budget fully, one needs to include measurements of the deep oceans in the OHC analysis. The surface layers, even down to 700 m, are not robust indicators of total OHC, although the OHC at 700 m is strongly correlated (though not perfectly) to the TOA imbalance.
The key questions here is whether this CO2-induced extra heat into the oceans is indicative of changes in the earth’s top-of-atmosphere radiation balance and thus an increase of the accumulation of energy in the earth system. Figure 2a shows the full depth ocean heat content anomaly (OHCA) of the centennial2005 (physics1 minus physics2) runs. This global average show a steady rise, to the end of the 100 yr, of the full depth ocean heat content for the centennial2005 (physics1) relative to the control run with fixed CO2 (physics2).
To assess the energy flux change on the top of the atmosphere due to different CO2 forcing conditions, we compute the top-of-atmosphere radiation flux imbalance implied by the ocean heat content difference in terms of its trend over the 100-yr period shown in Fig. 2b. The estimated trend in TOA for the 2005–2105 period is +0.3 W m−2 yr−1, represented by the red dashed line in Fig. 2b. In a recent determination of OHC for the 0–2000-m layer, Levitus et al. (2012) report an implied rate in TOA of +0.27 W m−2 per unit area of earth’s surface for the 1955–2010 period.
The radiative flux change on the top of the atmosphere due to different CO2 forcing conditions can be affected by changes in clouds radiative forcing (e.g., Senior 1999; Stephens 2005). Radiative feedbacks involving low-level clouds are a primary cause of uncertainty in global climate model projections (Solomon et al. 2007). To investigate the low cloud responses to climate perturbations, we computed the time series of global low cloud cover physics1 − physics2 runs during the 2005–2105 period. The increase in atmospheric CO2 concentration in BESM-OA2.3 produces a small change in low cloud cover, around 0.1% increase over 100 yr (figure not shown). The result is in accordance with other cloud responses to climate perturbations studies (Wyant et al. 2006; Williams et al. 2006).
The results in this subsection indicate that there are mechanisms operating in the BESM-OA2.3 climate model that are capable of redistributing substantial quantities of heat at all depths of the ocean on decadal time scale because of changes in atmospheric CO2 concentration.
b. The Atlantic meridional overturning circulation
The BESM-OA2.3 is able to reproduce the Atlantic meridional overturning circulation (AMOC), a fundamental process in the Atlantic Ocean. North Atlantic Deep Water (NADW) is formed on the north Atlantic and exported southward below 1000 m, while a warmer near surface flow closes the system determining a northward net heat transport (Talley 2003). Figure 3 shows the meridional spatial pattern of the AMOC as represented by the BESM-OA2.3 and obtained from the average along 25 yr (from years 95 to 120 of physics2 centennial run) of the vertical cumulative sum of the zonal total transport for each depth and longitude. In other words, each point shows the average meridional total transport through a section along that longitude, from that depth up to the surface. The upper 1000 m in the Atlantic basin is dominated by a northward transport up to 35°N. Between 1000 and 3500 m the cumulative transport reduces as a consequence of the southward NADW, while below 3500 m the Antarctic Bottom Water flows northward.
The BESM-OA2.3 not only reproduces the three-layer meridional overturning circulation system but also well determines its intensity. Probably the best estimate of the AMOC has been done by the Rapid Climate Change–Meridional Overturning Circulation (RAPID–MOC), an observing array system along the 26.5°N cross section operating since 2004, suggesting a transport of 18.5 ± 4.9 Sv (1 Sv ≡ 106 m3 s−1) (Johns et al. 2011). At the same section, the BESM-OA2.3 reproduces a maximum transport of approximately 20 Sv, slightly higher than but still within the uncertainty of the observations.
The most impressive point is that this is the preferred stable state of the model. The small frame in Fig. 4 shows the maximum AMOC transport throughout the 122 yr elapsed after the spinup: that is, once the coupled model runs free with physics2 configuration (constant atmospheric CO2 concentration). The transport is initially just below 10 Sv, nearly half of the observed value, but it continually increases until stabilizing at approximately 20 Sv, after 80 yr. There was no prescribed forcing driving to this result, except for the model internal dynamics. The adequate adjustment of the AMOC reproduced on the BESM-OA2.3, despite the bias from the spinup, is a strong argument that the model is capturing the main earth system processes and hence should be a good approximation of reality.
c. Model climatology and systematic errors
SST and precipitation annual biases
To evaluate the ability of BESM-OA2.3 to reproduce the observed climate, the annual biases of SST and precipitation were computed and contrasted with the biases for other three CMIP5 global models: namely, Hadley Centre Coupled Model version 3 (HadCM3), NCAR’s CCSM4 coupled model, and the Fourth Generation Canadian Coupled Global Climate Model (CanCM4). The optimal interpolation SST (OISST) dataset, combining in situ and satellite SST interpolated onto a 1° grid with weekly resolution (Reynolds et al. 2002), was used to assess the SST biases, while the Climate Prediction Center (CPC) Merged Analysis of Precipitation (CMAP), constructed from gauge observations, satellite estimates, and numerical model outputs (Xie and Arkin 1998), was used to evaluate the precipitation biases. Both datasets were averaged in time, resulting in global maps of climatological values.
The annual mean SST biases relative to OISST are presented in Fig. 5. The results indicate that BESM-OA2.3 has an extensive negative SST bias over the tropical and subtropical oceans (Fig. 5a), much like both the HadCM3 model (Fig. 5b) and the CanCM4 (Fig. 5d). The exception is on the eastern boundaries of the Pacific and Atlantic Ocean, where the models are unable to properly reproduce the tilt of the equatorial thermocline (not shown), resulting in a warm bias of approximately 2 K. The pattern of a colder tropical North Atlantic and a hotter tropical South Atlantic is related to wind stress biases and has also been observed in other models (Muñoz et al. 2012).
The most outstanding cold bias occurs in the simulation of the North Atlantic in all four CMIP5 models in Fig. 5, which is also similar to GFDL CM2.5 (Delworth et al. 2012), since both CM2.5 and BESM-OA2.3 coupled models share the same ocean component, and on the Southern Ocean, where warm biases are predominant much like the NCAR CCSM4 (Figs. 5a,c). The model overestimates SSTs by as much as 3 K in eastern boundary currents (e.g., Humboldt, California, and Benguela Currents). One possible explanation for the spatial extent and magnitude of the warm bias in the eastern currents is that it results from positive feedback between the processes involved in the formation of marine stratocumulus clouds over cold waters and their cloud shading effect reducing the net surface radiative forcing (Ma et al. 1996; Siqueira and Nobre 2006; Zheng et al. 2011).
The SST biases in BESM-OA2.3 show no significant variation when biases from the ensemble mean of individual experiments decadal1960, decadal1980, and decadal2005 are intercompared. The biases are also indifferent to CO2 values, displaying no significant difference between the physics1 and physics2 integrations (figures not shown). These results suggest that the observed biases in SST are essentially determined by the model formulation.
Annual mean precipitation biases are shown in Fig. 6, following the same approach used for SST. While BESM-OA2.3 is able to reproduce global observed values of precipitation, it displays biases that arise from the inability to position the intertropical convergence zone (ITCZ) correctly, a common feature of all the four CMIP5 models in Fig. 6. The tropical Pacific is characterized in BESM-OA2.3 by a double ITCZ and a missing South Pacific convergence zone (SPCZ) (Fig. 6a). The tropical Atlantic is better reproduced and although the ITCZ is shifted toward the south, the SACZ has an improved representation compared to its Pacific counterpart. These deficiencies, of the double ITCZ in the Pacific and the southward shift of the ITCZ in the Atlantic, are known problems of current global models (Lin 2007; Richter and Xie 2008).
The biases originated from not representing the ITCZ properly include an excess of precipitation over the South Atlantic Ocean, Indian Ocean, and most of the South Pacific Ocean (with the exception of the SPCZ region). The equatorial Pacific and northern tropical Atlantic are characterized by negative biases in precipitation, while most of the North Pacific and North Atlantic Oceans show positive biases. Over land, the strongest signal comes from the Amazon basin, where the model is characterized by negative biases, similar to the HadCM3 (Fig. 6b) and CanCM4 (Fig. 6d). These biases are caused by the reduced convection over the Amazon forest and are possibly connected to the SST biases via the Walker circulation over the Atlantic Ocean, through trade winds and the tilt of the equatorial thermocline (Braconnot et al. 2007; Breugem et al. 2008; Richter and Xie 2008; Tozuka et al. 2011). In earlier versions of the coupled ocean–atmosphere model (e.g., Nobre et al. 2009), base for the BESM-OA coupled model, the trade winds over the Atlantic Ocean had their direction reversed during part of the year, producing strong positive SST biases when the equatorial thermocline collapsed, reducing convection over land. This has been improved in the current version of the model, mostly due to the Grell ensemble cumulus parameterizations used (Grell and Devenyi 2002), although some bias still remains. Similar to the SST, precipitation biases are time and physics invariant, showing no significant differences between periods or CO2 physics.
Figure 7 shows the seasonal cycle in the deviations of SST at the equator about the annual mean for OISST and four current CMIP5 models (top panel). The variation in observed SST is dominated by the annual harmonic that peaks at about 100°W, with the warm phase in April slightly stronger than the cold phase in September–October (Fig. 7a). All four models fail to capture the correct start and duration of the warm phase (Figs. 7b–e). The main part of seasonal variation in equatorial SST is realistically located over the eastern Pacific in all four models, although deviations are concentrated too close to the coast for BESM-OA2.3 (Fig. 7b) or extend too far west as in CanCM4 (Fig. 7e). The seasonal cycle is quite weak in HadCM3 and CCSM4 (Figs. 7c,d) and the simulated SSTs are correctly dominated by the annual harmonic, except in BESM-OA2.3 for which the semiannual harmonic is the dominant feature (Fig. 7b).
To investigate to what extent convection is taking place south of the equator in the eastern Pacific and its seasonality, Hovmöller diagrams of precipitation for CMAP and four current CMIP5 models averaged from 100° to 150°W are also plotted in Fig. 7 (bottom). The rainfall from CMAP observations is largely confined to an ITCZ located between 5° and 15°N (Fig. 7f). The seasonal variation in latitude for the observations is small with the ITCZ reaching 15°N around September and its southernmost position in March–April, almost reaching the equator. The rainfall maximum simulated by BESM-OA2.3 has realistic magnitude, but the precipitation is severely overestimated south of the equator (much like CSSM4 and CanCM4 in Figs. 7i,j) and produces a rather persistent double ITCZ throughout the year (Fig. 7g).
d. Interannual variability
1) El Niño–Southern Oscillation
The spatial patterns of tropical interannual SST variability for BESM-OA2.3 and observations are shown in Fig. 8. BESM-OA2.3 exhibits weaker equatorial Pacific variability than that observed in recent decades, correctly placing the center of the Pacific SST variability relative to the observed pattern, although with too little variability along the coast of Peru. It is worth noting that the tropical Atlantic SST simulations also show weaker interannual SSTA fluctuations than the observations.
The El Niño–Southern Oscillation (ENSO) variability is frequently characterized by indexes of SST variability for specific regions over the equatorial Pacific; however, the SST and wind patterns may differ among models and observation. Therefore, performing comparisons can be difficult. Nobre and Shukla (1996), for example, discussed the tropical Atlantic ocean–atmosphere variability through empirical orthogonal function (EOF) analysis (Lorenz 1956), joining the fields of SST and wind stress anomalies normalized by the local standard deviation prior to the EOF calculations, a technique referred to as “joint EOF” (jEOF) analysis. Hence, it is more convenient to characterize ENSO variability in terms of joint EOF of the SST and wind stress fields so that the model’s ENSO can be studied. In this case, ENSO variability is considered to be represented by the leading jEOFs of anomalous monthly SST and wind stress over a region that nearly spans the equatorial Pacific and by the time series of the corresponding principal components. In this analysis, the jEOFs are normalized to have unit spatial variance, and hence the variances of each principal component (PC) reflect spatiotemporal variances of corresponding SST and wind stress anomalies, with variances of all the PCs summing to the total variance. The jEOF analysis is presented in Fig. 9, where Fig. 9a shows the first jEOF of observed SSTs (Reynolds et al. 2002) and wind stress fields estimated with observed winds from NCEP–NCAR reanalysis (Kalnay et al. 1996), from January 1950 to December 2010 (60 yr).
The main pattern of variability associated with ENSO on the outputs of BESM-OA2.3 physics2 centennial run with fixed atmospheric CO2 concentration can be determined through the first jEOF of monthly SSTs and wind stress and is plotted in Fig. 9b for a 100-yr period. The first jEOF (jEOF1) reflects the mature phases of the ENSO life cycle: that is, the mature phases of El Niños and La Niñas. Despite broadly resembling the observed variability (Fig. 9a), we detect a strong meridional confinement of the SST anomalies around the equator and the center of action of this mode is displaced eastward when compared to the observations. The meridional confinement can be caused by a narrow meridional span of zonal wind stress anomaly and has implications in the interannual variability (Kirtman 1997; Deser et al. 2006). The first jEOF explains only 19% of variance versus over 38% in the observations.
The first jEOF pattern of BESM-OA2.3 in the equatorial Pacific does not change significantly over the 100-yr period for a twin physics1 simulation (not shown). This suggests that the ENSO variability is only affected in its characteristic period by weak global temperature increase.
The power spectra of the first principal component (PC1) for observed SST (solid black) and for the twin 100-yr simulations starting at 2005 physics1 (solid red) and physics2 (solid blue) are shown in Fig. 9c. The seasonal cycle and linear trends were removed when computing the jEOFs and the respective PCs to better isolate the modes of variability. The observed ENSOs have a relatively broad spectrum spanning 3–5 yr. The model captures this feature qualitatively with peak ENSO period in BESM-OA2.3 at approximately 5 yr, which is on the higher end of the 3–5-yr band (toward lower frequencies). Changes in ENSO frequency under CO2 increase can be noticed in Fig. 9c (solid red) and an overall tendency for ENSO period to decrease is evident, in accordance with Merryfield (2006) multimodel ensemble analysis.
Overall the model has a robust ENSO signal with multidecadal fluctuations in amplitude (not shown), and an irregular period between 3 and 5 yr with peak ENSO period at approximately 5 yr. The main ENSO pattern of SST variability simulated by the model is shifted about 30° east of the observed patterns, and the transition pattern matches the observations well with some deficiencies mainly over the equator (not shown). Such problems may be linked to the model mean state biases: namely, the equatorial cold bias, the double ITCZ, and the overstratified surface waters near the South American coast.
2) Extratropical teleconnections
The extratropical teleconnections are recurrent atmospheric circulation patterns associated with variations in the strength and location of the polar jet stream, the subtropical jet stream, or midlatitude storm tracks (Chang and Fu 2002). The Pacific–North American (PNA) teleconnection pattern is one of the most prominent and influential modes of low-frequency variability in the Northern Hemisphere midlatitudes beyond the tropics (Wallace and Gutzler 1981). The PNA pattern features large-scale above-average surface pressure anomalies that arc from the Hawaiian Islands and over the intermountain region of North America and below-average surface pressure located south of the Aleutian Islands and over the southeastern United States (Fig. 10a).
The PNA pattern is associated with noticeable fluctuations in the strength and location of the East Asian jet stream. The positive phase is associated with an intensified East Asian jet stream and with an eastward shift in the jet exit region toward the western United States. The negative phase is associated with a westward shrink of the jet stream toward eastern Asia and blocking activity over the high latitudes of the North Pacific.
Although the PNA and ENSO are distinct natural internal modes of climate variability, the PNA has been found to be strongly influenced by the ENSO phase (Straus and Shukla 2002). The positive phase of the PNA pattern tends to be associated to El Niños; during La Niñas, the PNA tends to be in the negative phase.
BESM-OA2.3 matches well the PNA pattern and its amplitude (Fig. 10b) with the center of action over the intermountain region of North America shifted about 20° east of the observed pattern. This mode of climate variability explains 12% of variance in the model versus over 17% in the observations.
In the Southern Hemisphere (SH), there is also one prominent low-frequency mode in the extratropics that resembles the PNA pattern, which is characterized by a well-defined wave train extending from the central Pacific to Argentina, and is referred to as the Pacific–South American (PSA) pattern (Ghil and Mo 1991; Mo and Higgins 1998; Lau et al. 1994). The center of action of the PSA teleconnection pattern is located in the west of the Antarctic Peninsula in the observations and comprises a wave train of alternate positive and negative geopotential height (or pressure) anomalies across the South Pacific emanating from the tropics (Fig. 11a).
It has been suggested that this teleconnection pattern that acts over South America impacts the South Atlantic convergence zone, which results in a regional seesaw pattern of alternating dry and wet conditions (Nogués-Paegle and Mo 1997). The PSA has also been found to be influenced by ENSO conditions (Mo and Nogués-Paegle 2001).
The model simulations show two main centers of action instead of one, a negative center located in the west of the Antarctic Peninsula similar to observations and another with similar amplitude over the Antarctic Peninsula. This mode of climate variability in the SH explains 18% of variance in the model versus over 14% in the observations.
3) Tropical Atlantic variability
The tropical Atlantic variability (TAV) is characterized by two main modes of variability, an ENSO-like equatorial mode (Zebiak 1993) and an interhemispheric gradient mode (e.g., Nobre and Shukla 1996). This last mode, also referred as the “Atlantic dipole mode” or the “Atlantic meridional mode,” is related to the cross-equatorial gradient of the SST anomalies and resultant migrations of the intertropical convergence zone. Consequently, it has an impact on rainfall over the tropical Atlantic, including the continental areas of the northeast Brazil and African’s Sahel region, at interannual and decadal climate scales (Doi et al. 2010, and references therein).
Several studies have employed the technique of EOF analysis to study covariant patterns of variability over the Atlantic Ocean (Venegas et al. 1997; Weare 1977). According to Nobre and Shukla (1996), from 1963 to 1987 the first jEOF mode spatial pattern of SST and wind stress anomalies showed a zonal, dipole-like structure in both SST and wind stress anomalies over the tropical Atlantic, suggesting that warmer SSTs are connected to weaker trade winds and colder SSTs are connected to stronger trade winds (Dommenget and Latif 2000). A similar analysis is presented in Fig. 12a, considering observed SSTs (Reynolds et al. 2002) and wind stress fields of the NCEP–NCAR reanalysis, from January 1960 to December 2010 (51 yr), where all linear trends and seasonal cycles were removed prior to the analysis. The spatial pattern of the first jEOF, accounting for 12.5% of the total variance (Fig. 12a), shows that the structure of the tropical Atlantic variability described by Nobre and Shukla (1996) was sustained until 2010, with the dipole-like feature in the SST pattern and analogous connections SST–trade winds. The power spectrum analysis of the amplitude time series of this mode (Fig. 12c) reveals a low frequency oscillation, which is represented as the main peak at 11.75 yr, confirming the decadal time scale behavior of this pattern as noticed by Nobre and Shukla (1996).
The occurrence of SST biases in the tropical Atlantic is a common issue of coupled models and, as discussed by Richter and Xie (2008) and Doi et al. (2010), almost all CMIP3 models show the zonal SST gradient along the equator opposite to what is observed. This deviation is possibly attributed to the lack of the cold tongue in the equatorial Atlantic in many coupled model simulations. Even recent coupled models, with higher horizontal resolution and updated parameterizations present SST biases in the eastern equatorial Atlantic and the Angola–Benguela area (Doi et al. 2012). BESM-OA2.3 also has the SST bias in the tropical Atlantic as discussed before (Fig. 5). Notwithstanding, it is able to simulate the Atlantic dipole mode quite well. Figure 12b shows the jEOF analysis results of the physics2 centennial run simulation, starting from December 2005 and with fixed CO2 concentration. The spatial patterns of SST and wind stress anomalies are very similar to the ones for observed data (Fig. 12a), except for the region near the Gulf of Guinea, where the negative SST patterns are stronger and the wind stress patterns are preferably from the southeast instead of from the south–southwest, reflecting the SST biases and the consequent discrepancies in the trade wind anomalies over this region. The main period of this simulated Atlantic dipole mode is 15 yr (Fig. 12c). For the twin centennial simulation with increasing CO2 of physics1, jEOF analysis of the resultant simulation also shows a dipole structure (figure not shown) similar to Fig. 12b, but with the period of maximum spectral energy density at 6.43 yr (Fig. 12c). Apparently, the variations of the atmospheric radiative balance caused by an increasing atmospheric CO2 concentration act to increase the Atlantic dipole frequency. Yet, the mechanisms by which it occurs still have to be further investigated.
4. Summary and conclusions
This article has explored the global features and regional aspects of the BESM-OA2.3 coupled model for climate change research. The numerical experiment design was set accordingly to the CMIP5 protocol and was aimed at explicitly revealing the model sensitivity to atmospheric CO2 concentration increase, as a representation of the buildup of greenhouse gases in the atmosphere. Extended runs with over 2000 yr of ensemble members were generated to scrutinize the model behavior for decadal climate variability and change. It was demonstrated that BESM-OA2.3 responds to increasing atmospheric CO2 concentrations in a consistent manner. The model simulations showed changes in the earth’s top-of-atmosphere radiation balance, an increase of the accumulation of energy in the earth system, and associated global sea surface temperature increase and marine ice reduction.
In terms of annual biases, the simulations show some of the same characteristics of other global coupled ocean–atmosphere models: namely, the double ITCZ and an almost absent SPCZ over the Pacific Ocean, in addition to the southward displacement of the Atlantic ITCZ. Concurrent with most global coupled ocean–atmosphere models, a warm SST bias is present along the western coasts of South America and Africa. Regarding precipitation, the model shows excess rainfall over the oceans and deficiencies over the continents, particularly over the Amazon basin.
The model biases showed no significant sensitivity to atmospheric CO2 concentrations, suggesting that such biases are a consequence of more basic model physics deficiencies. The excess rainfall over the oceans was diagnosed after the fact, and a new formulation for water vapor surface fluxes over the ocean has been developed at CPTEC and will be tested in a future experiment. Yet, the BESM-OA2.3 simulations presented an improved representation of the SACZ relative to previous versions of its base coupled ocean–atmosphere model, which we interpret as a direct consequence of the improvements on deep cumulus parameterizations used for this research.
On the interannual spectrum of climate variability, the simulations have shown a robust ENSO signal over the equatorial Pacific Ocean, with peak energy in the 3–5-yr range and a dipole-like pattern of SST and surface winds with peak energy on the decadal scale in the Atlantic Ocean, in agreement with observations. However, for the centennial model runs under increasing atmospheric CO2 concentrations the coupled ocean–atmosphere modes of variability over both the Pacific (e.g., ENSO) and the tropical Atlantic (e.g., the meridional mode) presented signals of enhanced periodicity, which may be an indication that, under growing atmospheric CO2 concentrations, interannual climate variability may grow. The model also represented atmospheric teleconnections in close agreement with observations, which is an important benchmark for its use for climate variability and change research.
BESM-OA2.3 has been used for operational seasonal predictions at CPTEC. It constitutes a building block on the Brazilian framework for climate change research [Brazilian Network on Global Climate Change Research (Rede CLIMA)].
The authors acknowledge Dr. Felipe M. Pimenta for his contributions to the development of a previous version of BESM-OA. This research was partially funded by FAPESP (2009/50528-6), MCTI/FINEP (01.08.0405.00), and CNPq (550990/2011-9; 573747/2008-0). The model integrations were done on Cray EX6 supercomputer at INPE, funded by the Brazilian Network on Global Climate Change Research–Rede CLIMA and FAPESP Research Program on Global Climate Change.