The distribution of radiocarbon (14C) in the ocean and atmosphere has fluctuated on time scales ranging from seasons to millennia. It is thought that these fluctuations partly reflect variability in the climate system, offering a rich potential source of information to help understand mechanisms of past climate change. Here, a long simulation with a new, coupled model is used to explore the mechanisms that redistribute 14C within the earth system on interannual to centennial time scales. The model, the Geophysical Fluid Dynamics Laboratory Climate Model version 2 (GFDL CM2) with Modular Ocean Model version 4p1(MOM4p1) at coarse-resolution (CM2Mc), is a lower-resolution version of the Geophysical Fluid Dynamics Laboratory’s CM2M model, uses no flux adjustments, and is run here with a simple prognostic ocean biogeochemistry model including 14C. The atmospheric 14C and radiative boundary conditions are held constant so that the oceanic distribution of 14C is only a function of internal climate variability. The simulation displays previously described relationships between tropical sea surface 14C and the model equivalents of the El Niño–Southern Oscillation and Indonesian Throughflow. Sea surface 14C variability also arises from fluctuations in the circulations of the subarctic Pacific and Southern Ocean, including North Pacific decadal variability and episodic ventilation events in the Weddell Sea that are reminiscent of the Weddell Polynya of 1974–76. Interannual variability in the air–sea balance of 14C is dominated by exchange within the belt of intense “Southern Westerly” winds, rather than at the convective locations where the surface 14C is most variable. Despite significant interannual variability, the simulated impact on air–sea exchange is an order of magnitude smaller than the recorded atmospheric 14C variability of the past millennium. This result partly reflects the importance of variability in the production rate of 14C in determining atmospheric 14C but may also reflect an underestimate of natural climate variability, particularly in the Southern Westerly winds.
One of the great challenges in understanding the ocean’s role in climate variability is the scarcity of direct measurements, particularly prior to the last few decades. Such data limitations seriously hinder progress in understanding modes of interannual variability. These modes include the El Niño–Southern Oscillation (ENSO), for which continuous observations of many centuries are required for accurate characterization (Wittenberg 2009; Vecchi and Wittenberg 2010), and changes in the Southern Ocean overturning and Southern Annular Mode (SAM), which vary on time scales of multiple decades (Visbeck 2009). Proxy data can help to fill this gap, but many proxy measurements are confounded by spurious influences and are difficult to interpret.
Radiocarbon is one such proxy, directly recorded with high fidelity by marine microplankton providing decadal–centennial resolution (Hughen et al. 1998) and in both corals and tree rings at annual or even subannual resolution. Thousands of individual tree-ring records have been spliced together to reconstruct the radiocarbon history of the troposphere in both Northern and Southern Hemispheres (Fig. 1). Meanwhile, long coral records of radiocarbon have been recovered from the tropics (Druffel 1997; Guilderson et al. 1998; Druffel et al. 2004; Fallon and Guilderson 2008; Guilderson et al. 2009), revealing interannual variability in local seawater radiocarbon well in excess of the analytical error associated with coral measurements (~3‰ measurement error, Guilderson et al. 2009). Looking even further back in time, large changes in the radiocarbon activities of the ocean and atmosphere accompanied the warming that occurred at the end of the last ice age (Shackleton et al. 1988; Hughen et al. 1998; Galbraith et al. 2007; Marchitto et al. 2007; Skinner et al. 2010). Unfortunately, the manifold influences on radiocarbon activities can confuse the translation of these geochemical archives to useful climate information.
Radiocarbon is produced in the stratosphere and upper troposphere by cosmic rays, at a rate dependent on the strengths of the geomagnetic field and the solar wind (Masarik and Beer 1999) and decays with a 5730-yr half-life. Radiocarbon is quite well mixed in the atmosphere, with a small interhemispheric gradient (Fig. 1). A small portion of radiocarbon is lost through decay in the terrestrial biosphere, while the majority dissolves in the ocean, at a rate dependent on the ocean state (Siegenthaler et al. 1980). The fact that 14C production takes place exclusively in the atmosphere, coupled with the long air–sea equilibration time scale of carbon, caused the ocean to be everywhere undersaturated with 14C under preindustrial conditions, but to varying degrees. The slow pace of deep ocean circulation allows radiogenic decay to significantly deplete the 14C activity of abyssal waters so that strong undersaturation of 14C is a clear tracer of these “old” waters where they come to the surface (Toggweiler et al. 1991; Toggweiler and Samuels 1993; Rodgers et al. 2000; Gnanadesikan et al. 2004).
The complexity of interpreting radiocarbon records arises from the fact that both its natural rate of production and its distribution within the earth system could have changed simultaneously (e.g., Broecker and Barker 2007). Variability in the relatively small atmospheric 14C reservoir over the last millennium (Fig. 1), including “Suess wiggles” with a characteristic time scale of ~200 yr (Suess 1968), are frequently argued to reflect changes in the production rate of radiocarbon (Siegenthaler et al. 1980; Stuiver and Quay 1980; Masarik and Beer 1999; Muscheler et al. 2007). At the same time, interannual modes of climate variability are thought to drive redistribution of radiocarbon within the ocean (Rodgers et al. 1997; Rodgers et al. 2004; Druffel et al. 2007), with some impact on air–sea exchange. Meanwhile, sustained changes in ocean circulation, such as those thought to have occurred during the end of the last ice age, are argued to have caused very large changes in the redistribution of radiocarbon between the ocean and atmosphere (Siegenthaler et al. 1980; Hughen et al. 1998; Broecker and Barker 2007). These processes have been variously explored with box models (Siegenthaler et al. 1980), ocean GCMs (Toggweiler et al. 1989; Toggweiler et al. 1991; Rodgers et al. 1997), and ocean models coupled to simplified atmospheres (Meissner et al. 2003; Schmittner 2003). Through such works, radiocarbon has proven a useful diagnostic of the mean ocean circulation (Duffy et al. 1997; Guilderson et al. 2000; Gnanadesikan et al. 2004; Matsumoto et al. 2004; Grumet et al. 2005). However, the relationship between coupled ocean–atmosphere variability and the global distribution of natural radiocarbon on centennial time scales has received little attention.
In this paper we take a broad view of how oceanic radiocarbon is linked to natural, unforced variability in a global coupled climate model with no flux adjustments, but with a fixed atmospheric radiocarbon activity. The model, climate model version 2 with the Modular Ocean Model version 4p1 at coarse resolution (CM2Mc), includes more complex oceanic and atmospheric physics than those used in prior coupled model studies with radiocarbon (Marchal et al. 1999; Meissner et al. 2003; Schmittner 2003; Stocker and Wright 1996). CM2Mc is a new 3° configuration of the Geophysical Fluid Dynamics Laboratory (GFDL) coupled atmosphere–ocean GCM, using code similar to that of CM2.1 (Delworth et al. 2006) but with new parameterizations for subgrid-scale ocean physics and a new ocean biogeochemistry module. The coarse resolution enables us to simulate hundreds of years with many different tracers using reasonable computational resources, while maintaining a relatively realistic solution with robust modes of unforced ocean–atmosphere variability on decadal time scales. Section 2 describes the radiocarbon model and the mean state simulation produced under preindustrial radiative forcing. Section 3 examines the relationship between simulated modes of climate variability and oceanic radiocarbon at the sea surface. Section 4 discusses how the simulated variability of ocean circulation would modulate the air–sea balance of radiocarbon, applying the disequilibrium flux concept of Rodgers et al. (2011). Section 5 concludes the paper. Two appendices detail the model formulation and the control simulation.
2. Radiocarbon simulation
Radiocarbon is simulated as part of the idealized ocean biogeochemistry (iBGC) module, integrated simultaneously with the ocean component of CM2Mc (see appendix A for a complete model description). Dissolved inorganic carbon (DIC), that is, the sum of dissolved CO2, H2CO3, , is carried as a prognostic tracer. Air–sea gas exchange follows the Ocean Carbon Model Intercomparison Project (OCMIP2) protocol (Najjar and Orr 1998) and the Wanninkhof (1992) quadratic wind-speed-dependent piston velocity, using salinity to estimate alkalinity (i.e., with no carbonate cycling). Organic carbon, dissolved and particulate, is not carried as an explicit tracer but is linked to the cycling of organic phosphorus through a strict stoichiometric ratio of 106 C:P (Redfield et al. 1963). Dissolved inorganic radiocarbon (DI14C) is also carried as a prognostic tracer, with a concentration equal to its true concentration multiplied by , the standard “preanthropogenic” atmospheric 12C/14C ratio, that is, the ratio preceding contamination by fossil fuel burning and above-ground nuclear weapon testing, defined relative to tree rings from circa AD 1890 (Stuiver and Polach 1977). This multiplication causes the concentrations to be similar to those of DIC, avoiding numerical difficulties. DI14C undergoes the same air–sea exchange processes as DIC but, in contrast to total carbon, its cycling also includes uptake by plankton and export as particulate organic 14C, as well as partitioning of uptake into a dissolved organic 14C tracer, DO14C. The Δ14C (see Stuiver and Polach 1977 for the definition of this quantity) in equilibrium with the atmosphere is prescribed (equal to 0‰ in the simulations shown here), and the half-life for decay is 5730 yr.
The physical model was initialized from observed temperature and salinity, with the ocean at rest and integrated for 1000 years. Following this, the biogeochemical model was initialized from observed phosphate, DIC and DI14C, and integrated for 1000 years using preindustrial boundary conditions. The latter 500 years of this run are presented here. Note that radiocarbon is not fully equilibrated: the global average Δ14C gradually increases in the ocean throughout this interval (by 1.1‰ century−1), including a surface Δ14C increase of 0.3‰ century−1. Linear trends are removed from the 500-yr time series in all statistical analyses.
The sea surface radiocarbon simulation is difficult to evaluate precisely, given that almost all available observations were made after massive contamination of surface waters by above-ground nuclear weapons testing in the 1950s and early 1960s. Chemistry-based regression estimates of the prebomb radiocarbon activities are available (Rubin and Key 2002), but these are less reliable in surface waters and do not correct for the invasion of 14C-depleted CO2 emissions from fossil fuel burning (Suess 1953). A comparison of the model simulation to the Global Ocean Data Analysis Project (GLODAP) regression-based product shows similar large-scale surface Δ14C patterns, with simulated Δ14C values generally higher in the model (Fig. 2). The higher model Δ14C is largely attributable to the lack of correction for fossil-fuel carbon but may also reflect air–sea exchange errors. At the same time, the Δ14C in tropical upwelling regions is lower than expected from the coral-based preindustrial estimates, summarized in Toggweiler and Samuels (1993), which could reflect excessive vertical mixing across the tropical thermocline.
The distribution of radiocarbon in the ocean interior, isolated from the fossil fuel effect, agrees better with estimates of natural radiocarbon but highlights circulation error biases also apparent in the phosphate simulation (see Fig. B5). Within the Atlantic, the water mass structure is well reproduced, while in the Pacific, weak stratification of the upper ocean allows excessive ventilation from above (Fig. 2). Both the overactive North Pacific ventilation and excessive vertical mixing across the thermocline would contribute to the general overestimate of radiocarbon activities in the deep Pacific Ocean. Basinwide averages of deep ocean Δ14C, calculated according to Matsumoto et al. (2004) for the final 20 years of the simulation, lie within the error of the available observational data, as shown in Table B1. However, there is an ongoing drift to higher (more radiogenic) values, consistent with excessive ventilation of the deep ocean, that would likely cause these errors to be larger when equilibrated.
3. Variability of sea surface radiocarbon
We analyze the temporal variability of Δ14C at the ocean surface by dividing the globe into three regions: tropics, northern extratropics, and southern extratropics. Empirical orthogonal functions (EOFs) are calculated from deseasonalized monthly anomalies of surface Δ14C for the last 100 years of the simulation. We examine the mechanisms by which ocean circulation and mixing affect the leading modes of Δ14C variability at the sea surface and links to modes of climate variability.
The first EOF of tropical surface seawater Δ14C shows a dipole between the eastern equatorial Indian Ocean and the western equatorial Pacific (Fig. 3). This EOF accounts for 69% of the variability and is linked to the dynamics of the west Pacific warm pool.
Modeled interannual variability of SST in the tropics is dominated by a structure similar to ENSO. The SST anomaly in the Niño-3 region (Fig. 4, top left), has an annual peak that is somewhat larger than observed, while the interannual peak is roughly the right size but shifted toward slightly lower frequency. This small degree of disagreement between ~50 years of data and 500 years of model output may not be significant, given the pronounced heterogeneity of ENSO over interdecadal and centennial time scales in models and, presumably, the real world (Wittenberg 2009). The structure of the warm anomaly is also very similar to observations (Fig. 4, top right) with a correlation coefficient of 0.92, though the zero line is shifted to the west slightly. The wind stress anomaly (Fig. 4, bottom left) resulting from this temperature anomaly has about the right magnitude but is also shifted to the west while the wind stress curl anomaly is also about the right magnitude though shifted slightly to the north. In addition, the lagged regression of global SST monthly anomalies onto the Niño-3 index is structurally very similar between the model and observations (the first and second columns of Fig. 5).
The structural similarity suggests that the unforced model ENSO variability is operating via similar large-scale mechanisms to those underlying true ENSO variability, which can be monitored by SST within the Niño-3 region. The third column of Fig. 5 shows the regression coefficients of the monthly Δ14C anomaly onto the Niño-3 SST anomaly. The El Niño phase of ENSO, characterized by a high Niño-3 SST anomaly (Fig. 5, zero lag), is well correlated with positive anomalies of surface Δ14C in the western equatorial Pacific Ocean. The correlation arises from a weakening of the easterly winds during El Niño events, which reduces the upwelling of 14C-depleted subsurface waters in the central and eastern equatorial Pacific (Toggweiler et al. 1991). As a result of the reduced zonal gradient of surface Δ14C combined with a weakening of the equatorial counter currents, eastward advection of 14C-rich waters becomes reduced in the western equatorial Pacific. Consequently, the western pool of 14C-enriched surface waters intensifies and expands eastward, resulting in a net accumulation of high-Δ14C waters in the western and central equatorial Pacific Ocean, as previously shown by Rodgers et al. (1997).
Meanwhile, the transport of Pacific waters by the Indonesian Throughflow has a large impact on the surface Δ14C of the eastern Indian Ocean. Illustrating this, the linear correlation of Δ14C near Lombok Strait (11°S, 115°W) with the volume transport of the Indonesian Throughflow is r = 0.69 with a slope of ~2‰ Sv−1 (Sv ≡ 106 m3 s−1). A large part of the Indonesian Throughflow variability is related to reduced transport during El Niño events (~10 Sv rather than ~14 Sv) which restricts the supply of 14C-rich waters from the Pacific warm pool and results in the negative correlation shown in Fig. 5 for the eastern Indian Ocean (right column). Reduced Indonesian Throughflow transport during El Niño events has also been well documented by observations (Tillinger and Gordon 2009). However, the correlation of the Δ14C near Lombok Strait with Niño-3 SST (r = 0.55) is less than with the Indonesian Throughflow itself since the throughflow responds to more than just ENSO. Therefore, the Δ14C of the eastern Indian Ocean also reflects these additional sources of Indonesian Throughflow variability, as well as local changes in circulation.
We note that the response pattern of Δ14C to ENSO is different from that of SST to ENSO in two distinctive ways (Fig. 5, compare middle and right columns). While the SST response is characterized by an out of phase relationship between the tropical and extratropical Pacific, the Δ14C response is dominated by a zonal contrast between the tropical Pacific and Indian Oceans, part of which reflects the change in the strength of the Indonesian Throughflow. In addition, the Δ14C response persists with very little alteration at a 6-month lag, unlike the SST signal which has largely decayed by this point. This persistence reflects the long air–sea equilibration time of Δ14C in seawater (~10 yr), allowing it to maintain a more distinct memory of ocean circulation changes: as a result, the ENSO signal persists throughout an entire year, providing ample opportunity for it to be recorded in corals (Druffel 1997).
b. Northern Hemisphere
Interannual variability of surface seawater Δ14C in the Northern Hemisphere is dominated by the subarctic Pacific, which displays two distinct patterns (Fig. 6). The first EOF is centered on the Aleutian Islands and accounts for 36% of the variance. The second EOF reveals a dipole, with opposing poles in the Gulf of Alaska and western subarctic Pacific, accounting for 19% of the variance. Both EOF patterns are evident in subsurface waters to a depth of ~300 m.
The North Pacific has a circulation bias that includes an erroneously low-latitude position of the Oyashio extension, a weak Alaskan stream, an underexpression of the western subarctic gyre, and excessive ventilation to intermediate depths. Given these biases in the simulated mean state, the simulated variability must be viewed with caution.
The two leading EOFs of surface Δ14C reflect a basinwide propagating mode. The first EOF leads the second by ~4 yr (correlation coefficient of 0.55) while the second EOF leads the first by ~13 yr (correlation coefficient of 0.46). Basinwide decreases of surface Δ14C (the first EOF) occur with weaker subpolar gyre circulation. This EOF correlates with anomalous eastward flow of the upper-ocean currents north of 40°N and reduced exchange with the subtropical gyre (Fig. 7). This results in anomalous advection of 14C-depleted waters from the northwest Pacific basin to the northeast Pacific, reducing the zonal gradient. In addition, increased mixed layer depths occur in the northwestern edges of the Bering Sea. Entrainment of 14C-depleted subsurface waters into the mixed layer contributes to the basinwide reduction of surface Δ14C.
In contrast, the second EOF of the surface Δ14C anomaly is well correlated with a strengthening of westerly winds over the Kuroshio Extension region. Enhanced wind stress curl over the subpolar gyre intensifies the gyre circulation (Fig. 7), resulting in anomalous southward flow of 14C-depeleted surface waters in the western subarctic gyre and anomalous northward flow of 14C-rich subtropical water on the eastern side. This increases the zonal contrast of surface Δ14C in the subpolar North Pacific.
The time series for the second leading EOF of surface Δ14C anomaly correlates weakly (r = 0.14, 95% significance) with the Arctic Oscillation (AO) index (the simulation of which is shown in Fig. B4) owing to its dependence on the westerly winds. Meanwhile a less weak correlation (r = −0.35, 95% significance) occurs between the second EOF and the Pacific decadal oscillation time series (Mantua et al. 1997). The first EOF is similarly correlated with the Pacific decadal oscillation, but with the opposite sign (r = 0.33, 95% significance), while it shows no relationship with the Arctic Oscillation, given its less direct dependence on the westerly winds.
c. Southern Hemisphere
The first EOF of sea surface Δ14C in the Southern Hemisphere extratropics has the largest amplitude of anywhere in the ocean (Fig. 8). This variability is focused within the Weddell Sea. The second EOF is likely part of a propagating mode, of which a main standing mode is represented by the first EOF. The third EOF is much weaker.
In assessing the causes underlying this sea surface Δ14C variability, we consider first the Southern Annular Mode (SAM), a dominant pattern of climate variability over the Southern Ocean on interannual to decadal time scales. The SAM is characterized by meridional fluctuations of air masses south of 20°S, associated with changes in sea level pressure and the location and intensity of the Southern Hemisphere westerlies (e.g., Hall and Visbeck 2002; Marshall 2002). Figure 9 shows the first EOF of the monthly sea level pressure anomaly south of 20°S as obtained from the model and observations. Although the zonal structure is somewhat different than observed, the model captures the gross features of the observed leading EOF pattern with a correlation coefficient of 0.96. The positive SAM correlates well with small negative anomalies of Δ14C (~1‰) exclusively in the Pacific sector of the Southern Ocean, reflecting the third EOF. This association can be attributed to a strengthening of the westerlies during the positive phase of the SAM, inducing northward Ekman transport and bringing more 14C-depleted water from south of the Antarctic Circumpolar Current (ACC) to the north.
However, much larger amplitude variability occurs in the Weddell Sea when periodic convective events exhume 14C-depleted waters from the deep sea; similar, but weaker events occur in the Ross Sea as well as in the Indian sector of the Southern Ocean. These events arise from disruptions of the near-surface stratification. It is entirely possible that this mode of convection is an artifact of the coarse model resolution, particularly the lack of dense water production on shelves and weak downslope flows (Toggweiler et al. 2006). Nonetheless, the convective events are reminiscent of the Weddell Sea polynya that occurred between 1974 and 1976 (Carsey 1980), when vigorous convection occurred to a depth of 3000 m (Gordon et al. 2007). Given the similarity to what may be a natural mode of variability in the Weddell Sea and the large impact these events have on radiocarbon, we discuss the convective events in some detail.
Examination of the average annual-mean temperature stratification in the Weddell Sea (Fig. 10a) shows that, whereas surface temperatures are higher during surface 14C-depletion events, the waters at 300 m are actually slightly colder during convection, as wintertime cooling stirs cold water to depth. The high surface temperatures occur because convection brings up enough warm Upper Circumpolar Deep Water during the wintertime to maintain an ice-free polynya, facilitating more absorption of solar radiation and much warmer springtime temperatures. Since these changes result in increased stratification, they cannot by themselves explain the increase in convection. Instead, the cause is a breakdown in the salinity stratification, with convection exposing low Δ14C water at the surface when the salinity gradient between the surface and 300 m collapses (Fig. 10b).
We can use times when there is an anomalous decrease in surface radiocarbon [beyond one standard deviation (std dev)] to detect convective onset and times when there is an anomalous increase in surface radiocarbon to detect convective shutoff. As shown in Fig. 10c, five years before the onset of convection, salinities are lower than average to the west of the Antarctic Peninsula and higher to the east of 40°W. At the time of convective onset (Fig. 10d), the saline anomaly has strengthened and propagated into the central gyre. This is similar to the preconditioning mechanism for the Weddell Polynya suggested by Gordon et al. (2007), who related extended periods of a negative SAM to a buildup of salinity at the surface, destabilizing the upper water column. In the model, the simulated high-salinity anomaly reaches a maximum at the height of convective activity, partly due to mixing with the saline waters below (Fig. 10e). Finally, as fresher water is advected from the east along the southern flank of the gyre, the high salinity, low Δ14C anomaly dwindles and the convection shuts off (Fig. 10f). Convective events are generally followed by a multidecadal period of greater near-surface stability in the Weddell Sea, with higher Δ14C at the surface.
4. Variability in the air–sea flux of radiocarbon
The prior section discussed redistributions of radiocarbon within the ocean due to changes in ocean mixing and advection. However, oceanic radiocarbon is also altered by changes in the air–sea balance of radiocarbon that, in turn, alter the radiocarbon activity of the relatively small atmospheric reservoir. Although we do not explicitly model atmospheric 14C here, we can infer its behavior by analyzing the air–sea fluxes. Below, we illustrate how the sea surface Δ14C distribution and the wind-speed-dependent air–sea exchange both vary the air–sea flux of radiocarbon.
a. The disequilibrium flux
The air–sea fluxes of 12C and 14C are highly correlated in space and time, as are their anomalies (Fig. 11). This is because, despite globally lower saturation states for 14C, most of the spatial and temporal variability in the saturation state of both isotopes is driven by changes in water temperature and the biological pump of carbon. To highlight differences between the fluxes, which are responsible for changing the Δ14C of the atmosphere, we use the concept of a “disequilibrium 14C flux,” as introduced by Rodgers et al. (2011). This construct simply normalizes the air–sea flux of 14C relative to that of 12C, to avoid the unwieldy fact that their concentrations differ by approximately 12 orders of magnitude, allowing a straightforward comparison. The disequilibrium flux is defined as
where is the air–sea 12C flux and the air–sea 14C flux. The disequilibrium flux is given as 12C-equivalent units of 14C and can be integrated globally to provide the net effect of air–sea exchange on atmospheric Δ14C. According to the convention used here, a positive disequilibrium flux indicates a net transfer of 14C into the ocean such that the atmospheric Δ14C would decrease.
Figure 11 shows that the disequilibrium flux is a relatively invariant residual between the much more variable 14C and 12C fluxes. Moreover, the temporal structure of the disequilibrium flux is very different from that of 14C and 12C. For example, four intense multiyear anomalies, associated with the Weddell Sea convection events discussed above, show strong, highly correlated decreases in the flux of both 12C and 14C into the ocean owing to the high dissolved inorganic carbon (DIC) concentration of the exposed deep waters (Fig. 11). The high degree of correlation causes changes in the disequilibrium flux to be relatively small.
Figure 12 provides an alternate perspective, showing the interannual variability of the zonally integrated disequilibrium flux, compared to the interannual variability of surface Δ14C. The surface Δ14C variability is large at high latitudes of the North Pacific and Southern Ocean, as expected from the analysis in section 3. However, the strongest disequilibrium flux variability is centered at 58°S. As discussed by Rodgers et al. (2011), this portion of the Southern Ocean also dominates the magnitude of the integrated global disequilibrium flux due to the low 14C content of deep waters drawn up to the south and the large area in which strong winds drive rapid air–sea exchange.
b. Mechanisms of variability
To understand the mechanisms underlying the the air–sea balance of radiocarbon, it is instructive to recall that, in the model, the air–sea exchange flux scales with
where U10 is the wind speed at 10-m elevation, C(surf) is the sum of dissolved CO2 and H2CO3, and C(surf)sat is C(surf) at saturation (largely a function of temperature and salinity). Because 14C(surf)sat is equal to C(surf)sat , where is the atmospheric 14C/12C, the air–sea flux of 14CO2, , scales with
In the special case where atmospheric Δ14C is 0‰, that is, , we can define 14C(surf)* as 14C(surf) (similar to our disequilibrium flux convention) so that, by substitution of Eqs. (2) and (3) in Eq. (1), the disequilibrium flux φdiseq scales with
Figure 13a shows [C(surf) − 14C(surf)*] at the surface layer of the Southern Hemisphere during May, a time of vigorous exchange. High values occur close to Antarctica owing to the mixing with waters below, north of which they decrease; this is very nearly a mirror image of the surface Δ14C. Figure 13b shows the time-averaged disequilibrium flux φdiseq. The disequilibrium flux is clearly focused within the southern portion of the westerly wind belt, indicated here as the region between the two black contour lines (time-averaged surface winds between these contours exceed 9 m s−1). The high disequilibrium flux per unit area within this belt is due to the coincidence of large [C(surf) − 14C(surf)*] and strong winds. In addition to the weaker winds, air–sea exchange at higher latitudes is impeded by winter sea ice (not shown).
The final panel of Fig. 13 shows the correlation coefficient (r) between the local disequilbrium flux and the global disequilibrium flux (both smoothed with a 40-yr boxcar filter). The relatively high correlation within the westerly wind belt is due to the strong role that this zone plays in controlling the global air–sea balance of radiocarbon. Notably, the local disequilibrium flux north of ~45°S is often anticorrelated with the global disequilibrium flux owing to the fact that winds in this region tend to be weaker when winds intensify within the westerly belt, that is, during positive phases of the SAM; however, this more northerly region has little impact on the global integral because of the small magnitude of the disequilibrium flux there (Fig. 13b). Instead, the southern part of the wind belt dominates the air–sea balance variability.
The control by the winds is locally very important on short time scales, with very strong local correlation (r > 0.9) at most points in the ocean, because of the quadratic dependence of the disequilibrium flux on wind speed (Wanninkhof 1992) and the fact that (C(surf) − 14C(surf)*) is positive everywhere in the surface ocean. However, on longer time scales, the globally integrated flux is dominated by the air–sea gradient. This dominance is illustrated by the temporal correlation of the global disequilibrium flux and the average Δ14C of the global surface ocean on long time scales (r = −0.87, boxcar 20 yr smooth, rather than r = −0.54 for interannual variations). Because the westerlies drive the upwelling of 14C-depleted waters to the Southern Ocean surface, in addition to their role in accelerating local air–sea exchange within the wind belt, they play a central part in setting the air–sea balance of radiocarbon (Rodgers et al. 2011); however, the high frequency of the model SAM limits its ability to drive decadal–centennial changes in upwelling in the simulations shown here.
c. Simulated changes in atmospheric radiocarbon
Despite the demonstrated potential for the Southern Ocean to drive changes in the air–sea balance of 14C, the variability exhibited by the simulation shown here would have little impact on the resultant atmospheric Δ14C. We estimate the expected atmospheric response by removing the detrended average flux of φdiseq into the ocean, which balances the decay of radiocarbon within the global ocean and model drift, to estimate the effect of variable air–sea exchange on the atmospheric 14C and 12C inventories. Starting from 0‰, we integrate the air–sea fluxes, assuming an atmospheric C inventory of 600 Pg, a terrestrial biosphere of 2200 Pg, and an exchange between them of 60 Pg yr−1 (Hughen et al. 2004). Note that this calculation will overestimate true changes in atmospheric Δ14C since, over time, an increase (decrease) in atmospheric Δ14C will drive an increased (decreased) disequilibrium flux back into the ocean, ignored by our use of a constant atmospheric Δ14C of 0‰. Figure 14 shows the prediction of atmospheric Δ14C variability that results from this calculation over the last 400 yr of unforced variability in the model. Comparison to Fig. 1 clearly shows that the magnitude of atmospheric Δ14C variability estimated from model φdiseq is much smaller than that reconstructed from tree-ring records of atmospheric Δ14C: less than ±1‰, an order of magnitude smaller than the reconstructed changes.
Two possible, nonexclusive conclusions could be drawn from the small amplitude of the simulated variability. The first possibility is that, despite the apparent fidelity of the simulated modes of interannual variability, unforced variability in the model is significantly smaller than that of the real world on multidecadal time scales. Indeed, the lack of variability in the radiative boundary conditions, such as fluctuations in solar input, volcanic activity, and stratospheric chemistry, could be expected to produce an insufficiently variable simulation. In particular, if southern westerlies and/or stratification of the polar haloclines have undergone significant low-frequency variability, it could have resulted in larger disequilibrium flux variability than simulated here. This may be testable against observations in future, given well-constrained 14C production rates and through comparison to the interhemispheric gradient of atmospheric Δ14C derived from tree rings (Rodgers et al. 2011).
The second possibility is that natural, decadal to centennial time-scale changes in air–sea exchange had little impact on atmospheric Δ14C during the preindustrial, postdeglacial period. This possibility would require the production rate of 14C by cosmic rays to have been the primary driver of observed atmospheric Δ14C variability during the Holocene. Indeed, this argument was put forth by previous workers based on scaling arguments, including Stuiver and Quay (1980) and Siegenthaler et al. (1980). In addition, a production-rate control on atmospheric 14C is supported by radionuclide production records derived from 10Be abundances in ice cores (Bard et al. 2000; Muscheler et al. 2007) and is commonly assumed by those interpreting the atmospheric radiocarbon record as an indicator of solar output over the Holocene (Knudsen et al. 2009). If true—that is, if the simulated disequilibrium flux variability is not dramatically less than the true variability—this would simplify the interpretation of both atmospheric and oceanic records of Δ14C over the Holocene.
The CM2Mc model simulation explored here exhibits modes of interannual variability similar to ENSO, the AO, and the SAM, as characterized by observations. Overall, the CM2Mc tropical and ENSO biases are qualitatively similar in character to those seen in the higher-resolution model, CM2.1 (Wittenberg 2009), though somewhat larger in amplitude. This suggests that CM2Mc may provide a cost-effective new tool to explore the causes of these biases, which are shared among many state-of-the-art CGCMs (Guilyardi et al. 2009) and complicate future projections of tropical climate and ENSO (Collins et al. 2010; Vecchi and Wittenberg 2010). In addition, the model is well suited for long (millennial time scale) simulations.
The distribution of radiocarbon within the model was well correlated with these modes of interannual variability, with distinct relationships linking tropical 14C to Niño-3 SST, as previously shown with forced ocean models (Rodgers et al. 1997). The modeled surface water Δ14C at the Galapagos site (Fig. 14) is remarkably stable in contrast to the observed variability (Fig. 1). This disagreement suggests that the unforced model underestimates true variability at this site, perhaps due to difficulties in capturing the fine details of the equatorial current structure of Galapagos and/or resolving tropical instability waves, possibilities that we leave to further investigation. In contrast, in the vicinity of the Indonesian Throughflow, the modeled variability (1 std dev = 3.8‰) is of similar magnitude to that observed in corals (1 std dev = 4.2 ‰, Guilderson et al. 2009). Although there are certainly problems with the model’s ability to resolve the complex Indonesian Throughflow region as well, it does support this location as a particularly sensitive location in which to use coral Δ14C as a monitor of circulation changes (Fig. 14) owing to the coincidence of strong Δ14C gradients with variable water transport.
Meanwhile, the variability of 14C at the surface of the Southern Ocean was dominated by changes in deep convection on a multidecadal time scale. In the model, these convection events are not obviously related to the relatively high frequency SAM but are associated with dramatic changes in sea surface salinity. In the Weddell Sea, salty anomalies originate near the prime meridian due to atmospheric forcing and advect southwestward toward the gyre center where they set off vigorous convection. The convection brings warm, 14C-depleted, deep water to the surface until the halocline reestablishes itself (Fig. 14). Given the difficulty of representing convection accurately in a global, level-coordinate ocean model, the details of this mechanism must be regarded with some caution. Indeed, they are reminiscent of “flushing” events, noted in coarse-resolution models with mixed boundary conditions (Winton and Sarachik 1993) and may reflect the inability of the model to produce dense waters on shelves, as occurs in nature (Toggweiler et al. 2006). Nonetheless, the Weddell Polynya of 1974–76 appears to have been very similar to our simulated convective events (Gordon et al. 2007), suggesting that this is a real mode of variability, even if perhaps exaggerated by the model. If so, radiocarbon records may help to document the occurrence of Weddell Sea polynyas that predated the satellite era (e.g., Domack et al. 2001).
Although the simulated variability of Δ14C in the Weddell Sea was dramatic, the resulting local disequilibrium flux anomalies did not dominate the global integral, given their restricted areal extent. Instead, the global air–sea disequilibrium flux was dominated by exchange within the westerly wind belt over the open Southern Ocean. There, high wind speeds and ice-free winters conspire with the low Δ14C of upwelling Upper Circumpolar Deep Water to drive most of the global disequilibrium flux, as discussed by Rodgers et al. (2011). When the Δ14C of waters exposed in this region of vigorous air–sea exchange decreased, an increased disequilibrium flux of 14C into the ocean tended to occur.
Nonetheless, the magnitude of simulated air–sea disequilibrium anomalies was insufficient to alter atmospheric Δ14C by more than a few permil (Fig. 14), far less than the magnitude of the Suess wiggles (~20‰, Fig. 1). This discrepancy may indicate that natural variability in the climate system, produced by external forcings or internal feedbacks, significantly exceeds that of the model. In particular, large changes in the position or strength of the Southern Hemisphere westerlies may have caused sustained changes in the disequilibrium flux not captured here (Rodgers et al. 2011). However, previous authors have presented evidence that the Suess wiggles were dominated by changes in the production rate of 14C, caused by fluctuations in the geomagnetic field and solar activity (Stuiver and Quay 1980; Bard et al. 2000; Muscheler et al. 2007), and the results here do not provide any reason to disagree with this interpretation.
This paper has used an earth system model simulation to explore Δ14C variability at the sea surface on decadal to centennial time scales, given steady boundary conditions and constant atmospheric Δ14C, and linked this variability to modeled climate modes. Future modeling work could explore the effect of changes in radiative (or other) boundary conditions and how variability in the production rate of 14C propagates through the system. The development of new high-resolution radiocarbon records in regions with significant Δ14C variability could help to constrain climate modes on long time scales, identify model biases, and contribute to developing improved understanding of natural climate variability.
We thank Ellen Druffel for providing coral data and for comments on an early draft of the manuscript, and Michael Winton, Matthew Harrison, and Hyo-Seok Park for discussions. Sergey Malyshev assisted with the land model. Robbie Toggweiler and Stephanie Downes provided internal reviews. Karen Assmann and two anonymous reviewers provided thoughtful, constructive comments. The OAA_OI_SST_V2 data was provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, from their Web site at http://www.esrl.noaa.gov/psd/ and ECMWF ERA-40 data were provided by ECMWF and have been obtained from the ECMWF Data Server. This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), the Canadian Institute for Advanced Research (CIFAR), and the Carbon Mitigation Initiative (CMI) project at Princeton University, sponsored by BP and Ford Motor Company. This paper was prepared by Eric Galbraith under awards NA17RJ2612 and NA08OAR4320752 from the National Oceanic and Atmospheric Administration (NOAA), U.S. Department of Commerce. The statements, findings, conclusions, and recommendations are those of the authors and do not necessarily reflect the views of the National Oceanic and Atmospheric Administration or the U.S. Department of Commerce.
This paper presents the Coupled Model 2 with Modular Ocean Model version 4p1 (MOM4p1) at coarse resolution (CM2Mc). This model is similar to the Geophysical Fluid Dynamics Laboratory (GFDL) 1° coupled model (CM2.1) (Delworth et al. 2006), used in the Intergovernmental Panel on Climate Change Fourth Assessment Report. However, the resolution of both atmospheric and oceanic grids has been coarsened so as to decrease computational cost, and the model uses an updated version of the GFDL code, prepared in anticipation of a suite of simulations for the IPCC Fifth Assessment Report. This updated code includes the ocean model MOM4p1 (Griffies 2009), which includes a number of numerical and physical improvements, as well as the option of using pressure for the vertical coordinate instead of depth so that steric sea level is directly simulated. Modifications made to the CM2.1 atmosphere are minor, involving a reduction of spatial resolution and minimal parameter changes. More substantial changes were made to the ocean component, including improvements to the subgrid-scale parameterizations of ocean mixing, which is of critical importance to low-resolution climate models incapable of explicitly representing processes such as mesoscale eddies and gravity-driven bottom flows. CM2Mc requires approximately one tenth the computational cost of the higher resolution model CM2M. Additional model configuration details and diagnostics can be found on the Web site http://sites.google.com/site/cm2cmodel/.
a. Atmospheric model configuration
The atmosphere is very similar to the GFDL Atmospheric Model 2 (AM2), as described by the GFDL Global Atmospheric Model Development Team (Anderson et al. 2004), but using the finite volume dynamical core of Lin (2004), as implemented in CM2.1 (Delworth et al. 2006). The atmosphere used here employs the M30 grid, with a latitudinal resolution of 3° and a longitudinal resolution of 3.75°, instead of the M45 grid used in CM2.1 (2° by 2.5°), but uses the same vertical configuration with 24 levels. The decreased horizontal resolution results in 44% of the number of grid cells in the CM2.1 atmosphere, a relatively modest decrease in resolution. However, the atmospheric time steps are also increased, from 0.5 to 1.5 h for the tracer time step, from 6 to 9 min for the dynamical time step, and from 2 to 3 h for the radiative time step. The coupled model includes an explicit representation of the diurnal cycle of solar radiation.
The change in discretization required a small readjustment of two cloud parameters (the critical droplet threshold radius and the cloud erosion time scale) to maintain a net radiation balance close to zero under modern climate, which had negligible impact on atmospheric dynamics. In addition, the mountain gravity drag parameter G* (Anderson et al. 2004) was decreased from 1.0 to 0.5 to compensate for the coarser horizontal resolution. The final alteration was the replacement of the sea salt aerosol climatology of Haywood et al. (1999), used for the radiative conditions, with a more recent estimate from Ginoux et al. (2006). As in CM2.1, the Land Dynamics model of Milly and Shmakin (2002) is used, including a river routing scheme but no terrestrial ecosystem (see Anderson et al. 2004).
b. Ocean grid and bathymetry
The CM2Mc ocean uses a tripolar grid (Murray 1996) to avoid the grid singularity at the north pole, as in the GFDL CM2.0 and CM2.1 models (Griffies et al. 2005; Gnanadesikan et al. 2006). The horizontal dimensions of the ocean grid vary according to latitude, with the finest latitudinal resolution of 0.6° at the equator allowing an explicit physical representation of some equatorial currents. There are 28 vertical levels, the uppermost eight of which are each 10 dbar thick below which the levels gradually increase in thickness to a maximum of 506 dbar. The model employs partial bottom cells (Adcroft et al. 1997; Pacanowski and Gnanadesikan 1998) to allow a more realistic representation of the bathymetry.
The bathymetry was initially generated from averaging satellite-derived bathymetry onto the coarser model grid with a smoothing algorithm. The result was then carefully examined, and critical features destroyed in the regridding process were manually restored. To compensate for mixing between ocean basins and marginal seas where the coarse resolution of the topography does not resolve the narrow connections, the cross-land mixing scheme of Griffies et al. (2005) is employed, generating less than 2 Sv of total mixing at each of six locations. Finally, where the discharge locations for precipitation falling on Antarctica led to the margins of the Ross or Filchner–Rønne Ice Shelves, they were shifted north of the coast by a few degrees, so as to roughly account for the transport of freshwater by drifting icebergs.
c. Ocean–ice model and parameterizations
CM2Mc uses the MOM4p1 code with pressure as the vertical coordinate, thus allowing the model to solve the non-Boussinesq primitive equations. This is coupled to the GFDL thermodynamic–dynamic sea ice model (SIS) (Delworth et al. 2006). The vertically integrated mode is time stepped using an explicit bottom pressure solver, and the transfer of water mass across the ocean surface occurs using real freshwater fluxes. Tracer advection uses the Multidimensional Piecewise Parabolic Method (MDPPM) scheme for temperature and salinity (Adcroft et al. 2004), and the Sweby Multidimensional Flux Limited (MDFL) scheme for other tracers (Griffies et al. 2005). Light is absorbed by water and a chlorophyll field following the Manizza et al. (2005) algorithm. Over the initial millennium of simulation a smoothly varying, satellite-derived climatological chlorophyll field was used, after which the chlorophyll concentrations used for shortwave absorption were those predicted by a coupled biogeochemical model [Biogeochemistry with Light, Iron, Nutrients and Gases (BLING), described by Galbraith et al. (2010), not shown here]. The gross metrics of the simulation did not change significantly with the change to prognostic chlorophyll.
Most subgrid-scale parameterizations for mixing are similar to those used in the CM2 series (Griffies et al. 2005). The lateral friction uses an isotropic Smagorinsky viscosity in midlatitudes, while within 20° of the equator the anisotropic National Center for Atmospheric Research viscosity is used (Large et al. 2001; Smith and McWilliams 2003). Lateral diffusion and skew diffusion of tracers along isopycnals is represented using the parameterization of Gent and McWilliams (1990) with a spatially varying diffusion coefficient, AGM, as determined in CM2.1 (Griffies et al. 2005). The coefficient depends on the horizontal shear between 100 and 2000 m. A minimum coefficient of 200 m2 s−1 and a maximum coefficient of 1400 m2 s−1 are imposed. The slope-dependent thickness transport also stops increasing at a value of AGMSmax, where the maximum isopycnal slope Smax is set to 0.01, whereas it is set to 0.002 in CM2.1 (see Gnanadesikan et al. 2007 and Farnetti et al. 2010 for discussions of this parameter). Overall, these changes generate a more active eddy mixing parameterization than CM2.1, with a greater responsiveness to changes in baroclinicity, as suggested by recent observations in the Southern Ocean (Böning et al. 2008). Meanwhile, the neutral diffusion coefficient for tracers, AI, is set to a constant value of 800 m2 s−1.
Within the boundary layer, we use the K-profile parameterization of Large et al. (1994). Away from the boundary layer, a background diffusivity of 0.1 × 10−4 m2 s−1 and a background viscosity of 1 × 10−4 m2 s−1 are used. Below 300 m, these background coefficients are punctuated by locally intense vertical mixing, as shown by tracer studies (Naveira Garabato et al. 2004). This local mixing occurs as part of the K-profile parameterization, where convection is driven by along-isopycnal mixing and where the shear is large relative to the stratification, as well as where gravity waves are generated over regions of rough topography and break as they propagate upward into the water column. The latter process is parameterized with the tidally dependent mixing scheme of Simmons et al. (2004), which depends on local stratification as well as prescribed bottom roughness and tidal amplitude. In addition, tidal mixing in coastal areas is simulated by the mixing scheme of Lee et al. (2006).
d. Marine biogeochemical model
We use a new idealized ocean biogeochemistry (iBGC) model to simulate radiocarbon distributions, based on simple representations of nutrient uptake and remineralization. The effect of organic matter cycling on the radiocarbon distribution is small, but we describe the model briefly here for completeness.
Conceptually, this model is very similar to the pioneering biogeochemical model of Bacastow and Maier-Reimer (1990). In brief, it carries a single prognostic inorganic macronutrient tracer, which we call phosphate (PO4), and represents the net activity of phytoplankton through an uptake function in all ocean grid cells as a function of temperature T and available light I. Once taken up, the macronutrient is instantaneously partitioned between a dissolved organic phosphorus (DOP) tracer and a flux of sinking particulate organic phosphorus (POP), according to a fixed fraction φ. The sinking of POP is not achieved by using a prognostic tracer but, instead, by instantaneously redistributing remineralized PO4 throughout the water column beneath the layer in which uptake occurred, similar to the OCMIP2 scheme (Najjar and Orr 1998). Meanwhile, DOP decays to PO4 according to a temperature-dependent first-order rate, , where kr gives the temperature dependency of remineralization and T is temperature (°C). The total mass of phosphorus (PO4 + DOP) is conserved within the ocean. The overall equation for PO4 is
where circ indicates the net effect of ocean advection and mixing, Vmax is the maximum uptake rate, k is a constant determining the temperature sensitivity of uptake, Ik is a constant determining the tendency for light limitation, remiPOP is the remineralization source of phosphorus from sinking POP, and kPO4 is a constant determining the tendency for nutrient limitation.
All parameters used in iBGC are given in Table A1.
a. Atmosphere-only simulation
As a preliminary test, the atmospheric model was integrated for twenty years with prescribed SST and constant “1990” radiative conditions, following the Atmospheric Model Intercomparison Project method described by Anderson et al. (2004). In general, this simulation shows a very similar pattern of biases to the 2° configuration of the atmosphere, with slightly greater errors for all variables, and a radiative balance of +0.45 W m−2. The atmosphere was also tested with preindustrial radiative conditions and modern SST, in which case the radiative balance was −1.67 W m−2.
b. Coupled simulation
When coupled to the ocean model, the atmospheric simulation exhibits some pronounced differences in comparison to the prescribed-SST run, generally consistent with the effect of coupling the 2° M45 atmosphere to the 1° ocean of CM2.1 (Delworth et al. 2006; Wittenberg et al. 2006). Key differences include the appearance of a “split” intertropical convergence zone, that is, excessive rainfall in the Pacific just south of the equator (Fig. B1), and an increase in the excess heating above eastern boundary upwelling regions and the Southern Ocean (Fig. B1).
When forced with preindustrial radiative conditions, the global energy balance of the coupled model is significantly more positive than in the prescribed-SST case. The net radiation imbalance quickly rises as the ocean surface adjusts, stabilizing at +0.68 W m−2 over the latter half of the first century, then gradually declines as the deep ocean warms, reaching ~0.4 W m−2 by the end of the simulation (approximately 2000 years). Note that the atmospheric model does not conserve energy perfectly, consistently leaking ~0.4 W m−2 owing to the fact that, unlike the ocean and land models, the atmosphere does not account for the heat content of water (Delworth et al. 2006). A time series of ocean temperature shows that the model is very close to thermal balance at the end of the simulation (not shown).
The most significant temperature and salinity errors of the ocean simulation (Fig. B2) can be attributed to the errors of the atmospheric simulation. Foremost among these are 1) a cold bias in the North Pacific, particularly during boreal summer; 2) a warm bias in the Southern Ocean, particularly during austral summer; 3) warm biases in all eastern boundary upwelling regions; and 4) a low salinity bias in the South Pacific and high salinity bias in the North Pacific (cf. Figs. B1 and B2). The first, second, and third errors can be attributed to the regional energy balances, apparently due to inaccuracies in the representation of clouds that produce very similar errors in absorbed shortwave radiation in both the prescribed SST and coupled simulations. Meanwhile, the fourth, a result of the split ITCZ, arises from the coupling between the ocean and atmosphere in the tropical Pacific, and is a common symptom of coupled models (de Szoeke and Xie 2008). The warm bias of the Southern Ocean, combined with a tendency for North Atlantic Deep Water (NADW) to originate as excessively warm and salty, leads to a globally warm bias of the whole ocean temperature, 4.6°C (versus 4.3°C for CM2.1 and 3.0°C for observations, Conkright et al. 2002). Gross metrics of the ocean simulation are given in Table B1.
Compared to CM2.1, the overturning circulation is significantly different, most importantly in the circulation of deep waters formed in the Southern Ocean. CM2.1 ventilates a continuous water column throughout the Weddell Sea, evident in the ideal age and salinity sections (see Fig. 7 of Gnanadesikan et al. (2006)) and in disagreement with observations. In contrast, CM2Mc ventilates at a lower rate, closer to the Antarctic continent, evident in the overturning cell at σ2 ~ 1037 (cf. Fig. B3, top right, of this paper with Fig. 6d of Gnanadesikan et al. 2006). Despite the lower ventilation rate, the bottom waters produced by CM2Mc circulate more vigorously throughout the global abyss (cf. Fig. B3, top panels, with Figs. 6c,d of Gnanadesikan et al. (2006)), consistent with observations, such as the strong bottom water flows recently observed near the Kerguelen Plateau (Fukamachi et al. 2010). The contrast in Southern Ocean ventilation appears to be due to a combination of the eddy parameterization (which enhances the eddy flux to high latitudes in CM2Mc, thereby abetting deep sinking), the altered freshwater discharge map (which allows higher nearshore salinities in CM2Mc), and the model grid. The supply of NADW to the deep ocean is very similar between the two models.
The Arctic Oscillation, based on the first leading mode of sea level pressure anomaly north of 20°N, is shown in Figure B4 for both the model and observations. Plots for the Southern Annular Mode and ENSO are given in the main text.
c. Marine biogeochemical simulation
The marine biogeochemical tracers were integrated synchronously with the physical simulation for the latter half of the run (approximately 1000 years). The mean-state description here, and in section 2, pertains to the final 100 years of this run, while the discussion of variability includes the latter 500 years of the run.
The simulated PO4 field is fairly close to observations, especially considering the simplicity of the biogeochemical model and errors in the simulated ocean circulation. The principal features of the surface PO4 field are well captured (Fig. B5), including the very elevated Southern Ocean concentrations, low oligotrophic gyre concentrations, and elevated concentrations in low-latitude upwelling regions. We note that the high nutrient surface regions occur in the absence of explicit iron limitation.
The largest discrepancy versus observations occurs in the North Pacific where the surface nutrient concentration is dramatically lower than observed. This can be attributed to errors in the physical ocean circulation, such that nutrients are not maintained at intermediate depths where they can be resupplied to the surface (Fig. B5). A combination of the atmospheric simulation bias, poor resolution of topography near the Oyashio–Kuroshio confluence, and incorrect vertical mixing may be to blame. In contrast, the Atlantic meridional section shows PO4 concentrations that are close to the observations (Fig. B5).