Multidecadal accumulation of anthropogenic and remineralized dissolved inorganic carbon along the Extended Ellett Line in the northeast Atlantic Ocean

Marine carbonate chemistry measurements have been carried out annually since 2009 during UK research cruises along the Extended Ellett Line (EEL), a hydrographic transect in the northeast Atlantic Ocean. The EEL intersects several water masses that are key to the global thermohaline circulation, and therefore the cruises sample a region in which it is critical to monitor secular physical and biogeochemical changes. We have combined results from these EEL cruises with existing quality‐controlled observational data syntheses to produce a hydrographic time series for the EEL from 1981 to 2013. This reveals multidecadal increases in dissolved inorganic carbon (DIC) throughout the water column, with a near‐surface maximum rate of 1.80 ± 0.45 µmol kg−1 yr−1. Anthropogenic CO2 accumulation was assessed, using simultaneous changes in apparent oxygen utilization (AOU) and total alkalinity (TA) as proxies for the biogeochemical processes that influence DIC. The stable carbon isotope composition of DIC (δ13CDIC) was also determined and used as an independent test of our method. We calculated a volume‐integrated anthropogenic CO2 accumulation rate of 2.8 ± 0.4 mg C m−3 yr−1 along the EEL, which is about double the global mean. The anthropogenic CO2 component accounts for only 31 ± 6% of the total DIC increase. The remainder is derived from increased organic matter remineralization, which we attribute to the lateral redistribution of water masses that accompanies subpolar gyre contraction. Output from a general circulation ecosystem model demonstrates that spatiotemporal heterogeneity in the observations has not significantly biased our multidecadal rate of change calculations and indicates that the EEL observations have been tracking distal changes in the surrounding North Atlantic and Nordic Seas.


Introduction
Emissions of carbon dioxide (CO 2 ) from human activities have increased the atmospheric partial pressure of CO 2 (pCO 2 ), in particular during the last 200 years [Ahn et al., 2012]. This increase, and its well-documented implications for global climate [IPCC, 2013], would have been significantly greater without CO 2 uptake by the ocean, which currently sequesters about a quarter of anthropogenic CO 2 emissions each year [Le Quéré et al., 2009]. Ocean CO 2 uptake also induces a decline in ocean pH, commonly known as ocean acidification, which will persist for centuries after anthropogenic CO 2 emissions cease [Caldeira and Wickett, 2003]. Ocean acidification will have repercussions for marine ecosystems and biogeochemistry that we have only recently begun to understand [Doney et al., 2009;Gaylord et al., 2015].
Open ocean time-series sites that monitor marine carbonate chemistry provide essential observational data to quantify long-term trends in anthropogenic CO 2 uptake and acidification [e.g. Dore et al., 2009;Olafsson et al., 2009;González-Dávila et al., 2010;Bates et al., 2012]. The time-series data are also used to validate output from global coupled ocean-atmosphere models [Le Quéré et al., 2010]. However, only a handful of these sites exist globally [Bates et al., 2014]. We present a new time-series of marine carbonate chemistry measurements for the Extended Ellett Line (EEL), an open ocean transect in the northeast Atlantic Ocean. The EEL runs from Iceland to Scotland via the Rockall Plateau ( Figure 1), and repeated physical measurements have been carried out on parts of it since 1975 [Holliday and Cunningham, 2013]. The transect captures the flow of warm, salty water from the North Atlantic into the Nordic Seas, and around half of the returning deep, cold overflow current. The remaining overflow returns south via the west of Iceland [Hansen and Østerhus, 2000]. The North Atlantic is an important region to monitor because of its global ©2016 American Geophysical Union. All rights reserved.
importance for oceanic uptake and accumulation of anthropogenic CO 2 , accounting for about 23 % of global oceanic anthropogenic CO 2 storage despite covering only 15 % of the global ocean surface area Khatiwala et al., 2009Khatiwala et al., , 2013. As the EEL will continue to be surveyed by UK research vessels, our analysis provides a baseline within this critical region that can be extended in future years.
Our time-series consists of measurements carried out during annual EEL cruises from 2009 to 2013, augmented by hydrographic data from several quality-controlled compilations Schmittner et al., 2013]. We have quantified the rate of change of dissolved inorganic carbon (DIC) throughout the water column along the EEL, and used simultaneous changes in apparent oxygen utilization (AOU) and total alkalinity (TA) to quantify its anthropogenic (DIC anth ) and biogeochemical components. The approach that we have taken to partition the changes in total DIC into these components is based on the same principles as established back-calculation methods for estimating DIC anth [Brewer, 1978;Chen and Millero, 1979;Gruber et al., 1996]. However, we have determined the relative accumulation of DIC anth during the observational time period (i.e. 1981 to 2013), rather than seeking to evaluate the total DIC anth increase since pre-industrial times. This avoids the necessity to estimate pre-industrial fields for DIC, AOU, TA and other variables, which is a key source of uncertainty in back-calculation methods [Matsumoto and Gruber, 2005;Sabine and Tanhua, 2010]. The extended multi-linear regression (eMLR) technique is an alternative way to evaluate relative changes in DIC anth [e.g. Friis et al., 2005;Tanhua et al., 2007], but it does not provide information about non-anthropogenic changes in DIC, which we find to be significant at the EEL. The eMLR technique may also be inappropriate due to our study's multi-decadal duration [Goodkin et al., 2011]. We have instead determined multi-decadal rates of change of the relevant hydrographic variables with linear regressions using all of the data, thereby producing results that are not strongly biased by any individual cruise. These ©2016 American Geophysical Union. All rights reserved. regressions were carried out on constant potential density surfaces, in order to track water masses between cruises [Pérez et al., 2010;Wanninkhof et al., 2010]. Observations of the stable carbon isotopic composition of DIC (δ 13 C DIC ) provided independent support for our deconvolution of the changes in DIC, again through regression-derived rates of change and a process-based approach. Finally, using output from a coupled ocean general circulationbiogeochemical model [Yool et al., 2013a], we demonstrated that the calculated rates of change along the EEL do not appear to be biased by the spatiotemporal heterogeneity of the observations and track equivalent changes in the surrounding North Atlantic and Nordic Seas.
The model study emphasises the importance of the EEL for continued future monitoring.
©2016 American Geophysical Union. All rights reserved.

Recent EEL cruises
The EEL was occupied annually from 2009 to 2013 by RRS Discovery cruises D340 [Sherwin, 2009], D351 [Read, 2010], D365 [Read, 2011] and D379 [Griffiths, 2012], and RRS James Cook cruise JC086 [Griffiths and Holliday, 2013]. During these cruises, seawater samples for DIC and TA were collected and measured as described in the Supporting Information (Text S1), and CTD sensor measurements of temperature, practical salinity, and dissolved oxygen (DO) were carried out and calibrated using manual measurements of discrete samples (Text S2). During cruise D379, samples were also collected to measure δ 13 C DIC (Text S3), as detailed by Humphreys et al. [2015]. The distributions of these variables from cruise D379 are illustrated in the Supporting Information ( Figure S1).

Syntheses
Data from the GLobal Ocean Data Analysis Project (GLODAP)  and CARbon dioxide IN the Atlantic Ocean (CARINA)  syntheses were combined with the measurements undertaken during the recent EEL occupations. Carbonate chemistry data in GLODAP from the Transient Tracers in the Ocean -North Atlantic Study (TTO-NAS) [Brewer et al., 1985] were adjusted following Tanhua and Wallace [2005].
The δ 13 C DIC data in GLODAP and CARINA  have not undergone a secondary quality-control process, so we instead used the compilation prepared by Schmittner et al. [2013]. This consists of a high-quality subset of the GLODAP and CARINA results, augmented by data from additional cruises [Gruber et al., 1999]. The estimated accuracy of these δ 13 C DIC values is between 0.1 and 0.2 ‰ [Schmittner et al., 2013]. ©2016 American Geophysical Union. All rights reserved.

Bathymetry
Bathymetric data from the GEBCO_2014 30 arc-second grid (version 20141103, http://www.gebco.net) were obtained for the EEL and its immediate surrounding area. The bathymetry of the idealized EEL route (Table S1) was derived from this data by linear interpolation of depth from the GEBCO_2014 latitude and longitude grids.

Model output
We obtained the output of a simulation described in detail by Yool et al. [2013a] (referred to as the 'anthropogenic simulation'). This had been run from the year 1860 to 2100, and consisted of the size-based intermediate complexity ecosystem model MEDUSA-2.0 [Yool et al., 2013b] coupled to the physical model version 3.2 of the Nucleus for European Modelling of the Ocean (NEMO) [Madec, 2008]. The horizontal resolution is approximately 1° × 1° (with 292 × 362 grid points), and vertical space is divided into 64 levels that increase in thickness from about 6 m at the surface to 250 m at a depth of 6 km. We refer to each vertical column of grid points as a 'model station'. Surface forcing of NEMO used output from the HadGEM2-ES Earth-System model [Collins et al., 2011], and the DIC and TA fields in MEDUSA-2.0 were initialized using the GLODAP climatology .
Atmospheric pCO 2 followed historical data from 1860 through 2005, and then switched to Representative Concentration Pathway 8.5 [Riahi et al., 2011] for the rest of the simulation.
We also obtained output from a second 'control simulation', which had the same set-up except that the mean atmospheric pCO 2 was held constant at a pre-industrial value of 286 µatm throughout.
©2016 American Geophysical Union. All rights reserved.

Observations
An 'observational' dataset (Table S2) was created using all data in GLODAP, CARINA and from the recent EEL cruises from within 167 km of an idealized EEL route (Figure 1). The route runs in straight lines (great circles, Figure 1) through the waypoints listed in Table S1.
Many of the cruises in GLODAP and CARINA passing through the EEL region did not follow the EEL route, so 167 km was chosen as the optimal radius to satisfy the trade-off between capturing sufficient historical data to perform an effective analysis whilst remaining local to the EEL. Essentially, 167 km was the smallest possible distance that included enough historical data from the earliest time points (in 1981) to perform a robust analysis. The same processing was carried out separately for the Schmittner et al. [2013] dataset plus the D379 δ 13 C DIC measurements [Humphreys et al., 2015], to create the 'isotopes' dataset (Table S3).

Model output
The model outputs were subsampled into several different datasets (Table 1). Firstly, monthly mean fields from both simulations (anthropogenic and control) were subsampled to match the spatiotemporal distribution of observational data as closely as possible, using a nearestneighbour approach. These datasets are hereafter referred to as 'Subsampled Anthropogenic Monthly' and 'Subsampled Control Monthly' (SAM and SCM respectively). The SAM dataset is therefore the model equivalent of the real EEL observations, and SCM is the same but with no anthropogenic CO 2 . Secondly, annual mean fields from the anthropogenic simulation, from 1981 to 2013 inclusive and within the region from 25°N to 75°N and 70°W to 10°E, were extracted to form the 'Full Anthropogenic Annual' dataset (FAA). Finally, all output at the closest model locations to the idealized EEL transect route (Figure 1) was ©2016 American Geophysical Union. All rights reserved. selected from FAA to form the 'Transect Anthropogenic Annual' dataset (TAA). This TAA dataset can be considered to represent the EEL as if it had been sampled perfectly throughout the study period. The rates of change of variables calculated using TAA are therefore the standard against which the quality of the other model datasets are judged.

Derived variables
The potential density anomaly at pressure (P) = 0 dbar (σ 0 ), in situ density (ρ), and potential temperature (θ) were calculated from T, S and P using the Gibbs-SeaWater Oceanographic Toolbox for MATLAB® (MathWorks®, USA) [McDougall and Barker, 2011]. Apparent oxygen utilization (AOU) was calculated from θ, S and DO using the combined fit coefficients of García and Gordon [1992]. The Revelle factor was calculated from DIC, TA, T, S and P in the surface ocean, assuming negligible silicate and phosphate concentrations, using version 1.1 of the CO 2 SYS program for MATLAB [van Heuven et al., 2011], the carbonic acid dissociation constants of Lueker et al. [2000], and the boron:chlorinity of Lee et al. [2010].

Interpolations
At each sampling station in each dataset, DIC, TA, AOU, S, δ 13 C DIC and depth were interpolated to σ 0 values ascending in units of 0.001 from 26 to 28 kg m -3 (called 'σ 0 levels'), using piecewise cubic Hermite interpolating polynomial (PCHIP) fits to the observations (e.g. Figure 2(a)-(c)) [Fritsch and Carlson, 1980;Kahaner et al., 1988]. This was in order to compare these variables between cruises. Potential density is a better interpolant than depth in this context, because it tracks any vertical movements of water masses in the time between successive observations. A small number of stations had fewer than the 4 unique measurements required to carry out the interpolation, so the measured values were instead ©2016 American Geophysical Union. All rights reserved. assigned directly to their closest σ 0 levels. No extrapolations were performed beyond the measured σ 0 range at any station.

Rates of change
For the observational, SAM, SCM and TAA datasets, ordinary least-squares regressions between each variable and the sampling date aacross all sampling stations were used to determine the rate of change at each σ 0 level (e.g. Figure 2(d)-(f)). The mean value of each variable was also calculated for each σ 0 , again across all sampling stations. We report the rate of change of a variable X at any given σ 0 as dX/dt = R ± U, where R is the rate of change of X, and U is its 1σ uncertainty taking into account any autocorrelation in X. These calculations were carried out using the MATLAB® function 'regress2' written by I. Eisenmann to 030°W for the North Atlantic, and 66 to 72°N and 012°W to 001°E for the Nordic Seas.
©2016 American Geophysical Union. All rights reserved.

Using TA and AOU
We can use changes in other observed variables to deconvolve the total DIC change into its component driversthe carbonate pump (DIC carb ), soft tissue pump (DIC soft ), and the solubility pump (DIC sol ) [Gruber et al., 1996]: (1) The 'carbonate pump' is the formation and dissolution of calcium carbonate (CaCO 3 ).
Increasing its rate of dissolution relative to formation at a given point drives an increase in DIC, accompanied by an increase in TA of double the magnitude [Wolf-Gladrow et al., 2007]. We can therefore determine its contribution to the total dDIC/dt as: ( 2) where the dAOU/dt term corrects for changes in TA driven by nitrate release during organic matter remineralization, and R N/O2 is -0.0941 ± 0.0081 [Anderson and Sarmiento, 1994].
Biological activity, concentrated near the ocean surface, converts dissolved inorganic nutrients to particulate organic matter (POM), some of which sinks and remineralizes at depth, returning the nutrients to solution: the soft tissue pump. Remineralization also takes up DO, thereby increasing AOU. The component of dDIC/dt caused by changes in organic matter remineralization can therefore be predicted from dAOU/dt: where R C/O2 is the increase in DIC as a fraction of DO consumption during this process, which we assume takes a constant value of -0.688 ± 0.092 [Anderson and Sarmiento, 1994].
The remaining DIC increase is attributed to increases in air-to-sea CO 2 transfer at the surface ©2016 American Geophysical Union. All rights reserved.
outcrop regions for these σ 0 levels. Assuming no significant long-term trend in air-sea pCO 2 disequilibrium from 1981 to 2013 in the ventilation regions (i.e. dDIC diseq /dt = 0), this increase represents the accumulation of anthropogenic DIC (DIC anth ): (4) 3.5.2. Using DIC stable isotopes We can relate dδ 13 C DIC /dt to changes in the other variables in order to independently test our attribution of dDIC/dt to its components. At each σ 0 level, the total change in δ 13 C DIC from 1981 to 2013 (Δδ 13 C DIC ) is the sum of the changes caused by the same components that drove the total changes in DIC: Formation and dissolution of CaCO 3 minerals does not cause any significant fractionation of carbon isotopes [Romanek et al., 1992;Lynch-Stieglitz et al., 1995;Gruber et al., 1999], so the Δδ carb term is set to 0. The Δδ anth can be calculated from the total changes in DIC (ΔDIC) and AOU (ΔAOU), and the ratio between anthropogenic changes in δ 13 C DIC and DIC (called ΔRC following e.g. McNeil et al. [2001]): The remineralization component Δδ soft also depends on the initial DIC and δ 13 C DIC (DIC i and δ i respectively), and the isotopic composition of POC (δ 13 C POC ): ©2016 American Geophysical Union. All rights reserved.
The values of the rate-of-change regression lines for DIC and δ 13 C DIC at the mid-point of the year 1981 were used for DIC i and δ i . For any variable X, conversion between ΔX and its rate of change is straightforward for the 32-year observational period: Combining (6) and (7) into (5) and rearranging, we obtained the following relationship for each σ 0 level: Considering the terms Δδ 13 C DIC , δ i , DIC i , ΔDIC and ΔAOU to represent column vectors in which each row corresponds to a different σ 0 level, and ΔRC and δ 13 C POC as unknown scalar constants, (9) was further rearranged and rewritten in matrix form: where the open circle symbols denote the Hadamard (element-wise) product. For the system of linear equations thus generated, the least-squares best solution for ΔRC and δ 13 C POC across all σ 0 levels was then determined. The 95 % bootstrap confidence intervals for ΔRC and δ 13 C POC were calculated using the bias corrected and accelerated percentile method with 10 4 bootstrap resamples.
©2016 American Geophysical Union. All rights reserved.

Column inventories
To volume-integrate the rates of change for each variable, we first calculated the average depth of each σ 0 level and determined its lateral extent using the GEBCO_2014 bathymetric grid. The column inventories were then determined from the σ 0 depths, widths and rates of change. The approach is described in detail in the Supplementary Information (Text S5).

Water column changes in the observations
We observe increases in DIC from 1981 to 2013 throughout the water column. The maximum dDIC/dt of 1.80 ± 0.45 µmol kg -1 yr -1 is found in the upper 30 m of the water column (σ 0 ≈ 26.7 kg m -3 , Figures 3(a) and 4(a)). This corresponds to an increase in seawater pCO 2 of about 3.6 μatm yr -1 (calculated from the mean Revelle factor for this σ 0 layer). This is greater than the atmospheric pCO 2 increase rate of about 1.6 μatm yr -1 [Tjiputra et al., 2014], which supports some previous studies that suggest that the oceanic sink for atmospheric CO 2 has been decreasing in this region [e.g. . Below the surface layer, dDIC/dt decreases to a deep minimum of 0.02 ± 0.10 µmol kg -1 yr -1 (σ 0 ≈ 27.9 kg m -3 , c. 2 km and deeper, Figures 3(a) and 4(a)). The non-zero dDIC/dt means that the carbonate, soft tissue and solubility pump processes controlling DIC [Gruber et al., 1996] are not operating in a steady state.
Changes in TA in the observations are very small: dTA/dt is between 0.23 ± 0.26 and -0.19 ± 0.10 µmol kg -1 yr -1 , and consequently dDIC carb /dt is in the range between 0.16 ± 0.13 and -0.09 ± 0.05 µmol kg -1 yr -1 (Figures 3(a) and 4(b)). Therefore the observed DIC increase was not significantly driven by changes in the carbonate pump.
©2016 American Geophysical Union. All rights reserved.
Changes in AOU are also virtually zero near the surface, but there have been increases in AOU deeper in the water column. As a result, dDIC soft /dt closely tracks dDIC/dt for σ 0 > 27.45 kg m -3 (Figures 3(a) and 4(c)). This component of dDIC/dt can be attributed to a multi-decadal increase in the amount of remineralized organic matter at these σ 0 levels. Two possible mechanisms could explain this phenomenon.
Firstly, there could have been an increase in export and remineralization of POM at the EEL itself. The indirect evidence does not support this hypothesis. Although there have not been sufficient observations to directly confirm the presence or absence of a multi-decadal trend in POM export and remineralization in the EEL region, POM export is unlikely to be increasing fast enough to cause the observed pattern in dDIC soft /dt (as derived from dAOU/dt). Export rates can be estimated as a function of surface chlorophyll-a, with higher concentrations accompanying higher export rates [e.g. Dunne et al., 2007]. However, satellite observations have detected a small decline in chlorophyll-a from 1998 to 2012 for the northern North Atlantic [Gregg and Rousseaux, 2014], which is inconsistent with increasing in situ POM export.
Alternatively, there could have been changes in the lateral distribution of water masses along isopycnals, bringing waters with higher AOU (more remineralized POM) into the EEL region. At σ 0 levels lighter than about 27.7 kg m -3 , the EEL samples a mixture of waters from the subtropical and subpolar gyres. Contraction of the subpolar gyre increases the subtropical component, while expansion decreases it. The subpolar gyre index (SPGI) metric can be interpreted as a measure of subpolar gyre contraction, with lower values indicating a more contracted gyre [Hátún et al., 2005]. Overall, there has been a decrease in the SPGI during the period of our study, especially since the early 1990s [Häkkinen and Rhines, 2004;Hátún et al., 2005;Hughes et al., 2012]. This phenomenon has separately been shown to drive declining macronutrient concentrations in the Rockall Trough [Johnson et al., 2013]. The ©2016 American Geophysical Union. All rights reserved.
subtropical waters influencing the EEL in this way are a combination of Eastern North Atlantic Water, formed in the Bay of Biscay [McGrath et al., 2012a], and highly saline Mediterranean Water [Burkholder and Lozier, 2011;McGrath et al., 2012a]. Data from GLODAP and CARINA  indicate that, to first order, DIC and AOU increase to the south of the EEL region along σ 0 levels in these water masses. An increasing southern influence on the water at the EEL would therefore be expected to increase DIC and AOU at the EEL, like we observe.
The remaining DIC increase is interpreted as anthropogenic, and is confined in and above the thermocline (σ 0 < 27.5 kg m -3 , Figures 3(a) and 4(d)). This matches previous global-scale studies of the DIC anth and anthropogenic tracer distributions [e.g. Sabine et al., 2004]. In this upper part of the water column, our calculated dDIC anth /dt is consistent with similar analyses in the nearby or overlapping regions of the Iceland Basin [Pérez et al., 2010] and southern Rockall Trough [McGrath et al., 2012b]. At greater depths, dDIC anth /dt is virtually zero.
Between a σ 0 of 27.70 and 27.85 kg m -3 , the EEL samples Labrador Sea Water (LSW). In the EEL region, the properties of LSW are highly variable both spatially and temporally. This is because LSW undergoes extensive mixing with other water masses, including recirculating LSW ventilated in earlier years, during its transport from the Labrador Sea to the northeast Atlantic [Yashayaev et al., 2007]. So, although we do not observe an increase in DIC anth in the LSW during the study period, any anthropogenic signal from the source region could have been suppressed by mixing. However, we observed a small positive dDIC/dt in LSW (c. 0.3 μmol kg -1 yr -1 ). It is unlikely that this was caused by increased in situ remineralization, not only for the reasons already described for waters nearer the surface, but also because virtually all POM that is generated in the EEL region is remineralized within the mesopelagic zone, which is shallower than the LSW [Henson et al., 2012]. Other studies have identified increases in LSW DIC anth during time periods and within regions that were similar, but ©2016 American Geophysical Union. All rights reserved.
The dDIC anth /dt discrepancy is likely to be an artefact of differences between studies in the distribution of observations of this highly variable water mass.
We do not observe any DIC anth accumulation in the water masses below the LSW, in the deepest part of the water column (σ 0 > 27.85 kg m -3 ), although there are small increases in DIC and DIC soft (Figure 3(a)). These waters, like the LSW, have been significantly altered by mixing since their formation in the Nordic Seas, and subsequent flow through one of several narrow channels over the Greenland-Scotland ridge and into the EEL region [Hansen and Østerhus, 2000]. This prevents us from directly attributing the DIC and AOU increases to a specific driver.

Choice of POM stoichiometry
The choice of POM stoichiometry controls the values of R C/O2 and R N/O2 , which clearly influences our partitioning of dDIC/dt into its carbonate (2) and soft tissue components (3).
We selected values for R C/O2 and R N/O2 of -0.688 ± 0.092 and -0.0941 ± 0.0081 respectively, which are based on global macronutrient measurements [Anderson and Sarmiento, 1994].
These feature a higher O 2 coefficient than the 'original' stoichiometry of Redfield et al. [1963], which gives -0.768 for R C/O2 and -0.116 for R N/O2 . This higher O 2 coefficient is supported by considerations of the composition of several groups of algal biomolecules [Anderson, 1995]. Although switching between these stoichiometries does create a systematic offset in the results, the size of this offset is no larger than the random uncertainty inherent in the calculations. Switching creates a difference of 0.05 µmol kg -1 yr -1 in the mean dDIC anth /dt across all σ 0 levels, and a difference in the column inventory changes for the entire EEL of about 8 %, with the newer stoichiometry [Anderson and Sarmiento, 1994] giving a higher DIC anth inventory. For comparison, the random uncertainty propagated into this inventory ©2016 American Geophysical Union. All rights reserved. estimate from the rates of change themselves is about 9 % of the mean value. Furthermore, the original value for R C/O2 [Redfield et al., 1963], which has a much greater influence on the DIC anth calculation than R N/O2 does, falls within the stated uncertainty of the more recent result [Anderson and Sarmiento, 1994]. Consequently, we do not consider this choice to be a particularly important source of uncertainty in our final results.

Column inventories
The global ocean anthropogenic CO 2 sink was about 2 Pg-C yr -1 for the period from 1981 to 2013 [Le Quéré et al., 2009], which corresponds to a global mean DIC anth accumulation rate of about 1.5 mg-C yr -1 m -3 [Eakins and Sharman, 2010]. For the idealized EEL route, the column inventory C(dDIC anth /dt) is 2.8 ± 0.4 mg-C m -3 yr -1 , which is about double the global average value. The equivalent C(dDIC/dt) is 9.0 ± 1.0 mg-C m -3 yr -1 , so the DIC anth increase accounts for only 31 ± 6 % of the total DIC accumulation. Virtually all of the remainder is a result of the increased remineralized organic matter, contained in the increased supply of southern-sourced waters that have been brought into the region by contraction of the subpolar gyre, as discussed in Section 4.1.1.
The EEL region is part of the largest source of DIC into the Nordic Seas through advection of Atlantic waters over the Greenland-Scotland ridge [Jeansson et al., 2011]. Presently, the Nordic Seas are an important sink for anthropogenic CO 2 , convectively transporting it from the surface layer into the interior and then returning it back into the deep North Atlantic [Jutterström et al., 2008;. Increasing DIC concentrations in the North Atlantic waters prior to their transport over the ridge might therefore hinder the efficiency of the Nordic Seas CO 2 sink, by inhibiting further uptake of atmospheric CO 2 across the air-sea interface. However, the impact of this effect may be limited, as much of the DIC anth ©2016 American Geophysical Union. All rights reserved.
transported into the ocean interior in the Nordic Seas arrives in the surface ocean through advection, and is not taken up locally by air-sea exchange [Olsen et al., 2006].

Measurements of chlorofluorocarbon (CFC) inventory changes in the North Atlantic have
demonstrated that CFC column inventory variability can be dominantly controlled by changes in volume of different water masses at a given location, rather than changes in the CFC concentration within each water mass [Kieke et al., 2007;Steinfeldt et al., 2009]. A similar conclusion has also been suggested for anthropogenic DIC [Pérez et al., 2010].
However, the timescale of volumetric variability in these studies is sub-decadal. For our longer, multi-decadal study period from 1981 to 2013, we find that changes in σ 0 layer thicknesses are negligible, and their inclusion in the inventory calculations would change the final result by an order of magnitude less than its uncertainty.

Stable isotopes of DIC
There are far fewer δ 13 C DIC observations than for the other variables, and the length of the time-series is shorter, running only from 1993 to 2012. Nevertheless, changes in δ 13 C DIC can be used as an independent test of our attribution of the changes in DIC to its anthropogenic and soft-tissue components, as these are the two main processes influencing δ 13 C DIC in the interior ocean. First, uptake of anthropogenic CO 2 results in a decrease in δ 13 C DIC , known as the Suess effect, because fossil-fuel carbon is isotopically light relative to modern seawater [Keeling, 1979]. Second, particulate organic carbon (POC) remineralization also decreases δ 13 C DIC , because POC has a lighter isotopic signature than typical seawater. According to a meta-analysis study, surface ocean δ 13 C POC is between -20 and -30 ‰ at the latitude and sea surface T of the EEL [Goericke and Fry, 1994], compared with typical seawater DIC values near 0 ‰ [Olsen and Ninnemann, 2010;Schmittner et al., 2013;Humphreys et al., 2015].
The carbonate pump does not significantly affect δ 13 C DIC , because marine carbonate mineral ©2016 American Geophysical Union. All rights reserved. formation (calcification) does not significantly fractionate carbon. Carbonate minerals usually have a δ 13 C composition similar to that of the surrounding seawater [Romanek et al., 1992;Lynch-Stieglitz et al., 1995;Gruber et al., 1999].
We observed negative dδ 13 C DIC /dt values throughout the water column. The magnitude of dδ 13 C DIC /dt decreased from a maximum of -0.038 ± 0.026 ‰ yr -1 at near-surface σ 0 levels to a minimum of -0.002 ± 0.006 ‰ yr -1 at depth ( Figure 5). Quay et al. [2007] used a multi-linear regression approach to identify a mean dδ 13 C DIC /dt of -0.018 ± 0.002 ‰ yr -1 for the entire Atlantic Ocean surface mixed layer from 1981 to 2003. They also found that dδ 13 C DIC /dt increased to between -0.04 and -0.05 ‰ yr -1 in the subpolar region between 40°N and 60°N, due to a combination of changes in water mass properties and anthropogenic CO 2 uptake.
Their findings are consistent with our near-surface results for the EEL.
We deconvolved the dδ 13 C DIC /dt distribution into anthropogenic and remineralized components, which are controlled by the variables ΔRC and δ 13 C POC . The least-squares best fit solutions of (10) for ΔRC and δ 13 C POC across all σ 0 levels were -0.0166 ± 0.0003 ‰ (µmol kg -1 ) -1 and -27.0 ± 0.5 ‰, respectively. To visualize the results, dδ 13 C DIC /dt was predicted using (9), with the best fit values of ΔRC and δ 13 C POC , and the observed rates of dDIC/dt and dAOU/dt ( Figure 5). It is inevitable that mean value of the best fit dδ 13 C DIC /dt profile will match that of the observations, because of how the ΔRC and δ 13 C POC were determined. However, if there were elements of the observed dδ 13 C DIC /dt profile that were not driven by DIC anth or DIC soft inputs (e.g. driven by DIC carb ), then we would expect the shape of the predicted profile to deviate from the observations in the relevant σ 0 range. This does not occur and hence we conclude that DIC anth and DIC soft inputs are indeed the dominant drivers of the observed interior δ 13 C DIC changes. ©2016 American Geophysical Union. All rights reserved.
It was originally proposed that ΔRC might take a relatively globally-uniform value between -0.016 and -0.019 ‰ (μmol kg -1 ) -1 [Heimann and Maier-Reimer, 1996]. More recently, it has been demonstrated that ΔRC can deviate from this global average to exhibit significant spatial variation in certain regions [e.g. McNeil et al., 2001;Olsen et al., 2006], because the air-sea equilibration time is an order of magnitude faster for DIC than for δ 13 C DIC [Lynch-Stieglitz et al., 1995]. Körtzinger et al. [2003] calculated ΔRC throughout the North Atlantic, reporting a value of -0.022 ± 0.002 ‰ (μmol kg -1 ) -1 for the σ 0 levels observed in this study, while Olsen et al. [2006] found a wide range of ΔRC between about 0.00 and -0.03 ‰ (μmol kg -1 ) -1 in the upper 100 m of the nearby Nordic Seas. Our best fit value for ΔRC, -0.0166 ± 0.0003 ‰ (µmol kg -1 ) -1 , is at the lower end of, but consistent with, these published results.
We next evaluate our value for δ 13 C POC relative to previous estimates in a similar way. Using the linear regression between sea surface T and δ 13 C POC proposed by Goericke and Fry [1994] (for the northern hemisphere and T > 5 °C), and a value of 11 °C for T (the mean of all EEL T observations for which P ≤ 10 dbar), we would predict δ 13 C POC = -23 ± 4 ‰ for the EEL (we have estimated the uncertainty by eye from the figures presented by Goericke and Fry [1994]). Congruently, a more recent global compilation of δ 13 C POC results reports δ 13 C POC in the approximate range from -18 to -27 ‰ for the latitude range of the EEL [Young et al., 2013]. Thus our least-squares solution for δ 13 C POC of -27.0 ± 0.5 ‰ is concordant with these and other published values [e.g. Rau et al., 1997].
We conclude that the dδ 13 C DIC /dt observations provide independent support for our quantitative attribution of dDIC/dt to anthropogenic and remineralization components.
©2016 American Geophysical Union. All rights reserved.

Subsampled model output and observational data
Before discussing the rates of change calculated from the model datasets, we first assess how well the distributions of the absolute values of the modelled variables agreed with the observations. For this, we will use the SAM and SCM datasets, which have been subsampled to match the spatiotemporal distribution of the observational data. It is not necessary for these absolute value distributions to be identical in order to compare rates of change between the different datasets. However, if the model distributions were to diverge significantly from the observations for reasons that could not be explained, then the utility of the model as an analogue to the real world would be severely limited. An important caveat is that only one model has been used here, and others might result in different outcomes.
There was no significant systematic offset between the latitude, longitude and date of the observations and their matching points in the SAM and SCM datasets (Figure 6 unrealistic northward penetration of Antarctic Bottom Water (AABW) in the NEMO run used here [Yool et al., 2013a], combined with the model tendency to underestimate the density of this AABW [Heuzé et al., 2013]. The SAM and SCM AOU fields are also very similar to each other, and both represent their matching observations relatively well (Figure 6(f)). In SAM, DIC takes consistently high values relative to SCM (Figure 6(g)), so as expected some anthropogenic CO 2 should be detectable in the SAM dataset. The DIC in SAM is also consistently high relative to the observations, but the offset is fairly consistent across the entire DIC range, with the fit quality otherwise similar as for AOU. The TA fields from both SAM and SCM diverge considerably from reality, covering a much wider and higher range of values ( Figure 6(h)). However, as they have similar distributions to each other, this should not adversely affect identification of the anthropogenic CO 2 signal in SAM.
In terms of rates of change, SAM does represent the pattern for the observations in the upper water column reasonably well (Figure 3(b)). Its dDIC anth /dt decreases away from surface, where it takes values close to 1.0 µmol kg -1 yr -1 , to effectively 0 at a σ 0 of about 27.2 kg m -3 , in agreement with the observations. The rates dDIC/dt and dDIC soft /dt are similarly wellmatched. Deeper in the water column, between mean depths of about 1 and 2 km, there is a small increase of about 0.5 µmol kg -1 yr -1 in DIC anth in SAM. This is absent from the observations, we find a similar pattern in the same depth range for DIC anth in SCM ( Figure   3(d)), indicating that it may be due to model drift in the absence of a long spin-up period.
Alternatively, it may be associated with the high northward AABW penetration that we identified as a possible cause of relatively low σ 0 in the models [Hieronymus and Nycander, 2013;Heuzé et al., 2014]. Otherwise, SCM does not show significant changes in any of the tested variables, and the patterns with depth appear mostly random. The dDIC carb /dt, which is mostly dependent upon dTA/dt (2), exhibits changes in SAM that are absent from the observations throughout the water column. However, like for DIC anth in the deeper part of the ©2016 American Geophysical Union. All rights reserved.
water column, we find a similar pattern in SCM, again suggesting that it may result from model drift.

Spatiotemporal sampling heterogeneity
The model datasets can be used to estimate the uncertainty introduced into the observational rates of change from spatiotemporal heterogeneity in the data distribution. To do this, we compared the rates of change calculated for the anthropogenic simulation subsampled to match the observations (SAM, Figure 3(b)) with those from the same simulation but with no missing values (TAA, Figure 3(c)). The mean ± standard deviation of the differences in rates of change between SAM and TAA across all σ 0 are: -0.01 ± 0.14 for DIC, 0.25 ± 0.36 for TA, -0.27 ± 0.18 for AOU and +0.06 ± 0.17 for DIC anth , all in µmol kg -1 yr -1 . These differences, particularly for DIC anth , are up to an order of magnitude smaller than typical uncertainties in the rates of change themselves, indicating that the spatiotemporal heterogeneity of the observational data distribution has not adversely affected the calculated rates of change for these variables.

Sub-decadal variability
It has been separately shown using observational data that multi-decadal trends in DO (and therefore AOU) can be identified despite substantial short-term interannual variability in a shorter, 19-year time-series transect close to the EEL, which samples several water masses also present at the EEL [Stendardo et al., 2015]. However, difficulties are presented over shorter timescales, as described below. It has been suggested that a higher rate of DIC anth accumulation can be identified in the Iceland Basin during the high North Atlantic Oscillation (NAO) index period from 1991 to 1998, compared with the lower NAO index (NAOI) period from 1997 to 2006 [Pérez et al., 2010]. The NAOI can be defined in several different ways, all associated with the atmospheric pressure difference between Iceland and the Azores, with ©2016 American Geophysical Union. All rights reserved. a more positive NAOI indicating a greater difference in pressure [Hurrell et al., 2003]. This pressure difference affects the local atmospheric circulation and surface wind speeds, and consequently can influence surface ocean currents and air-sea gas exchange [Thomas et al., 2008;Gruber, 2009].
To test for any NAOI signal in our data, we calculated dDIC anth /dt using the same methods as for the observational SAM and TAA datasets, but restricted to these two shorter date ranges (1991 to 1998 and 1997 to 2006). We find greater dDIC anth /dt for the latter, low-NAOI period in the observations, an opposite result to Pérez et al. [2010], although part of their DIC anth increase was due to changing σ 0 layer volumes rather than changes within σ 0 layers. More importantly, our calculated rates are barely distinguishable from uncertainties, because fewer data are available for shorter time periods, so the statistical significance of any apparent nonzero trends is much lower. For the full observational dataset (1981 to 2013), the median uncertainty in dDIC anth /dt across all σ 0 is 0.33 μmol kg -1 yr -1 , while the equivalent figures for 1991 to 1998 and 1997 to 2006 are 1.94 and 1.99 μmol kg -1 yr -1 respectively.
The atmospheric forcing used in the model simulations does not necessarily contain an NAOlike phenomenon, and even if there was one present it would not be expected to vary simultaneously with the real NAO. This is because the atmospheric data is entirely modelgenerated, rather than being from an atmospheric reanalysis. Consequently, an NAO effect cannot be directly observed in the model datasets. However, the model outputs can be used to indicate the unreliability of rates calculated using the EEL time-series data for these shorter time periods. The root-mean-square (RMS) difference between dDIC anth /dt in the SCM and TAA datasets is 0.19 μmol kg -1 yr -1 for 1981 to 2013, but it increases 0.93 μmol kg -1 yr -1 for 1997 to 2006, and 2.39 μmol kg -1 yr -1 from 1991 to 1998. The spatiotemporal heterogeneity of the observations therefore does adversely affect the calculated rates of change on these shorter timescales.
©2016 American Geophysical Union. All rights reserved.
That is not to say that the NAO does not influence dDIC/dt and its components. Indeed, we attribute the positive dDIC soft /dt to contraction of the subpolar gyre. This was probably itself driven by the NAO, as the atmospheric weather regimes associated with a positive NAOI phase tend to cause northward extension of the North Atlantic subtropical gyre [Gruber, 2009;Barrier et al., 2014]. Unfortunately, the significant increase in uncertainty that we find in the rates of change calculated over shorter timescales prevents these data from being used to support a direct link between the NAOI and the water column DIC at these shorter timescales.
In future studies, it may be instructive to investigate relationships between changes in water column DIC beneath the surface ocean mixed layer relative to the subpolar gyre index (SPGI), as well as the NAOI. At the sea surface, direct relationships between the NAOI and hydrographic properties might be expected [Thomas et al., 2008;Reverdin, 2010], as the NAOI is defined in terms of atmospheric conditions [Hurrell et al., 2003]. Despite the ability of the NAO to influence the subpolar gyre, the relationship between the NAOI and SPGI is non-linear and asymmetric (the response of the subpolar gyre to a negative NAOI phase is not simply the opposite of its response to a positive) [Lohmann et al., 2008]. As an oceanic property, the SPGI is perhaps more likely than the NAOI to directly correlate with changes in DIC, even if the ultimate driver of those changes is the NAO.

Applicability of EEL rates of change to wider area
The TAA dataset can be compared with FAA, in order to evaluate how changes observed at the EEL represent changes in the wider surrounding regions, in the model domain. This comparison suggests that changes observed in the EEL water column are representative of changes on a much larger spatial scale. However, the region that is most closely represented by the EEL varies with σ 0 . For illustrative purposes, we take the mean value of ΔdDIC/dt ©2016 American Geophysical Union. All rights reserved.
(and its standard deviation) and its components ((1) through (4)) at each σ 0 level across all model stations in the FAA dataset within the latitude range from 25 to 40°N and longitude range from 070 to 030°W to be representative of the North Atlantic, and from 66 to 72°N and 012°W to 001°E equivalently for the Nordic Seas (Figure 7). Therefore if the mean ΔdX/dt for any variable X in either region is close to 0, it means that the region's dX/dt is similar to that observed at the EEL, which is then considered to represent that region well. Positive ΔdX/dt indicates a faster increase (or slower decrease) in X at the station than at the EEL, and the opposite applies for negative values.
For most of the water column (σ 0 from 27.0 to 27.8 kg m -3 ), dDIC/dt is between 0.5 and 1.5 µmol kg -1 yr -1 higher in the Nordic Seas than in the North Atlantic, but the position of the EEL on this gradient shifts with depth (Figure 7(a)). In the upper part of the water column (σ 0 < 27.35 kg m -3 , e.g. Figure 8(a)), ΔdDIC/dt is close to 0 for the Nordic Seas, but at higher σ 0 (e.g. Figure 8(b)) the EEL rate more closely resembles the North Atlantic. The divide between upper and lower water column in this context, at σ 0 between 27.3 and 27.4 kg m -3 , corresponds to a depth of roughly 300 to 500 m at the EEL, and the bottom of the thermocline. Both the DIC soft (Figure 7(b)) and DIC carb (Figure 7(c)) components exhibit similar patterns as the total DIC changethat is, the EEL is changing more like the Nordic Seas in the upper water column, and North Atlantic lower down. However, because of how these components are combined to calculate the anthropogenic contribution to DIC change (1), the pattern is reversed for DIC anth . Its rate of accumulation in the upper water column is similar to the North Atlantic, while it matches the Nordic Seas at greater depths.
We can draw several conclusions from this part of the analysis. For a significant section of the water column, several variables are changing at the same rate at the EEL as they are throughout the wider surrounding regions in the model domain. However, different variables (and components of variables) at any given σ 0 level may not reflect changes in the same ©2016 American Geophysical Union. All rights reserved. adjacent region as each other. It is reasonable to expect that the base of the thermocline might be the σ 0 range where the EEL switches from representing one adjacent region to the other, as that is where the main currents change between travelling to the north and to the south at the EEL [Hansen and Østerhus, 2000]. Our analysis indicates that changes in DIC and its components are sufficiently spatially coherent that measurements of their changes along the EEL are representative of similar basin-wide changes.

Conclusions
Sufficient measurements have now been made along the Extended Ellett Line (EEL) transect to establish a time-series of data from which increases in DIC can be identified throughout the water column when combined with historical datasets. Most of the increase in DIC occurs in and above the thermocline. Anthropogenic CO 2 accumulation has driven 31 ± 6 % of the column inventory change in DIC, while the rest is due to increases in organic matter remineralization. The latter is likely driven by a net contraction of the North Atlantic subpolar gyre, and therefore an increasing influence of southern-sourced waters, during the time period considered in this study. The isotopic data provide independent supporting evidence for the attribution of DIC changes to these different driving processes. We found that values of -0.0166 ± 0.0003 ‰ (µmol kg -1 ) -1 for ΔRC and -27.0 ± 0.5 ‰ for δ 13 C POC were best able to explain the observed dδ 13 C DIC /dt profile.
Future EEL occupations will provide additional data to extend this time-series analysis. Once more consecutive years of data become available, it may also become possible to better assess the influence of processes operating on sub-decadal timescales, for example the North Atlantic Oscillation. A more robust quantification of interannual variability at the EEL will ©2016 American Geophysical Union. All rights reserved. be useful to better evaluate the confidence bounds on rate-of-change calculations, and timesof-emergence for long-term trends.
Combination of the observational data with output from model simulations demonstrated that the spatiotemporal heterogeneity in the distribution of the observations does not adversely affect the calculated rates of change. The model data have also provided insight into the relevance of the EEL in a larger regional setting, suggesting that changes observed locally may reflect much wider scale changes occurring in the North Atlantic Ocean and Nordic Seas.