The physical climate formulation and simulation characteristics of two new global coupled carbon–climate Earth System Models, ESM2M and ESM2G, are described. These models demonstrate similar climate fidelity as the Geophysical Fluid Dynamics Laboratory’s previous Climate Model version 2.1 (CM2.1) while incorporating explicit and consistent carbon dynamics. The two models differ exclusively in the physical ocean component; ESM2M uses Modular Ocean Model version 4p1 with vertical pressure layers while ESM2G uses Generalized Ocean Layer Dynamics with a bulk mixed layer and interior isopycnal layers. Differences in the ocean mean state include the thermocline depth being relatively deep in ESM2M and relatively shallow in ESM2G compared to observations. The crucial role of ocean dynamics on climate variability is highlighted in El Niño–Southern Oscillation being overly strong in ESM2M and overly weak in ESM2G relative to observations. Thus, while ESM2G might better represent climate changes relating to total heat content variability given its lack of long-term drift, gyre circulation, and ventilation in the North Pacific, tropical Atlantic, and Indian Oceans, and depth structure in the overturning and abyssal flows, ESM2M might better represent climate changes relating to surface circulation given its superior surface temperature, salinity, and height patterns, tropical Pacific circulation and variability, and Southern Ocean dynamics. The overall assessment is that neither model is fundamentally superior to the other, and that both models achieve sufficient fidelity to allow meaningful climate and earth system modeling applications. This affords the ability to assess the role of ocean configuration on earth system interactions in the context of two state-of-the-art coupled carbon–climate models.
We describe the physical formulation and simulation characteristics of two new global coupled carbon–climate Earth System Models (ESMs) developed over the last several years at the Geophysical Fluid Dynamics Laboratory (GFDL) of the National Oceanic and Atmospheric Administration (NOAA) to advance understanding of how the earth’s biogeochemical cycles, including human actions, interact with the climate system (Fig. 1). The models are the product of an effort to expand upon the capabilities of past GFDL models used to study climate on seasonal-to-centennial time scales (e.g., Manabe et al. 1991; Rosati et al. 1997; Delworth et al. 2002; Delworth et al. 2006). Our approach has been to develop two Earth System Models with different ocean dynamical/physical cores while keeping all other components the same. The atmosphere and sea ice components in the new ESMs are very similar to those in GFDL’s previous Climate Model version 2.1 (CM2.1; Delworth et al. 2006). The new Land Model version 3.0 (LM3.0) component includes new hydrology, physics, and terrestrial ecology components. The ocean dynamical/physical component in CM2.1 was replaced by two new ocean components using the same ocean ecology and biogeochemistry—one in which the vertical coordinate is based on depth (ESM2M) and another in which the vertical coordinate is based on density (ESM2G).
These independent frameworks for describing the ocean allow us to explore the sensitivity of the ocean response under anthropogenic CO2 forcing and climate change to ocean configuration. In traditional depth coordinates, the ocean is divided into boxes in which horizontally adjacent pressures interact, typically emphasizing resolution near the sea surface. While this framework holds many benefits, it leads to fundamental problems in representing the lateral connectivity between vertical boxes, given their extremely skewed horizontal/vertical aspect ratio, and results in spurious numerical mixing and difficulty in representing bottom-following flows, such as dense overflows down topography. An isopycnal framework resolves these problems by taking advantage of the relative ease of motion along isopycnal layers and explicitly characterizing diapycnal transfer between these layers, but adds its own representational challenges, mostly by excluding resolution from weakly stratified regions. Comparison between these formulations allows us to assess the relative fidelity of these approaches and reduce uncertainty of ocean sensitivity to model assumptions. Previous comparisons of ocean-only models with different vertical coordinates include Willebrand et al. (2001) and Chassignet et al. (2003). Recently, Megann et al. (2010) compared two coupled models with identical atmosphere and land components, where a z-coordinate ocean model was replaced with a hybrid-coordinate (largely isopycnic) model, with the result enhancing tropical surface temperature, shoaling the thermocline, and decreasing deep temperatures that they inferred was in response to reduced numerical mixing. This interpretation is also supported by analysis of spurious mixing in analysis of the ocean components discussed here (Ilicak et al. 2012).
We built ESM2M and ESM2G specifically to study carbon–climate interactions and feedbacks in the context of global climate change under the influence of increased greenhouse gases. Results obtained from these model integrations will be part of the Coupled Model Intercomparison Project version 5 (Taylor et al. 2012). This paper serves as the first of two companion papers. This work documents the physical climate formulation and baseline simulation characteristics of these two ESMs. It also discusses the initialization methods used to obtain stable climate states. Finally, it presents a few results from the preindustrial control integrations (nominally 1860). Dunne et al. (2012, manuscript submitted to J. Climate, hereafter Part II) document the ecological and biogeochemical formulation and baseline simulation characteristics of these models.
2. Model description
The component models pass fluxes across their interfaces using an exchange grid system. The exchange grid enforces energy, mass, and tracer conservation on the fluxes passed between the component models. The atmosphere, land, and sea ice radiative components are coupled every 30 min. The ocean tracer and atmosphere–ocean tracer-coupling time steps are 2 h.
The Atmospheric Model, version 2 (AM2) is virtually identical to that in CM2.1, incorporating only a few minor code updates that were found to not appreciably affect climate. AM2 uses a 2° latitude × 2.5° longitude horizontal grid with 24 vertical levels on a D grid using finite-volume advection (Lin 2004) with a 30-min dynamical time step and 3-h radiation time step. The atmospheric physical parameterizations are described in Anderson et al. (2004). The radiation parameterizations are also unchanged from CM2.1. Because the ESM2M ocean is very similar to the CM2.1 ocean, the fact that neither the ESM2M nor ESM2G atmosphere was retuned to improve initial radiative balances gives ESM2M a relative tuning advantage. The orbital parameters are held constant at 1860 values for the computation of incoming solar radiation. Radiatively active trace gases and aerosols are held constant at 1860 values. Three-dimensional, time mean distributions of aerosols (i.e., black carbon, organic carbon, sulfates, dust, sea salt, and volcanoes) are specified as input concentration fields (V. Naik and D. Schwarzkopf 2011, personal communication).
The land model in both ESMs (LM3.0) represents land water, energy, and carbon cycles. LM3.0 includes a multilayer snowpack above the soil; continuous vertical representation of soil water, spanning both unsaturated and saturated zones; frozen soil water; parameterization of water table height, saturated-area fraction, and groundwater discharge to streams derived from groundwater hydraulics and surface topography; finite-velocity horizontal transport of river runoff; lakes, lake ice, and lake-ice snowpacks that exchange mass and energy with both the atmosphere and rivers; and consistent, energy-conserving accounting of sensible heat content of water. Temperature is tracked in the vegetation canopy and leaves, and in multiple soil/snow layers. Vegetation radiation, hydrology, and carbon dynamics are tracked by an extension of Shevliakova et al. (2009). For radiative transport, the vegetation canopy is treated as a dispersed medium where leaf reflectance is dependent on leaf properties and intercepted snow. The canopy is underlain by a reflective, opaque surface of soil and/or snow, with snow extent depending on snow depth and sufficiently deep snow masking a fraction of the canopy. Bare soil and snow radiation parameters were specified from the Moderate Resolution Imaging Spectroradiometer (MODIS) Bidirectional Reflectance Distribution Function (BRDF)/Albedo Snow-Free Quality product parameters [MOD43C2 V004; see https://lpdaac.usgs.gov/lpdaac/products/modis_products_table/brdf_albedo_snow_free_quality/16_day_l3_global_0_05deg_cmg/mcd43c2; Land Processes Distributed Active Archive Center (LP DAAC), U.S. Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center (https://lpdaac.usgs.gov/)].
d. Sea ice
The sea ice component (Winton 2000) is very similar to that in CM2.1. It has full ice dynamics, three-layer thermodynamics (two ice and one snow layer), and five different ice thickness categories plus open water on the same tripolar grid as the B-grid ocean component (Murray 1996) of Modular Ocean Model, version 4p1 (MOM4p1). The ice albedos (snow on ice albedo = 0.85, ice albedo = 0.65 compared to snow on ice albedo = 0.80 and ice albedo = 0.5826 in CM2.1) and temperature range of melting (1°C compared to 10°C in CM2.1) were modified based on observations (Perovich et al. 2002). In addition, sea ice diffusion was added for narrow (one point wide) channels in ESM2G.
Icebergs are supplied through the ice-calving term in LM3.0. Whenever snow depth in LM3.0 exceeds a critical value, then excess snow is transported to the ocean–sea ice component in coastal cells via rivers. This rudimentary representation of land ice calving is accumulated and released as icebergs according to a prescribed size distribution. Within the iceberg model, clusters of icebergs are treated as Lagrangian particles with trajectories determined by drag by the atmosphere, ocean, and sea ice (Martin and Adcroft 2010). The icebergs impart weight on the ocean, deposit fresh and 0°C meltwater, but do not affect the ocean–ice–atmosphere radiation or momentum balance.
f. MOM4p1 (ESM2M ocean)
The ocean component of ESM2M employs the MOM4p1 code of Griffies (2009) configured with the same grid and bathymetry as the CM2.1 ocean component [Griffies et al. (2005) and Gnanadesikan et al. (2006); a 1° horizontal grid up to ⅓° meridionally at the equator and tripolar above 65°N, with 50 vertical levels]. ESM2M ocean uses a rescaled geopotential vertical coordinate (z*; Stacey et al. 1995; Adcroft and Campin 2004) for a more robust treatment of free surface undulations. For tracer advection, ESM2M uses the conservative, minimally diffusive, monotonic, multidimensional piecewise parabolic method (MDPPM) ported from the MIT General Circulation Mode (MITgcm; Marshall et al. 1997). Neutral (isopycnal) diffusion is based on Griffies et al. (1998), with constant diffusivity of 600 m2 s−1 and a slope-tapering scheme (Danabasoglu and McWilliams 1995) as in CM2.1, but uses a larger value of for the maximum slope (Gnanadesikan et al. 2007a) to enhance abyssal flow out of the Southern Ocean (Downes et al. 2011). For eddy parameterization, ESM2M uses the skew-flux approach of Griffies (1998), but computes the quasi-Stokes streamfunction via a boundary value problem extending across the full column after Ferrari et al. (2010), rather than locally (Gent and McWilliams 1990; Gent et al. 1995). As in CM2.1, horizontal variation of eddy diffusivity in ESM2M is based on local flow properties (Griffies et al. 2005), but with an expanded allowable range of 100–800 m2 s−1. The ESM2M implementation includes updates to the K-profile parameterization (Large et al. 1994) based on Danabasoglu et al. (2006), as well as model-predicted chlorophyll modulation of shortwave radiation penetration through the water column.
ESM2M also includes completely novel parameterizations relative to CM2.1, such as parameterization of submesoscale eddy-induced mixed layer restratification (Fox-Kemper et al. 2011). Instead of prescribed vertical diffusivity for interior mixing (Bryan and Lewis 1979), ESM2M employs the Simmons et al. (2004) scheme along with a background diffusivity of 1.0 × 10−5 m2 s−1 in the tropics and 1.5 × 10−5 m2 s−1 poleward of 30° latitude following a tanh curve. ESM2M implements geothermal heating following Adcroft et al. (2001). While CM2.1 used a horizontal anisotropic friction scheme (Large et al. 2001), ESM2M uses an isotropic Laplacian friction and western boundary–enhanced biharmonic friction with less frictional dissipation to allow more vigorous tropical instability wave activity at the expense of adding zonal grid noise, particularly in the tropics.
g. GOLD (ESM2G ocean)
The ocean component of ESM2G employs the Generalized Ocean Layer Dynamics (GOLD) model isopycnal code originally developed by Hallberg (1995), with a 1° horizontal grid up to ⅓° meridionally at the equator. The model includes two mixed layer layers, two buffer layers, and 59 interior layers to provide a reasonable resolution across both the narrow density range of the Southern Ocean and the broad dense-to-light water structure from bottom waters to low latitudes and marginal seas. A three-way barotropic, baroclinic, and diabatic/tracer split time-stepping scheme (Hallberg 1997; Hallberg and Adcroft 2009) with subcycling of shorter time steps ensures that the model can be integrated efficiently while conserving momentum; mass; tracers, including potential temperature and salinity; potential vorticity; and energy (Arakawa and Hsu 1990). Monotonic tracer advection is based on Easter (1993). The interior isopycnal layers track prescribed values of potential density with a 2000-dbar reference pressure (Sun et al. 1999); vertical advection in the diapycnal mixing scheme accounts for density drifts resulting from cabeling (Hallberg 2000). For all aspects of the dynamics, the full nonlinear equation of state is used. To avoid thermobaric instability arising from the traditional Montgomery potential form of the pressure gradient force (Hallberg 2005), ESM2G uses the analytically integrated finite-volume pressure gradient discretization of Adcroft et al. (2008).
Diapycnal diffusivity between interior layers in ESM2G combines multiple parameterizations. The background diapycnal diffusivity is prescribed using the latitude-dependent profile proposed by Henyey et al. (1986) and described in Harrison and Hallberg (2008) of 2 × 10−5 m2 s−1 at 30°N, decreasing to a minimum of 2 × 10−6 m2 s−1 at the equator. Turbulent entrainment in descending overflows and other regions of small resolved shear Richardson number is parameterized after Jackson et al. (2008). Additionally, the bottom boundary layer is mixed using the parameterization of Legg et al. (2006), in which 20% of the energy extracted by bottom drag is used to thicken a well-mixed bottom boundary layer. As in ESM2M, ESM2G uses the baroclinic tide mixing of Simmons et al. (2004) and geothermal heating after Adcroft et al. (2001). To maintain the weakly stratified abyss, ESM2G uses a floor on turbulent dissipation [ɛmin = 1.0 × 10−7 m2 s−3 + 6.0 × 10−4 m2 s−2 N, where N is buoyancy frequency (Gargett and Holloway 1984)] combined with an assumed flux Richardson number of FRi = 0.2N2/(N2 + Ω2) for a minimum diapycnal diffusivity of κMin = 0.2N2ɛMin/(N2 + Ω2).
ESM2G incorporates a bulk mixed layer based on Hallberg (2003) and Thompson et al. (2003) with two additional buffer layers for smooth water mass exchange with the isopycnal interior. The turbulent kinetic energy budget determines mixed layer depth (MLD) after Kraus and Turner (1967) divided into two layers that are well mixed with respect to tracers but not velocity. Submesoscale eddy-driven restratification of the mixed layer is parameterized after Fox-Kemper et al. (2011) to avoid mixed layer depths becoming excessively deep (Hallberg 2003). The new two–buffer layer scheme allows the mixed layer to retreat and advance without introducing unphysical surface property biases. Like ESM2M, ESM2G includes dynamic representation of model-predicted chlorophyll modulation of shortwave radiation penetration through the water column. ESM2G uses a thickness diffusion parameterization motivated by Gent and McWilliams (1990) with a local slope-dependent coefficient varying between 10 and 900 m2 s−1. As in ESM2M, momentum dissipation is implemented with Laplacian and biharmonic viscosities. The Laplacian coefficient is Δx multiplied by 1 cm s−1. The biharmonic viscosity is the larger of Δx3 multiplied by 5 cm s−1 and a Smagorinsky (1993) viscosity with a nondimensional coefficient of 0.06 (Griffies and Hallberg 2000).
In contrast to the smoothed and otherwise modified topography in ESM2M (Griffies et al. 2005), bathymetry in ESM2G is simply the areal depth average of high-resolution bathymetry adjusted to give the full depth for sills. While ESM2M uses “cross-land mixing” to diffusively connect the properties across 8 such straits (Griffies et al. 2005), ESM2G represents explicit exchanges across the open faces of 14 straits restricted to approximate their actual widths (e.g., the Gibraltar and Bosporus Straits were restricted to 12 and 2.5 km, respectively). While ESM2M uses a constant of 600 m2 s−1, ESM2G uses a spatially variable 50–900 m2 s−1 with the buoyancy frequency and isopycnal slope.
To obtain 1860 initial conditions for ESM control and perturbation integrations, a method similar to that described by Stouffer et al. (2004) is used. The ocean model is initialized with present-day ocean temperature (Locarnini et al. 2006) and salinity (Antonov et al. 2006) from World Ocean Atlas 2005 and run for 1 yr forced by atmospheric conditions from a CM2.1 1990 control run. The fully coupled models were then integrated over 1000 model years with 1860 solar and radiative forcing before declaring “quasi-steady-state equilibrium” and beginning the 1860 control and perturbation integrations. In addition to the many qualitative requirements, we define acceptable quasi-steady-state equilibrium with quantitative metrics: net top-of-the-atmosphere radiative fluxes less than 0.5 W m−2, surface temperature drifts less than 0.1°C century−1, stable Atlantic meridional overturning circulation (AMOC; Delworth et al. 1993) above 10 Sv (1 Sv ≡ 106 m3 s−1), local sea surface temperature (SST) biases less than ~9°C, global 70°S–70°N root-mean-square SST errors less than 1.9°C, and global net CO2 fluxes between the atmosphere and both land and ocean lower than 20 PgC century−1 (averaged over two centuries). For accelerated initialization of the terrestrial soil component, years 551–600 of the integration were used to project equilibrium values offline (Shevliakova et al. 2009).
Here we compare the model control integrations to observations that are mainly taken near the present day. We do this knowing that 1860 conditions were likely colder than the present day, but focus our analysis on those differences that are expected to be larger than the forced climate changes over the last 150 yr. We show results averaged over years 1001–1100 of the integrations unless otherwise noted.
a. Climate drift from present-day initialization
When initialized with present-day ocean observations and given 1860 radiative forcing, the climate in both models cools (Fig. 2). The net radiation at the top of the model atmosphere (TOA; Fig. 2a) and the oceanic heat flux (Fig. 2c) become strongly negative (implying cooling of the climate system and the ocean, respectively). In response, the surface air temperature (SAT; Fig. 2b) and SST (Fig. 2d) cool in both models by more than 1°C in the first century. In subsequent centuries, however, the SAT and SST warm about 0.6°C in ESM2M and 0.3°C in ESM2G. After detrending, the models have similar total interannual variability (σ = 0.13°C for both models), with SAT variability largely following SST variability (r2 = 0.9 for both models) and ESM2G partitioning more variability into the multidecadal range.
While the initial evolution of the surface temperature from the present-day to 1860 conditions in the two models is quite similar, the overall ocean heat flux and subsequent ocean volume average temperature responses are quite different. In ESM2M, the TOA and oceanic heat flux reverse sign and begin to warm the ocean as a whole after year 11. The long-term drift in volume-averaged temperature (Fig. 2e) is 0.60°C over the 1000-yr period. By year 1000, the rate of warming in ESM2M diminishes to 0.038°C century−1. This deep ocean warming phenomenon is common to earlier versions of this model (Delworth et al. 2006), as well as many other atmosphere–ocean models (Randall et al. 2007, see their supplemental material), and tends to occur as a typical z-coordinate model drift away from the observed initial condition rather than as a consequence of the particular radiative forcing, which is qualitatively similar to the result of Megann et al. (2010). Ilicak et al. (2012) quantified the spurious numerical mixing in both models, finding the much larger values in ESM2M consistent with its warming trend relative to ESM2G. In ESM2G, the TOA and ocean surface fluxes gradually approach zero over the first 50 yr, but the ocean heat flux remains slightly negative throughout the 1000-yr period with the volume-averaged temperature (Fig. 2e) cooling by 0.18°C. By year 1000, the rate of cooling in ESM2G diminishes to −0.010°C century−1.
Both models undergo a spinup phase of the AMOC over the first century followed by a slight drop over the second and third centuries before becoming stable (Fig. 2f). After spinup, AMOC in ESM2M is stronger (28.1 ± 1.8 Sv) than in ESM2G (23.6 ± 1.7 Sv). As was seen for SAT and SST, AMOC variability in ESM2M is characterized by a relatively strong, regular 18-yr period while that in ESM2G is longer (26 yr) and less regular. This highlights the important role of ocean formulation in determining modes of climate variability.
The vertical structure of global ocean temperature drift (Fig. 3a) is characterized by a thermocline deepening in ESM2M (red) as it warms over much of the water column relative to present-day observations (black) and by thermocline shoaling in ESM2G (green) as the upper 2000 m cools with warming below. Both models warm extensively in the deep Atlantic (Fig. 3c), while ESM2M also warms extensively in the Southern Ocean (Fig. 3e). In all oceans, biases and root-mean-square errors in ESM2G are found to be significantly lower than in ESM2M. The global thermocline volumes (>8°C) in ESM2M (173 × 1015 m3) and ESM2G (131 × 1015 m3) straddle the observations [145 × 1015 m3 for the present day, and 137 × 1015 m3 after accounting for 0.5°C warming over the historical period (Levitus et al. 2005)].
ESM2M biases and root-mean-square errors in global ocean salinity profile (Fig. 3b) are smaller than those found in ESM2G. Both models underestimate the shallow salinity maximum and redistribute salinity into the deep ocean through the creation of overly saline (Fig. 3d) and overly warm (Fig. 3c) North Atlantic Deep Water (NADW) with approximately the correct density. In the Southern Ocean, the observed upper water column inversion of the halocline is poorly simulated in ESM2M and entirely absent in ESM2G (Fig. 3f). These biases are dominated by atmospheric model precipitation biases and are common in phase 3 of the Coupled Model Intercomparison Project (CMIP3) class of models (Randall et al. 2007, see their supplemental material).
b. Shortwave radiation
Some of the major model accomplishments and continued challenges in the development of these models can be seen in model representation of the shortwave radiation budget. One of the most robustly interpretable satellite measurements is the shortwave albedo from the perspective of the top of the atmosphere (the flux off the earth normalized to the flux the earth receives by the sun; Fig. 4a), because it requires none of the complicating factors involved in atmospheric corrections. Beyond the excellent global-scale agreement between the present-day observational and preindustrial modeled fields (r2 = 0.93 in ESM2M and r2 = 0.92 in ESM2G), local biases (Figs. 4b,c) reveal large scale patterns of 1) overly high albedos over the Rocky Mountains; 2) overly high albedos in the Bering, Greenland, and Barents Seas associated with sea ice (which is appropriate given that the model is preindustrial); 3) overly low albedos in the Southern Ocean, particularly the Atlantic sector; 4) overly low albedo near the coasts of California, Peru, Chile, Angola, and Namibia, and low albedo in the eastern tropical ocean basins associated with excessive low cloud; and 5) farther offshore from these coasts, overly low albedo in the eastern tropical ocean basins, associated with too little low cloud.
As an indication of the role of cloud-driven differences, we show the downward radiation impinging on the land and ocean surface in Fig. 4d. The major biases (Figs. 4e,f) include the Amazon and India having far too little cloud cover, resulting in overly high shortwave radiation reaching the surface. Over the ocean, the eastern boundary current upwelling areas and Southern Ocean receive too much radiation. The Southern Ocean shortwave bias results in a lack of sea ice and low surface albedo (Fig. 4, bottom row), with major implications for the representation of the seasonal mixed layer dynamics. The north equatorial and western Pacific regions also receive too much radiation, and the south equatorial Pacific region receives too little associated with the difficulties in representing the intertropical convergence zone (ITCZ). On the land, the primary surface albedo biases are the northward extension of the African tropical region, low albedo all along the boreal polar region, and high albedos associated with excessive Northern Hemisphere sea ice.
c. Surface climate
The surface ocean patterns of temperature (SST) are preserved well in both ESM2M and ESM2G (Fig. 5, top row; r2 of 0.99 and 0.98, respectively). The major differences relative to present-day observations are a general cooling consistent with their being forced by 1860 radiative conditions, a general hemispheric bias of cold north subpolar Pacific and warm Southern Ocean, an equatorial cold Pacific bias, and eastern boundary condition warm biases. These problems are all similar to biases seen in GFDL CM2.1 and other models of this class (Delworth et al. 2006; Wittenberg et al. 2006; Randall et al. 2007, their supplemental material; Richter et al. 2012), with the biases in ESM2G (Fig. 5c) being slightly larger in each case. The smaller surface temperature biases in ESM2M may relate to its greater similarity to the ocean in CM2.1 for which the atmosphere was originally developed and radiatively tuned. Both models do a fairly good job at representing sea surface salinity (SSS) structure (Figs. 5d–f; r2 = 0.80 for ESM2M and r2 = 0.76 for ESM2G). SSS bias patterns are quite similar between the two models, with ESM2G (Fig. 5f) again manifesting larger biases than ESM2M (Fig. 5e). M. J. Harrison et al. (2011, unpublished manuscript) show that these surface salt biases in the western North Atlantic (Figs. 5d–f), which are driven by low biases in Atlantic-draining land precipitation and related weak Amazon outflow, in turn drive the deep Atlantic temperature and salinity biases (Fig. 3c,d).
As an overall indicator of the models’ abilities to capture large-scale circulation, sea surface height (SSH) structure is well preserved in both models (Figs. 5g–i; r2 = 0.96 for ESM2M and r2 = 0.92 for ESM2G), with the major bias common to both models being a lower North Atlantic as a consequence of its salinification. ESM2M and ESM2G exhibit opposing biases in meridional SSH gradient and Antarctic Circumpolar Current (ACC), with Drake Passage transports being low in ESM2G (108 Sv) and high in ESM2M (164 Sv) compared to the observational estimates (140 ± 6 Sv; Ganachaud 2003). This difference between the models is most strongly seen southeast of New Zealand, where the strength and path of the ACC meanders show opposite biases relative to observed SSH. These differences, as well as differences with CM2.1 (132 Sv; Griffies et al. 2005) can be explained in terms of differences in lateral mixing parameterizations.
Overall, both models reproduce present-day observed precipitation patterns (Figs. 6a–c), particularly well in Africa and Oceania. However, they suffer from a variety of biases, including the common “double ITCZ” problem of too much precipitation in the central and eastern Pacific south of the equator (Lin 2007) and the “dry Amazon” problem of too little precipitation in both the equatorial Atlantic and equatorial South America (Delworth et al. 2006). In addition, ESM2G has too little precipitation in the central equatorial Pacific. Precipitation in both models is low in midnorthern latitudes and high in midsouthern latitudes in accordance with the interhemispheric SST bias (Figs. 5a–c).
Various precipitation biases are further highlighted in precipitation variability (Figs. 6d–f). While the models do not represent the magnitude of observed Amazon precipitation, both show strong Amazon precipitation variability with heavy precipitation in boreal winter. While both models represent the patterns of precipitation magnitude in the Indian Ocean, ESM2G has much better precipitation variability than ESM2M in this region. Conversely, while both models are biased low in central and western equatorial Pacific precipitation, variability is too high in ESM2M and too low in ESM2G as a consequence of their differing equatorial cold biases (Figs. 5a–c) and El Niño–Southern Oscillation (ENSO) statistics shown via the Niño-3 index (SST in 5°S–5°N, 150°–90°W) in Fig. 7.
Like CM2.1 (Wittenberg et al. 2006; Wittenberg 2009), both models exhibit a broad spectrum of ENSO variability consistent with observations, with amplitude and frequency modulation on multidecadal time scales. Relative to observations, however, ESM2M partitions too much variability into interannual modes relative to the seasonal cycle (period = 1 yr), while ESM2G does not partition enough variability into interannual modes. While previous work has suggested a dominant role for the atmospheric component (Guilyardi et al. 2004, 2009), our results highlight the crucial role of ocean formulation on ENSO behavior. This difference is likely due to the difference in climatological mean state simulations of west Pacific SST (Figs. 5a–c) and surface mixing. Relative to ESM2M, the stronger west Pacific SST gradient in ESM2G prevents convection over Indonesia from spreading eastward during El Niño, thus weakening the coupling between SST (central Pacific warming) and wind response (westerly wind anomalies). This westward-shifted wind response in ESM2G also reduces the time for wind curl–induced ocean Rossby waves to impact the western boundary, thus shortening the time available for equatorial thermocline adjustment to these wind anomalies (Kirtman 1997; Wittenberg 2002; Capotondi et al. 2006), and further weakening the El Niño and shortening its period. While the climatological eastern equatorial Pacific SST in ESM2G is not very much higher than in ESM2M (Figs. 5b,c), the subsurface is much colder, allowing for ESM2G to tap into a stronger thermal contrast for the seasonal upwelling and stronger seasonal cycle of Niño-3 SST.
d. Ocean meridional mass and heat transport structure
Both ESM2M and ESM2G match observations of NADW formation and circulation by generating warm, salty Atlantic midlatitude water flowing northward near the surface that becomes much denser as it is cooled by the atmosphere in high northern latitudes, sinks to depth, and flows southward through and out of the basin (Figs. 8a–c). While ESM2M (35 Sv) and ESM2G (22 Sv) span the observational estimates for the Florida Strait transport (31 ± 1 Sv; Ganachaud 2003; Schott et al. 1988), ESM2G builds a robust Gulf Stream transport (51 Sv) comparable to that in ESM2M (49 Sv) by 27°N. Comparison with zonally integrated NADW profiles at 24°N show both models to have higher and shallower maximum overturning (22 Sv) than estimates from Talley (2003) and Ganachaud and Wunsch (2003; 17.2–19.8 Sv). Under historical radiative scenarios, however, both models fall to 19 Sv [not shown; see Stouffer et al. (2006) for the CM2.1 response]. While ESM2G agrees well with observations for the streamfunction depth scale, the southward flow is relatively shallow in ESM2M (Fig. 8d). The deeper and more realistic southward NADW flow in ESM2G is consistent with previous comparisons of isopycnal and z-coordinate models (Willebrand et al. 2001; Megann et al. 2010), relating to the isopycnal models’ ability to maintain realistic bottom flows (Legg et al. 2006). By 30°S (Fig. 8e), this distinction is largely removed. Antarctic Bottom Water also flows into the Atlantic and upwells to about 3-km depth before flowing southward back out of the basin. In the North Pacific at 24°N (Fig. 8f), ESM2G has a stronger deep circulation than ESM2M. Compared to observationally based estimates of the Indonesian Throughflow of 8–17 Sv (Figs. 8g,h), ESM2M gives a value of the high end (16 Sv), while the ESM2G value is excessive (22 Sv). ESM2G has 16 Sv of Southern Ocean transport into the deep Pacific, while ESM2M has only 6 Sv (Fig. 8g). In the Southern Indian Ocean (Fig. 8h), both models have smaller deep flows than those observed. Taking the Indo-Pacific as a whole (Fig. 8i), ESM2G behaves similarly to the stronger Talley (2003) estimate, while ESM2M behaves similarly to the weaker Ganachaud (2003) estimate.
Both ESM2M and ESM2G match observations of global ocean heat well. The global ocean heat transports from the two ESMs are much closer to each other than are the different observational estimates (Fig. 9). In the Northern Hemisphere, the model transports lie within the observational estimates, though the initiation of northward transport near the equator is shifted southward in both models by approximately 5°. Atlantic meridional heat transport at 26.5°N in ESM2G (1.14 PW) and ESM2M (1.07 PW) is slightly lower than the Rapid Climate Change program (RAPID) estimate (1.33 ± 0.40; Johns et al. 2011). We attribute this bias to a lack of northward extent of subtropical gyre circulation in these models. In the Southern Hemisphere, both models lie within the observational estimate of Ganachaud and Wunsch (2003), but have less southward heat transport than National Centers for Environmental Prediction (NCEP) and European Centre for Medium-Range Weather Forecasts (ECMWF) reanalyses. We suspect these errors are related to problems in simulation of the cloud distribution, height, and brightness in the Southern Ocean (e.g., Delworth et al. 2006).
e. Sea ice
Both models agree fairly well with the observed minimum in sea ice extent in both hemispheres (Fig. 10). ESM2G is too extensive in Northern Hemisphere winter while ESM2M has a low bias of similar magnitude in Southern Hemisphere winter. Relative to observations, both models have larger ice extent in the Northern Hemisphere. This is consistent with common climate model bias patterns of the warm Southern Hemisphere–cold Northern Hemisphere (Fig. 5; see Randall et al. 2007, Fig. 8.2) present in CM2.1 (Delworth et al. 2006). The larger influence of ocean formulation on winter ice cover is consistent with studies showing that ocean circulation strongly controls winter sea ice extent (e.g., Winton 2003).
f. Ocean mixed layer dynamics
One of the fundamental differences between ESM2M and ESM2G is their treatment of the ocean surface mixed layer. ESM2M uses the vertically resolved K-profile parameterization (Large et al. 1994) while ESM2G uses a bulk energetic parameterization (Kraus and Turner 1967). Resultant explicit diffusivities in ESM2M are about a factor of 3 higher than ESM2G in the upper 200 m, but fall to extremely low explicit values in the thermocline in both models. While ESM2M requires the upper three 10-m levels to interact, ESM2G represents mixed layers as shallow as 2 m. To compare these models with observations, we show climatologies for monthly minimum and maximum MLD using the Levitus (1982) 0.125 density criterion in Fig. 11. Comparison of minimum MLD highlights the relatively shallow summer Southern Ocean mixed layer values in these models associated with the surface shortwave bias (Figs. 4e,f) and warm SST bias (Figs. 5b,c), with ESM2G suffering more so than ESM2M. In addition, ESM2G exhibits a stronger zonal equatorial gradient relative to ESM2M. ESM2G’s extremely shallow values in many areas (e.g., 5–10 m near Peru) highlight the capacity for ESM2G to maintain average mixed layers that are much shallower than those of ESM2M.
Comparison of maximum MLD (Figs. 11d–f) illustrates deep mixing in the North Atlantic of both models consistent with their robust NADW formation (Fig. 8). In the Southern Ocean, however, both models exhibit a lack of deep mixing, with observed MLDs at 50°S reaching 717 m, while those in ESM2M reach only 653 m and those in ESM2G only 559 m. These relative bias intensities are reversed in the western North Pacific because observed MLDs do not exceed 232 m, while those in ESM2M reach 422 m and ESM2G only 358 m. Overall, ESM2G (r2 = 0.39) correlates better with observations than ESM2M (r2 = 0.24) for minimum MLD, while both correlate equally with maximum MLD (r2 = 0.41).
We further examine the different interior character of these two models through their propensities to create and maintain mode waters—minimally stratified zones within the main thermocline with distinct minima in vertical density gradient (also referred to as a thermostad or pycnostad). Mode waters form through Ekman pumping into deeply mixed water columns where they are isolated from the surface and advect equatorward within the main thermocline. This mechanism serves to complement deep water formation to ventilate the ocean interior at intermediate depths (e.g., Talley et al. 2003), and supply the majority of nutrients for tropical productivity (e.g., Sarmiento et al. 2004).
Maps of the depth gradient in potential density (dρ0/dz; kg m−3 km−1) are shown in Fig. 12 for 300 m, 600 m, and a meridional section along 160°W for World Ocean Atlas 2009 (Figs. 12a,d,g, respectively), and ESM2M (Figs. 12b,e,h, respectively) and ESM2G (Figs. 12c,f,i, respectively). At 300 m, both models do a reasonable job at representing regional variability in stratification, with ESM2G (Fig. 12c; r2 = 0.77 and bias = 0.35) having slightly higher covariance with the observations (Fig. 12a) but larger bias than ESM2M (Fig. 12b; r2 = 0.74 and bias = −0.068). Defining mode waters with a stratification cutoff of 0.6 kg m−3 km−1, ESM2M reproduces the observed areas of mode water in the Southern Ocean between 40° and 50°S while ESM2G has too little. By 600 m, the distinction seen between the two ocean models is largely erased because neither ESM2M (Fig. 12e) nor ESM2G (Fig. 12f) maintain the low stratification (mode water) regions observed throughout the Southern Ocean near 40°S, though ESM2M shows faint traces of it that boost the r2 from only 0.453 in ESM2G to 0.605 in ESM2M. Meridional sections along 160°W highlight the scope of these differences (Fig. 12). While the observed low stratification region between isopycnals σθ = 27.2 and σθ = 26.7 extends down to 900 m at 50°S, these areas are restricted to depths shallower than 400 m in ESM2M and are not well pronounced at all in ESM2G. Overall, the Southern Ocean Mode Water volume above σθ = 27.2 in the observations (19 × 1015 m3) is reproduced partially in ESM2M (14 × 1015 m3) and only a small amount in ESM2G (2 × 1015 m3).
This dissimilarity between these models has the opposite implications for relative fidelity when considering the central North Pacific because ESM2M builds an excessive pool of minimally stratified mode water in the central North Pacific above σθ = 26.7 (1.9 × 1015 m3) that is not expressed in ESM2G (Fig. 12i). North Atlantic Mode Water volume above σθ = 27.6 in ESM2G (2.2 × 1015 m3) is only slightly lower than that of the observations (2.5 × 1015 m3), while that of ESM2M (5.3 × 1015 m3) is approximately double. This propensity of GFDL’s MOM-based models to create unrealistically large volumes of mode water in the Northern Hemisphere was previously described by Gnanadesikan et al. (2007b) in the context of potential vorticity.
g. Ocean ideal age distributions
Ocean ideal age (Thiele and Sarmiento 1990) serves as a valuable tracer of ocean ventilation as a combined analog of chlorofluorocarbon, radiocarbon, and other tracers. After 1100 yr of integration, ESM2M and ESM2G have very similar basin-averaged age (Fig. 13) in both the Atlantic (230 yr for ESM2M and 224 yr for ESM2G) and Pacific (600 yr for ESM2M and 650 yr for ESM2G). The vertical age structure, however, strongly differs between the two models, with ESM2G having a much older upper water column but younger deep waters than ESM2M. The deep patterns are a direct consequence of differences in meridional overturning (Fig. 8). In the North Atlantic, relatively old bottom water persists in ESM2M, while ESM2G maintains bottom water properties to more actively ventilate the abyss. Similarly, in the Pacific, the larger supply of Indo-Pacific bottom waters in ESM2G propagates farther northward than in ESM2M. These differences are consistent with the isopycnal ESM2G’s more robust ability to maintain bottom flows.
Much of the upper water column age difference between these two models can be tracked to the shallower thermocline in ESM2G than ESM2M (Fig. 13, bottom row). Even though the thermocline volume (>8°C) is much larger in ESM2M than in ESM2G, their average age is similar (86 yr in ESM2G and 81 yr in ESM2M). As noted above, explicit vertical diffusion is much higher in ESM2M than in ESM2G both within and just below the mixed layer. In the tropical thermocline, however, ESM2M diffusivities are generally lower than ESM2G diffusivities. While the causality behind this different behavior is not completely understood, we suspect that it is related to the propensity of ESM2G to maintain coherence in water mass properties compared to ESM2M.
A combined view of Southern Ocean zonal mean age (color), meridional density (σθ; black lines), and zonally integrated meridional overturning (green) in ESM2M (Fig. 14a) and ESM2G (Fig. 14b) highlights how differently these two models respond under conditional stability in deeply convective areas. Contours of density illustrate ESM2M’s stronger meridionally sloping isopycnals between 45° and 60°S consistent with deeper maximum mixed layer depths in ESM2M than ESM2G (Fig. 11). The overturning streamfunction (contours) illustrates the much deeper southward flow in ESM2M, resulting in a much younger upper water column (color) than that in ESM2G. South of the deep mixing region, westerly winds drive year-round Ekman upwelling supplying relatively dense waters that maintain upper water column stratification. The transport of those dense waters to north of 50°S (25 Sv in both models) supports deep wintertime convection with associated intense heat loss. Because the meridional temperature gradient is strongest in the deep convective region and decreases to the north, horizontal mixing in the meridional direction provides a significant local cooling. In ESM2M, this effect overcomes the weak vertical stratification in the deep wintertime convection region to invert the vertical density gradient and induce further convection year-round. Because horizontally adjacent layers only interact within the mixed layer in ESM2G, and cabeling during isopycnal mixing does not induce significant cooling, convection occurs only in the mixed layer during wintertime with southward flow restricted to a fairly shallow depth (Fig. 14b).
In ESM2M, at least three mechanisms drive this enhanced lateral mixing convection and deepening southward circulation of the Deacon Cell relative to ESM2G. First, lateral mixing otherwise oriented along locally referenced potential density surfaces defaults to horizontal, with diapycnal mixing beyond the threshold slope of . Second, the vertical transition between ESM2M’s surface-oriented submesoscale scheme (Fox-Kemper et al. 2011) and deeper-scale thickness mixing parameterization (Gent and McWilliams 1990) drives a secondary subsurface overturning circulation stimulating convection. Finally, high shear at the edge of the Antarctic Circumpolar Current induces significant advective exchange across the meridional temperature gradient. Given these complex interactions, we suggest that further development of parameterizations relating to convection is warranted. Potential development foci include improved representation of the mesoscale and submesoscale schemes (e.g., Boning et al. 2008; Farneti et al. 2010). Southern Ocean biases would also benefit from amelioration of atmospheric biases in radiative transfer (Fig. 4) that drive overly warm SST biases (Fig. 5) and enhance summer surface stratification in these models.
In summary, we have developed two new Earth System Models (ESM2M and ESM2G) that differ only in ocean formulation and compared them via a suite of factors influencing climate. Climate simulation quality was found to be generally very similar to that of GFDL’s earlier CM2.1. Demonstrating the importance of ocean configuration on climate fidelity, the two models straddle observed estimates for a number of climate indices, including ENSO variability being overly strong in ESM2M and weak in ESM2G, the volume of the thermocline being too deep in ESM2M and slightly too shallow in ESM2G, deep Pacific ventilation being relatively weak in ESM2M (6 Sv) and strong in ESM2G (16 Sv), and North Pacific thermocline ventilation being too strong in ESM2M. While each model has its relative strengths and weaknesses, our overall assessment of the models is that their quasi-equilibrium simulations both achieve sufficient fidelity to allow meaningful perturbation studies. Like Megann et al. (2010), we find that the insertion of an isopycnal model under the same atmosphere used by a z-coordinate model results in interior ocean changes consistent with reduced mixing. While some of these differences may result from numerical artifacts, the present runs cannot prove causality because of numerous explicit configuration and parameterization differences, as well as coupled climate feedbacks (e.g., sea ice and albedo). We demonstrate that, when equilibrated with the same radiatively tuned atmosphere, ESM2G develops a climate with similar overall fidelity to that of ESM2M, but with a shallower and less-ventilated thermocline, more vigorous bottom flows, weaker ENSO, more extensive sea ice, shallower mixing and mode water formation, and absence of long-term warming drift. The relative biases in circulation between the two models suggest that ESM2G might better represent climate changes relating to total heat content variability, given its lack of long-term drift; gyre circulation and ventilation in the North Pacific, tropical Atlantic, and Indian Oceans; and depth structure in the AMOC and abyssal flows. Alternatively, ESM2M might better represent climate changes relating to surface circulation given its superior SST, SSS, and SSH patterns; tropical Pacific circulation and variability; and Southern Ocean dynamics. In the companion paper (Part II), we describe the implications of these differences for simulation of the preindustrial carbon cycle. In future papers, we will apply these models to evaluate the impact of human activities on past and future climate changes.
Vaishali Naik and Dan Schwarzkopf processed the radiative forcing. Tony Rosati, Anand Gnanadesikan, Jorge Sarmiento, Joellen Russell, Robbie Toggweiler, and Stephanie Downes provided valuable insights on Southern Ocean circulation that assisted in our interpretation.