Open Access Controls on the spatial distribution of oceanic δ 13 CDIC

We describe the design and evaluation of a large ensemble of coupled climate–carbon cycle simulations with the Earth system model of intermediate complexity GENIE. This ensemble has been designed for application to a range of carbon cycle questions, including the causes of late-Quaternary fluctuations in atmospheric CO 2 . Here we evaluate the ensemble by applying it to a transient experiment over the recent industrial era (1858 to 2008 AD). We employ singular vector decomposition and principal component emulation to investigate the spatial modes of ensemble variability of oceanic dissolved inorganic carbon (DIC) δ 13 C, considering both the spun-up pre-industrial state and the transient change. These analyses allow us to separate the natural (pre-industrial) and anthropogenic controls on the δ 13 C DIC distribution. We apply the same dimensionally-reduced emulation techniques to consider the drivers of the spatial uncertainty in anthropogenic DIC. We show that the sources of uncertainty related to the uptake of anthropogenic δ 13 C DIC and DIC are quite distinct. Uncertainty in anthropogenic δ 13 C uptake is controlled by air–sea gas exchange, which explains 63% of modelled variance. This mode of variability is largely absent from the ensemble variability in CO 2 uptake, which is rather driven by uncertainties in thermocline ventilation rates. Although the need to account for air–sea gas exchange is well known, these results suggest that, to leading order, uncertainties in the ocean uptake of anthropogenic 13 C and CO 2 are governed by very different processes. This illustrates the difficulties in reconstructing one from the other, and furthermore highlights the need for careful targeting of both δ 13 C DIC and DIC observations to better constrain the ocean sink of anthropogenic CO 2 .


Introduction
A substantial component of the uncertainty associated with both anthropogenic and natural climate change derives from uncertainties in the carbon cycle. Improved understanding is required to better quantify the response of the Earth system and changing strength of various feedbacks to ongoing carbon emissions, especially given recent observations that the oceanic carbon sink is weakening (Le Quéré et al., 2007). For instance, while coupled climate-carbon cycle models consistently predict a weakened efficiency of the ocean-terrestrial carbon sink under future-emissions scenarios, they do so with a highly uncertain magnitude with projections of year 2100 atmospheric CO 2 (SRES A2 forcing) ranging from 740 to 1030 ppm (Friedlingstein et al., 2006). Entraining observations in addition to that of atmospheric CO 2 is essential to better constrain such carbon cycle changes.
Anthropogenic emissions of CO 2 from the burning of fossil fuels and deforestation are strongly depleted in 13 C, reflecting the preferential uptake of 12 C during photosynthesis. As a consequence, the δ 13 C composition of the atmosphere, where δ 13 C = 1000 (( 13 C/ 12 C) sample /( 13 C/ 12 C) standard -1), has decreased from about −6.5 ‰ in pre-industrial times to −8.2 ‰ in the present day (Keeling et al., 2010) and is is known as the 13 C Suess effect (Keeling, 1979). It has also produced an imprint on the oceanic δ 13 C DIC distribution (Gruber et al., 1999), which can then be employed to help constrain estimates of fossil fuel emissions and their uptake (e.g. Joos and Bruno, 1998;Quay et al., 2003). However, the oceanic δ 13 C distribution reflects a complex interplay between marine productivity, water column remineralisation of organic matter, ocean circulation and mixing, and air-sea gas exchange (Gruber et al., 1999;Tagliabue and Bopp, 2008). A P. B. Holden et al.: Controls on the spatial distribution of oceanic δ 13 C DIC clear understanding of the factors that control the distribution of the 13 C/ 12 C ratio in the ocean is essential if observations of δ 13 C DIC are to be of use in constraining the carbon cycle.
Faced with large uncertainties in multiple processes and feedbacks in the global carbon cycle, Earth system models of intermediate complexity have become important tools for helping explore process sensitivities and quantifying uncertainty. Carefully designed ensembles of simulations can be used to derive data-constrained probability distributions and complement high-complexity (but small in number) simulations of fully coupled ocean/atmosphere-based Earth system models. We take this approach here, describing an ensemble of instances of a low-resolution Earth system model with a non-dynamical energy balance atmosphere, and associated analysis of the controls on anthropogenic CO 2 and δ 13 C uptake by the ocean.
In this paper Sects. 2 to 5 describe the development and evaluation of a perturbed-parameter ensemble of the Grid ENabled INtegrated Earth system model (GENIE). This ensemble has been applied to several experiments addressing the historical and future carbon cycle (Eby et al., 2012;Holden et al., 2013;Joos et al., 2013;Zickfeld et al., 2013). Sections 6 and 7 then describe the application of this ensemble to investigate the processes controlling the δ 13 C DIC distribution in both the pre-industrial (assumed equilibrium) and modern ocean states. In both cases we apply singular vector decomposition to identify the spatial patterns corresponding to the dominant modes of variability amongst the simulations. Emulation of the principal components, following , then enables us to evaluate the processes and parameters driving this simulated variability. It is important to note that this evaluation is not purely objective as it is dependent upon the prior distributions of parameter ranges. We further apply this analytical approach to the oceanic uptake of anthropogenic CO 2 in order to examine the relationship between the resulting changes to oceanic distributions of DIC and δ 13 C DIC .

The GENIE configuration
We utilise a coupled carbon cycle-climate configuration of GENIE (version 2.7.7). The physical model comprises the 3-D frictional geostrophic ocean model GOLDSTEIN (at 36 × 36 × 16 resolution) coupled to a 2-D energy moisture balance model of the atmosphere (EMBM) and a thermodynamic/dynamic sea ice model (Edwards and Marsh, 2005;Marsh et al., 2011). Adjustments to the parameterisations of outgoing long-wave radiation (OLR) and lapse rate effects are described in . Notable recent refinements to the ocean model are the addition of stratificationdependent mixing (Oliver and Edwards, 2008), the inclusion of the thermobaricity term in the equation of state (K. I. C. Oliver, personal communication, 2013) and the incorporation of alternative wind forcing fields which have improved the representation of global ocean circulation, especially in terms of the strength of the Antarctic Circumpolar Current (Marsh et al., 2011). The land surface module is the model of terrestrial carbon storage ENTS . Ocean biogeochemistry is modelled with BIOGEM, based on the model of Ridgwell et al. (2007a) including the cycling of iron described by Annan and Hargreaves (2010), except with biological uptake following Doney et al. (2006). For consistency with the wind stress forcing of the ocean model, the prescribed wind speed field used by Ridgwell et al. (2007a) to compute the squared wind speed coefficient in the air-sea gas exchange parameterisation was replaced with wind speeds derived from the annual mean wind stress climatology applied in the ocean model component. Sediments are modelled with SEDGEM at 36 × 36 resolution (Ridgwell and Hargreaves, 2007). The rock-weathering module ROKGEM (Colbourn, 2011) is included to redistribute prescribed weathering fluxes according to a fixed river-routing scheme.
For each ensemble member (see Sect. 3), the model was first spun up for 25 000 yr in order to bring the sediments close to equilibrium with pre-industrial boundary conditions. Atmospheric CO 2 and δ 13 C were relaxed to pre-industrial values of 278 ppm and −6.5 ‰ , respectively. During model spin-up, ocean biogeochemistry was simulated as a closed system (i.e. sedimentary fluxes, except iron, were returned at the bottom of the ocean) and weathering fluxes required to balance sediment burial were diagnosed at the end of the spin-up phase (Ridgwell and Hargreaves, 2007). In order to perform the industrial transient simulations (1858 to 2008 AD), the sediment system was opened (i.e. exchanges of sediment/ocean fluxes were applied as simulated by the sediment model), applying the constant weathering flux. Atmospheric CO 2 and δ 13 C were relaxed to (time-varying) observations.

Ensemble design
The philosophy for the design process has been described in detail elsewhere . In short, the approach attempts to vary key model parameters over the entire range of plausible input values and to accept parameter combinations which lead to climate states that cannot be uncontroversially ruled out as implausible . The approach represents an attempt to find plausible realisations of the model from all regions of the highdimensional input parameter space in order to capture the range of possible feedback strengths. In general, the "best" parameter values are not known as both model structural errors and observational errors limit the degree to which a parameter can be constrained. In our approach we apply broad parameter ranges, generally broader than ranges that are applied in model-tuning exercises. This is in part to enable us to fully quantify model behaviour over not-implausible P. B. Holden et al.: Controls on the spatial distribution of oceanic δ 13 C DIC parameter space and in part to improve the validity of the ensemble for application to diverse climate states, such as the Last Glacial Maximum. We refer to the resulting uncertainty as parametric uncertainty, but are cognisant of the fact that the broad ranges applied are in some cases deliberately larger than the "true" parametric uncertainty associated with the pre-industrial state.

Parameters
We varied 24 model parameters in the ensemble (see Table 1). Overall, the model configuration and choice of variable parameters was governed by consideration of the processes that are thought to contribute to variability of atmospheric CO 2 on glacial/interglacial timescales (Kohfeld and Ridgwell, 2009), and hence to which the distributions of carbon and carbon isotopes may be sensitive in general. The varied parameters are summarised as follows:

Ocean
The ocean frictional parameter (ODC) and three ocean diffusivity parameters (OHD, OVD and OP1) are varied. The uncertain impact of winds on stratification is captured through the wind stress scaling parameter (WSF). The ranges for ODC, OHD and WSF were taken from . The parameter OP1 defines the e-folding depth of the reference profile of diapycnal mixing (Oliver and Edwards, 2008). The ranges for OP1 and OVD were chosen to allow diapycnal diffusivities that vary from 10 −6 to 2 × 10 −4 m 2 s −1 in the surface ocean and from 2 × 10 −5 to 4 × 10 −3 m 2 s −1 in the deep ocean.

Atmosphere
The EMBM parameters were chosen as those that control global temperature (OL0, a perturbation to the clearskies outgoing long-wave radiation), Atlantic overturning strength (AMD, atmospheric moisture diffusivity, and APM, Atlantic-Pacific moisture flux correction) and Antarctic sea P. B. Holden et al.: Controls on the spatial distribution of oceanic δ 13 C DIC ice coverage (OL0 and AHD, the atmospheric heat diffusivity through its control on the equator-pole temperature gradient). These parameter choices, and ranges, were identified from previous ensembles (Fig. 2 of . Where applicable, ranges were reduced so as not to exceed the mean ±2σ of the  plausibility-filtered parameter distributions.

Sea ice
By varying sea ice diffusivity (SID) we attempt to represent uncertainty introduced by brine rejection on Antarctic Bottom Water (AABW). Brine rejection is a potential mechanism for increasing deep ocean stratification and drawing down CO 2 from the atmosphere in the glacial state (Bouttes et al., 2010). Varying SID provides a means to vary this feedback strength in a way that naturally connects the dynamics of brine rejection and sea ice. A wide range of Atlantic circulation states, ranging from the near-total dominance of AABW to the total dominance of North Atlantic Deep Water (NADW) can be achieved by spanning the ranges for SID in Table 1. While it is difficult to quantify the plausible range of input space, we rather choose an input range that spans the not-implausible range of output space, and then constrain output space accordingly (Sect. 3.2). We note that the SID range considered here extends to substantially higher values than have been considered in previous ensembles (e.g. Edwards and Marsh, 2005;Edwards et al., 2011). We note that the uncertainty in sea ice extent is controlled by atmospheric parameters rather than SID (Sect. 3.1.2). A change in sea ice extent can modify the rate of air-sea gas exchange in polar oceans. A wide range of sea ice coverage was allowed in the ensemble (Sect. 3.2).

Ocean biogeochemistry
Oceanic biological activity affects atmospheric CO 2 concentrations through changes in the relative rates of CaCO 3 and particulate organic carbon (POC) supply to the seafloor. Uncertainty is captured through parameters controlling the export ratio and remineralisation depths. Changes in Redfield ratios or remineralisation depths have the potential to cause significant changes in carbon isotope and DIC storage and may be linked to climate (Omta et al., 2006). While the dynamics of possible feedbacks may not be captured, the wide range of input values allows us to capture possible sensitivity to these processes. Key biogeochemistry parameters (PHS, PRP, PRD, RRS, TCP, PRC and CRD) were chosen, corresponding to seven of the eight parameters selected for the Ensemble Kalman Filter (EnKF) tuning of the biogeochemical model (Ridgwell et al., 2007). The ranges for these parameters are approximately centred on the EnKF tuning ranges, but are broadened to reflect the philosophy of the ensemble design. We note that an eighth parameter from Ridgwell et al. (2007), the maximum phosphate uptake, was neglected in error. The basic configuration of the biogeochemical model was changed late in the design process to incorporate the most up-to-date iron-cycling scheme. Instead of a parameter for the maximum phosphate uptake, the ironcycling scheme prescribes a phosphate uptake timescale but the configuration files were not updated to reflect this. The effect of not varying the phosphate uptake will be discussed below.
Productivity is iron limited in regions such as the Southern Ocean. Iron supply to the ocean is included, and varied, in the experiment by prescribing the pre-industrial dust fluxes of Mahowald et al. (2005) and varying the soluble iron fraction (FES) over ranges based on exploratory ensembles (Ridgwell and Death, 2013).
Air-sea gas exchange rates exert controls on the degree of ocean-atmosphere equilibration and 13 C fractionation. The globally uniform scaling constant of the air-sea gas exchange parameterisation ASG (Eq. (13), Ridgwell et al., 2007) was varied in the ensemble design. We note that the air-sea gas exchange parameterisation (Wanninkhof, 1992) is proportional to the square of the wind speed, but is applied to the long-term average wind speed, thus neglecting non-linear contributions of short-term variability which contribute substantially to the globally averaged air-sea gas transfer velocity (Wanninkhof, 1992). This introduces a systematic, spatially non-uniform bias in the disequilibrium of δ 13 C at the air-sea interface simulated by the model (typically underestimating the air-sea exchange at high latitudes). An additional spatial distortion is introduced close to continental boundaries because, in making an estimation of wind speed from wind stress, no attempt is made in the model to transform from edge-centred to tracer grid. This introduces an underestimation of the wind speed in some grid cells close to continental boundaries. The wide range over which the air-sea gas exchange was varied is likely to compensate for these shortcomings on the global scale. The ensemble average of the global mean gas exchange coefficient for CO 2 of 0.038 ± 0.014 mol m −2 yr −1 µatm −1 compares to the range of (0.052 ± 0.010) mol m −2 yr −1 µatm −1 , based on the estimate of the air-sea gas exchange velocity (and an estimate of the globally averaged solubility of CO 2 ) by Naegler (2009). Thus, while the ASG parameter (0.1 to 0.5) spans the value for long-term average wind speeds (0.39), given by Wanninkhof (1992), our central estimate for the global mean gas exchange coefficient is on the low side, which may reflect the information about fluctuating wind direction that is lost in calculating annually averaged wind stress (from which wind speed is subsequently derived).

Terrestrial carbon
Uncertainty in terrestrial vegetation is captured through five key parameters of the terrestrial vegetation model (VFC, VBP, VRA, LLR and SRT). These parameters exert the  Konkright et al. (2002) dominant controls on uncertainty in vegetative and soil carbon ( Fig. 2 of . The ranges of these five parameters were taken from  and, where applicable, these ranges were reduced so as not to exceed the mean ±2σ of the Holden et al. (2010) plausibility-filtered parameter distributions. We note that in the experiments described here, vegetation uncertainty is only relevant through its direct impact on climate (by way of changes in albedo, surface roughness and soil moisture storage) as atmospheric concentrations of CO 2 and carbon isotopes are relaxed to observational estimates. The vegetation parameters were varied in anticipation of future experiments with freely evolving CO 2 . We note that two parameters that are not varied in the ensemble do not take the default values of . These parameter values (baseline autotrophic respiration rate k 24 = 0.225 yr −1 and baseline heterotrophic respiration rate k 29 = 0.165 yr −1 ) were set equal to the mean values of the precalibration of  in an attempt to best centre the ensemble distribution of global terrestrial carbon pools on plausible ranges.

Sediments
When the distribution of oceanic dissolved carbon and/or alkalinity is changed, the stability of CaCO 3 in the sediments responds. No sediment parameters were varied in this study, but a wide range of sediment states (Sect. 3.2) is allowed in order to capture the uncertainty in this feedback.

Statistical design
An investigative ensemble was performed with a 500-member maximin Latin hypercube (MLH) design for 26 parameters. The 26th parameter was included as a dummy parameter but was not required in this analysis as the inactive maximum phosphate uptake parameter (Sect. 3.1) serves to fulfil this role. There are 24 active parameters, summarised in Table 1. A regression-based predictor, including linear and quadratic terms, was built for each of eight metrics from the output of the MLH ensemble . The eight metrics were chosen to provide global-scale constraints on atmosphere, ocean (NADW and Atlantic AABW), Antarctic sea ice coverage, terrestrial carbon (vegetation and soil), ocean biogeochemistry and ocean sediments. The metrics are summarised in Table 2. Parameters were then sampled randomly from a uniform distribution between plausible ranges (Table 1) and applied as input to these emulators. A parameter set was accepted as potentially plausible when the emulators predicted values within the "approximate Bayesian calibration" (ABC) ranges in Table 2 for all eight metrics. The ABC filtering ranges are centred on observations. They are reduced relative to the plausibility ranges in order to more efficiently generate plausible simulations (i.e. acknowledging that the emulations are imperfect) with a reduction chosen cognisant of the individual emulator errors. Note the poorly centred ABC range for Antarctic sea ice. The range was selected based on estimates of annually averaged area (∼ 10 million km 2 ), but was applied to an emulation of 31 December coverage (estimated to be ∼ 7 million km 2 , Cavalieri et al., 2003). Given the similarity of these metrics (given the wide range accepted as plausible and the emulator error), this did not result in a bias to the ensemble (see Table 2), although it did result in a ∼ 12 % decrease in the ABC success rate.
Each emulator-filtered parameter set was then applied as input to a further simulation. As simulations completed, the emulators were rebuilt four times including the additionally available data. This process progressively improved the success rate of the ABC selection (i.e. the percentage of ABC selected parameter sets that lead to plausible simulations) from 24 % to 65 %. In total, 1000 emulator-filtered parameter sets were applied to the simulator. These produced 885 completed simulations, of which 471 were plausible. This 471-member subset forms the emulator-filtered plausibilityconstrained (EFPC) ensemble. The parameter ranges of the EFPC ensemble ("Filtered range", Table 1) have, in general, been narrowed relative to the prior range (Table 1) by the plausibility filtering process.

Plausibility emulators and total effects
It is useful to apply the plausibility emulators (Sect. 3.2) in order to investigate the role played by the parameter values in determining the spun-up state, the principal motivation being to inform future ensemble design. To achieve this we calculate the "total effect" (Saltelli et al., 2000) of each parameter in each of the emulators. The total effect of a parameter is equal to the expectation of the variance that remains when we have learnt everything about the model besides the value of that parameter. See  for a more complete description of the method. Figure 1 illustrates the total effects for each of the eight plausibility emulators. For each emulator, the total effect is normalised to total 100 %, thus illustrating the relative importance for each parameter in determining the plausible modern state (as defined by the acceptable input ranges in Table 1 and output ranges in Table 2). This analysis was previously performed on the SAT, MAXA, SHSI, CVEG and CSOIL emulators and is discussed in some detail in . Three parameters (AMD, OL0 and VFC) appear less important in this revised analysis. However, this likely reflects the fact that the ranges for these parameters were narrowed for this ensemble (by 40 %, 25 % and 23 %, respectively), highlighting that total effects are not purely objective measures as they are dependent upon estimates of the ranges of input parameters ("Prior range", Table 1). We note further that the MAXA emulator is intrinsically different given the refinements to the ocean model and the increased depth resolution applied in this study (Sect. 2).
In addition to the tests used by , a test for minimum Atlantic overturning strength was introduced to enforce the constraint that plausible ensemble members have an overturning cell associated with AABW. As with all eight plausibility filters, no constraint was imposed upon its spatial distribution, just the maximum value of the stream function (at depths below 500 m and latitudes north of 30 • S). The strength of the AABW cell is controlled by atmospheric moisture transport, isopycnal diffusivity, ocean drag and sea ice diffusivity, each of which play comparably important roles. It is notable that high values of the SID parameter are not necessary to generate a strong AABW cell: for instance, low values of isopycnal diffusivity seem to achieve this objective equally efficiently. However, when several parameters are capable of achieving a similar result independently, it is useful to vary them all in order to span as much plausible input space as is possible. Moreover, the primary reason for including SID as a variable parameter in the ensemble was to span a range of circulation feedback strengths in response to changing sea ice coverage.
The percentage of CaCO 3 in the surface sediment is unsurprisingly controlled by the parameters that control the export of CaCO 3 , i.e. the two parameters (RRS, TCP) controlling the CaCO 3 :POC rain ratio at every grid point (see : Ridgwell et al., 2007a, b) and the fraction of exported CaCO 3 assumed to reach the sediment surface without being subject to water column dissolution (PRC). Also apparent is the importance of parameters controlling the flux of organic matter to the sediments (PRP, PRD) and thereby promoting CaCO 3 dissolution. Dissolved oxygen is mainly controlled by circulation strength and mixing (through APM, OVD, OP1, ODC, WSF), presumably through their role in redistributing nutrients (and oxygen) within the ocean, the POC remineralisation depth (PRD, via its control on respiration in the deep ocean) and iron solubility (FES), which exerts a strong control on productivity. The phosphate half-saturation concentration (PHS) does not appear as a strong control, suggesting that uncertainty in productivity is not dominated by the rate of phosphate uptake, but rather the rate of phosphate redistribution. The dominance of parameters that control the supply of phosphate, rather than the rate of phosphate uptake, suggest that the omission of the maximum rate of phosphate uptake has not resulted in a substantial underestimation of the uncertainty associated with phosphate limitation. In contrast, iron solubility does play a significant role through its influence in determining the extent of iron-limited regions.

Ensemble evaluation
The ensemble design intentionally avoids tuning the model in the conventional sense. The philosophy reflects the acceptance of the fact that, on account of the unknown dynamic behaviour of structural errors that are inevitable in any model, the simulation that best predicts unobserved variables or future change is not necessarily that which most closely fits observations. Consequently, the evaluation performed here is intended to demonstrate the degree to which the eight plausibility constraints are sufficient to produce simulations that we would not rule out as uncontroversially implausible, and to identify biases or structural errors in the ensemble-averaged spatial distributions. Figure 2 plots Taylor diagrams for the model data comparison between the EFPC ensemble and observation fields for six oceanic properties: temperature, salinity, dissolved phosphate, dissolved oxygen, alkalinity and dissolved inorganic carbon. The diagrams depict the relative performance of the 471 ensemble members as well as the outputs of an optimised biogeochemical configuration (Ridgwell and Death, 2013) for comparison. Figures 3 to 8 illustrate surface fields and latitude-depth transects through the Atlantic (25 • W) and Pacific (155 • W) for the same properties, in each case comparing observations with EFPC ensemble averages. Observed surface fields are depth averages over the vertical extent of the surface layer of the ocean model, approximately 81 m. Figure 9 plots ensemble-averaged Pacific and Atlantic cross sections of δ 13 C DIC , in both the pre-industrial and modern states. Figure 10 compares the respective Atlantic cross sections of δ 13 C DIC with the modern (1990s) GLODAP bottled observations (Key at al., 2004; WHP A16 section) and the pre-industrial reconstructions of Olsen and Ninnemann (2010).

Temperature
Temperature exhibits an ensemble-averaged correlation with observations , interpolated onto the GENIE grid) of +0.96, comparable to the optimised model. The spatial variability tends to be greater than observations, a feature which is shared with the optimised model. The ensemble-averaged SST field agrees very well with observations (Fig. 3), so that the overestimate of spatial variability is not reflected in this layer. The ensemble-averaged equator-topole SST temperature gradient reflects observations closely. This is important given the strong temperature dependence of carbon isotope fractionation during air-sea gas exchange (Mook, 1986). The overestimate of the spatial variance is rather related to the vertical temperature gradients, which are greater than observations. This is likely, at least in part, a consequence of the excessive strength of the AABW cell in the ensemble average (Sect. 5.3).

Salinity
The ensemble-averaged salinity correlation with observations  is +0.45 (cf. +0.77 in the optimised model). Poor salinity correlations are expected as atmospheric moisture transport in the simple EMBM is dominated by diffusion. In the Bern3D model (structurally very similar to the configuration of GENIE used here as both are based on the same progenitor ocean model), optimised salinity correlations of 0.65 were obtained (Ritz et al., 2011), compared to 0.9 when the ocean model was restored to climatological sea surface temperature and salinity fields (Müller et al., 2006). The salinity correlation across the EFPC ensemble (i.e. the correlation displayed separately by each of the 471 simulated fields with observations) exhibits a correlation of 0.54 with the atmospheric moisture diffusivity parameter (AMD). Low atmospheric diffusivity inhibits moisture transport out of low latitudes, restricting the development of the increased surface salinities that are observed there. EFPC-averaged AMD of 8.4 × 10 5 m 2 s −1 is somewhat low compared to optimised values ranging between ∼ 1 × 10 6 and 3 × 10 6 m 2 s −1 . We note that the sign of the vertical salinity gradient in the Pacific can be reversed due to this effect resulting in the very low correlations that are seen in some ensemble members.
An important decision therefore was whether to apply additional filters to the EFPC ensemble in order to eliminate those simulations with poorly correlated salinity fields. Outside of the Arctic, the spatial pattern of the ensemble mean surface salinity is reasonable, although it exhibits a bias of ∼ −1 psu (Fig. 4). We note that the ensemble variability of Atlantic (Pacific) surface salinity varies spatially from ∼ 0.2 to 0.4 psu (∼ 0.5 to 1 psu). The salinity latitude-depth cross sections compare reasonably with observations, though again they exhibit a low bias. The major disagreement is thus the understated surface salinities at low to mid latitudes, which may have important consequences for biological production via impacts on stratification-dependent mixing and the redistribution of nutrients back to the photic zone. To test for this, a salinity correlation filter was applied to leave only the 112 ensemble members with a salinity correlation > 0.6. The cross sections of 112-member ensemble-averaged phosphate, oxygen and δ 13 C DIC in the Atlantic, Pacific and Indian oceans were very similar those of the complete EFPC ensemble, suggesting that the impact of the freshened surface waters on biological production is not significant, especially given the philosophy of the ensemble design. To illustrate, a comparison of the δ 13 C DIC cross sections of the means of the two ensembles exhibits an R 2 of > 99 % in all three basins, although the salinity-filtered ensemble exhibits www.biogeosciences.net/10/1815/2013/ . We decided to retain all 471 of the EFPC members, anticipating that the benefits of a large ensemble likely outweigh any beneficial effects of applying this filter.

Dissolved phosphate
The ensemble-averaged correlation of dissolved phosphate with respect to observations  is +0.88, compared to the optimised model correlation of +0.91. The ensemble-averaged variability of the spatial distribution captures the variability of the observations more closely than the optimised model (the cloud is approximately centred on a normalised standard deviation of 1). These statistics together suggest that the neglect of the maximum phosphate uptake as a variable parameter has not biased the ability of the ensemble to capture a wide range of phosphate distributions that are consistent with observations. The ensemble-averaged surface field of dissolved phosphate compares reasonably with observations (Fig. 5), although the increased concentrations associated with upwelling regions are not well captured so that phosphate is near fully utilised throughout low-and mid-latitude oceans. Only 20 simulations exhibit surfacelayer-averaged phosphate concentrations that exceed observation estimates of 0.6 µmol kg −1 . The cross sections compare favourably with observations, although the volume and northerly extent of AABW are overestimated in the ensemble average. This may be a consequence of the very high sea ice diffusivities (SID) that have been allowed in the ensemble leading to increased brine rejection and AABW production. The EFPC distribution for SID is 52 000 ±26 000 m 2 s −1 (1σ uncertainties are quoted throughout). Previous studies have generally not considered values above 25 000 m 2 s −1 (e.g. Edwards et al., 2011).

Dissolved oxygen
Dissolved oxygen exhibits an ensemble-averaged correlation of +0.76 (cf. +0.83 in the optimised model) and the normalised standard deviation is approximately centred on unity, demonstrating a comparable spatial variability to observations .  equatorial upwelling are not captured by the model (Fig. 6). The dissolved oxygen cross sections compare well to observations, a conspicuous difference being a ∼ 20 % underestimate of oxygen in NADW which, in the absence of large differences in temperature (Sect. 5.1), is suggestive of excess productivity in high northern latitudes and consistent with the underestimate of surface phosphate in the Arctic (Fig. 5). It could also reflect a bias towards low rates of airsea gas exchange (see Sect. 3.1.4), preventing the North Atlantic surface from fully equilibrating with the atmosphere. The plausibility constraints require reasonable globally averaged oxygen concentrations (Table 2). It appears that the underestimate of equatorial productivity, which we interpret as a structural issue, may favour parameters that compensate by increasing productivity in high latitudes, notably high northern latitudes that are not iron limited.

Alkalinity
Alkalinity exhibits an ensemble-averaged correlation of +0.65 (cf. +0.80 in the optimised model). The normalised standard deviation is 1.9 (cf. 1.6 in the optimised model) demonstrating, in both the ensemble and the optimised model, greater spatial variability than observations (Key et al., 2004). The ensemble-averaged surface layer displays a low bias with respect to observations. Conversely, the 40 ensemble-averaged deep ocean exhibits a high bias with respect to observations. We note that the global average alkalinity is held fixed at 2363 µmol kg −1 during these "closed" pre-industrial spin-ups. The source of these biases is not related (as in the case of salinity) to atmospheric moisture diffusivity, but rather to parameters that drive ocean tracer distributions.

Dissolved inorganic carbon
Dissolved inorganic carbon (DIC) exhibits an ensembleaveraged correlation of +0.86 (cf. +0.88 in the optimised model) and the normalised standard deviation is approximately centred on unity (average 1.3), demonstrating a comparable spatial variability to observations (Key et. al., 2004). The ensemble-averaged spatial distribution of DIC displays similar characteristics to the distribution of ALK, notably a low bias in the surface ocean and a high bias in the deep ocean. The close relationship between DIC and ALK is further demonstrated by the correlation of +0.86 between their standard deviations (i.e. data plotted in Fig. 2e and f). The low surface bias may have implications for the evaluation of anthropogenic carbon uptake. In the Pacific, anthropogenic carbon penetrates to depths of ∼ 1000 m (Sect. 7). Ensemble-averaged pre-industrial DIC in the Pacific cross section to depths of ∼ 930 m is 2110 µmol kg −1 , compared to P. B. Holden et al.: Controls on the spatial distribution of oceanic δ 13 C DIC 41 observations of 2160 µmol kg −1 , suggesting this bias is likely to be modest. Figure 9 shows EFPC pre-industrial and modern (2008 AD) ensemble-averaged latitude-depth transects of δ 13 C DIC through the Atlantic (25 • W) and Pacific (155 • W). The discussion in the following paragraphs addresses the evaluation of the simulated present day δ 13 C DIC distribution with respect to the observational data of Kroopnick (1985, Figs. 3 and 5) and Gruber et al. (1999). In Fig. 10 a further comparison is provided, and illustrated, against Atlantic GLODAP modern observations (Key at al., 2004) and the pre-industrial reconstructions of Olsen and Ninnemann (2010).

Ocean δ 13 C DIC
In the Atlantic, ensemble-averaged δ 13 C DIC of AABW and AAIW masses varies in the range ∼ 0.3 ‰ to ∼ 0.7 ‰ and NADW varies in the range ∼ 0.7 ‰ to ∼ 1.2 ‰, in both cases consistent with observations (Kroopnick, 1985). The simulated 0.75 ‰ contour extends to ∼ 35 • S, again in good agreement with the observational data, but only reaches depths of ∼ 3000 m -in disagreement with the observations -which indicate that values in excess of 0.7 ‰ extend to the ocean floor north of the Equator. This is consistent with the excessive northerly penetration of AABW that is apparent in the phosphate distribution (Sect. 5.3). We note that ensemble-42 averaged oceanic-mean δ 13 C DIC = 0.38 ±0.22 ‰ is slightly lower than observational estimates of 0.5 ‰ (Quay et al., 2003), also consistent with too much simulated AABW.
In the Pacific, ensemble-averaged δ 13 C DIC at depths below 1000 m ranges from ∼ −0.5 ‰ at high northern latitudes to ∼ 0.4 ‰ at high southern latitudes. These are in reasonable agreement with observations, although observational data reach values as low as −1.0 ‰ at high northern latitudes. In the simulations, the 0 ‰ contour extends south to ∼ 10 • S and to depths of ∼ 3000 m, in good agreement with observations.
Simulated surface values of δ 13 C DIC (not shown) in the Pacific of ∼ 1.5 ‰ are ∼ 0.3 ‰ lower than observations of Gruber et al. (1999), consistent with a 13 C Suess effect of −0.018 ‰ yr −1 (Gruber et al., 1999). Simulated Atlantic surface δ 13 C DIC of ∼ 2 ‰ is comparable to Gruber et al. (1999), so that the lack of a "Suess offset" indicates a high bias (∼ 0.3 ‰) in simulations of surface Atlantic δ 13 C DIC . Meridional surface profiles display similar characteristics to observations of Gruber et al. (1999): low values (∼ 1 ‰ ) at high southern latitudes of the Atlantic, a maximum in both basins at ∼ 40 to 50 • S, a peak in the Pacific at ∼ 15 • S and a trough in the Atlantic at ∼ 30 • N (located at ∼ 20 • N in the simulations). These complex surface distributions can be explained through the interplay of the temperature-dependent fractionation of air-sea exchange, upwelling of 13 C depleted waters in the Southern Ocean and meridionally variable biological uptake rates and fractionation, greatest in regions of nutrient upwelling and lowest in the centres of subtropical gyres (Gruber et al., 1999). The Suess effect has been observed to homogenise the δ 13 C DIC distribution in the North Atlantic, largely erasing the observed tongue of high δ 13 C DIC waters at mid-depths in the reconstructed pre-industrial ocean (Olsen and Ninnemann, 2010). This effect is illustrated in Fig. 10. The meridional gradient of reconstructed pre-industrial δ 13 C DIC at depths of ∼ 2000 m is well captured in the ensemble average (Fig. 10b). This gradient is substantially weaker in the present-day ensemble average, consistent with observations (Fig. 10a). We note that the GLODAP observations (Key et al., 2004) were taken during cruises in the 1990s (cf. 2008 simulations). The principal weakness of the simulations apparent from these comparisons is, as with the comparison with Kroopnick (1985), the underestimation of δ 13 C DIC in the deep North Atlantic, likely reflecting the excessive northerly penetration of AABW in the ensemble mean. 44

Drivers of the spatial distribution of pre-industrial δ 13 C DIC
Aside from the excessive northerly intrusion of AABW and the ∼ 0.3 ‰ high bias in surface Atlantic δ 13 C DIC , the comparisons of Sect. 5.7 demonstrate that the ensemble-averaged δ 13 C DIC spatial distributions are in good agreement with observations. The following analysis seeks to explain the processes and parameters that determine these spatial distributions. We apply singular vector decomposition (SVD) to determine the dominant patterns (empirical orthogonal functions (EOFs)) of variability across the ensemble and then to emulate the principal components (PCs) as functions of model parameters. Any of the simulated fields can be constructed as a linear combination of EOFs, weighted by their respective PCs. As such, an emulation of a PC can be used to quantify the contribution of each parameter to the uncertainty associated with the respective EOF. We analyse both the preindustrial spun-up states and the transient anomaly (the difference between years 2008 and 1858) in order to separate equilibrium effects from anthropogenic 13 C uptake.  .

EOF decomposition
We consider two latitude-depth transects through the Atlantic (25 • W) and Pacific (155 • W). These transects run between latitudes 71 • S and 63 • N, thus including the Southern Ocean but not the Arctic Ocean. The EFPC ensemble was decomposed into EOFs by applying SVD to the ensemble of δ 13 C DIC fields. The analysis was applied by first forming a 1152-element vector for each ensemble member, constructed from the combined Atlantic and Pacific fields (16 depth cells × 36 latitude cells × 2 ocean basins). The ocean basins are combined in this way to derive a single set of EOFs that describes the variability in both basins self-consistently. The 471 ensemble members are then combined to form a 1152 × 471 centred matrix of δ 13 C DIC fields Y and singular vector decomposition is applied: where L is the (1152 × 471) matrix of left singular vectors ("EOFs"), D is the 471 × 471 diagonal matrix of eigenvalues and R is the 471 × 471 matrix of right singular vectors ("principal components"). Figure 11 illustrates the first three EOFs of the preindustrial distribution, plotted separately in the two basins. These three EOFs explain 42 %, 27 % and 11 % of the variance across the ensemble. Figure 11 additionally shows the emulator coefficients associated with each EOF, discussed in Sect. 6.2.
The first EOF, explaining 42 % of the ensemble variance, is of the same sign at every location across the two basins.  Therefore, processes controlling the first PC (i.e. the coefficients of the first EOF) may be interpreted as processes that control the partitioning of δ 13 C between the ocean and the atmosphere. Note that the experimental setup precludes partitioning of δ 13 C between the ocean and sediments (and between the atmosphere and terrestrial biosphere). Since atmospheric δ 13 C was relaxed to observations, the atmosphere/ocean partitioning is expressed in our ensemble as a change in global ocean mean δ 13 C. Had we instead held constant the ocean/atmosphere inventory of δ 13 C constant, this partitioning would have been expressed as a change in atmospheric δ 13 C of a similar magnitude, but the opposite sign. This change in δ 13 C DIC is not uniformly distributed throughout the global ocean, but is qualitatively similar to the spatial pattern of the ensemble mean in both basins (Fig. 9). The correlation between globally averaged δ 13 C DIC and the 1st PC is −0.93, confirming this interpretation. (Note that both the EOF and the correlation coefficient have a negative sign.) The change in globally averaged δ 13 C DIC is dominantly associated with changes in the δ 13 C DIC of AABW and AAIW. The Pacific EOF (average −0.035) has a greater magnitude than the Atlantic EOF (−0.026) so, volume effects aside, a global increase in δ 13 C DIC appears to be more pronounced in the Pacific Ocean, reflecting the additional influence of Antarctic water masses in this basin.
The second EOF does not exhibit a constant sign, but rather reflects a change in the spatial distribution of δ 13 C DIC . In the Pacific, the EOF clearly reflects the variability in the exchange of δ 13 C DIC values between surface and deep water (i.e. via the export of particulate organic carbon POC). This interpretation also appears to be valid for the Atlantic, especially given that these two fields are from a single EOF and so should reflect modes of variability that are closely related. The Atlantic EOF is complicated by the fact that changes in the surface δ 13 C DIC are transmitted to NADW, so that increased POC export results in elevated δ 13 C DIC in both surface waters and NADW. Whilst the basin-wide average of the second EOF in the Atlantic is of positive sign, in the Pacific it is negative, suggesting that some of the light carbon exported from the surface waters of the Atlantic ends up in the deep Pacific. The second PC exhibits a correlation of −0.30 with global POC export, so as POC export increases the vertical gradient of δ 13 C DIC is reduced. This is discussed further in Sect. 6.2. Furthermore, the second PC exhibits a weak correlation (+0.2) with the globally averaged DIC inventory, so that an increase in the DIC reservoir is associated with an increase in the vertical gradient of δ 13 C DIC , consistent with the application of this metric as a proxy for oceanic carbon storage and atmospheric CO 2 (Shackleton et al., 1983).
The third EOF takes generally positive values throughout the Atlantic and generally negative values throughout the Pacific, suggesting it describes a mechanism for inter-basin exchange of δ 13 C DIC . The greatest variability in the Atlantic is apparent at the boundary between Atlantic and Antarctic source waters, suggesting that this EOF reflects the relative influence of these water masses. Increased influence of NADW leads to increased values of deep Atlantic δ 13 C DIC . This may be associated with a depletion of 13 C DIC in the Arctic Ocean that possibly explains the lowered δ 13 C DIC in the surface waters of the northern Pacific (connected, in this model configuration, via diffusive transport through the shallow Bering Strait).

Principal component emulation
In order to better understand which parameters, and hence which processes, are driving the spatial variability in simulated δ 13 C DIC , simple linear emulators of the first three PCs corresponding to the EOFs described in Sect. 6.1 were derived.
where each PC ζ is expressed in terms of the 25 parameters θ i and the scalar coefficients a and b i . Note that not all possible terms are included as the function is required to satisfy the Bayes information criterion. Linear terms were found sufficient to produce models with an adequate fit (model R 2 ∼ 85 % in each of the three PC models) and are more straightforward to interpret than more complex (e.g. quadratic) functions. Strong correlations exist between many parameter pairs; twelve parameter pairs exhibit a correlation P. B. Holden et al.: Controls on the spatial distribution of oceanic δ 13 C DIC (positive or negative) of > 0.3. These correlations were introduced by the ensemble design (which constrains for plausible pre-industrial states) and their presence makes cross-terms difficult to reliably interpret.
The parameter coefficients b i for each of the PC emulators are plotted in Fig. 11. The relative magnitudes of the coefficients reflect the relative importance of the parameters in determining each PC. However, it is important to note that, as with total effects (Sect. 4) and indeed with the apportionment of variance to the EOFs, these are not purely objective measures as they are dependent upon the ranges over which parameters were varied. We further note that linear terms are likely to be complicated by the inter-parameter correlations as discussed above, so that a degree of caution is required with respect to their interpretation. In order to validate the robustness of parameter coefficients, we additionally applied the techniques of EOF decomposition and PC emulation to the MLH ensemble (Sect. 3.2). The first three EOFs of the MLH ensemble exhibit correlations of 0.83, 0.96 and 0.87 with the respective EOFs of the EFPC ensemble. The two ensembles thus display similar high-order modes of variability, so that an emulation of the MLH ensemble (with uncorrelated parameter inputs) provides a useful validation of the robustness of the EFPC emulators.

1st EOF of pre-industrial δ 13 C DIC
In Sect. 6.1, the 1st EOF was identified as corresponding to changes in globally averaged δ 13 C DIC . A variety of processes control the atmosphere/ocean partitioning of δ 13 C. Firstly, fractionation during air-sea gas exchange depends on surface temperature, with warming acting to decrease δ 13 C DIC for a given δ 13 C atm . In the ensemble, increasing OL0 warms the surface ocean by decreasing outgoing long-wave radiation, decreasing global δ 13 C DIC via temperature fractionation during air-sea gas exchange. The EFPC emulator coefficient for OL0 (AMD) is increased (decreased) relative to the MLH coefficient. These differences may, at least in part, be a consequence of the correlation between OL0 and AMD (−0.41) in the EFPC ensemble that is introduced through their similar effects on Antarctic sea ice area in the plausibility filtering. However, ocean mean δ 13 C DIC across the ensemble satisfies the relationship δ 13 C DIC = −0.091T + 2.159 (R 2 = 21 %), where T is the globally average sea surface temperature. This is consistent with the temperature variation of oceanic δ 13 C in thermodynamic equilibrium with atmospheric δ 13 C of CO 2 of −0.1 ‰ • C −1 (Mook, 1986), suggesting that the interpretation of a temperature dependence is robust.
Secondly, if air-sea exchange is sufficiently slow compared with vertical exchange within the ocean, surface δ 13 C DIC can be prevented from approaching thermodynamic equilibrium with the atmosphere. This effect is most important in the high-latitude Southern Ocean, where the vertical exchange is rapid and sea ice cover impedes air-sea gas exchange. The heat flux out of the ocean at these latitudes acts to increase δ 13 C DIC , due to temperature dependent fractionation, so that a decrease in air-sea gas exchange would be expected to reduce δ 13 C DIC . This is consistent with the negative dependence of the first PC on air-sea gas exchange and atmospheric heat diffusivity (the dominant control on sea ice cover).
Neither of these two processes can explain the amplification of the first EOF in subsurface waters, especially Antarctic Intermediate Water, and other lower pycnocline waters. This is controlled by ventilation of these water masses, which acts to erode the vertical gradient in δ 13 C DIC (Fig. 9) within the pycnocline that is created by the downward particulate organic flux of isotopically light carbon. Increasing wind stress increases ageostrophic flow in the near-surface layers, enhancing this ventilation and transporting isotopically heavy surface waters into the pycnocline. By reducing surface δ 13 C DIC , which must be compensated by air-sea exchange, this controls atmosphere/ocean partitioning of δ 13 C.
The effect of export of isotopically light organic carbon is apparent through the weak dependencies on the parameters that control the soft-tissue pump (PHS, PRP and PRD), with coefficients of opposite sign to those of the parameters that control (weakly fractionated) CaCO 3 export (RRS, TCP and PRC). Although the EFPC coefficients appear, surprisingly, to suggest that the inorganic pump dominates over the soft-tissue pump, this behaviour is not apparent in the MLH emulator. There are significant correlations between these three parameters (RRS:TCP 0.53 and TCP:PRC 0.42), suggesting that these emulator coefficients may not be robust in the EFPC emulator. A separate emulator, built without the three inorganic pump parameters exhibited a standard error of 0.0213, compared to 0.0185 in the full emulator, apparently confirming this interpretation; the inorganic pump EFPC emulator terms are largely cancelling through the autocorrelations.

2nd EOF of pre-industrial δ 13 C DIC
The second EOF describes the partitioning of δ 13 C between the upper ocean (plus NADW), and the deep ocean. This is somewhat analogous to the partitioning of δ 13 C DIC between the surface ocean and the lower pycnocline described within the 1st EOF, but different parameters control ventilation of the deep ocean within the model. Ocean parameters control the strength of the energy source for diapycnal mixing and hence the strength of the AABW overturning cell that erodes the vertical gradient in δ 13 C DIC . On the other hand, a greater particulate flux to the deep ocean enhances this gradient. The apparent dependency on inorganic pump parameters in EFPC is, again, probably an artefact of inter-parameter correlations (see Sect. 6.2.1). If the emulator is rebuilt neglecting these three parameters, the POC remineralisation depth (PRD) dominates the resulting emulator (consistent with the MLH emulator). Increasing the POC remineralisation depth leads to the transfer of nutrients to greater depths from where they cannot be re-mixed or upwelled as readily into the euphotic zone. This enhanced export of nutrients to depth decreases productivity and export from the mixed layer, but the increased export from the pycnocline to the deep ocean dominates the overall effect so that an increased in PRD results in an increased near surface-to-deep gradient of δ 13 C DIC . Since the second EOF contains the deep expression of both Southern Ocean surface δ 13 C DIC disequilibrium and the vertical gradient in δ 13 C DIC , for which the near-surface expressions are in the first EOF, it has a signature in the ocean/atmosphere partitioning of δ 13 C DIC . This is an example how principal component analysis, being a purely statistical technique, may not cleanly distinguish the key physical processes.

3rd EOF of pre-industrial δ 13 C DIC
The third PC is determined by OHD and APM (in both the EFPC and MLH emulators). The isopycnal diffusivity OHD correlates strongly with the strength of the AABW cell (+0.37), so that high values of OHD are associated with strengthened AABW circulation (within which DIC is depleted in δ 13 C DIC relative to NADW). Conversely, high values of APM lead to strengthening of NADW (through a reduction in Atlantic surface salinity) and relative enrichment of δ 13 C DIC in the region where the boundaries between NADW and AABW are least clearly defined across the ensemble. As noted in 6.1, this effect appears to dominantly control the inter-basin contrast of δ 13 C DIC , although the third PC is weakly correlated with the global average δ 13 C (+0.09), suggesting it is also associated with an effect on the oceanic-mean of δ 13 C DIC . The inter-basin contrast may arise from the controls exerted by APM and OHD on the depth (and therefore the δ 13 C DIC ) at which water flows from the Atlantic through controls on the degree to which winddriven circulation is balanced by differing contributions of surface waters and NADW.
Note that the three PCs together explain 99.4 % of the variance in the ocean mean δ 13 C DIC , (86.0 %, 12.5 % and 0.9 % respectively) so that almost all of the ensemble variability in ocean mean pre-industrial δ 13 C DIC is described by these three modes of variability.

Decomposition and emulation of the anthropogenic 13 C and DIC imprints
The method presented in Sect. 6 is applied again here to explore the changes in the δ 13 C DIC and DIC distributions from pre-industrial to modern times. The EOF analysis and PC emulation is now applied in order to investigate the processes and parameters that control the transient response to changes in the atmospheric carbon dioxide concentrations and its isotopic composition (Fig. 12). For δ 13 C DIC , we restrict our analysis to the first two EOFs, which describe 63 % and 18 % of the variance of the transient response, respectively. For DIC, we consider only the first EOF, which describes 54 % of the ensemble variance.

1st EOF of anthropogenic δ 13 C DIC
The first EOF is of constant sign throughout both basins and so represents a change in global δ 13 C DIC . This can be interpreted as simply due to the uptake of relatively light carbon from the atmosphere (in which the 13 C isotopic ratio has moved to lower values since pre-industrial times as a result of the burning of isotopically light fossil fuels). This change is restricted to depths up to ∼ 1000 m in the Pacific, but the changes are redistributed to greater depth by NADW in the Atlantic. Unsurprisingly, the first PC emulator is dominated by the air-sea gas exchange parameter, which exhibits a correlation of −0.89 with the first PC (a sufficient demonstration of emulator robustness). Globally increasing air-sea gas exchange strength acts to reduce the equilibration timescale of surface δ 13 C DIC with atmospheric δ 13 C. (The negative correlation is consistent with a lowering of the δ 13 C DIC as the EOF is of positive sign.) The equilibration timescale for δ 13 C DIC at the sea surface is of the order of a decade, hence there are large disequilibria at the air-sea interface (see e.g. Gruber et al., 1999).

2nd EOF of anthropogenic δ 13 C DIC
The second EOF describes the partitioning of δ 13 C DIC between surface and intermediate waters, indicative of the degree to which the isotopically lightened surface carbon is mixed with lower pycnocline and intermediate waters. The linear emulator of the second PC is dominated by the global wind stress scaling factor and the ocean drag coefficient. High values of wind stress scaling increase mixing and upwelling, resulting in preferential lowering of δ 13 C DIC in intermediate waters relative to surface waters. WSF exhibits a correlation of 0.99 with the maximum of the overturning stream function in near-surface waters, supporting the interpretation that controls on the ventilation of intermediate waters dominate this mode of variability. The specific role played by ocean drag ODC is less clear. Drag retards and modifies the circulation at all depths, but its overall effect is complex and difficult to interpret (Killworth, 2003). We note that there is significant correlation (−0.39) between ODC and WSF. However, the emulator error increases from 0.0223 to 0.0274 when the EFPC emulator is rebuilt neglecting the ODC term, suggesting that the presence of ODC in the emulator does have a physical basis (i.e. which cannot be represented through inter-parameter correlations).

1st EOF of anthropogenic DIC
A similar analysis was performed for the anthropogenic change in DIC. The first EOF and the emulator coefficients for the change in DIC are plotted in Fig. 12c. The dominant mode of DIC variability (54 % of the ensemble variance) is strikingly similar to the second EOF of the industrial change in δ 13 C DIC (although it differs in the sense that it represents an inventory change, i.e. of constant sign everywhere, rather than surface-intermediate exchange). The corresponding PC emulators ( Fig. 12b and c) exhibit similar dependencies on model parameters, negatively correlated with the ocean drag (ODC) and wind scaling (WSF), suggesting a close physical relationship between these modes of variability, both driven by physical processes of advection and mixing. The anthropogenic change in the global DIC inventory exhibits a correlation of −0.70 with respect to the second PC of anthropogenic δ 13 C DIC , providing further support for this close relationship. A significant, though weaker, correlation (−0.36) is also apparent between anthropogenic DIC and the first PC of anthropogenic δ 13 C DIC , suggesting that air-sea gas exchange plays a relatively minor role in the uncertainty associated with the uptake of CO 2 . This is consistent with the relative lack of importance of ASG in the anthropogenic DIC emulator (Fig. 12c). The minor role for air-sea gas exchange reflects the fact that the characteristic air-mixed layer exchange time for CO 2 is short (∼ 1 yr) compared to the mean residence time of surface waters (Gruber et al., 1999). This is in contrast to the previously noted importance of the air-sea gas exchange rate for δ 13 C DIC (Fig. 12a), arising from carbonate buffering that leads to air-mixed layer equilibration times of ∼ 10 yr. We note that, despite their close relationship, the first PC of anthropogenic DIC is more difficult to emulate than the second PC of anthropogenic δ 13 C DIC (R 2 of 48 % and 77 %, respectively). This may reflect the fact that the DIC EOF is weakly affected by air-sea gas exchange, a process that is orthogonal to the dominant driver of thermocline mixing.
For completeness, we note that the Southern Ocean carbon sink is absent from the first EOF of anthropogenic DIC uptake. This signal is rather apparent in the third EOF (not shown) explaining 9 % of the ensemble variance. The third PC exhibits a correlation of −0.41 with the anthropogenic change in the global DIC reservoir. It is driven by ocean parameters (OHD, ODC) and sea ice diffusivity (SID), the latter presumably through its control on high-latitude deep water formation via brine rejection (Sect. 3.1.3).
In summary, the dominant mode of anthropogenic 13 C uptake variability (arising from uncertainties in air-sea gas exchange) is largely absent from DIC variability (instead arising from uncertainties in thermocline ventilation rates). The dominant processes controlling the distributions of industrial 13 C and CO 2 thus appear to be quite distinct. However, δ 13 C DIC in lower pycnocline and intermediate waters do bear the signature of anthropogenic DIC, as illustrated by the similarity between the second anthropogenic δ 13 C DIC EOF and the first anthropogenic DIC EOF.

Summary and Conclusions
We have described the development and evaluation of a large ensemble of precalibrated parameterisations of GENIE for application to coupled carbon cycle problems. These parameterisations have already been applied to a range of experiments (Eby et al., 2012;Holden et al., 2013;Joos et al., 2013;Zickfeld et al., 2013). Work is presently in progress to apply them to a transient ensemble of simulations over the most recent glacial cycle, an experiment that has the potential to be constrained by a rich compilation of watercolumn δ 13 C DIC reconstructions from marine sediment data (e.g. Oliver et al., 2010).
We have applied singular vector decomposition and principal component emulation  to investigate the modes of ensemble variability of the oceanic distribution of δ 13 C DIC , considering both the pre-industrial state and the industrial change. These analyses allow us to separate the natural and anthropogenic controls on the δ 13 C DIC distribution. The two analyses are connected by the possibility for compensating errors in the two distributions when modelling the modern ocean. For example, the second mode of variability in the pre-industrial state (Fig. 11b) has a similar spatial distribution to the dominant mode of variability of anthropogenic uptake (Fig. 12a), although they are controlled by quite distinct processes (surface-deep ocean exchange and air-sea gas exchange, respectively).
The analysis provides a quantitative (albeit partially subjective) demonstration of the well-known balance between physical and biogeochemical processes that control the vertical gradient of δ 13 C DIC (Tagliabue and Bopp, 2008). These same controls exert an influence on the global DIC inventory in the equilibrium state (as distinct from changes in the inventory due to anthropogenic forcing), supporting the use of the surface-deep gradient as a proxy for past ocean productivitydriven changes in atmospheric CO 2 (Shackleton et al., 1983).
Uncertainty in the distribution of δ 13 C DIC related to the Suess effect is dominated by air-sea gas exchange. The first EOF describes 63 % of the variability across the ensemble and is controlled almost entirely by the air-sea gas exchange parameter. The second EOF, describing 18 % of the variance, is driven by uncertainties in thermocline ventilation rates. In contrast, the decomposition of the distribution of anthropogenic DIC demonstrates that the air-sea gas exchange rate does not contribute substantially to uncertainty in the ocean carbon sink. However, the dominant mode of variability in DIC (54 %) is strikingly similar to the second mode of variability of anthropogenic 13 C uptake. Furthermore, these two modes of variability are controlled by the same parameter relationships and hence, presumably, by the same physical processes (controls on thermocline ventilation rates).
This comparison highlights well-known complications that arise when applying δ 13 C DIC as a constraint on the ocean sink and the need to separate the effects of uncertainty due to air-sea gas exchange. However, the similarity of the modes of variability of DIC and δ 13 C DIC after removal of the airsea gas imprint strongly supports the value of observationally based constraints imposed by δ 13 C DIC on DIC uptake, e.g. Quay et al. (2003). (Note these methodologies correct for the air-sea gas exchange signal.) The need for a precise separation of the influence of air-sea gas exchange suggests careful targeting of δ 13 C DIC observations, with emphasis on lower pycnocline and intermediate waters is essential for improved quantification of the ocean sink.