Mechanisms governing interannual variability in upper-ocean inorganic carbon system and air–sea CO2 fluxes: Physical climate and atmospheric dust

,


Introduction
The ocean exhibits variability in physical circulation on subannual to decadal and longer time-scales that in turn drives substantial changes in regional to basin-scale biogeochemistry and air-sea CO 2 fluxes (Chavez et al., 1999;Le Que ´re ´et al., 2000, 2003;Gruber et al., 2003;Dore et al., 2003;Bates et al., 2003;Corbiere et al. 2007).Variations in aeolian dust deposition may generate additional biogeochemical variability in iron-limited ocean regions (Jickells et al., 2005;Patra et al., 2007;Cassar et al., 2007;Tagliabue et al., 2008;Aumont et al., 2008).Such variability offers unique, natural perturbation ''experiments'' for probing underlying mechanisms and characterizing potential responses to and feedbacks on natural and anthropogenic climate change (Boyd and Doney, 2002).Variability in air-sea CO 2 fluxes also imprints on the overlying atmosphere (Nevison et al., 2008), and thus directly impacts on so-called ''top-down'' efforts to reconstruct oceanic and terrestrial carbon sources and sinks by inverting the time/space evolution of atmospheric CO 2 (e.g., Bousquet et al., 2000;Gurney et al., 2002;Peylin et al., 2005).
The amplitude, patterns and phasing of ocean interannual variability are modulated primarily by the major atmosphere and atmosphere-ocean climate modes (Wang and Schimel, 2003).Globally, the largest marine physical and biogeochemical signals arise in the tropical Pacific, modulated by the El Nino-Southern Oscillation (ENSO).The mechanisms of the ENSO biogeochemical response are reasonably well understood (Chavez et al., 1999;Behrenfeld et al., 2006).Regional changes in atmospheric convection and trade winds affect upwelling of subsurface water containing high dissolved inorganic carbon (DIC) and carbon dioxide partial pressure (pCO 2 ) while remotely forced Kelvin waves influence the depth of the thermocline and thus the biogeochemical concentrations of the source waters (Feely et al., 1999).ENSO-related variability also extends over much of the globe because of ocean wave propagation from the tropical Pacific and atmospheric teleconnections (Wang and Schimel, 2003).
Repeat spatial surveys of the air-sea CO 2 partial pressure difference, DpCO 2 , document a strong correlation between depressed CO 2 efflux to the atmosphere and the onset of El Nino events; observations of La Nina events show a corresponding enhancement of CO 2 efflux (Feely et al., 1999(Feely et al., , 2002)).Field-based estimates of air-sea CO 2 flux interannual variability are about 70.2 Pg C yr À1 (Feely et al., 2002).Ocean models suggest comparable, or somewhat lower, estimates of the interannual variability of 70.13 to 70.3 Pg C yr À1 , depending in part on the spatial and temporal analysis domains (Le Que ´re ´et al., 2000; Obata and Kitamura, 2003;McKinley et al., 2004;Wetzel et al., 2005).
In the extratropics three major climate modes cause interannual variability, the North Atlantic Oscillation (NAO) (Visbeck et al., 2001;Gruber et al., 2003;Bates et al., 2003), Pacific Decadal Oscillation (PDO) (Takahashi et al., 2003;Feely et al., 2006), and the Southern Annular Mode (SAM) (Lenton and Matear, 2007;Lovenduski et al., 2007;Le Que ´re ´et al., 2007;Verdy et al., 2007).All three modes involve atmospheric pressure oscillations that drive substantial changes in the strength and location of the surface winds, ocean upwelling, ocean convection patterns, seasurface temperature (SST), and air-sea heat and freshwater fluxes on regional scales.The impact on air-sea CO 2 fluxes depends upon the interaction of several, often competing, climatic factors such as thermal solubility (SSTs), biological drawdown of DIC, upwelling/mixing of nutrient-and DIC-rich waters, net surface freshwater fluxes (through dilution of DIC and alkalinity, Alk), and wind speed.
Using ocean time-series data from Bermuda, Gruber et al. (2003) suggest a correlation of negative NAO index years with deeper mixed layers, lower SSTs, increased entrainment and biological production, and enhanced CO 2 uptake.Extrapolating from Bermuda to the entire subtropical gyre leads to an interannual air-sea CO 2 flux variability of 70.2 Pg C yr À1 ; assuming that the subtropic and subpolar air-sea CO 2 flux variability is in phase, as done by Gruber et al. (2003), increases the variability of the whole North Atlantic to 70.3 Pg C yr À1 .From an extended Bermuda time-series, Bates (2007) estimates subtropical interannual air-sea CO 2 flux variability of 0.2-0.3Pg C yr À1 .However, as discussed below there is reason to question the assumptions that the North Atlantic basin as a whole or even just the subtropics vary coherently, suggesting that these extrapolations from the Bermuda record may overestimate basin-integrated interannual variability.
Model simulations indicate considerably weaker Northern Hemisphere extratropical variability (Le Que ´re ´et al., 2000, 2003).Obata and Kitamura (2003), for example, predict that the interannual variability in North Atlantic air-sea CO 2 flux is only 70.04 Pg C yr À1 in the subtropics and 70.03 Pg C yr À1 in the subpolar gyre.McKinley et al. (2006) analyze tropical and North Pacific CO 2 variability from seven biogeochemical models.They report relatively weak interannual variability in North Pacific air-sea CO 2 flux of 70.03 to 70.11 Pg C yr À1 for the basin, reflecting in part cancellation of out-of-phase regional air-sea flux anomalies.Although the model air-sea CO 2 fluxes are correlated with the Pacific Decadal Oscillation, the total projection on the PDO is small because PDO-generated surface temperature, DIC and alkalinity anomalies have opposing effects on surface pCO 2 .
The Southern Annular Mode is the dominant climate mode over the Southern Ocean and is approximately zonally symmetric and circumpolar about the Antarctic continent.A positive SAM phase is associated with increased surface wind stress, enhanced upwelling of carbon rich, subsurface circumpolar deep water, and efflux of ocean CO 2 to the atmosphere (Lenton and Matear, 2007;Lovenduski et al., 2007;Le Que ´re ´et al., 2007).Simulated air-sea CO 2 flux variability in the Southern Ocean ranges from 70.07 (Obata and Kitamura, 2003) to 70.19 Pg C yr À1 (Lovenduski et al., 2007).
Historical ocean carbon data are too sparse, except for a few regions, to fully resolve upper-ocean carbon system and air-sea CO 2 flux variability on the required regional and monthly scales (Bender et al., 2002).This will likely remain true in the near-term on a global-scale, even with the recent growth in instrumented biogeochemical moorings and volunteer observing ship (VOS) pCO 2 transects (e.g., Metzl et al., 2007;Doney et al., 2009b).Further, while critical, the data from observational networks often cannot differentiate clearly among different potential factors driving the observed variability.
Numerical simulations provide important tools to aid the interpretation of observational data, extrapolate/interpolate airsea CO 2 fluxes to the required regional and monthly scales in a dynamically consistent manner, and characterize the underlying physical-biogeochemical processes (e.g., Le Que ´re ´et al., 2000; Obata and Kitamura, 2003;McKinley et al., 2004McKinley et al., , 2006;;Wetzel et al., 2005).Here we present a globally consistent analysis of upper-ocean biogeochemical interannual variability from a numerical hindcast  that exhibits good skill relative to observations.The analysis includes both air-sea CO 2 fluxes and surface-water chemistry.We also quantify and partition the underlying forcing factors including atmospheric physical forcing, dust deposition, and ocean circulation and biology.

Model formulation
The Community Climate System Model (CCSM-3) ocean Biogeochemical Elemental Cycle (BEC) model consists of an upper-ocean ecological module (Moore et al., 2004) and a fulldepth ocean biogeochemistry module (Doney et al., 2006) both embedded in a three-dimensional (3-D) global physical ocean general circulation model (Smith and Gent, 2004;Collins et al., 2006).The model is forced with physical climate forcing from atmospheric reanalysis and satellite data products (Doney et al., 1998(Doney et al., , 2007) ) and time-varying dust deposition (Mahowald et al., 2003).The CCSM BEC model is cast as a set of 3-D, time-varying advection diffusion equations for a suite of tracers K: where the second and third terms on the left-hand side are the advective and diffusive divergence terms, respectively, and r is the 3-D del operator.The physical transport is partitioned into resolved advection and parameterized eddy mixing terms; all of the ecological-biogeochemical source/sink terms and surface and sediment fluxes are grouped into the right-hand-side term RHS bio .

Ecosystem-biogeochemistry modules
The ecosystem module builds on traditional phytoplanktonzooplankton-detritus-nutrient food web models and incorporates multi-nutrient limitation (N, P, Si, Fe) on phytoplankton growth and specific phytoplankton functional groups (Moore et al., 2004).The ecosystem module is coupled with an ocean biogeochemistry module with full carbonate system thermodynamics and air-sea CO 2 and O 2 fluxes (Doney et al., 2004(Doney et al., , 2006)), nitrogen fixation and denitrification (Moore and Doney, 2007), and a dynamic iron cycle with atmospheric dust deposition, water-column scavenging and a continental sediment source.There are 14 main compartments: pico/nano-plankton, diatoms, and diazotrophs; zooplankton; suspended and sinking particulate detritus; and dissolved nitrate, ammonia, phosphorus, iron, silicate, oxygen, dissolved inorganic carbon (DIC), and alkalinity (Alk).
The pico/nano-plankton size class is designed to replicate the rapid and highly efficient nutrient recycling found in many subtropical, oligotrophic (low nutrient) environments.Diatoms model a larger, bloom-forming size class.Phytoplankton growth rates are determined by available light and nutrients using a dynamic growth model.Photoadaptation is parameterized with dynamically adaptive chl/C ratios.The diazotrophs fix all required nitrogen from N 2 gas, and calcification is parameterized as a fraction of the pico/nano-plankton production as a function of temperature and nutrients adapted for coccolithophores.The model has one adaptive zooplankton class that grazes on phytoplankton and zooplankton.Size-structure effects are included by varying key zooplankton (e.g., partitioning of fecal pellets between suspended and sinking detritus) depending on the food source.Many of the biotic and detrital compartments contain multiple elemental pools, in addition to carbon, to track flows through the ecosystem.
The biogeochemistry module (Doney et al., 2006) is based on an expanded version of the Ocean Carbon Model Intercomparison Project (OCMIP) biotic model (Najjar et al., 2007).The module includes full carbonate system thermodynamics and air-sea CO 2 and O 2 fluxes.A dynamic iron cycle is incorporated with atmospheric dust deposition, water-column scavenging, and a continental sediment source (Moore et al., 2006).Denitrification is simulated in oxygen minimum zones following Moore and Doney (2007), and subsurface particle remineralization is parameterized incorporating the mineral ballast arguments of Armstrong et al. (2001).
The model equations are identical to those reported for the 3-D implementation of Moore et al. (2004) with two important modifications as documented in more detail in Moore et al. (2006).First, water-column denitrification has been added to the model in order to close the global nitrogen cycle.Second, a number of the parameters associated with the model iron dynamics and scavenging have been adjusted to improve the simulated dissolved iron fields (see Table 1 of Moore et al., 2006).The ecosystem-biogeochemistry simulation is in generally good agreement with bulk metrics (e.g., total biomass; productivity; nutrients; export; pCO 2 , air-sea CO 2 flux) across diverse regions that include both macro-nutrient-and iron-limited conditions (Moore et al., 2004;Doney et al., 2009a).

Atmospheric dust deposition
Time-varying mineral aerosol deposition to the ocean is simulated using a 3-D atmospheric chemical transport model (Mahowald et al., 2003;Luo et al., 2003) based on National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis (Kistler et al., 2001).The dust source and deposition scheme is based on the Dust Entrainment and Deposition (DEAD) scheme (Zender et al., 2003), and the chemical transport model is the Model of Atmospheric Transport and Chemistry (MATCH) (Rasch et al., 1997), which has been developed specifically to be used with reanalysis winds (Mahowald et al., 1997).The dust source areas are defined as dry, poorly vegetated regions that have easily erodible sources, using topographic lows as preferential sources areas (Ginoux et al., 2001).Dust is removed by wet deposition during precipitation events, and by dry deposition from gravitational settling and turbulent processes.
The MATCH/DEAD model captures fairly well the observed mean seasonal cycle and distribution of the aerosols (e.g., in-situ concentration, in-situ optical depth, deposition data and satellite optical depth), although with somewhat reduced skill in the Southern Hemisphere (Luo et al., 2003).Some skill is also shown in capturing interannual variability in dust concentrations downwind of the source regions (Mahowald et al., 2003), most of which is driven by atmospheric transport (or transport/source correlations) and not by source interannual variability alone (Tegen and Miller, 1998;Mahowald et al., 2003).Model monthly averages exhibit statistically significant correlations against observations at the few available ($10) in-situ dust concentration stations.For ocean regions where mineral aerosols dominate the aerosol optical depth, the model appears to be able to capture much, but not all, of the ocean interannual variability.

Atmospheric physical forcing and ocean hindcasts
The Parallel Ocean Program (POP) is a z-level, hydrostatic, primitive equation model integrated here with a resolution of 3.61 in longitude, 0.8-1.81 in latitude, and 25 vertical levels (Yeager et al., 2006).The effects of mesoscale eddy transport are parameterized according to Gent and McWilliams (1990).The Large et al. (1994) K-Profile Parameterization is implemented in the vertical to capture surface boundary-layer dynamics and interior diapycnal mixing.The historical simulation ) (Doney et al., 2009a) is integrated with air-sea heat, freshwater, and momentum fluxes derived from a bulk flux forcing method that combines 6-h atmospheric surface fields (temperature, humidity, winds) from the NCEP reanalysis (Kistler et al., 2001) with satellite-and in-situ-derived clouds, precipitation, runoff and sea-ice fraction (Large and Yeager, 2004).The ocean physical solutions replicate skillfully observed interannual variability of temperature, sea-surface height, and circulation, as demonstrated in Doney et al. (2007) for an earlier variant of the physical model component.
The initial conditions for the nutrient and inorganic carbon variables are prescribed from data-based climatologies (e.g., Key et al., 2004).The ecological-biogeochemical simulation is spunup for several hundred years, prior to initiating the interannual varying forcing, using a repeat annual cycle of physical forcing and dust deposition.Model ecosystem components converge within a few years.There remains a slow drift in the subsurface nutrient and inorganic carbonate fields; the global net air-sea CO 2 uptake flux is 0.150 Pg C yr À1 (mean areal flux 0.025 mol C m À2 yr À1 ) but the change in the drift over the 26-yr integration ) is only À0.010 Pg C yr (mean À0.002 mol C m À2 yr À1 ), much smaller than the simulated interannual variability.We focus our analysis primarily on a ''pre-industrial'' model simulation with fixed atmospheric CO 2 mole fraction (280 ppmv); in a companion ''anthropogenic CO 2 '' simulation, atmospheric CO 2 is increased following observations beginning in model year 1870.The full interannual variability in physics and dust forcing is initiated in model calander year 1979.

Model analysis
We focus our analysis on three key carbon metrics, air-sea CO 2 flux F CO 2 , surface-water CO 2 partial pressure pCO 2 , and upperocean dissolved inorganic carbon inventory I DIC .Surface-water pCO 2 is computed from model prognostic variables DIC, Alk, temperature (T) and salinity (S) using a full seawater carbonate thermodynamics code.The air-sea CO 2 flux: depends on CO 2 solubility a and a gas transfer coefficient k derived using a quadratic wind-speed formulation (Wanninkhof, 1992) and NCEP 10-m winds.In regions of sea-ice, the transfer velocity k is multiplied by (1Àf), where f is the fractional sea-ice cover (0À1).The climatological seasonal patterns compare reasonably well against observations for surface-water pCO 2 (correlation 0.78; rms error 18.6 matm) and air-sea CO 2 flux (correlation 0.61; rms error 1.53 mol CO 2 m À2 yr À1 ) (Fig. 1; see also Doney et al., 2009a and Gruber et al., in press).To explore interannual variability signals, we construct spatial maps of the simulated monthly anomalies, denoted with a prime as Y 0 ¼ Y À ȲM , by removing a model mean monthly climatology ȲM (i.e.mean January, mean February, etc.) computed at each grid point.
To diagnose the mechanisms governing interannual biogeochemical variability, we construct linear Taylor expansions of the F 0 CO 2 , pCO 0 2 , and I 0 CO 2 anomalies at each grid point in terms of the time-varying component of the underlying driving factors X 0 i : where qY=qX i is the partial derivative with respect to forcing factor i holding all other forcing factors constant.For two forcing components X i and X j , the expansion resembles the Reynold's decomposition for turbulent flow: where the final term includes second-order cross-terms.For air-sea CO 2 flux anomaly (Eq.( 2)), the corresponding decomposition (Eq.( 4)) is where the right-hand-side terms (RHS) reflect, respectively, the contributions to the variability in air-sea flux from gas transfer/ solubility anomalies, DpCO 2 anomalies, and CO 2 and gas transfer/ solubility anomaly cross-term corrected for the, typically small, mean cross-term of the anomalies.The temperature dependence of k and a approximately cancel, and the (ka) 0 term primarily reflects variations in wind speed and sea-ice fraction.

ARTICLE IN PRESS
The interannual surface-water anomalies pCO 0 2 also can be expanded linearly into thermal, freshwater, and circulation/ biology components (Takahashi et al., 1993;Lovenduski et al., 2007): where the second-order terms are neglected.Precipitation and evaporation of freshwater drive correlated variations in surfacewater DIC 0 and Alk 0 , which have opposing effects on pCO 2 .To remove this effect we use a modified form of Eq. ( 6) (Lovenduski et al., 2007): where nDIC and nAlk are salinity normalized, and the partial derivative with respect to salinity includes the effects of freshwater changes on DIC and Alk concentrations.Similar linear pCO 2 decompositions are used, for example, by Wetzel et al. (2005).As the size of the forcing anomalies X 0 grows, deviations from the linear approximations in Eqs. ( 6) and ( 7) begin to appear because of non-linearities in carbonate system thermodynamics.To address this, we compute numerically the partial derivatives in Eq. ( 7) for the full range of the anomaly: where the pCO 2 values are calculated with full non-linear seawater carbonate thermodynamics.
In our ocean hindcast simulations, interannual variability in the DpCO 0 2 term is controlled primarily by surface-seawater conditions because the atmospheric mole fraction xCO atm 2 is held fixed globally.The only variations in the pCO atm 2 term (Eq.( 2)) in the model are driven by factional change fluctuations in total atmospheric pressure: However, the interannual variability in atmospheric xCO atm 2 over the ocean is small relative to surface-water variations, with rmsðxCO atm 2 Þo0:5 À 1:0 ppmv (Nevison et al., 2008).As shown in the results (Section 4.2) below, temporal variations in surface-water nDIC 0 play an important role in interannual pCO 0 2 anomalies (Eqs.( 6) and ( 7)).We therefore conduct a further decomposition on the model DIC field to explore the contributions of different physical, chemical and biological processes to variations in the water-column inventory of upperocean (0-100 m) dissolved inorganic carbon (mol C m À2 ): Similar to an analysis of upper-ocean heat content presented in Doney et al. (2007), we partition the net annual change of the DIC inventory anomaly into specific components: where each of the RHS terms are interannual anomalies integrated in time over a year Dt and, if appropriate, integrated over depth (0-100 m).Annual, rather than monthly, time-steps are used to because the CCSM-3 BEC model diagnostics are stored for monthly means, which introduces some residual noise in constructing DI 0 DIC on monthly time-scales because of interpolation issues.The RHS terms are given as follows: the air-sea CO 2 flux F CO 2 ; the surface virtual flux of DIC in the volume conserving model due to freshwater fluxes V DIC ; the vertical integral of net biological release of inorganic carbon B DIC ; the vertical integral of the convergence of the resolved advective DIC transport: where ṽres is the resolved model velocity; and the vertical integral of the convergence of the eddy-parameterized DIC: where ṽbolus is the bolus velocity from the mesoscale eddy parameterization, and k is the parameterized model diffusivity.
All terms are defined such that positive values result in a net increase in the upper-ocean DIC inventory.Following Doney et al. (2007), we also adopt a novel approach to display the results of our analysis for the case when we have more than two forcing components.Multivariate regressions and simple anomaly rms plots, by themselves, fall short because of non-trivial correlations between the different forcing terms.To include the interconnection among the different terms, we examine instead the slopes b from linearly regressing each RHS term for pCO 0 2 (Eq.( 7)) and DI 0 DIC (Eq.( 11)) individually onto the temporal anomalies at each grid point.For the general case: Note that the intercept is approximately zero because the average of the anomalies is zero.
Our treatment differs somewhat from traditional multiple linear regression where the focus is often deconvolving the effect of cross-correlated forcing variables on the predicted variable.This is often done by estimating unknown partial regression slopes and correlation coefficients for individual forcing variables on the predicted variable, holding all other forcing variables constant.For the hindcast model solutions, the partial regression slopes qY=qX i for the RHS terms are already known: from the model gas exchange equations for F 0 CO 2 (Eq.( 2)); from carbonate thermodynamics for pCO 0 2 (Eq.( 7)); and by construction (and equal to 1.0) for the DI 0 DIC rate of change (Eq.( 11)).Our objective, in contrast, is to highlight the relative dominance of terms compared to overall variance in Y 0 , for which we need to include the interconnection or cross-correlation among the different forcing or budget terms.The regression slopes b i (Eq.( 14)) reflect this and differ from 1.0 either because the variance of X 0 variable is substantially smaller than other forcing terms or because X 0 is out of phase or anti-correlated with more dominant terms.
A slope b i near 1 indicates that a particular term produces in phase anomalies of comparable magnitude to Y 0 .A slope greater than 1 indicates that the term is producing even larger anomalies compensated by other terms.A negative slope indicates such a compensating term, and the more negative the slope the more effective the compensation.A slope near zero indicates that the term is not important in generating the flux anomalies.To highlight the regions with substantial variability, we normalize the slopes with the total rms of anomalies from a model monthly climatology, computed at each grid cell, by displaying spatial maps of b i Â rms(Y 0 ) for each of the RHS terms in the various linear expansions.

Air-sea CO 2 flux variability
Fig. 2 illustrates a spatial map of the monthly anomalies in air-sea CO 2 flux F 0 CO 2 from a particular month in the model hindcast.The amplitude of typical anomalies are on the order of 74.5 mol C m À2 yr À1 in high variability regions in the Northern Hemisphere mid-to high latitudes, tropical Pacific, and Southern Ocean.The air-sea flux anomalies tend to be elongated along the axis of major current systems and have spatial extents of 1000-5000 km.The coarse horizontal resolution of the simulation does not allow for the formation of mesoscale eddies and intense coastal upwelling regimes, and in the real ocean mesoscale variability would layer on top of the climate-scale variability patterns in Fig. 2. Within a particular ocean region, large anomalies of opposite sign occur commonly, sometimes in dipole patterns.This can result in cancellation and damping of the variability signal when the fluxes are averaged over sub-basin and basin scales.This is especially true in the extratropics, explaining in part the relatively small interannual variability in regionally integrated air-sea CO 2 flux found in the CCSM simulation and other models; it also argues strongly against the extrapolation of point time-series measurements to gyre or basin scale.
A spatial map of the temporal, root-mean-squared (rms) variability in the model interannual air-sea CO 2 flux anomalies rmsðF 0 CO 2 Þ is displayed in Fig. 3 (top left panel).The local (gridpoint) values of rmsðF 0 CO 2 Þ range from 0.1 to 0.2 mol C m À2 yr À1 in parts of the subtropics to 0.8-1.8mol C m À2 yr À1 in temperate to high latitudes.Maxima in rmsðF 0 CO 2 Þ are found in the tropical Pacific, northern North Atlantic and North Pacific, and Southern Ocean.
Time-series of spatially averaged, monthly air-sea CO 2 flux anomalies hF 0 CO 2 i reg are shown in Fig. 4 for different sub-basin regions.The sub-plots are displayed uniformly in terms of flux per unit area (mol m À2 yr À1 ) and on a common y-axis scale to facilitate inter-comparison, independent of the size of the regions.The amplitude of the model temporal rms variability of the regional, spatially averaged time-series rmsðhF 0 CO 2 i reg Þ varies from 0.08 to 0.12 mol C m À2 yr À1 in low-variability regions (e.g., North Indian Ocean, north subtropical Atlantic Ocean) to as high as 0.28-0.43mol C m À2 yr À1 in high-variability regions (e.g., south Southern Ocean and western and eastern Equatorial Pacific).The low variability regions have short temporal decorrelation timescales, on the order a few months, while the correlation scales are 1-2 yr in the south Southern Ocean to about 5 yr in the Equatorial Pacific.The rms temporal variability in the globally integrated, monthly air-sea CO 2 flux anomaly rmsð R globe F 0 CO 2 Þ is 0.335 Pg C yr À1 , similar to that found by previous regional and global model studies (e.g., Le Que ´re ´et al., 2000; Obata and Kitamura, 2003;McKinley et al., 2004McKinley et al., , 2006;;Wetzel et al., 2005).The main contributions to this variability come from the tropical Pacific (0.181 Pg C yr À1 ), northern North Atlantic (0.035 Pg C yr À1 ) and North Pacific (0.050 Pg C yr À1 ), and Southern Ocean (0.190 Pg C yr À1 ), where the numerical values represent rmsð R reg F 0 CO 2 Þ.The large, ENSO-generated simulated interannual variability averaged over the equatorial Pacific hF 0 CO 2 i reg (Fig. 4) is similar in both magnitude and phasing to that found in a recent data-based synthesis (Feely et al., 2006) (Fig. 5).The correlation coefficient for the monthly model and observed values is 0.55 for 1982-2004 inclusive.The correlation coefficient increases to 0.64 when the data and simulation are smoothed to remove subannual variability and the model is shifted back in time by a few months.The model-data agreement for the equatorial air-sea CO 2 flux variability compares favorably against previous simulations.

ARTICLE IN PRESS
Historical data coverage is considerably sparser in time and space for the other, extratropical high variability regions, making the verification of the model predicted variability signal more difficult.The interannual air-sea CO 2 flux variability in the CCSM-3 ocean BEC model is at the upper range reported for other models for the tropical Pacific Ocean models (70.13 to 70.3 Pg C yr À1 ) and about twice that estimated for the Southern Ocean (70.1 Pg C yr À1 ) (Le Que ´re ´et al., 2000;Obata and Kitamura, 2003;McKinley et al., 2004;Wetzel et al., 2005).

Diagnosing ocean carbon variability mechanisms
The mechanisms driving the interannual variability in simulated air-sea CO 2 flux are displayed in the other panels of Fig. 3 following the linear decomposition of F 0 CO 2 (Eq.( 5)).Model variability is dominated in most regions of substantial rmsðF 0 CO2 Þ by surface-water pCO 2 anomalies modulated by mean transfer velocity ðkaÞDpCO 0 2 .Noticeable contributions from gas transfer variability ðkaÞ 0 DpCO 2 occur in the tropical Pacific and extratropical storm tracks, regions with both large wind speed variability and large mean air-sea pCO 2 differences.Gas transfer velocity-driven flux variability also arises in high-latitude, seasonal sea-ice zones, but in this case the mechanism is through variations in ice cover fraction f rather than wind speed.The anomaly cross-term ðkaÞ 0 DpCO 0 2 À ðkaÞ 0 DpCO 0 2 is minor.The ocean processes governing variability in surface-water pCO 2 are examined in Fig. 6 using the linear decomposition of pCO 0 2 in Eq. 7. The values of rms(pCO 0 2 ) are 2.5-7.5 matm in the subtropics increasing to 10-30 matm in the tropical Pacific and high latitudes.

ARTICLE IN PRESS
Variations in surface-water pCO 2 arise from different mechanisms in different regions.Examples of the linear decomposition are shown in Fig. 7 for two grid points, one in the tropical Pacific and one in midlatitude North Pacific.The right column displays the model time series of pCO 0 2 anomalies ðqpCO 2 =qXÞX 0 from forcing components X 0 due to S 0 FW , T 0 , nAlk 0 , and nDIC 0 (Eq.( 7)); the left column shows property-property plots of X 0 i qpCO 2 =qX i (y-axis) versus the total pCO 0 2 anomaly (x-axis).Linear regression slopes b i (Eq.( 14)): are calculated for each forcing component.
The analysis at the North Pacific site is straightforward to understand (Fig. 7, bottom panel).The largest pCO 0 2 anomalies are created by variations in nDIC 0 , and the slope b nDIC is positive and slightly greater than 1.0 indicating that nDIC 0 is the dominant forcing factor.The pCO 0 2 anomalies from T 0 are smaller in amplitude and anti-correlated with those from nDIC 0 ; the corresponding slope is negative, showing that the temperaturedriven pCO 0 2 anomalies partially damp the pCO 0 2 anomalies due to inorganic carbon variations.The analysis is consistent with variations in the supply of cold (lower pCO 0 2 ), DIC-rich (higher pCO 0 2 ) subsurface water.The nAlk 0 -and S 0 FW -driven pCO 0 2 anomalies are small.Note that the anomalies in the property-property plot are not always perfectly linear and suggest a somewhat different slope for positive versus negative pCO 0 2 anomalies.Such behavior could arise from interdependencies (cross-correlations) with other thermodynamic factors driving pCO 0 2 anomalies.At the tropical Pacific site, substantial variations are generated by interannual anomalies in all four forcing components (Fig. 7, top panel).Again nDIC 0 -driven pCO 0 2 anomalies dominate the total pCO 0 2 variability with partial temperature compensation, at least for the large positive pCO 0 2 anomalies associated with La Nina conditions.The pCO 0 2 anomalies from nAlk 0 and S 0 FW are smaller in amplitude and positively correlated with those from nDIC 0 ; the corresponding slopes are smaller but still positive.The analysis is consistent with correlated upwelling and surface freshwater signals modulated by ENSO.

ARTICLE IN PRESS
Spatial maps are shown in Fig. 6 of the contributions of different terms to rms(pCO 0 2 ).The spatial patterns are generally consistent with the analysis presented for the two example sites in Fig. 7. Thermally generated interannual pCO 0 2 anomalies dominate in the subtropics and tropical Atlantic.In tropical Pacific, interannual variability in surface-water pCO 2 primarily reflects anomalies in surface-water nDIC 0 ; thermal pCO 0 2 anomalies partially cancel the nDIC 0 -driven variations.Salinity (surface freshwater flux) variations contribute to pCO 0 2 , mainly through concentrating or diluting surface-water DIC and Alk, in the western tropical Pacific and near the mouths of the major river systems (e.g., Amazon, Orinoco, Zaire, Ganges/Bhramaputra).The contribution to rms(pCO 0 2 ) from surface-water nAlk 0 anomalies is small and not shown.
At present, there is no interannual variability in river inputs in the ocean-only version of the model, and therefore the variability reflects internal ocean circulation acting on the river inflow.In the CCSM-3 BEC model, rivers input pure freshwater and thus act like precipitation; their contribution is correctly diagnosed using Eq. ( 7).However, for a more realistic model treatment of river geochemistry, the analysis framework would have to be modified for these near-coastal regions to account for the non-zero DIC and Alk end-member concentrations from different rivers.
A spatial map is displayed in Fig. 8 (top left panel) of the interannual variations in the annual rate of change upper-ocean DIC water-column inventory rms(DI 0 DIC ).The spatial patterns are similar to those of the surface-water nDIC 0 contributions to rms(pCO 0 2 ), with maxima of 1-5 mol C m À2 yr À1 in the Indian and Pacific tropics and 1-2 mol C m À2 yr À1 in temperate and polar regions; the subtropics show considerably weaker variability.
The variability in surface-water DIC and I DIC are not identical in regions of shallow mixing such as in the tropics, where vertical heaving of shallow isopycnals combined with sharp vertical gradients can induce variability in water-column inventories that is more weakly expressed at the surface (Doney et al., 2007).But an analysis of only the model surface layer folds in additional complications because of strong vertical mixing in the mixed layer and fact that subsurface biological carbon uptake/release can be transported by vertical mixing/advection into the surface layer.An integration depth of 100 m was chosen as a compromise, as this is below the annual mean mixed layer in the model for a large fraction of the globe and is the approximate depth of net biological carbon uptake in the model (e.g., Jin et al., 2007).

ARTICLE IN PRESS
Using the linear expansion in Eq. ( 11), we diagnose the contributions to rms(DI 0 DIC ) from biological and physical circulation terms (horizontal and vertical divergence, surface flux, biological uptake and remineralization) (Fig. 8).The strongest contributors are the convergence terms from DIC transport due to resolved advection A 0 DIC and eddy-parameterized transport and mixing E 0 DIC .The two terms are spatially distinct to a degree, with the resolved advection term dominating in the tropics and high latitudes and the eddy-term contributing most in regions of steep isopyncal depth variations (e.g., western boundary currents, edges of the Antarctic Circumpolar Current) and in deep winter mixing zones.Air-sea CO 2 flux variations have a relatively small impact on interannual water-column DIC anomalies, at least on an annual time-scale and comparing to an inventory integrated to 100 m.Most of the variability is internal to the ocean involving either lateral or vertical reorganizations in the distribution of DIC.The virtual flux of DIC does contribute to variability in upper-ocean DIC inventory in the tropics associated with ENSO-driven precipitation anomalies.

ARTICLE IN PRESS
Net biological uptake is out of phase with DIC carbon inventory changes in the tropics.This anti-correlation can occur by the following example.Consider an ENSO-generated increase in subsurface upwelling (La Nina event) increases both upper-ocean DIC and biological productivity; the magnitude of the negative anomaly from larger export production anomaly B 0 o0 is smaller than the positive input from the circulation, resulting in a net positive increase in DI 0 DIC 40.The impact of biological uptake is small, relative to other terms, over most of the rest of the ocean.The regions where there is a substantial positive contribution in the Southern Ocean appear to be dominated by atmospheric dust-driven variability, which is decoupled with the circulation-driven variations in subsurface DIC and nutrient supply (see Section 4.3).

Dust-driven interannual variability in ocean carbon
The interannual variability in model air-sea CO 2 flux and ocean carbon system dynamics discussed to this point reflects the combined biogeochemical responses to time variations in both physical climate and dust deposition forcing.To partition these two effects, we conducted sensitivity experiments using interannually varying physics and repeat seasonal dust deposition (Fig. 9, middle and bottom rows, left column) and repeat seasonal physics and interannually varying dust deposition (Fig. 9, middle and bottom rows, right column).Since dust variations do not alter gas transfer rates, sea-surface temperature, or freshwater, the dust-driven rmsðF 0 CO 2 Þ patterns are isolated to variations in DpCO 0 2 from altered biological uptake of surface-water nDIC 0 and nAlk 0 anomalies (middle row, Fig. 9).The rms variability from physical climate-driven and dust-driven export production (sinking particle flux) anomalies rms(EP 0 ) are shown in the bottom row of Fig. 9.
Atmospheric dust deposition to the ocean surface varies by several orders of magnitude between areas of high deposition downwind of desert source regions and the most remote locations with very low deposition (Fig. 9, upper left panel).Interannual variability in the deseasonalized dust deposition anomalies also varies dramatically and scales to a degree with deposition rates.Therefore, we normalize dust deposition variability by the annual mean deposition, rmsðF 0 dust Þ=F dust (i.e.coefficient of variation; Fig. 9 upper right panel).The coefficient of variation for dust deposition has values of 0.5-1.0 over much of the globe.But it exhibits maxima with values 41 in some areas of moderate to low deposition, often in areas with steep spatial gradients in annual mean deposition.
The strongest impacts of dust deposition variability on F 0 CO 2 , net biological uptake B 0 , and export production EP 0 are in highnitrate, low chlorophyll regions in the western North Pacific, Equatorial Pacific, and Southern Ocean downwind of dust source regions (Fig. 9, right column).In those regions, dust-driven rmsðF 0 CO 2 Þ and rms(EP 0 ) values are comparable to those generated by physical climate variability.Using similarly formulated numerical experiments, Aumont et al. (2008) recently reported substantial dust-induced surface ocean chlorophyll variability in approximately the same areas.Dust and physical climate forcing are almost entirely independent on a local scale in our simulations; although dust sources on land and transport in the atmosphere reflect regional climate variability modes (e.g., ENSO) that also effect ocean circulation and mixing, there is not a strong correlation between either air-sea flux or export production between the climate-only and dust-only simulations (not shown).
On a global scale, the temporal variability in spatially integrated air-sea CO 2 flux rmsð R globe F 0 CO 2 Þ is 0.29 Pg C yr À1 for physical climate alone and 0.14 Pg C yr À1 for dust deposition alone.The corresponding values of rmsð R globe EP 0 Þ are 0.19 Pg C yr À1 for physical climate alone and 0.18 Pg C yr À1 for dust deposition alone.The large dust-driven interannual variability levels arises from a global shift toward lower dust input, reduced export production, and net efflux of CO 2 to the atmosphere (F 0 CO 2 o0) in the late 1990s.The air-sea CO 2 flux and export production anomalies in the late 1990s are approximately À0.2 Pg C yr À1 .The ratio of 1:1 is in agreement with other modeling studies with the CCSM modeling showing that iron-induced changes in ocean biology are effective at altering net ocean carbon storage at least on decadal scales (Moore et al., 2006;Jin et al., 2008).
The dust deposition anomalies in the Mahowald et al. ( 2003) simulation exhibit relatively short spatial correlations, and the sharp reductions in the global integrals reflect the superposition of somewhat more variable, and not fully coherent regional patterns.These are exhibited in regional time-series of deseasonalized dust deposition and ocean biogeochemical anomalies shown in Fig. 11 for the main regions of dust-induced variability (western and eastern Southern Ocean, subtropical South Pacific, North Pacific, and eastern equatorial Pacific).The dust model predicts trends toward reduced dust deposition in the 1990s for all five regions.
The resulting air-sea CO 2 efflux anomalies in the Southern Ocean are consistent with positive net CO 2 efflux anomalies of 0.10-0.15Pg C yr À1 derived from atmospheric CO 2 inversions and ocean models (Le Que ´re ´, et al., 2007;Lovenduski et al., 2007Lovenduski et al., , 2008)).These anomalies have been ascribed to increased upwelling of CO 2 -rich subsurface waters because of a more positive Southern Annular Mode.Our simulations suggest that dust deposition also may be playing a role.
In the Southern Hemisphere, the dust variations largely come from the South American source.There are not long term aerosol records available for these regions that would show the relatively small variations in aerosol amounts from dust (compared to sea salts) or in-situ datasets covering this time period.The only data we have is from in-situ visibility datasets from meteorological stations.These data are consistent with higher dust during the early 1990s, falling off in concentration at the end of the 1990s (Fig. 16, Mahowald et al., 2007), and individual visibility plots of WMO stations 847820, 852300 (near Atacama desert) and 876230 (in Argentina-other stations in this region have too few visibility events to see trends).This timing of a reduction in dust in the source areas is consistent with the late 1990s reduction in modeled dust deposition to the oceans and resulting reduction in CO 2 fluxes, which can be seen in the global average (Fig. 10).This is also consistent with the increased precipitation seen in the sources regions during the late 1990s (Mahowald et al., 2007, Fig. 16).Some caveats should be kept in mind regarding the potential magnitude of atmospheric dust-induced variability.First, the atmospheric dust model predicts too high dust concentrations (and therefore likely predicts too high deposition rates) in the Southern Hemisphere (Wagener et al., 2008), although most of this will be offset by the higher than average iron solubilities seen in the Southern Ocean (Baker et al., 2006).Second, while the CCSM simulations reported here include a continental shelf iron source, some studies argue for a larger shelf source (Moore and Braucher, 2008); this would tend to reduce the sensitivity of the ocean model biogeochemistry to dust deposition variations.

Summary and conclusions
Dynamical systems are characterized by their equilibrium states (both stable and unstable), internally driven oscillations, limit cycles and chaotic behavior, and transient responses to external perturbations.In the case of ocean ecology and biogeochemistry, natural atmosphere-ocean climate modes stimulate interannual variability in ocean physics, which in turn drives regional variations in the ocean carbon system.Using a numerical hindcast of ocean biogeochemistry, we present a systematic global analysis of the magnitude and processes governing carbon system variability on subannual to decadal time-scales.The deseasonalized anomalies in air-sea CO 2 flux, surface-water pCO 2 , and upper-ocean inorganic carbon inventory are decomposed into the underlying processes using a linear Taylor expansion and partial derivatives of the variable of interest to individual forcing terms.Specifically, we diagnose the atmospheric, biological and ocean physical circulation terms (horizontal and vertical divergence, surface flux, biological uptake and remineralization) contributions.
Based on our model analysis, the mechanisms governing interannual variability in the upper-ocean carbon system and air-sea CO 2 flux differ markedly with region.As has been identified previously, the major regions of variability occur in the Southern Ocean, tropical Indo-Pacific, and Northern Hemisphere temperate and subpolar latitudes.Ocean circulation is the dominant underlying factor driving biogeochemical variability over much of ocean.The physical circulation acts by modulating the surface concentration of salinity normalized dissolved inorganic carbon, which in turn effects surface-water pCO 2 and air-sea CO 2 flux.Gas-transfer-driven variability in the air-sea CO 2 flux occurs in the equatorial Pacific, extratropical storm tracks and the seasonal sea-ice zones.
Biological export and thermal solubility effects on pCO 2 act to damp partially the circulation-driven biogeochemical variability in the tropics.In the subtropics, thermal solubility contributes in a positive sense to surface-water pCO 2 and air-sea CO 2 flux variability.Net biological carbon uptake and export are important factors in upper-ocean carbon system variability in the temperate Southern Ocean, a combination of ocean physics and atmospheric dust driving.The model interannual surface-water pCO 2 variability patterns with thermal forcing in the subtropics and circulation/biology forcing in the tropics and subpolar gyres mirror those found by Takahashi et al. (2002) for the surface ocean pCO 2 seasonal cycle.Net freshwater inputs from precipitation/ evaporation and rivers induce variability in the tropical Indo-Pacific and near major river mouths, respectively.
Variations in atmospheric dust deposition generate substantial variability in ocean export production and air-sea CO 2 flux in HNLC in the Southern Ocean, equatorial Pacific and subpolar North Pacific.Globally, dust deposition variations induce air-sea CO 2 flux anomalies of about 0.2 Pg C yr À1 .Beginning in the mid-1990s, reduced global dust deposition generates increased air-sea CO 2 outgassing in the Southern Ocean, consistent with trends derived from atmospheric CO 2 inversions that have previously been ascribed to changes in ocean circulation.

Fig. 1 .
Fig.1.Annual-mean air-sea CO 2 flux (mol m À2 yr À1 ) in the CCSM ocean model (top panel) and an in-situ, data-based climatology(Takahashi et al., 2009).Positive fluxes are into the ocean, and negative fluxes are out of the ocean.

Fig. 2 .
Fig. 2. Monthly air-sea CO 2 flux anomalies (mol m À2 yr À1 ) from a particular month (January, 1995) in the CCSM-3 ocean BEC model.Positive fluxes are into the ocean, and negative fluxes are out of the ocean.

Fig. 3 .
Fig. 3. Partitioning of the mechanisms driving interannual variability in air-sea CO 2 flux (mol m À2 yr À1 ) in the CCSM ocean BEC model.The panels show the root mean square (rms) of the model deseasonalized CO 2 flux anomalies (1979-2004) (top left) and the contributions from gas transfer velocity (wind speed and ice cover) (top right), surface-water DpCO 2 (lower left), and the cross-correlation of gas transfer and pCO 2 anomalies (lower right).The black contour lines mark the seasonal maximum and minimum sea ice extent.

Fig. 7 .
Fig. 7. Factors driving interannual variability in monthly air-sea pCO 2 anomalies (matm) in the CCSM ocean BEC model at two locations, equatorial Pacific (top row) and mid-latitude North Pacific (bottom row).The right column displays the contributions of different processes: temperature (green), salinity/freshwater (black), salinity normalized dissolved inorganic carbon (blue) and alkalinity (red) to monthly DpCO 2 anomalies versus time.The left panel shows regressions of the individual forcing factor contributions (y-axis) against the full pCO 2 anomalies.

Fig. 8 .
Fig. 8. Partitioning of the mechanisms driving interannual variability in the upper-ocean dissolved inorganic carbon inventory (0-100 m) time rate of change DI DIC /Dt (mol m À2 yr À1 ) in the CCSM ocean model.The panels show the rms of the model annual DI DIC /Dt anomalies (1979-2004) (top left) and the contributions from gas exchange (top right), particulate export (middle left), resolved horizontal and vertical advection (middle right), horizontal and vertical diffusion (bottom left), and surface freshwater fluxes (bottom right).

Fig. 9 .
Fig. 9. Comparison of interannual ocean carbon cycle variability driven by ocean physics versus atmospheric dust deposition.The top panel shows a global spatial map of the mean dust deposition field (g m À2 yr À1 ; note non-linear color bar) and the coefficient of variation (standard deviation/mean; unitless) in the deseasonalized monthly dust deposition anomalies.The middle and bottom panels show the rms of the monthly deseasonalized anomalies in air-sea CO 2 flux (mol m À2 yr À1 ) (middle row) and export production (mol m À2 yr À1 ) (bottom row) from the CCSM ocean BEC model forced by interannually varying physics and repeat annual dust (left column) and repeat annual physics and interannually varying dust (right column).

Fig. 10 .
Fig. 10.Time-series of globally integrated, monthly deseasonalized anomalies in air-sea CO 2 flux (Pg C yr À1 ) (top row) and export production (Pg C yr À1 ) in the CCSM ocean BEC model forced by interannually varying physics and repeat annual dust (left column) and repeat annual physics and interannually varying dust (right column).

Fig. 11 .
Fig. 11.Regional time-series of the monthly deseasonalized anomalies from the CCSM-3 ocean BEC model for the repeat annual physics and interannually varying dust experiment: atmospheric dust flux (g m À2 yr À1 ) (left column), particle export production (mol m À2 yr À1 ) (middle column), and air-sea CO 2 flux (mol m À2 yr À1 ) (right column).Note that the scale for the dust deposition differs from panel to panel.Positive fluxes are into the ocean, and negative fluxes are out of the ocean.A 2-yr running mean is overlaid in red.The rms of the deaseasonalized anomalies is reported for each region.