The authors describe carbon system formulation and simulation characteristics of two new global coupled carbon–climate Earth System Models (ESM), ESM2M and ESM2G. These models demonstrate good climate fidelity as described in part I of this study while incorporating explicit and consistent carbon dynamics. The two models differ almost exclusively in the physical ocean component; ESM2M uses the Modular Ocean Model version 4.1 with vertical pressure layers, whereas ESM2G uses generalized ocean layer dynamics with a bulk mixed layer and interior isopycnal layers. On land, both ESMs include a revised land model to simulate competitive vegetation distributions and functioning, including carbon cycling among vegetation, soil, and atmosphere. In the ocean, both models include new biogeochemical algorithms including phytoplankton functional group dynamics with flexible stoichiometry. Preindustrial simulations are spun up to give stable, realistic carbon cycle means and variability. Significant differences in simulation characteristics of these two models are described. Because of differences in oceanic ventilation rates, ESM2M has a stronger biological carbon pump but weaker northward implied atmospheric CO2 transport than ESM2G. The major advantages of ESM2G over ESM2M are improved representation of surface chlorophyll in the Atlantic and Indian Oceans and thermocline nutrients and oxygen in the North Pacific. Improved tree mortality parameters in ESM2G produced more realistic carbon accumulation in vegetation pools. The major advantages of ESM2M over ESM2G are reduced nutrient and oxygen biases in the southern and tropical oceans.
We describe carbon system formulation and simulation characteristics of two new global coupled climate–carbon Earth System Models (ESMs; Fung et al. 2000) developed at the National Oceanic and Atmospheric Administration (NOAA)/Geophysical Fluid Dynamics Laboratory (GFDL). Like global box models (e.g., Siegenthaler and Sarmiento 1993; Fig. 1a), ESMs represent major carbon reservoirs and fluxes (Figs. 1b,c). As ESMs are extensions of climate models, they are based on mechanistic geophysical understanding with geographically explicit atmosphere, ocean, land, and sea ice dynamics. To physical climate, ESMs add interactive carbon dynamics and associated chemistry and ecology to explore the earth’s system behavior at both equilibrium and in transient. ESMs resolve coupled climate–carbon responses to diverse anthropogenic perturbations such as fossil fuel emissions, agriculture and forestry, and aerosol chemistry within a single self-consistent system.
Carbon in both the land and ocean has been implicated as having potentially strong responses and feedbacks to anthropogenic forcing. In a climate model with collapsing North Atlantic overturning (Manabe and Stouffer 1993), Sarmiento and Le Quéré (1996) added simple ocean biogeochemistry to suggest ocean carbon uptake and biological feedbacks exhibit high sensitivity and uncertainty. Friedlingstein et al. (2003) compared two early ESMs to show CO2 feedbacks were highly sensitive to representation of both the Southern Ocean and land vegetation and soil responses. Friedlingstein et al. (2006) compared a suite of early ESMs to suggest land and ocean carbon responses each accounted for >100 PgC in emissions uncertainty at an atmospheric CO2 concentration of 700 ppm. Since then, much effort has gone into developing fully coupled ESMs with long-term stability in both climate and carbon (e.g., Doney et al. 2006). One fundamental development has been the simulation of climate without flux adjustment (e.g., Delworth et al. 2006; Reichler and Kim 2008). Considerable development has also taken place in both ocean biogeochemistry (e.g., Moore et al. 2004) and land carbon (e.g., Thornton et al. 2007) models. Initial comparison of second generation ESMs has shown large variability (e.g., Steinacher et al. 2010). More extensive comparison is expected as part of phase 5 of the Coupled Model Intercomparison Project (CMIP5; Taylor et al. 2011).
Our goal was to develop two ESMs with different ocean vertical coordinate and dynamical—physical core while adhering to the climate fidelity of the GFDL Climate Model version 2.1 (CM2.1; Delworth et al. 2006). The primary benefit of our approach is to explore the sensitivity of ocean carbon and heat uptake to anthropogenic forcing under fundamentally different frameworks describing ocean dynamics. While the z-coordinate framework holds many benefits, the extreme anisotropy of mixing and advection leads to spurious numerical mixing (Griffies et al. 2000; Ilicak et al. 2012). While the isopycnal framework takes advantage of the relative ease of motion between isopycnal layers, it loses resolution in the unstratified ocean. Comparison between these formulations thus allows assessment of the relative and absolute fidelity of both approaches in representing ocean physics, climate, and biogeochemistry.
Dunne et al. (2012a, hereafter Part I), described the physical climate simulations of both the z-coordinate (ESM2M) and isopycnal-coordinate (ESM2G) models, finding them of similar overall fidelity but having distinct differences. Among these, ESM2M warms with more thermocline ventilation, while ESM2G cools with more bottom water ventilation. While ESM2M has an overly vigorous El Niño–Southern Oscillation, ESM2G’s is weak.
We first describe the ESM carbon cycle components. We then discuss initialization to stable carbon states with 1860 radiative forcing and potential vegetation (i.e., no land use) over >1000 yr. Finally, we present preindustrial control integrations with CO2 interacting freely between the ocean–atmosphere–land system for another 600 yr.
2. Model description
The physical components for the ESMs presented here are described in Part I. Only their carbon system components are described below.
a. Ocean ecology and biogeochemistry
The current GFDL ocean biogeochemical and ecological component is Tracers of Ocean Phytoplankton with Allometric Zooplankton code version 2.0 (TOPAZ2). A full technical description of TOPAZ2 is available in the supplement. A technical description of TOPAZ version 0 is available in Dunne et al. (2010). TOPAZ version 1 was discussed in Henson et al. (2009). TOPAZ2 includes 30 tracers to describe the cycles of carbon, nitrogen, phosphorus, silicon, iron, oxygen, alkalinity, and lithogenic material as well as surface sediment calcite (Dunne et al. 2012b). TOPAZ2 considers three explicit phytoplankton groups (“small,” “large,” and diazotrophic) utilizing modified growth physiology after Geider et al. (1997) with Eppley (1972) temperature functionality and iron colimitation with luxury uptake after Sunda and Huntsman (1997). The small group represents mostly prokaryotic picoplankton and nanoplankton, whereas the large group represents diatoms, greens, and other large eukaryotes. Diazotrophs represent facultative nitrogen fixers with nitrogen fixation inhibited by nitrate (NO3), ammonia (NH4) and oxygen (O2). CO2:NO3:O2 stoichiometry is 106:16:150 after Anderson (1995). Phytoplankton N:P utilizes optimal allocation after Klausmeier et al. (2004).
Phytoplankton loss and production of sinking detritus utilize the size-based relationship of Dunne et al. (2005) with mineral-driven penetration of sinking detritus (Klaas and Archer 2002; Dunne et al. 2007). TOPAZ2 diagnoses plankton mineral formation of opal, calcite, and aragonite. TOPAZ2 includes seasonal time-scale dissolved organic material and heterotrophic biomass with fixed N:P and multiannual dissolved organic material with variable N:P. Gas exchange of O2 and CO2 follows Najjar and Orr (1998). Nitrification is inhibited by light after Ward et al. (1982). TOPAZ2 includes second-order iron scavenging with ligand kinetics, lithogenic particle scavenging, water column denitrification under suboxia, and sediment denitrification after Middelburg et al. (1996). In the absence of both NO3 and O2, a respiration deficit is accumulated as negative O2. TOPAZ2 includes external inputs of atmospheric nitrogen deposition (Horowitz et al. 2003); lithogenic dust and soluble iron (Fan et al. 2006); river nitrogen (Seitzinger et al. 2005); and river inputs of dissolved inorganic carbon, alkalinity, and lithogenic material set to balance Holocene burial of calcite and lithogenic material (Dunne et al. 2007).
b. Land vegetation and carbon cycle
The current GFDL land model (LM3.0) represents vegetation as five dynamically competing vegetation types: deciduous, coniferous, and tropical trees and warm and cold grasses. They combine three characteristics: C3 versus C4 leaf physiology; leaf longevity (i.e., temperate cold deciduous, tropical broadleaf, and coniferous evergreen); and allocation ratios among stems, roots, and leaves (i.e., trees versus grass). Tropical trees behave as evergreen or deciduous depending on the drought conditions. Vegetation height varies as a function of plant biomass across a continuum from shrubs through trees. All vegetation types have five carbon pools: leaves, fine roots, sapwood, heartwood, and labile. The sizes of all pools are modified daily through allocation rules based on carbon accumulated in the biophysics module governing exchanges of water and CO2 on fast time scales (30 min). LM3.0’s photosynthesis is based on Farquhar et al. (1980) and Collatz et al. (1991, 1992). The relationship between stomatal conductance and net photosynthesis is based on Leuning (1995).
LM3.0 simulates changes in vegetation carbon pools through phenology (e.g., leaf drop and emergence), natural mortality, and fire. Carbon loss from vegetation pools is deposited into two soil pools decomposing at seasonal and multiyear time scales with rates dependent on carbon amount, temperature, and water. Annual fire loss is proportional to fuel available and number of drought months. Since the ESM2M runs were conducted first, we were able to make use of them to ameliorate the ESM2M temperate biomass biases in the later ESM2G runs. To do so, ESM2G mortality was increased for evergreen coniferous trees from 0.015 to 0.0275 yr−1 and for temperate deciduous trees from 0.015 to 0.025 yr−1 while their branch-wood turnover rate was doubled.
One unique feature of LM3.0 is land-use heterogeneity in which each grid cell is described as a combination of four land-use categories (tiles): undisturbed lands (i.e., “primary” or “potential”), croplands, pastures, and lands either previously harvested or used in agriculture (i.e., “secondary,” not used in the present experiments). Each tile has its own carbon and physical state, above and below ground. In the absence of human disturbances, all grid cells on land represent potential vegetation. These simulations were performed for the potential vegetation state and be more aptly considered “prehuman” rather than “preindustrial” from a carbon perspective because of their lack of incorporation of land use history up to 1860. Nonetheless, we retain the preindustrial moniker for consistency with previous work.
c. Atmospheric CO2
Atmospheric CO2 was treated as an explicit prognostic tracer exchanging CO2 with land and ocean every 30 min and 2 h, respectively. For the initial phase of the long spinup (~1000 yr) described here, this atmospheric CO2 tracer was subjected to global, linear restoring to the 1860 reference value of 286 ppm dry volume mixing ratio throughout the atmosphere with a restoring time scale of 1 yr. This approach provides considerable advantages over prescribing a fixed surface CO2 concentration in enabling realistic diurnal and seasonal time-scale CO2 variability over land, which in turn influences plant growth and carbon inventory. Furthermore, it eliminated atmospheric CO2 drift in the later phase of the spinup for “restored” perturbation and “free” CO2 scenarios (i.e., where atmospheric CO2 exchanges without any restoring in both 1860 control and emissions perturbation scenarios). The radiation calculation uses the global average atmospheric CO2 concentration.
To obtain the 1860 initial conditions for the ESM2M and ESM2G control integrations, a method similar to Stouffer et al. (2004) was used (and see Part I). Biogeochemical tracers were initialized from World Ocean Atlas 2005 observations for nitrate, phosphate, silicate, and oxygen (Garcia et al. 2006a,b; Collier and Durack 2006) and the Global Ocean Data Analysis Project (Key et al. 2004) for alkalinity and preindustrial dissolved inorganic carbon (DIC). Initial sediment calcite was derived from bottom water conditions and fluxes after Dunne et al. (2012b). All 1860 radiative forcings were included simultaneously at spinup. The land model was initialized offline with CM2.1 climate forcing to obtain initial distribution of vegetation and soil carbon pools. After 500 yr, the last 50 yr was used for offline soil equilibration after Shevliakova et al. (2009) and the soil values reset. ESM2M and ESM2G were integrated over 1000 model years to achieve quasi equilibrium defined by net atmosphere–land and atmosphere–ocean fluxes being <0.1 PgC yr−1 on the multicentennial average. The 1860 control integrations and perturbation scenarios were then integrated another 600 yr.
a. Initial carbon cycle drift
While the ESM2M and ESM2G land carbon and sediment calcite components are spatially static and thus amenable to offline initialization, their ocean physical and biogeochemical components are spatially dynamic and require long, online equilibration to achieve our multicentennial 0.1 PgC yr−1 quasi-equilibrium goal. Characterization of the physical climate drifts are included in Part I. Initial ALK and DIC inventory differences arise from differences in total ocean volume (ESM2G has 0.5% more volume than ESM2M) with ESM2M’s relatively smooth bathymetry limited to 5500 m while ESM2G extends to 6000 m. Upon initialization with preindustrial DIC and present-day alkalinity (ALK) and nutrient distributions, both ESMs exhibit initial outgassing before ingassing CO2 over the first few centuries (Figs. 2a,b). ESM2G takes up about 10 Pg more carbon than ESM2M and comes into equilibrium more quickly.
The rapid uptake of CO2 in both these models over the first 300 yr and long-term DIC inventory response relates to a complex interplay of factors including temperature drift, accumulation of remineralized DIC in the tropics, loss of DIC in the deep Southern Ocean, and long-term equilibration of the calcite and nitrogen cycles. While these models include externally fixed concentrations of DIC and ALK in river fluxes to the ocean to match the long-term burial of calcite in surface sediments, they do not simulate river export of organic carbon or DIC from land. Because of the long-term equilibration time scales with respect to sediment calcite, the DIC inventory (Fig. 2c) does not follow only the integrated gas exchange flux (Fig. 2b), but it also follows the long-term changes in ALK inventory (Fig. 2d) associated with drifting bottom water saturation state. The ocean nitrogen inventory in both models (Fig. 2e) also undergoes a long-term loss under imbalanced nitrogen supply from N2 fixation, rivers, and atmospheric deposition versus loss from water column and sediment denitrification. As we describe below, this nitrogen loss is due to the common overexpression of tropical pelagic suboxia in this class of models (Najjar et al. 2007).
b. Surface ocean biogeochemistry
ESM2M and ESM2G share similar overall surface biogeochemical fidelity but also contrast in important ways (Fig. 3). While both models capture most regional variability in surface NO3 (r2 = 0.87 for ESM2M and r2 = 0.91 for ESM2G), ESM2M has a surplus of 1.6 μM compared to the observed average of 5.1 μM 70°S–70°N expressed in the equatorial Pacific upwelling region, the Subpolar Front near 40°S, and subpolar North Atlantic. ESM2G, in contrast, has a muted surplus in the Subpolar Front near 40°S and minimal global-scale bias (−0.062 μM). Both models underestimate surface NO3 in the subpolar North Pacific: ESM2G more severely.
ESM2M and ESM2G capture slightly more than half of log(Chl) spatial variability (r2 = 0.54 for ESM2M; r2 = 0.52 for ESM2G; Figs. 3d–f) as both models capture the broad transitions between oligotrophic subtropical gyres to mesotrophic, high nitrate regions, though neither model captures the high values in the north polar, coastal upwelling, and shelf regions. ESM2M better represents the tropical west Pacific, while ESM2G better represents the tropical Atlantic and Indian Oceans and subpolar regions of all oceans. We attribute these differences in surface NO3 and Chl fidelity to differences in mixed layer dynamics (Part I) with ESM2G’s capacity to represent shallow mixed layers allows phytoplankton maintaining high growth rates with lower chlorophyll than ESM2M. Strong shallow biases in minimum Southern Ocean mixed layer depths in both models (Part I) lead to locally low biases in surface NO3 (ESM2G more severely than ESM2M) and high bias in Chl (ESM2M more severely than ESM2G).
Both models capture slightly more than half of spatial ΔpCO2 variability (r2 = 0.56 for ESM2M and r2 = 0.55 for ESM2G; Figs. 3g–i). The primary bias patterns include over expression of high ΔpCO2 in the northerly and westward expressions of equatorial Pacific upwelling, more so for ESM2G than ESM2M. The ΔpCO2 in both models is too low in the northwest Pacific near 40°N, particularly in ESM2G, and too high in the Southern Ocean, particularly east of Australia and Argentina. While the global patterns in PO4 (Figs. 3j–l) bear much similarity with observations (r2 = 0.87 for both models), ESM2M has a significant low bias in PO4 (−0.087 μM) whereas ESM2G does not (0.0063 μM), opposite from their relative NO3 biases.
c. Interior ocean biogeochemistry
Comparisons of biogeochemical distributions in the ocean interior highlight the limits of these models’ abilities to represent large-scale ventilation patterns at 500 m (Fig. 4). While ESM2M reproduces 72% the overall interior O2 regional variability, it is overly ventilated in both the North Pacific and south of 50°S and lacks ventilation of the eastern equatorial Pacific and Atlantic (Fig. 4b). While ESM2G represents the lack of oxygen ventilation in the North Pacific much better than ESM2M, ESM2G reproduces slightly less (69%) spatial variation overall as it exacerbates the ESM2M bias in the tropical and southern subtropical oceans and Southern Ocean north of 60°S. In addition, ESM2G is overly ventilated south of 60°S. These patterns are mirrored in the PO4 patterns at 500 m (Figs. 4d–f). Overall, ESM2M has minimal average bias in O2 and PO4 at 500 m (−5.2 and 0.005 μM, respectively), while ESM2G has a much larger negative bias in O2 (−32 μm) and positive bias in PO4 (0.31 μM). Given their generally good surface nutrient representation (Figs. 3j–l), we infer these O2 biases as due to tropical underventilation leading to excess accumulation of remineralization. Comparison of NO3 at 500 m (Figs. 4g–i) highlights the severe implications of these ventilation biases for the nitrogen cycle as overexpression of suboxia leads to excess denitrification and subsequent NO3 deficiency in the eastern tropical Pacific of ESM2M and even more so in ESM2G. This leads to both significant global NO3 drift (Fig. 2e) and to ESM2G’s relative surplus in surface PO4 (Fig. 3l).
Model fidelity for both DIC and excess alkalinity (ALKEXCESS = ALK − DIC; an approximation of the ocean buffering potential with respect to addition of CO2 and other acid) is demonstrated in Fig. 5 through regional comparison of observed and modeled profiles. In the North Atlantic (Fig. 5a), distributions of DIC and ALKEXCESS have intermediate values throughout the water column, with ESM2G DIC values slightly elevated above 3000 m, and ESM2M DIC values slightly depressed below 700 m. ALKEXCESS in both models, however, agrees well with observations in the upper water column but is slightly elevated below 700 m. In all tropical areas (Fig. 5b), both models have elevated DIC at the base of the thermocline down to 2000 m consistent with their elevated phosphate and over expression of suboxia in the east Pacific (Fig. 4), with ESM2G more strongly and more into the thermocline than ESM2M. The ALKEXCESS is correspondingly depressed, though the observed structure is well preserved. In the North Pacific (Fig. 5c), ESM2G reproduces the DIC profile down to 700 m but has strongly depressed values below, while ESM2M underestimates DIC down to 2000 m with both models overestimating ALKEXCESS through much of the water column. The differences in these deep patterns is due to the much higher Antarctic Bottom Water transport and lower ideal age distribution in ESM2G than ESM2M discussed in Part I. In the Southern Ocean (Fig. 5d), DIC patterns similarly reflect these models’ differing ventilation, with ESM2M reproducing the DIC profile extremely well down to 1200 m but having depressed values below, while ESM2G has elevated values in the 300–1400-m range and depressed values below. Both models exhibit slightly depressed ALKEXCESS in this region above 1200 m and elevated values below.
d. Marine ecosystem functioning
The similar atmospheric coupling and identical biogeochemical algorithms but differing thermocline ventilation in these two models (Part I) results in roughly scaled marine ecological functioning. Global flux integrals of relevant carbon cycle parameters for both models are given in Fig. 1. Both models broadly center primary production on the equator with midlatitude secondary maxima (Figs. 6a,b) due to the combination of nutrient supply and equator-to-pole gradient in light and temperature. As a result of its more vigorous thermocline ventilation, ESM2M has higher total primary production than ESM2G in most ocean areas (Fig. 6c). Patterns in large phytoplankton production (Figs. 6d–f) are similar. Exceptions to this general pattern include ESM2G’s 1) high productivity regions more locally restricted to the site of ventilation, 2) Southern Ocean high productivity ribbons shifted slightly southward, 3) smaller western equatorial Pacific warm pool (Part I), and 4) California Current upwelling drawing from a shallower nutricline (Fig. 5c). Nitrogen fixation patterns (Figs. 6g–i) highlight the differing N:P stoichiometry in these models (Fig. 4) with ESM2G (Fig. 6g) having much higher nitrogen fixation in the North Pacific than ESM2M (Fig. 6h) and 19% higher nitrogen fixation globally (ESM2G: 173 TgN yr−1; ESM2M: 146 TgN yr−1).
Comparison of zonally integrated metrics of the organic carbon pump (Fig. 7a) illustrates both similarities and differences between ESM2M (solid) and ESM2G (dashed). In both models, sinking particle export (black) follows large phytoplankton production (green) and NO3-based production (light blue). Net community production (total primary production minus respiration; blue) and net phytoplankton production (total primary production minus grazing; red) exhibit more focused maxima in the high productivity regions. While the patterns of the two models are similar because of their similar atmospheric forcing, the ESM2M organic carbon pump is generally stronger than the ESM2G pump because of ESM2M’s more vigorous thermocline ventilation (Part I). The similar biogeochemical response to atmospheric forcing is highlighted in net phytoplankton production (red) as upwelling-driven biomass accumulation and lateral transport near the equator and convection-driven vertical transport at midlatitudes. South of 40°S, while NO3-based production (light blue) is nearly identical in these models, net community production (blue) and particle export (black) are considerably lower in ESM2G than ESM2M. This differing behavior in the Southern Ocean is due to their differences in depth-resolution of light-inhibited nitrification. While ESM2M always has 10-m euphotic zone resolution, ESM2G changes euphotic zone resolution based on stratification. Particularly in the Southern Ocean, ESM2G layers directly below the mixed layer can thicken enough to relieve light inhibition of nitrification and drive down ammonia.
The mineral pumps (Fig. 7b) each highlight different controls dominating in different regions. The silicon export (green) dominates in high productivity areas, the calcite pump (blue) dominates in the high temperature and saturation state tropics, the aragonite pump (black) dominates in high productivity tropical areas, and the lithogenic pump (red) is focused in the Northern Hemisphere. The similarity of the lithogenic pump between the two models is because each was given the same external supply from atmospheric deposition and river runoff. While the silicon pump is strongest in terms of surface export, the assumed mineral protection efficiency is higher for calcite (0.070), aragonite (0.070), and lithogenic (0.065) mineral relative to silicon (0.026; Klaas and Archer 2002), and calcite and lithogenic mineral penetrate the deepest through the water column.
e. Land vegetation and carbon
Both ESMs dynamically simulate vegetation characteristics as a function of climate and atmospheric CO2 concentration. To assess ESM2M and ESM2G fidelity with respect to land dynamics, we compare three bulk characteristics critical to land–climate interactions: the distribution of vegetation types, total carbon storage, and the seasonal extremes of partitioning among different pools into leaves. We have combined information on the ESMs simulated vegetation type distribution with other vegetation characteristics (e.g., biomass, seasonality) to classify each of the land model grid cells into one of the 10 aggregated Olson ecotypes. Overall ecotypes distributions (Olson et al. 1985) in both ESM2M (Fig. 8b) and ESM2G (Fig. 8c) compare well to observational estimates (Gibbs 2006; Fig. 8a) in representing the general distinction between warm grassland, tropical forest, cold grassland, hot sandy desert, deciduous forest, savanna, and coniferous forest, while tundra is underrepresented in both models and permanently iced areas such as Greenland and Antarctica are specified. The models both over extend the areal coverage of tropical forest in South America and Africa at the expense of savanna and overextend the areal coverage of hot sandy desert in Africa and South Asia at the expense of warm grassland. While much of this bias is related to precipitation biases, it may also be related to the simple treatment of postfire succession in LM3.0, a topic of ongoing work (Magi et al. 2012). The models also exaggerate the westward extent of deciduous forest in North America and East Asia. Also in comparison to observational estimates (Gibbs 2006; Fig. 8d), both models (Figs. 8e,f) overestimate both tropical and deciduous forest biomass, while because of the tuning differences coniferous biomass is better represented in ESM2G than ESM2M.
As an essential ingredient for hydrological, radiative, and carbon exchanges in the ESMs, we compare the ESM partitioning of biomass to leaves against the normalized difference vegetation index (NDVI) to evaluate whether the models make the correct amount of leaves at their seasonal extremes. Comparisons of NDVI against the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite (Figs. 8g–l) illustrates strong agreement as expected from the spectral calibration used in LM3.0 (Part I). The climatological maximum comparison (Figs. 8g–i) illustrates the underrepresentation of permanent deserts and semiarid regions are underrepresented in these models. This is primarily because the low precipitation extremes are not as extreme as in the satellite estimate. The climatological minimum comparison (Figs. 8j–l) illustrates the excellent observational agreement achieved with these models in both the Northern Hemisphere snow line extent and transition between coniferous and deciduous forests. However, the models appear to have overly high minimum NDVI in tropical forest indicating a possible bias in lack of vegetation sensitivity to climate in these regions.
As both ESMs share the same atmosphere and land model except for tree mortality parameters, the two models share nearly identical features. Three main issues complicate evaluation of their vegetation characteristics: 1) most land vegetation and soil datasets represent present-day conditions and include a substantial anthropogenic footprint; 2) most global datasets are derived from highly uncertain, extrapolated point measurements; and 3) as with many categorical variables modeled vegetation types cannot be exactly corresponded with the different vegetation classifications available. We choose Olson et al. (1985) for such comparison because it reflects carbon content of different ecotypes as well as degree of human management.
f. Atmospheric CO2 cycling
The land and ocean exhibit fundamentally different atmospheric CO2 flux characters in these ESMs (Fig. 9). In the long-term mean (Figs. 9a,b) only ocean fluxes are apparent as land biomass is near equilibrium on the centennial scale. While ESM2M and ESM2G share similar overall CO2 flux patterns, ESM2M has its strongest ingassing north of Norway whereas ESM2G’s occurs in the Sea of Japan and Kuroshio Extension. ESM2M also shows strong ingassing along the coast of Antarctica while ESM2G does not. Focusing on the monthly CO2 flux variability (Figs. 9c,d) nearly completely obscures the ocean due to the much higher land signal (note logarithmic scale), with both models showing nearly identical patterns of intensely high fluxes in the tropical South American and African forests and savannas and secondary zones of flux in tropical Asia, eastern North America, and Western Europe. Over the ocean, minimum variability is exhibited near the equator, and lower fluxes in ESM2G than ESM2M in the Southern Ocean. As expected because of the primary driver of land variability fluxes by precipitation, the partitioning of variability between seasonal versus interannual modes of flux (Figs. 9e,f) shows a strong dominance of seasonal variability except for the equatorial Pacific and Indian Oceans as a consequence of the El Niño–Southern Oscillation (ENSO), especially for ESM2M where the ENSO signal extends to the central Pacific, Indian, and Atlantic Oceans. Additionally, ESM2G shows a strong interannual signal in the Southern Ocean while interannual signals in ESM2M for this region are moderate. ESM2G also shows a more interannual component in the subpolar North Pacific.
The local factors driving this variability are explored via the monthly air to land–sea CO2 flux (FCO2) correlation with monthly surface temperature (TSurf) in Figs. 9g,h. Strong positive FCO2–TSurf correlations are seen throughout the land but most particularly in South America; Australia; and southeast areas of Asia, Africa, and the United States as a function of the seasonal cycle in irradiance. In the ocean, these positive correlations are restricted somewhat to the subtropics and Ross and Weddell Sea regions where seasonal heating drive solubility-induced fluxes, whereas negative correlations are seen in the tropical and Southern Ocean upwelling areas where remineralized CO2 drives outgassing during upwelling of colder subsurface waters. This latter response is seen more prominently in ESM2G than ESM2M because of its stronger pycnocline.
Comparison of time-averaged, zonally integrated, preindustrial air–sea fluxes (Fig. 10a) of heat downward [PW (° latitude)−1 into the ocean; dashed] and CO2 upward [PgC (° latitude)−1 yr−1; solid] and implied ocean transport (Fig. 10b) of heat northward (PW; dashed) and CO2 southward (PgC yr−1; solid) serves to highlight both the differing dynamics controlling heat versus CO2 fluxes as well as the differences between ESM2M and ESM2G. Relative to heat, CO2 flux is strong south of 45°S and meridionally broad in the tropics. As the Southern Ocean release is roughly compensated by uptake south of 15°S, the tropical release is largely balanced northern midlatitude uptake. This divergence of CO2 from heat is highlighted in the implied atmospheric CO2 transport being northward at all latitudes north of 60°S, whereas the implied ocean heat transport is southward south of the equator (Fig. 10b). While the two models have very similar fluxes and transport for heat (ESM2G’s northward heat transport being only slightly stronger), ESM2G has a 28% stronger maximum CO2 transport (0.97 PgC yr−1) than ESM2M (0.76 PgC yr−1) resulting from their differing extents of thermocline ventilation discussed in Part I. In ESM2M, CO2 locally remineralized in the subsurface tends to ventilate back to the atmosphere during winter convection at northern midlatitudes while ESM2G maintains a stronger biological pump and more uptake at subpolar latitudes. Because ESM2M and ESM2G represent only the component of river fluxes necessary to match long-term sediment burial of CaCO3 (i.e., fixed Alk = 0.48 mmol kg−1 and DIC = 0.24 mmol kg−1 in river water), other terms driving total estimated river DIC fluxes of 0.8 PgC yr−1 (Siegenthaler and Sarmiento 1993; Fig. 1) do not appear in model air–sea CO2 fluxes as they would in the real ocean.
g. Atmospheric CO2 variability
Much of our ability to infer present-day carbon fluxes in land and ocean depends on our ability to interpret patterns in atmospheric CO2 variability and the changes in this variability from the preindustrial conditions, which are unfortunately not well constrained from ice cores because of problems in intercalibration (Barnola et al. 1995). It is thus important to be able to constrain patterns and trends in atmospheric CO2 variability through models to provide the preindustrial reference. The preindustrial atmospheric signatures of regional land and ocean CO2 fluxes are presented in Fig. 11. At the regional scale, surface atmospheric CO2 anomalies from the global mean (Figs. 11a,b) are driven by the dominance of land variability over ocean variability (Fig. 9) and exhibit broad maxima centered on the productive land areas spreading over the marine areas, particularly in the tropics. These features are a consequence of the well-characterized atmospheric CO2 rectification effect (Denning et al. 1995) such that, as thick daytime boundary layers are capped to thin nighttime boundary layers when land undergoes net respiration, CO2 increases dramatically. At the global scale, these anomalies also demonstrate a general marine gradient of decreasing CO2 from south to north as a consequence of the general southward transport of CO2 in the ocean (Fig. 10). Temporal changes in global mean surface CO2 anomalies are highly correlated with temporal changes at regional scales (Figs. 11a,b), except for directly above the regions of high land productivity. Because of this globally synoptic behavior, local CO2 variability in such distant places as Antarctica (black in Figs. 11e,f) and Greenland (blue in Figs. 11e,f) is largely coincident (r2 = 0.81 for ESM2M and r2 = 0.84 in ESM2G) with concentrations at Antarctica approximately 1.2 ppm higher than at Greenland in ESM2M and 1.3 ppm higher in ESM2G.
We compare the seasonal surface CO2 cycle for select locations along the latitudinal gradient in ESM2M (red) and ESM2G (green) with observations from the NOAA/Global Monitoring Division (black) in Fig. 11. At high northern latitudes (Barrow, Alaska; Fig. 11a), the seasonal cycle is strongly damped and a month or two early. In middle northern latitudes (Mauna Loa, Hawaii; Fig. 11b), the spring maximum is well represented but the fall minimum is underestimated. In middle southern latitudes (Samoa; Fig. 11b), the models exhibit a strong winter maximum and a summer minimum is not observed and extends all the way to the South Pole (Fig. 11d).
Broad comparison of seasonal surface CO2 amplitude across the available suite of marine stations (Fig. 12) illustrates the scope of the damped meridional gradient in these models. This damped gradient has long been recognized as a deficit in models of the atmospheric CO2 cycle (Fung et al. 1987). ESM2G suffers slightly more from this bias than ESM2M due to ESM2G’s stronger seasonal cycle in the tropics (Fig. 7 of Part I). The causality behind this damped meridional gradient is threefold: 1) underpredicted strength and timing of the northern boreal seasonal cycle (e.g., Fig. 11a), 2) overpredicted tropical seasonal cycle, and 3) lack of fossil fuels and biomass burning terms in these preindustrial model configurations. Randerson et al. (1997) estimated the anthropogenic terms as 6%–10% of the seasonal cycle at high northern latitudes and 16%–33% in the tropics. In northern boreal regions, our models initiate CO2 drawdown approximately 1 month earlier than observations and return CO2 to higher, with winter values approximately 2 months early. While both models underestimate maximum NDVI in Siberian forests (Figs. 8j–l), the net primary productivity cycle matches satellite-based estimates (Randerson et al. 1997). Net ecosystem production, however, is damped in ESM2M and ESM2G compared to Randerson et al. (1997), pointing to suspect a problem in the seasonality of respiration. The tropics exhibit common Amazonian biases in overall low precipitation, relative intensification of the dry season, lack of intense cloud cover in the wet season, and related soil hydrology deficiencies intensifying water stress in the dry season (Fig. 6 of Part I). The modeled Amazonian seasonal CO2 flux cycle is thus reversed and synchronous with the African and Oceania rain forests, giving strong CO2 seasonality extending to the South Pole with net ecosystem production driven by precipitation.
To describe atmospheric CO2 variability time scales, we show as insets in Fig. 13 the temporal variability spectrum at Barrow, Mauna Loa, and the South Pole (detrended; black) compared to ESM2M (red) and ESM2G (green). While at high northern latitudes both models underestimate the seasonal cycle as in Fig. 12a, they have approximately correct interannual variability. While at middle northern latitudes both models overestimate interannual variability, they underestimate the seasonal cycle. In the Southern Hemisphere, both models overestimate variability at both annual and interannual time scales, with ESM2M exhibiting more bias on interannual time scales and ESM2G exhibiting more bias on annual time scales.
h. Overall land and ocean biosphere representation in GFDL’s ESM2 class of models
Overall, the models demonstrate the ability to represent a wide range of coupling between climate and the carbon cycle. Having described all of the various components of these models, we now turn back to the comparison of carbon inventory and flux estimates in Fig. 1. Both models demonstrate excellent agreement with previous estimates in terms of the partitioning of carbon between the various reservoirs, though the land vegetation is quite a bit higher and the ocean biota quite a bit lower than those in the observational synthesis and box model analysis of Siegenthaler and Sarmiento (1993), which is used as a quasi-consensus estimate used in the Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report (Randall et al. 2007). These estimates remain highly uncertain, however. While the Sabine et al. (2004) land gross primary production estimate (120 PgC yr−1) given in Fig. 1 is higher than the original Siegenthaler and Sarmiento (1993) estimate (100 PgC a−1) and similar to the Beer et al. (2010) estimate of 123 PgC a−1, it is well below both the latest estimate of 150–175 PgC yr−1 from the recent observational estimates (Welp et al. 2011) and range of 127–166 PgC yr−1 from the recent atmospheric CO2 data assimilation efforts (Koffi et al. 2012). Both ESM2M (141 PgC yr−1) and ESM2G (143 PgC yr−1) estimates fall between these more recent studies. The global biomass in both models (1044 PgC in ESM2M and 838 PgC in ESM2G) compares well to preindustrial estimates of 924–1080 PgC (Bazilevich and Rodin 1971; Adams et al. 1990). The Siegenthaler and Sarmiento (1993) biomass incorporates significant historical land use and is expected to be lower than preindustrial estimates. The global soil carbon estimates from ESM2M (1339 PgC) and ESM2G (1481 PgC) are within the range of the observational estimates for soil carbon (top 1 m) 1220–1502 PgC (Sombroek et al. 1993; Jobbagy and Jackson 2000).
Ocean primary productivity is higher than the Siegenthaler and Sarmiento (1993) value (50 PgC yr−1) in both ESM2M (82 PgC yr−1) and ESM2G (68 PgC yr−1) and at the upper end of the range of global productivity estimates (36–78 PgC yr−1) given in the comparative analysis of Carr et al. (2006) and much higher than the model range given in Steinacher et al. (2010) (24–49 PgC yr−1). However, ESM2M (8.0 PgC yr−1) and ESM2G (5.4 PgC yr−1) are within the range of modeled particulate organic carbon export (5.0–9.1 PgC yr−1) in Steinacher et al. (2010), indicating that the primary production differences are largely a function of differing recycling efficiency between these studies. We have split the ESM2M and ESM2G carbon export across 100 m in Fig. 1 into three terms: dissolved transport, which is much smaller than in Siegenthaler and Sarmiento (1993); sinking particulate organic carbon export (bold arrow), which is higher than Siegenthaler and Sarmiento for ESM2M but of similar magnitude in ESM2G; and combined organic particle (phytoplankton, zooplankton, and bacteria) transport and aragonite and calcite sinking. Combined, these export fluxes give an inferred ventilation of the upper 100 m in ESM2M (98 PgC yr−1) in excellent agreement with Siegenthaler and Sarmiento (1993) (100 PgC yr−1) while that in ESM2G (80 PgC yr−1) is moderately lower. Since these ESMs do not include an explicit representation of organic carbon cycling in rivers, our net atmospheric flux into land is zero, and the river runoff flux of DIC and ALK is specified with a 1:2 ratio to balance calcite burial at a steady state. This steady state is yet unachieved in our control runs as ESM2M with 642 PgC in surface sediment calcite is gaining 0.05 PgC yr−1 and ESM2G with 469 PgC in surface sediment calcite is gaining 0.12 PgC yr−1. Assessment of the sensitivity of the sediment calcite inventory to ocean acidification will thus require careful attention to the behavior of the control runs for comparison.
Representation of the carbon cycle in GFDL’s climate models has exposed a wide range of sensitivities beyond those normally focused upon in analysis of climate. As discussed above, many of the biases relative to observationally based estimates are common to both models. On land, the models overestimate the seasonal flux of CO2 to the atmosphere in the tropics and underestimate it in boreal high latitudes (e.g., Fig. 12). In the tropics, we attribute this bias to a combination of the physical model’s inability to represent the actual radiative and hydrological cycle in the Amazon (not shown) and the enhanced northward extent of seasonal productivity in the Sahel (Figs. 8g–i). In the ocean, both models overestimate surface chlorophyll (Figs. 3d–f) and the extent of interior suboxia (Figs. 4a–c). The relative equatorward restriction of the subtropical gyre oligotrophic region is a cause of concern for the representation of ecological biomes and overestimation of midlatitude nutrient and oxygen ventilation, which are both better represented in ESM2G than ESM2M. While not directly connected to the surface biogeochemistry, the oxygen bias has severe implications for the ability to represent hypoxia and denitrification in the tropics (Fig. 5b) and to achieve a steady-state budget for global alkalinity (Fig. 2d) and nitrogen (Figs. 2e, 4g–i). Efforts to address these biases are ongoing. Overall, however, the models’ very good representation of surface nutrients and ΔpCO2 (Fig. 3) and interior excess alkalinity (dashed lines in Fig. 5), in combination with the general agreement with physical metrics of transport and ventilation discussed in Part I, give us confidence both models should have good fidelity in representing ocean carbon uptake and acidification sensitivity under anthropogenic emissions.
The primary biogeochemical differences between ESM2M and ESM2G relate to their differing physical representations of thermocline ventilation and the magnitude of Antarctic Bottom Water (AABW) penetration into the deep Pacific as discussed in Part I. In ESM2G, lower thermocline ventilation leads to a lower productivity (Figs. 1, 6), higher Southern Ocean and tropical thermocline DIC (Figs. 5b,d), and larger northward atmospheric CO2 transport (Fig. 10b) than in ESM2M. In the North Pacific, higher AABW ventilation of the deep North Pacific in ESM2G leads to deficit of DIC relative to ESM2M and observations in deep waters, while higher thermocline ventilation in ESM2M leads to deficit of DIC relative to ESM2G and observations in the upper water column (Fig. 5c).
In summary, we find ESM2M and ESM2G of similar overall fidelity, both capable of representing the major carbon inventories and fluxes in a prognostic earth system context. The models are capable of simulating realistic levels of biogeochemical function (Fig. 1) with minimal drift (Fig. 2). Regionally, the models are capable of representing factors controlling productivity and net CO2 exchanges between the atmosphere and both land and ocean (Figs. 9, 10), though regional biases in variability do exist. The major advantages of ESM2G over ESM2M are improved representation of the subtropical, oligotrophic gyre structure in general and specifically lack of nutrient oxygen ventilation in the northeast Pacific. The major advantages of ESM2M over ESM2G are in reduced tropical nutrient trapping and suboxia biases: also due to the enhanced thermocline ventilation in ESM2M. While some ocean interior and sediment calcite drifts continue through our control runs, the surface fluxes and interior budgets for both models are balanced to within a 0.1 PgC yr−1 tolerance in order to provide a strong signal to noise ratio under conditions of anthropogenic CO2 emissions with IPCC Fourth Assessment Report consensus estimate for ocean CO2 uptake in the 1990s at 2.2 ±0.4 PgC yr−1 (Table 7.1 of Denman et al. 2007). We expect the model differences described here will drive much of the ocean carbon uptake and acidification sensitivity under anthropogenic CO2 emissions between ESM2M and ESM2G. Such simulations and analyses are ongoing as part of GFDL’s contribution to the IPCC Fifth Assessment Report.
Discussions with Inez Fung, Scott Doney and Chris Jones helped inform the experimental design of these simulations. Charles Stock, Robbie Toggweiler, Chris Jones, and two anonymous reviewers provided extremely helpful reviews of the manuscript.