Impact of physical and biological processes on temporal variations of the ocean carbon sink in the mid-latitude North Atlantic (2002–2016)

Abstract The ocean is currently a significant net sink for anthropogenically remobilised CO2, taking up around 24% of global emissions. Numerical models predict a diversity of responses of the ocean carbon sink to increased atmospheric concentrations in a warmer world. Here, we tested the hypothesis that increased atmospheric forcing is causing a change in the ocean carbon sink using a high frequency observational dataset derived from underway pCO2 (carbon dioxide partial pressure) instruments on ships of opportunity (SOO) and a fixed-point mooring between 2002 and 2016. We calculated an average carbon flux of 0.013 Pg yr−1 into the ocean at the Porcupine Abyssal Plain (PAP) site, consistent with past estimates. In spite of the increase in atmospheric pCO2, monthly average seawater pCO2 did not show a statistically significant increasing trend, but a higher annual variability, likely due to the decreasing buffer capacity of the system. The increasing Δ pCO 2 led to an increasing trend in the estimated CO2 flux into the ocean of 0.19 ± 0.03 mmol m−2 day−1 per year across the entire 15 year time series, making the study area a stronger carbon sink. Seawater pCO2 variability is mostly influenced by temperature, alkalinity and dissolved inorganic carbon (DIC) changes, with 77% of the annual seawater pCO2 changes explained by these terms. DIC is in turn influenced by gas exchange and biological production. In an average year, the DIC drawdown by biological production, as determined from nitrate uptake, was higher than the DIC increase due to atmospheric CO2 dissolution into the surface ocean. This effect was enhanced in years with high nutrient input or shallow mixed layers. Using the rate of change of DIC and nitrate, we observed Redfieldian carbon consumption during the spring bloom at a C:N ratio of 6.2 ± 1.6. A comparison between SOO and PAP sustained observatory data revealed a strong agreement for pCO2 and DIC. This work demonstrates that the study area has continued to absorb atmospheric CO2 in recent years with this sink enhancing over time. Furthermore, the change in pCO2 per unit nitrate became larger as surface buffer capacity changed.


Introduction
The cumulate emissions of carbon by human activities since the industrial revolution have caused atmospheric CO 2 concentrations to approach 410 ppm in 2018 (Le Quéré et al., 2018a), a number not seen on Earth for the past 3 million years (Seki et al., 2010). The atmospheric concentration would be even higher and climate change would proceed much faster had the ocean not taken up approximately a quarter of the cumulative anthropogenic carbon emitted between 1870 and 2016 (Le Quéré et al., 2018b). Although the global ocean is an overall carbon sink, the distribution of the sources/sinks varies across its surface (Takahashi et al., 2009). Of the sink regions, the North Atlantic is the largest per surface area, taking up a total of 0.7 ± 0.1 Pg C annually (Gruber et al., 2002). The size and evolution of the sink is a matter of considerable societal interest.
The ocean's ability to take up carbon from the atmosphere is determined by the Revelle factor, the ratio of change in carbon dioxide to the change in total DIC (Broecker et al., 1979;Egleston et al., 2010). Low Revelle factors correspond to high buffer capacity in regions such as the North Atlantic, and reflect an enhanced ability to absorb https://doi.org/10.1016/j.pocean.2019.102223 Received 17 April 2019; Received in revised form 12 September 2019; Accepted 13 November 2019 atmospheric CO 2 (Bates et al., 2014). The anthropogenic increase in atmospheric CO 2 and corresponding ocean uptake will drive an increase in oceanic pCO 2 and a consequential decrease in the oceanic buffer capacity Sabine et al., 2004). In the modern ocean, the Revelle factor is increasing in response to carbon accumulation (Fassbender et al., 2017). In addition, absorption of CO 2 leads to a reduction in pH and lowering of calcium carbonate saturation state in the surface ocean. This process, known as ocean acidification, potentially leads to decreased calcification and growth rates in some calcifying species (Doney et al., 2009). At the moment, natural variability is masking some of the anthropogenic effect on ocean sinks (McKinley et al., 2016), but studies from Le Quéré et al. (2010) -for the 1981 to 2007 period or Schuster and Watson (2007) -for 1994 showed that regionally, current sinks can become weaker. Long term studies found that between 1990 and 2006 ) and between 1993(McKinley et al., 2011, the mid and high-latitude North Atlantic seawater pCO 2 increased at faster rates than the atmospheric rate. In fact, due to the decrease in buffer capacity of the northern North Sea, the North Sea basin shifted from being an overall sink to an overall source of atmospheric CO 2 in 2011 (Clargo et al., 2015). Changes such as this one emphasise the importance of long-term observations of the surface ocean. The existing global network of fixed point ocean observatories has produced important time series of biogeochemical variables, but the effective footprints of these observatories are estimated to only reflect around 9-15% of the ocean surface (Henson et al., 2016). Ships of Opportunity (SOO) equipped with CO 2 measuring facilities are an additional way of monitoring by sampling different regions and having a large observational capacity. This is particularly important since better understanding of the role of oceans in the global carbon budget has been identified as a major challenge (Steinberg et al., 2001).
In this study we make use of observations of carbon and nutrients collected on a SOO route between 2002 and 2016 combined with ancillary modelled and observational oceanographic and atmospheric data to investigate (1) variability, (2) trends, and (3) physical and biological drivers of the mid-latitude North Atlantic Ocean carbonate system. The SOO route passes in the vicinity of the Porcupine Abyssal Plain sustained observatory (PAP-SO; 49°N, 16.5°W), an area where the variability of the carbon parameters has previously been studied . This allows an inter-comparison to be performed. The PAP-SO is found within the NA-STSS biome (Fay and McKinley, 2014), where data-based methods suggested a relatively stable air-sea flux between 1990 and 2010 (Rödenbeck et al., 2015), but neural network methods found a decrease in the carbon sink between 1998 and 2007 (Landschützer et al., 2013). There are still significant uncertainties regarding the interannual to decadal variability of the ocean carbon sink in this region.

Data acquisition and processing
The main source of data for this study derives from the analysis of water samples collected aboard commercial ships chartered by Geest Line Ltd. (Fareham, UK): MV Santa Maria, MV Santa Lucia and MV Benguela Stream. These ships transited between the UK and the Caribbean nearly every month over the study period ). An automated underway CO 2 system provides continuous measurements of surface seawater temperature (SST), salinity (SSS) and fCO 2 (carbon dioxide fugacity) -converted to pCO 2 (carbon dioxide partial pressure) . These data are available from the Surface Ocean CO 2 Atlas (SOCAT version 6, SOCAT flags A-E, WOCE flag 2; www.socat.info; Bakker et al., 2018). While SSS is not quality controlled in SOCAT, we are confident with the data. On the UK-Caribbean line, salinity bottle samples are collected regularly, analysed at the National Oceanography Centre, UK, and used to calibrate the sensor. In addition, water samples were collected every 4 h in 15 mL plastic tubes and frozen immediately for later analysis of inorganic nutrients (nitrate + nitrite -hereby nitrate, silicate, phosphate) in the laboratory on shore. Nutrient data are available from the British Oceanographic Data Centre (doi: 10.5285/7a7b497a-c47b-0a69-e053-6c86abc07ef0). We selected data within the 'footprint' of the PAP-SO. This footprint was defined by Henson et al. (2016) as the spatial area where the outputs from Earth System Models have similar mean and variability to the station time series. We used the footprint that they calculated with pH as a parameter since this present study deals with the carbonate system and pH is the only carbonate system parameter considered in their study. Thus, we selected data within the boundaries of this area. In addition, we excluded from the analysis any measurements taken in the shelf areas to avoid the influence of coastal processes ( Fig. 1). Our study area also lies within the North Atlantic Drift (NADR) biogeographic province as defined by Longhurst (2006). The SOO used in the present study have occupied this route since 2002, but other cruises were conducted in this area before, making it one of the most regularly observed ocean regions and allowing a long term assessment of change in the surface mid-latitude North Atlantic Ocean. Monthly averages were first prepared for all parameters and used in subsequent calculations. For months when surface pCO 2 data were not available on the usual route, the gaps were filled with values measured by other ships in the same area, also available from SOCAT. When nitrate measurements were not available for a particular month, the long-term climatological value (obtained from the results of this study) was used as a replacement. The same climatological gap-filling method was not used for pCO 2 since one of the objectives of this study was precisely to investigate interannual variability and differences from mean annual pCO 2 .
Available SST and SSS data were used to estimate total alkalinity using the empirical equation derived by Lee et al. (2006) for the 30°N -80°N North Atlantic region. We assume the empirical relationship between SSS and alkalinity is valid for the entire time series. With two of the carbonate system parameters now determined (pCO 2 and total showing all pCO 2 measurement locations in red symbols, the PAP footprint as the black contour, the excluded area as the grey contour, the location of the mooring at the PAP-SO site as the blue diamond and Mace Head atmospheric observatory as the green triangle. The orange box shows the extent of the NPP data selected. Produced with Ocean Data View (Schlitzer, 2016). alkalinity), all others can be mathematically derived. We used the Matlab CO2SYS toolbox (Heuven et al., 2011) with the Lueker et al. (2000) K1 and K2 dissociation constants and the Dickson et al. (1990) K SO 4 constant to obtain dissolved inorganic carbon (DIC) concentrations. In order to eliminate potential freshwater influences, the DIC was normalised (nDIC) to the area average salinity (35.6).
Atmospheric xCO 2 (dry mole fraction) data were downloaded from the World Data Centre for Greenhouse Gases from the Mace Head station (WDCGG, 2018). This was chosen due to the proximity to our study area.
Monthly mean atmospheric pressure, dew point temperature and wind speed were obtained from the ERA-Interim (Dee et al., 2011), a reanalysis that combines model derived data and observations to produce a best estimate of global atmospheric and oceanographic parameters.
We used mixed layer depth (MLD) from Argo floats with data hosted by the Scripps Institution of Oceanography, UC San Diego and available every month of our time series (Holte et al., 2017), defined as the depth at which a density difference of 0.03 kg m −3 relative to the surface was found.
Since no measurements on the UK-Caribbean route were taken north of 50°N, we further restricted the area part of the PAP-SO footprint that was investigated in this study (Fig. 1). Extra SOCAT data and ancillary variables were only selected from this defined area.
Satellite-derived net primary productivity was downloaded from the Ocean Productivity website (Ocean Productivity, 2016) calculated according to Behrenfeld and Falkowski (1997). These data were selected from within the area shown in Fig. 1.

Calculations
The dry mole fraction of atmospheric CO 2 is converted to the partial pressure of atmospheric CO 2 by (1)  )] 10 v dp dp 5 (2) where T dp is the dew point temperature [°C] (Alduchov and Eskridge, 1996). With both atmospheric and seawater pCO 2 known, sea-air CO 2 flux [mmol m −2 day −1 ] was calculated following Wanninkhof (2014): where k is the gas transfer velocity [m day −1 ] and is the solubility of CO 2 [mmol m −3 atm −1 ] (Weiss, 1974). The gas transfer velocity is calculated using the Schmidt number and the wind speed, while is a function of SST and SSS. Negative values indicate flux from the atmosphere into the ocean (Wanninkhof et al., 2013). We calculated air-sea CO 2 fluxes using code freely available online (https://github.com/ mvdh7/co2flux) -see Humphreys et al. (2018). Water density used for unit conversions was calculated using the Seawater Oceanographic Toolbox (McDougall and Barker, 2011). With the information generated as above, the contributions of physical and biological processes to monthly changes in nDIC can be calculated. We used the same 'model' employed by Jiang et al. (2013). Firstly, we determined the gas exchange influence using where DIC Gas is the gas exchange term, F is the sea-air CO 2 flux converted to μmol m −2 month −1 and MLD i ( ) is the average mixed layer depth for that particular month. Then, we determined the biological production influence using where DIC BP is the biological production term, NO 3 is the average nitrate + nitrite concentration, 'i' is the month, is the average mixed layer depth of two consecutive months and the C:N ratio is the carbon to nitrogen ratio of biological production/remineralization: 6.6 for this study Henson et al., 2003). When surface nitrate is replenished in winter, DIC is also brought up, but this mixing term is not given an individual decomposition term due to the high uncertainties.
The difference between the observed changes and the sum of the individual contributions gives a residual term, which includes vertical mixing.
where Temp i ( ) and Temp i ( 1) are the mean SST of the current and previous months respectively.
As with DIC, the difference between the observed changes and the sum of the calculated components was assigned to a pCO Residual ( ) 2 term. Finally, the individual monthly contributions to the carbonate system changes were combined to produce annual cycles. These show the relative importance of the different processes and also help assess the temporal variability within the PAP footprint when comparing to individual years.

Uncertainties
Given the multiple data sets used, with respective associated uncertainties and given the complex nature of the CO 2 toolboxes employed, we assess uncertainty of our calculations via a Monte Carlo simulation, similar to past studies (Fröb et al., 2018;Palevsky and Quay, 2017). We used the 6.4 μmol kg −1 uncertainty defined by Lee et al. (2006) for their North Atlantic alkalinity calculation and the 10 μatm error range for flags A-E in the SOCAT database . This is the largest uncertainty allowed in the database but a great majority of the measurements used were more accurate. The results of this simulation are therefore an upper-bound error estimation. We used a 0.1 μatm uncertainty for atmospheric pCO 2 measured at Mace Head (Derwent et al., 2002). Since most calculations in this study are performed using monthly averages, the SSS, SST and wind speed standard deviation of the monthly averages were used as their respective uncertainties. These uncertainty intervals are shown in Table 1.
A total of 1000 simulations were run with the input terms at each of 10 random locations using the CO2SYS toolbox to determine DIC and the CO2flux toolbox to determine CO 2 flux. The input parameters (Table 1) were allowed to vary randomly within the ± uncertainty interval. The DIC standard deviation as a percentage of the mean value of the simulations was found to vary between 0.27 and 0.30%. The CO 2 flux standard deviation as a percentage of the mean value of the simulations varied between 16 and 34%, with higher uncertainties for fluxes closer to 0. These large values are mainly driven by the large range in the wind speed estimations inside the PAP footprint. When all other input terms were kept fixed and only wind speed allowed to vary, the resulting error was of a similar range. The significant effect of the choice of wind product and the use of monthly mean wind speed on interannual variability of sea-air CO 2 exchange has been studied before (Olsen et al., 2005). In order to improve the precision of the estimates, wind speed measured on the sampling ship would ideally be used. Since such data were unavailable, reanalysis data, with their associated uncertainty, have been used in this study.

The average annual cycle
The annual cycles features ( Fig. 2) are typical of a mid-latitude northern-hemisphere oceanic area. These were obtained by averaging available data over the restricted PAP-SO footprint. Comparing with the annual cycles derived by Frigstad et al. (2015) from the PAP surface mooring, the SST, MLD and nitrate are very similar both in terms of the range and the annual evolution.
The double peak in seawater pCO 2 and CO 2 flux is not however observed in the PAP-derived annual cycles. It is also not observed in the North Atlantic climatology derived by Takahashi et al. (2009). However, this feature was observed in a study in the Bay of Biscay (Jiang et al., 2013). Further decomposition by latitude reveals that the summer peak is mostly driven by measurements in the southern half of the study area, which give an annual cycle closer to a typical subtropical one, while the measurements in the northern half give an annual cycle closer to a typical subpolar record, resembling a sinusoid. The PAP-SO site, where the double peak is not observed, is situated in the northern half.
Nitrate concentrations are highest in the months of January, February, March and April. These winter values are higher than 4 μmol L −1 , while the summer values are close to the detection limit, suggesting that nitrate is a limiting nutrient for phytoplankton growth following the spring bloom in this region.
The mixed layer deepens in the winter up to a maximum of over 200 m, while in the summer, it is shallower than 20 m in July and August. The strongest average winds (over 10 m s −1 ) are recorded in the months of December, January and February and summer values drop to below 7 m s −1 .
Features such as a double peak in seawater pCO 2 and CO 2 flux are observed. Carbon uptake is weak at the end of the winter season (February-March) and end of the summer season (August-September). Seawater pCO 2 has the highest average values in the months of February and March at over 370 μatm with a secondary peak in August at just over 360 μatm, while CO 2 flux into the ocean has two peaks in November and May. On average, over 5 mmol m −2 day −1 of CO 2 enter the surface ocean in the months of November and December.
In the following section, in order to avoid issues related to autocorrelation of the variables, monthly anomalies (seasonally detrended monthly averages) were calculated and used to determine long-term trends.

Long term variability
Atmospheric pCO 2 calculated from the measured xCO 2 at Mace Head Observatory in Ireland had a mean increase of 2.09 ± 0.03 μatm yr −1 (p < 0.01) from 2002 to 2016 (Fig. 3a), similar to the rate observed during the last decade at the Mauna Loa Observatory (Earth System Research Laboratory, 2018).
Seawater pCO 2 in our study area however does not follow the same trend (Fig. 3b). There is a small, not statistically significant (p = 0.09) rise of 0.37 ± 0.22 μatm yr −1 but the amplitude of the seasonal cycles becomes larger in the second half of the time series, with monthlyaveraged values in 2003 ranging from 343 to 368 μatm, while the range in 2016 is 336 to 392 μatm. This observation differs from those of Bates et al. (2012) and Takahashi et al. (2009), who found an increasing trend of 1.8 μatm yr −1 near Bermuda and of 1.66 to 2.13 μatm yr −1 in the NADR, respectively. In a similar study area, McKinley et al. (2011) and Schuster et al. (2009) found seawater pCO 2 trends exceeding the atmospheric one between 1993 and 2005 and 1990 and 2006 respectively. However, it is consistent with other studies where no trend in seawater pCO 2 was observed, such as in the Eastern Tropical Atlantic (Lefèvre et al., 2016). Although in a different study area, where different mechanisms are at play, this study shows that regionally, seawater pCO 2 may not reflect the rise in atmospheric levels. In fact, a recent neural network mapping-based study found that in our study area, observations indicate a lower seawater pCO 2 than predicted by the climatology (Fig. 6 in Landschützer et al., 2014).
There is no statistically significant trend in calculated DIC either (Fig. 3h). Monthly-averaged values range between around 2040 and 2140 μmol kg −1 , with a slight increase in the variability in the second half of the time series. Since there is a statistically significant decrease in salinity anomalies with time (−0.018 ± 0.005 yr −1 , p < 0.01, Fig. 3d), salinity-normalized DIC anomalies (nDIC) show a positive trend (0.68 ± 0.18 μmol kg −1 yr −1 , p < 0.01). This trend is lower than all seven of the carbon time series evaluated by Bates et al. (2014), spanning over various ocean provinces, which show positive trends in nDIC, but with faster increases with time. Their analysis period is also much longer than the present study. The recent freshening of the North Atlantic has been attributed to the interaction of increased freshwater fluxes and subpolar gyre circulation (Tesdal et al., 2018) and can also be a result of the global water cycle intensification (Durack et al., 2012).
Since atmospheric pCO 2 has been increasing and seawater pCO 2 has not, the flux of CO 2 into the ocean has increased by 1.19 ± 0.03 mmol m −2 day −1 per year (p < 0.01, Fig. 3k). Data show the study area acted as a CO 2 source only on a few occasions during the 15 year-long time series (September 2005, September 2006, August 2007, February 2010and August 2010. Apart from one, all of these instances occurred in late summer, early autumn, at the end of the productive season when atmospheric pCO 2 is at the annual minimum. Data from the nearby PAP observatory show this area is an overall carbon sink (Körtzinger et al., 2008). Consistently, we found an average ± standard deviation CO 2 flux of 1.33 ± 0.71 mol m −2 year −1 into the ocean, weaker than the flux calculated by Körtzinger et al. (2008) at the PAP-SO site. The PAP-SO sensor had however recorded much lower values of seawater pCO 2 than the SOO. The resulting disequilibrium between the ocean and the atmosphere pCO 2 leads to a strong flux into the ocean. This mismatch between the datasets in 2004 is discussed further in the following section.
Integrated over the 7.9 × 10 5 km 2 area of the restricted PAP footprint, the 2002-2016 SOO-derived flux yields an uptake of 0.013 Pg C yr −1 . This uptake is stronger than the area-weighted estimation of Table 1 Range of standard deviations as percentage of mean values for the output of the CO2SYS and CO2flux toolboxes obtained by running 1000 Monte Carlo simulations for 10 random locations. V.A. Macovei, et al. Progress in Oceanography 180 (2020) 102223 Takahashi et al. (2009) for the mid-latitude North Atlantic in the year 2000, but weaker than than the 1985-2010 area-weighted uptake estimated for the 'North Atlantic Subtropical Seasonally Stratified' biome by Rödenbeck et al. (2015).
Since the salinity time series shows no values below 35.4 measured on the UK-Caribbean route before 2011, but multiple in the second part of the time series, including 11 consecutive months below this value, trends were calculated for the 2002-2010 and 2011-2016 periods and synthesised in Table 2. The rise in atmospheric pCO 2 has accelerated with time. On the other hand, the seawater pCO 2 trend was only statistically significant before 2010. As a consequence, the CO 2 flux into the ocean trend only became significant after 2011.
While there is no statistically significant trend in seawater pCO 2 when all the monthly averages are analysed, only using the annual maximum value results in a positive trend of 1.63 ± 0.44 μatm yr −1 . This is closer to the atmospheric pCO 2 trend. Several studies use winter observations to detect pCO 2 trends in order to avoid summertime biological activity biases (Fröb et al., 2019;Metzl et al., 2010). At the same time, the average Revelle factor is increasing and thus the buffer factor of the system is decreasing. The long term increase in the Revelle factor is 0.019 ± 0.011, consistent with the time series investigated by Bates et al. (2014), albeit not statistically significant (p = 0.08). This means that the amount of pCO 2 change expected to occur for a given change in DIC will be larger (Egleston et al., 2010). Towards the end of the time series, there is a larger difference between the annual seawater pCO 2 maximum (at the end of the winter) and minimum (after the spring bloom). Nitrate concentrations are not higher at this time, but the spring bloom, which to a first order would be expected to draw down DIC in proportion to the amount of nitrate removed, is having an increased influence on the resulting summer pCO 2 decrease. This is consistent with the increased seasonality observed in the pCO 2 time series and visually represented in Fig. 4. The ratio of pCO 2 (μatm) -here used to signify the seasonal winter to summer difference -to Nitrate (μmol kg −1 ) increased with time by 0.31 per year (p < 0.05, R 2 =0.35) when one outlier was removed. Recently, Landschützer et al. (2018) found observations that confirmed the predictions of increased atmospheric CO 2 impacting the seasonal cycle of ocean carbon. Consistent with this study, they found an increasing difference between the winter and summer seawater pCO 2 .

Comparison with PAP-SO
A comparison of the footprint-wide data with the ones from the PAP-SO buoy sensors at 1 and 30 m is shown in Fig. 5. Although one of the sensors is located at 30 m, this depth is within the mixed layer for approximately 8 months on any given year. The MLD is shallower than 30 m only between June and September on average. These months are indicated on the comparison plots.
There   , wind speed (f), surface nitrate (g), DIC (h), nDIC (i), alkalinity (j) and CO 2 flux (k). The seawater pCO 2 , SST, SSS, DIC and alkalinity blue symbols represent data from the SOO usual route and the red symbols show months when data was filled in from other ships transiting our study area. Statistically significant trends in the monthly anomalies are shown in magenta lines. 40 μatm, but many of the measurements were taken when the sensor was outside the shallow summer mixed layer. The PAP sensor data for the 2003-2004 interval is deemed correct for the PAP location as it has been calibrated against pCO 2 derived from direct DIC and total alkalinity measurements (Körtzinger et al., 2008). Perhaps the seasonal variability has increased in the second part of the time series because the SOO sensor was not capturing the same signal as the 30 m PAP sensor in the early years. The lowest number of SOCAT observations in our study area are recorded in this 2003-2004 interval, which might be an additional reason for the mismatch. Even with a mismatch in the pCO 2 comparison, the nDIC from the ship and the PAP site show a strong correspondence. DIC is calculated from pCO 2 and alkalinity. Similar results with different pCO 2 inputs highlight the greater weighting factor of alkalinity compared to pCO 2 in this calculation. Although being derived through different methods, the SOO and PAP site nDIC are very similar. We investigated the influence of the differences between the two data sets on our CO 2 flux calculations. Monthly averages were computed for the PAP 1 m sensor data in order to directly compare and compute fluxes (as described earlier). For the 18 months when data was available from both sources (excluding the 2014 period affected by the storm), the average CO 2 flux into the ocean was 4.54 and 5.90 mmol m −2 day −1 when calculated with the SOO and the PAP data respectively. The relevance of this difference depends on the questions being addressed and the required accuracy. For our study, we find the agreement satisfactory in order to strengthen our conclusion that the area around the PAP site has remained a strong carbon sink region.

Discussion
Despite displaying a larger seasonal variability in the second part of the time series, seawater pCO 2 does not follow the rising trend of atmospheric pCO 2 and instead does not show a statistically significant trend. This finding has implications for the ability of this area of the ocean to remain a carbon sink, as it has been previously identified, in the context of increasing atmospheric forcing. In order to determine the reason behind this behaviour, we looked at the individual factors exerting control on seawater pCO 2 . Seawater pCO 2 is mainly a function of temperature, salinity, DIC concentration and alkalinity . Tjiputra et al. (2014) V.A. Macovei, et al. Progress in Oceanography 180 (2020) 102223 driver of pCO 2 changes. We therefore first investigate what causes changes in surface nDIC.

Processes affecting nDIC
DIC concentration is determined by physical and biological processes. Gas exchange affects it through the dissolution of carbon dioxide in seawater (Park, 1969). Influx of CO 2 into the sea increases the DIC concentration. DIC is consumed by phytoplankton during the productive season and replenished during winter through mixing/convection (Shadwick et al., 2011). Following calculations explained in Section 2.2, in Fig. 6a and b, we present monthly contributions of variables to nDIC.
Our assessment reveals DIC changes due to gas exchange increased with time (Fig. 6a), corresponding to the increasing pCO 2 . There is a statistically significant (p < 0.01) trend of 1.8 ± 0.3 μmol L −1 per year in the yearly contribution between 2003 and 2016, the complete data years. This is more than 7% of the average yearly sum of 25.2 μmol L −1 -the contribution of gas exchange to the nDIC concentrations. DIC Gas exchange is largest in May, June and July (Fig. 6d). These are the months that on average have the lowest seawater pCO 2 values (Fig. 2a) resulting in large pCO 2 . In addition, the same months also have shallow MLD (Fig. 2g) so, according to Eq. (4), the gas exchange contribution will be greater.
The DIC BP terms show a trend of 2.4 ± 1.5 μmol L −1 per year, albeit not statistically significant -more than 3% of the − 79.6 μmol L −1 average yearly sum of the DIC BP contributions. The monthly nitrate concentrations and average MLD do not show a statistically significant trend during our study period, so the DIC BP changes are likely caused by a change in the timing of the spring bloom. Palevsky and Quay (2017) found that the seasonal timing of net community production strongly influences the effect on air-sea CO 2 flux. Earlier onsets of the decrease in nitrate, during the months when MLD is deeper, lead to greater integrated DIC BP values according to Eq. (5). The contribution of the DIC BP is largest in 2008, 2011 and 2015. While nitrate concentrations were high in the winter season before the productive season in 2011 and 2014, the large value in 2008 is driven by the relatively shallow MLD of that year. For 2011 and 2015 this is coherent with independent net primary productivity (NPP) derived from satellite. The strongest uptake of DIC is observed during the spring bloom in the months of April, May and June with average DIC changes of around 20 μmol L −1 per month. It is not uncommon for biological production to have an influence on the carbonate system. In the Eastern North Pacific, DIC removal by net community production more than fully offsets the DIC increase due to air-sea gas exchange (Palevsky and Quay, 2017). The increase in the importance of the biological term is consistent with the increase in NPP observed in this region (black symbols in Fig. 6b -note the reversed axis).
The DIC Residual term annual contribution (Fig. 6c) is positive every Fig. 6. Annual sum of the monthly contributions of gas exchange (a), biological production (b) and the residuals (c) to the nDIC changes (all bars in μmol L −1 ); the long-term average nDIC monthly contributions with the average monthly nDIC change plotted in black symbols on top (d); the nDIC monthly contributions predicted from the typical annual cycle in Fig. 2 (e). Subfigure (b) also shows the annual average NPP (mg m −2 day −1 ) on a reversed y-axis.
year since it includes the influence of wintertime vertical mixing, which brings up DIC from the deep waters. The reventilated CO 2 produced by subsurface respiration was found to negate more than one third of the carbon drawdown by net community production at the PAP site (Körtzinger et al., 2008). Other contributions to the DIC Residual are likely due to horizontal advection as well as the propagation of the combined uncertainty of the terms used in the calculations. Fig. 6d and e show the average annual cycle of the contributions from the individual monthly contributions and from the mean annual cycle presented in Fig. 2. While mostly similar, some differences are seen in the magnitude of the wintertime DIC Residual terms and the absence of wintertime DIC BP term when the mean annual cycle is used. The rare occurrences of nitrate decreasing between two consecutive autumn/winter months and causing a DIC BP contribution are eliminated when using the typical annual cycle.
In our study, during an average year, the DIC Residual term accounts for 68.9% of the total changes in the nDIC concentration. This is a large value, but much of the contribution is given by winter mixing of DICrich water.

Processes affecting pCO 2
As discussed above, the calculated changes in nDIC were used to determine their respective influence on pCO 2 . The terms that cause an increase (a decrease) in surface nDIC concentration also cause a proportional increase (decrease) in pCO 2 . Decomposing the nDIC contributions in the previous section allows us to decompose the pCO 2 contributions in more detail than just using temperature, alkalinity and DIC changes.
In addition to the terms discussed above, the temperature influence is introduced as a new term ( pCO Temp ( ) 2 ). The individual monthly contributions as well as the typical annual cycle are shown in Fig. 7. The warm summer temperatures cause a rise in pCO 2 between the months of April and August and the reverse is true for the colder temperatures between September and March. The sum of the average temperature contribution during the April-August period would cause an increase of surface pCO 2 of 112 μatm. This is however balanced by the 106 μatm opposing change forced by biological production during the same period. Fay and McKinley (2017) found the pCO Temp ( ) 2 term is more important for the pCO 2 seasonality in the subtropics, while the pCO DIC ( ) 2 BP becomes more important in the subpolar areas. The PAP footprint lies in the mid-latitude ocean, between the subtropical and subpolar areas and it appears that these two terms are of similar magnitude, which is consistent with the conclusions of Takahashi et al. (2002). For our entire study period and study area there is no significant change in annual average SST. Between 2005 and 2014 however, there is a small decrease of − 0.07 ± 0.03°C yr −1 , consistent with the findings of Robson et al. (2016), who identified a cooling and freshening of the upper ocean over the 2005-2014 period in the northeast Atlantic, where our study area is situated. It has been shown however that the effect of large temperature differences between winter and summer is more important than the maximum temperature reached during the summer. Strong warming in the transition from the winter to the summer seasons can indeed outweigh the effect of intense phytoplankton blooms (Dumousseaud et al., 2010).
The influence of changes in alkalinity only produces a strong signal in years when there are rapid changes in alkalinity between two consecutive months. Because our derived alkalinity is largely related to SSS, these coincide with changes in salinity.
There is a statistically significant increasing trend in the annual sum of the pCO DIC ( ) Gas Exchange 2 term of 3.1 ± 0.5 μatm yr −1 . This is due to the increasing difference between the atmospheric and seawater pCO 2 and not due to increasing wind speed, which shows no significant trend in our analysis. Although future projections estimate an increase in global wind speed, only the magnitude, not the direction of the flux will be affected (Wanninkhof and Triñanes, 2017). Both sink and source areas are estimated to become stronger.
The difference between the observed changes in pCO 2 and the sum of the contributions of temperature, alkalinity and various DIC changes has been assigned to the pCO Residual ( ) 2 term. In years such as 2010, when some data are missing from the main UK-Caribbean route, this term is larger since the sum of the calculated terms does not match some of the abrupt changes observed between consecutive months in the SOCAT database. Although for an average year, this residual term only accounts for 23.4% of the total contribution, it is negative in all but one of the years of the time series (i.e. 2009), which implies there exists a sum of unknown processes lowering pCO 2 . It appears that part of the reason the seawater pCO 2 does not show an increment during our study period is due to this additional term.

Carbon to nitrogen stoichiometry and nitrate-pCO 2 relationship
To further investigate the effect of new production on the carbonate system variables, we focus on the spring bloom months. The PAP footprint area exhibits annual cycling typical of a mid-latitude Northern Hemisphere oceanic region (Wong et al., 2002). We work under the assumptions that the April, May and June months represent the bloom period since nitrate concentrations decrease abruptly. The nitrate concentrations are on average below 0.7 μmol L −1 during July, August, September and October. The average SST is above 16°C during these months, which were grouped together as the 'post-bloom' period. The remaining months were defined as the 'pre-bloom' period. We thus used the same notation as Jiang et al. (2013) in grouping the months around the timing of the bloom.
In this study we used a C:N ratio of 6.6 (i.e. Redfield ratio) to determine the influence of biological production on nDIC concentrations. There are studies however that suggest carbon overconsumption occurs in the North Atlantic (Koeve, 2004;Sambrotto et al., 1993;Toggweiler, 1993). The gas-corrected nDIC ( DIC Gas exchange excluded) and nitrate concentrations in our study are linearly correlated (R 2 =0.37, p < 0.01) during the 'bloom' months (blue symbols in Fig. 8a). The slope of this line, obtained by ignoring nitrate values < 0.5 μmol L −1 similarly to the method used by Jiang et al. (2013), is 6.2 ± 1.6, similar to the canonical Redfield ratio. While we found no evidence for carbon overconsumption when using the data within the restricted PAP footprint, performing the same calculations with data from the whole NADR (a rectangular box between 42 and 55°N and 10 and 42°W) revealed a C:N ratio of 9.6 ± 1.4, a value similar to that found in a study in the Iceland Sea . In both cases it seems likely remineralisation occurs within the mixed layer since the 'prebloom' (red symbols in Fig. 8a) have a similar trendline, indicating a comparable C:N ratio during the wintertime mixing.
In Fig. 2a we showed that the annual seawater pCO 2 displays a double peak distribution with high values in winter as well as late summer. This behaviour has also been reported by Jiang et al. (2013) in the Bay of Biscay. While subtropical regions display a distribution with a single summer maximum due to the temperature control (Bates et al., 1996), the DIC variations are the main driving factor for pCO 2 leading to a single winter peak in subpolar regions (Olsen et al., 2008). The double peak in our study area suggests that both controls are important, consistent with past studies (Takahashi et al., 2002). Winter mixing of DIC and nitrate-rich water coincides with the increase in pCO 2 , which is then removed by spring biological production. Warming in the summer season then increases pCO 2 without a corresponding increase in nitrate. The temperature-driven range of pCO 2 values during the 'post-bloom' summer period in Fig. 8b is similar but slightly smaller compared to the biological production and remineralization-driven range of the 'prebloom' and 'bloom' values. While the maximum pCO 2 value reached due to summer warming is smaller than 380 μatm, the winter upwelling of DIC-rich water makes the surface seawater pCO 2 reach almost 400 μatm. These results are consistent with a recent study that found the subpolar North Atlantic pCO 2 to be predominantly correlated with biological processes and the subtropical gyre pCO 2 with temperature on seasonal time scales (Henson et al., 2018).

The residual terms
Around 23% of the annual average seawater pCO 2 observed changes are not explained by the individual terms calculated in this study. Part of the difference could be due to using atmospheric reanalysis products and estimated (rather than measured) alkalinity. To allow inter-comparison, all the terms are converted into monthly averages which might introduce further error if the study area is not homogeneous over its surface. This study does not calculate the influence of advection, which is known to change carbon concentrations in the water column (Fröb et al., 2018). Lateral advection in the subpolar gyre can increase DIC concentrations on a similar magnitude to air-sea gas exchange (Zunino et al., 2015). However, advection can also increase carbon storage if the water that recently absorbed atmospheric CO 2 is moved to a mode water formation site (Andersson et al., 2013). The influence of calcification is also not included in this study. Calcification consumes DIC, but is a CO 2 source (Frankignoulle and Canon, 1994;Shutler et al., 2013). However, with no direct measurements, the influence of calcification would need to be determined through proxies such as using the particulate inorganic carbon to particulate organic carbon ratio (PIC:POC). This estimation has a high uncertainty and a small relative influence on the seawater pCO 2 (Palevsky and Quay, 2017). A recent study found that the presence in the surface water of biological surfactants can reduce gas exchange by up to 32% due to the influence on the gas transfer velocity (Pereira et al., 2018). In our calculations of CO 2 flux, the Fig. 7. The annual sum of the monthly contributions to seawater pCO 2 (a); the long-term average pCO 2 monthly contributions with the average monthly pCO 2 change plotted in black symbols on top (b); the pCO 2 monthly contributions predicted from the typical annual cycle in Fig. 2 (c). surfactant influence on k is not accounted for in Eq. (3). While a decrease of the gas flux into the ocean would move the system in the right direction in order to explain some of the pCO Residual ( ) 2 term, the pCO DIC ( ) Gas Exchange 2 contributes on average only 4.7% to the observed changes so a small change to it would not be enough. Whichever way we decompose the reasons for the variability, the key finding that our observations show us is the relatively stable annual mean seawater pCO 2 . A more sophisticated decomposition model, which would include terms not addressed in the current work, would be dedicated to another study.

Conclusion
In this study we used SOO observational data to show that the mid latitude North Atlantic Ocean around the PAP site remained a sink for atmospheric CO 2 in recent years. There is good agreement between SOO and fixed-point observatory data for the 15 year time series presented here. The study expands on previous work and constrains the relative contribution of physical and biological factors to observed changes in seawater DIC and pCO 2 . The annual average surface seawater pCO 2 showed no statistically significant increasing trend in the area around the PAP-SO site between 2002 and 2016 in spite of the observed increase in atmospheric pCO 2 and nDIC. This meant that the carbon sink capacity was maintained, especially due to the summer time pCO 2 deficit driven by biological processes. Consistent with other recent studies, the difference between winter and summer seawater pCO 2 values became larger during the study period, without a similar change for nitrate. Gas exchange, biological production, temperature and alkalinity changes explained over 76% of the annual pCO 2 variability in the study area.
We are now observing the ocean at a higher resolution than ever before. Using floats and unmanned autonomous vehicles has the potential to add the depth dimension to our observational capacity and also to reach areas where meteorological conditions hinder vessel-based observations. In addition, observations have now extended over climatologically significant time scales. In some cases however, even two decades of observations were deemed insufficient to investigate a relationship with decadal variability (Midorikawa et al., 2006). In much of the ocean, long term signals require a much longer period of observations than what is available in this study (McKinley et al., 2016). Continuous monitoring will improve our understanding of the carbonate system and help guide our future predictions.
By the end of the 21st century, the ocean carbon sink is expected to weaken overall since the effect of nutrient limitation will outweigh the effect of a longer growing season (Henson et al., 2013). The amount of CO 2 uptaken will determine the driving mechanisms for the North Atlantic carbon system. It has been modelled that in a low anthropogenic carbon uptake world, SST will be dominating the annual cycle of seawater pCO 2 , while mixing and biological production will be dominant in a high anthropogenic uptake world (Goris et al., 2018). It   Fig. 8. Scatter plots of gas-corrected nDIC versus nitrate (a) and nitrate versus pCO 2 (b). The red, blue and black symbols correspond to values measured or calculated for the 'pre-bloom', 'bloom' and 'post-bloom' months respectively. The blue line is the C:N regression for the 'bloom' period. Subfigure (b) contains a conceptual visualisation of the inter-seasonal changes.
is clear that without a solid understanding of the constrains to the marine carbon cycle, our ability to confidently predict the status of the global ocean carbon sink under increased atmospheric forcing will be limited.

Declaration of Competing Interest
The authors declare no conflict of interest. Sources of funding have been acknowledged.