Dynamics of air–sea CO 2 fluxes in the northwestern European shelf based on voluntary observing ship and satellite observations

Abstract. From January 2011 to December 2013, we constructed a comprehensive pCO2 data set based on voluntary observing ship (VOS) measurements in the western English Channel (WEC). We subsequently estimated surface pCO2 and air–sea CO2 fluxes in northwestern European continental shelf waters using multiple linear regressions (MLRs) from remotely sensed sea surface temperature (SST), chlorophyll a concentration (Chl a), wind speed (WND), photosynthetically active radiation (PAR) and modeled mixed layer depth (MLD). We developed specific MLRs for the seasonally stratified northern WEC (nWEC) and the permanently well-mixed southern WEC (sWEC) and calculated surface pCO2 with uncertainties of 17 and 16 μatm, respectively. We extrapolated the relationships obtained for the WEC based on the 2011–2013 data set (1) temporally over a decade and (2) spatially in the adjacent Celtic and Irish seas (CS and IS), two regions which exhibit hydrographical and biogeochemical characteristics similar to those of WEC waters. We validated these extrapolations with pCO2 data from the SOCAT and LDEO databases and obtained good agreement between modeled and observed data. On an annual scale, seasonally stratified systems acted as a sink of CO2 from the atmosphere of −0.6 ± 0.3, −0.9 ± 0.3 and −0.5 ± 0.3 mol C m−2 yr−1 in the northern Celtic Sea, southern Celtic sea and nWEC, respectively, whereas permanently well-mixed systems acted as source of CO2 to the atmosphere of 0.2 ± 0.2 and 0.3 ± 0.2 mol C m−2 yr−1 in the sWEC and IS, respectively. Air–sea CO2 fluxes showed important inter-annual variability resulting in significant differences in the intensity and/or direction of annual fluxes. We scaled the mean annual fluxes over these provinces for the last decade and obtained the first annual average uptake of −1.11 ± 0.32 Tg C yr−1 for this part of the northwestern European continental shelf. Our study showed that combining VOS data with satellite observations can be a powerful tool to estimate and extrapolate air–sea CO2 fluxes in sparsely sampled area.


Introduction
Continental shelf seas form a complex interplay between the land, ocean and atmosphere, hosting a multitude of biogeochemical processes (Walsh, 1991;Liu et al., 2010) and play a key role in the global carbon cycle (Walsh et al., 1981;Muller-Karger et al., 2005;Bauer et al., 2013). Even though marginal seas occupy only 7 % of global oceanic area, they host enhanced biological activity, which accounts for 15 to 30 % of global oceanic primary production (Gattuso et al., 1998). These productive regions are characterized by enhanced air-sea CO 2 fluxes compared to open oceans (Tsunogai et al., 1999;Thomas et al., 2004) and are particularly vulnerable to anthropogenic forcings such as eutrophication and ocean acidification (Borges and Gypens, 2010;Borges et al., 2010a;Wallace et al., 2014). In a context of climate change, with rising anthropogenic CO 2 levels in the atmosphere and the oceans (IPCC, 2013), it is essential to better constrain carbon cycle dynamics and particularly air-sea CO 2 fluxes. Given the large diversity and heterogeneity of coastal ecosystems, this goal remains challenging. pansion of partial pressure of CO 2 (pCO 2 ) observations over the past decade have allowed the first assessments of the contribution of coastal ecosystems in terms of global air-sea CO 2 fluxes Cai et al., 2006;Chen and Borges, 2009;Cai, 2011). However, extrapolation from local to global estimates still involves large uncertainties and many continental shelf seas remain under-sampled.
Accurate estimates of air-sea CO 2 fluxes in continental shelf seas still suffer from lack of sufficient spatial and temporal coverage. Surveys based on seasonal sampling during oceanographic campaigns and time series at fixed locations are limited due to the large temporal and spatial variability of these systems. The use of voluntary observing ships (VOSs) can improve the coverage of coastal areas at a lesser cost. The recent advances made in this field (Schneider et al., 2006(Schneider et al., , 2014Padin et al., 2007;Omar et al., 2010;Marrec et al., 2014) can be combined with other new approaches. Since the 2000s, pCO 2 predictions based on remote sensing techniques have been successfully developed for open ocean areas (Lefèvre et al., 2002;Ono et al., 2004;Olsen et al., 2004;Rangama et al. 2005;Gledhill et al., 2008;Padin et al., 2009;Chierici et al., 2009Chierici et al., , 2012. These estimates were based on the use of multiple linear regressions (MLRs) to relate surface ocean pCO 2 to sea surface temperature (SST), chlorophyll a concentration (Chl a) and occasionally also mixed layer depth (MLD), sea surface salinity (SSS) or geographical position (latitude and longitude). More complex neural networks techniques using self-organizing maps have also given promising results (Lefèvre et al., 2005;Telszewski et al. 2009;Friedrich and Oschlies, 2009). In continental shelf seas the development of remotely sensed approaches is more challenging because of higher temporal and spatial variability of biogeochemical processes. The complex optical properties of these systems can also impede computations based on satellite ocean-color data. These techniques have nevertheless been used to conduct successful assessments of pCO 2 variability in coastal areas (Lohrenz and Cai, 2006;Salisbury et al., 2008;Borges et al., 2010b;Shadwick et al., 2010;Hales et al., 2012;Jo et al., 2012;Signorini et al., 2013).
To efficiently constrain surface pCO 2 in dynamic shelf seas from remotely sensed data, a comprehensive pCO 2 data set with sufficient spatial and temporal resolution is essential. In addition to a robust intra-annual temporal resolution, acquisition of pCO 2 measurements over several years is necessary in order to take into consideration the important inter-annual variability of biogeochemical processes in coastal seas. From 2011 to 2013, we collected an extensive pCO 2 data set based on VOS observations in the western English Channel (WEC), which is part of the northwestern European continental shelf (NWES). Marrec et al. (2013Marrec et al. ( , 2014 assessed for the first time the seasonal and latitudinal dynamics of pCO 2 and air-sea CO 2 fluxes in the WEC, whereas previous studies of the CO 2 system in the WEC were either based on longitudinal transects (Borges and Frankignoulle, 2003;Padin et al., 2007;Dumousseaud et al., 2010) or a fixed station approach (Kitidis et al., 2012). Air-sea CO 2 flux and pCO 2 dynamics in the North Sea have been well constrained (Thomas et al., 2004Bozec et al., 2006;Schiettecatte et al., 2007;Prowe et al., 2009;Omar et al., 2010), but ocean carbon variability in other part of the NWES as the Celtic Sea (CS) and the Irish Sea (IS) is still poorly documented (Frankignoulle and Borges, 2001).
In this study, we will focus on the WEC, CS and IS provinces of the NWES in order to describe the dynamics of pCO 2 and air-sea CO 2 fluxes in these areas over a decade. We used MLR in WEC to develop algorithms to predict surface pCO 2 and air-sea CO 2 fluxes from remotely sensed SST, chlorophyll a concentrations (Chl a), wind speeds (WNDs), photosynthetically active radiation (PAR) and from modeled mixed layer depth (MLD). We extrapolated the relationships obtained in the WEC based on the 2011-2013 data set (1) temporally over a decade and (2) spatially in the adjacent CS and IS, two regions where pCO 2 data are very sparse. Based on the reconstructed decadal data set, we investigated the variability of pCO 2 and air-sea CO 2 fluxes over the shelf.

Study area
The WEC forms part of the northwestern European continental shelf, one of the world's largest margins. We studied this area from January 2011 with a VOS (Fig. 1) equipped with an autonomous ocean observing system, called Ferry-Box, featuring several sensors (Sect. 3.1, Marrec et al., 2013Marrec et al., , 2014. This area is characterized by relatively shallow depths and intense tidal streams with maximum speeds ranging from 0.5 to 2.5 m s −1 (Pingree, 1980;Reid et al., 1993). Along the French coast (southern WEC -sWEC), where the tidal currents are the strongest, the water column remains vertically mixed (Wafar et al., 1983;L'Helguen et al., 1996), whereas near the English coast (northern WEC -nWEC), where tidal streams are less intense, seasonal stratification occurs . Between these two distinct structures, a frontal zone oscillates, separating well-mixed and stratified waters (Pingree et al., 1975). In this complex hydrographical context, high-frequency measurements from Ferry-Box data allowed us to precisely locate this thermal front and to accurately identify the real extent of each hydrographical province (Marrec et al, 2014).
Satellite SST data (Fig. 2, Sect. 3.2) combined with Ferry-Box measurements allowed us to further define the different hydrographical provinces of the northwestern European continental shelf. Water column characteristics similar to those in the WEC are also observed in adjacent seas, i.e., the Irish Sea (IS) and the Celtic Sea (CS; Pingree and Griffiths, 1978;Pingree, 1980;Holligan, 1981;Simpson et al., 1981;Hill et al., 2008). column is well-mixed and the warmest SST, areas with seasonal stratification. The Ushant front (Pingree et al., 1975;Morin, 1984, Sournia et al., 1990 separates the seasonally stratified southern Celtic Sea (sCS) and nWEC from the permanently well-mixed sWEC. Such a frontal structure is also observed off the Penwith Peninsula (west of Cornwall, UK), around the Land's End (LE) waters; hereafter we refer to these well-mixed waters as LE. St. George's Channel front separates permanently well-mixed southern IS (sIS) waters from the seasonally stratified northern CS (nCS) waters. In addition to the similar hydrographical properties, the WEC, CS and IS also exhibited similar seasonal dynamics and biogeochemical processes Pemberton et al., 2004;Smyth et al., 2009). Based on these observations, we defined six key hydrographical provinces ( Fig. 2) with fixed boundaries, which represent the shifting area of thermal fronts. The use of fixed boundaries allows a direct comparison between the representative provinces. We did not include coastal areas strongly influenced by riverine inputs (Fig. 2) such as the Bristol Channel, coastal Irish waters, surface waters in vicinity of Plymouth and the eastern part of the sIS (which is also seasonally stratified). We chose to study only the southern part of the IS because of the complexity of the northern IS, which has successive stratified, frontal and mixed systems (Simpson and Hunter, 1974) and is influenced by freshwater inputs (Gowen and Stewart, 1995). The study of the permanently well-mixed part of the IS allowed us to apply our algorithm developed for the sWEC to estimate for the first time air-sea CO 2 fluxes in the IS. In the southwest corner of our study area, at the shelf break, internal tides and turbulence favor vertical mixing which sustains biological activity by supplying nutrients to the photic zone (Pingree et al., 1981;Joint et al., 2001;Sharples and Tweddle, 2007). Because the internal tides at the shelf break induce specific biogeochemical properties and our algorithms are not intended to predict surface pCO 2 in this province, we excluded the shelf break region (Fig. 1, southwestern area) from our study area.

FerryBox data sets
From January 2011 to January 2014, a FerryBox system was installed on the Voluntary Observing Ship (VOS) Armorique (Brittany Ferries). This vessel crossed the English Channel between Roscoff (France, 48 • 43 38 N, 3 • 59 03 E) and Plymouth (United Kingdom, 50 • 22 12 N, 4 • 08 31 E; Fig. 1) up to three times a day. The FerryBox continuously measured sea surface temperature (SST), salinity and partial pressure of CO 2 (pCO 2 , from April 2012) along the ferry track with more than 600 crossings with pCO 2 acquisition. Between January 2011 and January 2014, discrete sampling was performed on 57 return crossings between Roscoff and Plymouth with a total of 1026 sampling locations in the WEC. During each cruise, 18 water samples were taken from the FerryBox seawater circuit for the determination of dissolved inorganic carbon (DIC), total alkalinity (TA) and associated  (Marrec et al, 2013). Seawater pCO 2 values were calculated from TA, DIC, temperature, salinity and nutrient concentrations with the CO2SYS program (Pierrot et al., 2006) using the equilibrium constants of CO 2 proposed by Mehrbach et al. (1973), refitted by Dickson and Millero (1987) on the seawater pH scale, as recommended by Dickson et al. (2007). The methods used for the analytical determinations of DIC and TA are described in detail in Marrec et al. (2014) and gave accuracies of ±2 and 3 µmol kg −1 , respectively. Thus, the computed values of pCO 2 from DIC and TA have uncertainties at the lower end of ±6 µatm (Zeebe and Wolf-Galdrow, 2001). Sensors were calibrated and/or adjusted based on these bimonthly discrete measurements as described in Marrec et al. (2014). Based on the comparison between high-frequency pCO 2 data obtained with a Contros HydroC/CO 2 FT sensor and bimonthly pCO 2 data calculated from DIC/TA, we estimated high-frequency pCO 2 measurements uncertainties at the lower end ±6 µatm (Marrec et al., 2014), in the same range as computed values of pCO 2 from DIC and TA. We built a composite monthly data set of in situ SST and pCO 2 data over 3 years based on both high-frequency and bimonthly (twice a month) measurements. We used bimonthly discrete pCO 2 data between January 2011 and April 2012 and high-frequency pCO 2 data from April 2012 to January 2014. The monthly pCO 2 data were then adjusted to reference month July 2012, for MLR computation, using an atmospheric growth rate of 0.15 µatm month −1 (1.8 µatm yr −1 ) and assuming that the surface ocean pCO 2 is growing at the same rate as the atmosphere. This growth rate was calculated from atmospheric pCO 2 at Mace Head (Ireland, see Sect. 3.5 for further details) from January 2003 to December 2013 and is then representative of our study area.

Satellite and other environmental data
Satellite-derived Chl a concentrations (µg L −1 ) were acquired from the Moderate Resolution Imaging Spectroradiometer (MODIS) aboard the Aqua satellite. Daily images were provided by the Natural Environment Research Council (NERC) Earth Observation Data Acquisition and Analysis Service (NEODAAS) at a spatial resolution of 1.1 km. Monthly mean Chl a estimates were computed from January 2003 to December 2013 from these individual images over our study area (Fig. 1). WEC, CS and IS waters are optically complex shelf waters (Joint and Groom, 2000;Darecki et al., 2003;McKee et al., 2007). These shelf seas present both Case 1 and Case 2 optical water types (Morel and Prieur, 1977;Morel et al., 2006) depending on their hydrographical properties (seasonally stratified or homogeneous), the proximity to the coast, and the period of the year. In Case 1 waters, the optical properties are dominated by chlorophyll and associated degradation products as in open ocean waters. In coastal waters, classified as Case 2, suspended particulate sediments and yellow substances of terrestrial origin induce important biases on chlorophyll a concentration estimates and special algorithms have been developed for these waters (Gohin et al., 2002). As shown by Groom et al. (2009), who explain how a coastal station in the nWEC (L4) can be considered as Case 1 or Case 2 depending on various parameters, it is difficult to label our studied provinces as Case 1 or Case 2 waters. However sWEC, LE and IS present more similarities with Case 2 waters, especially during winter, whereas nWEC and CS are closer to Case 1 waters. The NEODAAS provided satellite Chl a estimates based on the OC3 algorithm, more specific to Case 1 waters, and on the OC5 algorithm (Gohin et al., 2002), developed in coastal waters of the eastern English Channel and the Bay of Biscay affected riverine input (Seine, Loire, Gironde). Chl a estimates based on the OC3 algorithm show enhanced Chl a concentrations during winter, particularly in near-coast and in well-mixed provinces, whereas Chl a estimates from the OC5 algorithm tend to underestimate the Chl a concentrations especially during spring and summer (data not shown). We chose to use the OC3 algorithm in this study, which seemed more suitable and more representative of the biological activity dynamics, and we gridded monthly 1.1 km satellite data into 0.05 • × 0.05 • grid cells over our study area. We extracted monthly mean Chl a values along the ship track from January 2011 to December 2013 ( Fig. 3b) to predict pCO 2 based on MLRs (see below).
Satellite-based SST ( • C) data were acquired from the Advanced Very High Resolution Radiometer (AVHRR) instrument. Monthly mean SST estimates were computed from January 2003 to January 2014 from individual images with a spatial resolution of 1.1 km by the NEODAAS. A validation between monthly in situ SST and associated satellite SST showed a robust correlation (R 2 = 0.97, N = 448, p < 0.001 and RMSE = 0.43 • C). We gridded 1.1 km resolution satellite SST into 0.05 • × 0.05 • cells as with all other remotely sensed and modeled parameters.
Photosynthetically active radiation (PAR, in E m −2 d −1 ) data were retrieved from the Ocean Biology Processing Group (McClain, 2009; http://oceancolor.gsfc.nasa.gov). We used the Level 3 monthly merged PAR product from MODIS Aqua. PAR was used as a variable in the MLRs as an indicator of the amount of light available for phytoplankton, which presented inter-annual variation over our study period (Fig. 3c). Based on the observations of L' Helguen et al. (1996), Marrec et al. (2014) suggested that light availability might be an important factor responsible for the strong inter-annual variability of phytoplankton blooms in the sWEC.
Mixed layer depth (MLD), which was one of the variables used in algorithm development for the seasonally stratified nWEC and in the spatial extrapolation of this algorithm in the adjacent CS, was computed from the MARS3D model (Lazure and Dumas, 2008;Berger et al., 2014) developed in the PREVIMER project (Charria and Repecaud, 2014). MLD was defined as the shallowest depth corresponding to a tem- perature or density difference with the surface water higher than δT = 0.5 • C or δDens = 0.125 (Monterey and Levitus, 1997). We compared the model outputs with MLD calculated from the temperature and salinity profile at the fixed station E1 off Plymouth (50.03 • N, 4.37 • W; depth 75 m) from January 2006 to January 2014 ( Fig. S1 in the Supplement). Measurements were undertaken every 2 weeks by the Western Channel Observatory (NERC National Capability of the Plymouth Marine Laboratory and Marine Biological Association, www.westernchannelobservatory.org.uk). Profiles were obtained by a Seabird SBE 19+ with precision for temperature and computed salinity of 0.005 • C and 0.002, respectively. We also compared the modeled MLD in the CS and at the E1 fixed station with Armor-3D L4 Analysis observation products provided by the Copernicus Ma-rine Environment Monitoring Service (ex-MyOcean, http: //marine.copernicus.eu). The latter are combined products from satellite observations (sea level anomalies, mean dynamic topography and sea surface temperature) and in situ (temperature and salinity profiles) data on a 1/2 • regular grid in our study area. Modeled and in situ (from the Western Channel Observatory) MLD at the E1 station showed a robust correlation (R 2 = 0.82, N = 89, Fig. S1). Comparing modeled and Armor-3D L4 Analysis observation products MLD of CS ( Fig. S1) clearly showed the robust approximation of MLD by the model, particularly concerning the start and the end of stratification, despite a small overestimation of the modeled MLD. These comparisons validated the use of modeled MLD in our computations. Modeled MLDs were gridded in the 0.05 • × 0.05 • grid in seasonally stratified provinces and were extracted along the ship track in the nWEC to be included in the pCO 2 algorithms. We chose to use the MLD over depth ratio (MLDr) in the MLR computation instead of MLD. During winter in seasonally stratified areas, the whole water column is mixed. However, depths are not homogeneous (ranging from −20 to −200 m); thus the use of MLD winter values, which corresponded approximately to the bathymetry, would lead to bias in MLR computation. MLD, in our algorithms, was only an indicator of the presence or absence of stratification of the water column, particularly concerning the start and the end of stratification. Figure 3e shows the monthly MLDr ratio in the nWEC between Roscoff and Plymouth.

Biogeosciences
Monthly wind speed data (WND, in m s −1 ) corrected to 10 m height were obtained from four-times-daily wind speed products from the ERA-interim re-analysis project (Dee et al., 2011) produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). We extracted the 0.125 • latitude by 0.125 • longitude global grid wind speed values over the study area, and we gridded these data into our 0.05 • × 0.05 • grid. We used WND data in the computation of the gas transfer velocity of CO 2 (k) as an indicator of wind stress for the calculation of air-sea CO 2 fluxes (Sect. 3.5). Figure 3d shows the monthly WND values used in the algorithm development along the ferry route from 2011 to 2013.

Development of pCO 2 algorithms
We developed two specific algorithms to estimate surface seawater pCO 2 in each of the hydrographical provinces of the WEC (seasonally stratified nWEC and permanently wellmixed sWEC) in order to apply them on a larger spatial and temporal scale in the adjacent Celtic and Irish seas. We used MLRs to predict pCO 2 in each province based on monthly mean values of Chl a, SST, PAR, WND (for the sWEC), MLD (for the nWEC) and from a time variable TI (Eqs. 1 and 2) representative of the seasonality (Friedrich and Os cording to where pCO 2, MLR is the predicted pCO 2 , a 0 is the intercept of the MLR and a i is the coefficient related to each variable p i . In Eq. (2), "Day" is the 15th day of each month (Julian day) and α a value between 0 and 365 chosen by iteration to optimize the seasonal phasing until the minimum standard deviation on residuals and the best correlation coefficient R 2 are obtained by the MLR. All of these parameters were gridded in 0.05 • latitude intervals ( Fig. 3 and 5) between 48.80 • N (off Roscoff) and 50.20 • N (off Plymouth).
The northern latitude limit of 50.20 • N is relatively far from Plymouth in order to exclude effects of freshwater inputs from the Tamar and Plym rivers, which influence the biogeochemical properties of the area  and are not representative of nWEC waters. The WEC is divided into sWEC and nWEC at 49.40 • N from the average position of the thermal front separating the two hydrographical provinces during the period of study ( Fig. 2 and Marrec et al. 2014). MLRs were applied on these gridded monthly values in each province using the "regress" Matlab ® function. The performance of regional algorithms was evaluated by the correlation coefficient R 2 , the adjusted R 2 , the root-meansquare error (RMSE) and the p values (for each of the parameters and for the regression). These parameters represent the capacity (the R 2 and the adjusted R 2 ) and uncertainty (RMSE) of the algorithms to predict pCO 2 . The coefficient of determination R 2 indicates the amount of total variability explained by the regression mode. The adjusted R 2 is the coefficient of determination of the MLR adjusted to the degree of freedom, which depends on the number of variables used.
In each MLR presented in the study, the adjusted R 2 and R 2 were similar; thus only R 2 is presented. MLR coefficients were calculated based on our 3-year data set, and the goal of the study is to apply the algorithms over a decade (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) over the study area (Fig. 1). Based on the atmospheric pCO 2 calculated from Mace Head atmospheric xCO 2 , we calculated a regional anthropogenic increase in atmospheric CO 2 of 1.8 µatm yr −1 , a similar rate as mentioned by Thomas et al. (2008) and Le Quéré et al. (2010). We assumed that the ocean surface pCO 2 increase is trending at the same pace as the atmospheric pCO 2 . Thus we consider a regional increase of seawater pCO 2 of 1.8 µatm yr −1 representative of our study region. When we computed the algorithms, we considered this factor in the computations by adding a correction term X (Eq. 3) on the right term of Eq. (1) (Shadwick et al., 2010;Signorini et al., 2013) with where m (month) is equal to the number of months since July 2012, the middle of our study period (2011)(2012)(2013). For example, in January 2013, X would be equal to (1.8/12) × (+6) and in January 2012 X would be (1.8/12) × (−6). The same reference month (i.e., July 2012) was used to extrapolate the algorithms from January 2003. We normalized each of the variables p i (Eq. 4) using the mean (p i,m ) and standard deviation (p i, SD ) of p i over the study period. The normalized coefficients, which are directly comparable and dimensionless, allowed us to evaluate the relative contribution, or weight, of each of the independent variables (i.e., SST, Chl a, TI, wind speed, PAR and MLD) in the prediction of the dependent variable (i.e., pCO 2, MLR ).

SOCAT and LDEO data
The  (Table 1), mainly from the same southwest-northeast route (Fig. 4) operated principally by three voluntary observing ships (Lüger et al., 2004;Steinhoff, 2010;Schuster et al., 2013;Lefèvre et al., 2014; more details available on the SOCAT website) which crossed these provinces up to twice per month, almost every month from 2003 to 2011. In nCS and IS, the data coverage was sparser and in LE very few data were available. We binned all of these data into the study grid on a monthly basis. We binned the SOCAT and LDEO data into 0.05 • × 0.05 • grid and computed the mean monthly value in each grid cell. The performance of the model was obtained by comparing the mean observed and predicted monthly value per grid cell (see Fig. S2 and 2). For each province, the observed and predicted monthly mean based on this 0.05 • × 0.05 • grid in each province are plotted in Fig. 8.

Calculation of air-sea CO 2 fluxes
The fluxes of CO 2 across the air-sea interface (F ) were computed from the pCO 2 air-sea gradient ( pCO 2 =pCO 2 water -pCO 2 air , µatm) according to  where k is the gas transfer velocity (m s −1 ) and K 0 is the solubility coefficient of CO 2 (mol atm −1 m −3 ) calculated after Weiss (1970). The exchange coefficient k (Eq. 6) was computed as a function of wind speed with the algorithm given by Nightingale et al. (2000) established in the southern bight of the North Sea (SBNS). The SBNS and the WEC present similar environmental characteristics: these two shallow continental shelves are both close to land with high tidal currents controlling the physical structure of the water column.
where u 10 is the wind speed data at 10 m height (m s −1 ), Sc the Schmidt number at in situ SST and C 2 the nonlinearity coefficients for quadratic terms of the gas transfer relationships according to Wanninkhof et al. (2002): where u j is the high-frequency wind speed (m s −1 ), u mean is the monthly mean wind speed (m s −1 ) and n the number of available wind speed data in the month. The empirical relationships of gas transfer with wind speed are nonlinear; then, by using monthly mean wind speeds, the temporal distribution of the wind speeds over the month will affect the gas transfer velocity. We apply the same method reported by Jiang et al. (2008) to account for the intrinsic variability of monthly wind speed data. We also computed gas transfer velocity with the Wanninkhof (1992) and with the Wanninkhof and McGillis (1999) parameterizations for long-term winds to give a range of computed air-sea CO 2 fluxes and we used the nonlinearity coefficients for quadratic and cubic terms of the gas transfer relationship (Wanninkhof et al., 2002;Jiang et al., 2008). Wind speeds along the ferry track and over the study area were extracted from four-times-daily wind speed data corrected at 10 m height from the ERA-interim reanalysis project (Sect. 3.2). Atmospheric pCO 2 (pCO 2 air ) was calculated from the CO 2 molar fraction (xCO 2 ) at the Mace Head site (53 • 33 N, 9 • 00 W; southern Ireland) of the RAMCES network (Observatory Network for Greenhouse gases) and from the water vapor pressure (pH 2 O) using the Weiss and Price (1980) equation. Atmospheric pressure (P atm ) over the study area was obtained from the ERA-interim re-analysis project (Dee et al., 2011).
To estimate the uncertainties of the monthly air-sea CO 2 flux values, we follow the error propagation method exposed in Omar et al. (2007) and used in Lauvset et al. (2013). Uncertainties were expressed by the random error associated with the monthly flux values (σ F ) according to where σ K0 , σ k and σ pCO 2 are the uncertainties associated with K 0 , k and pCO 2 , and Cov xy denotes the covariance between any two parameters x and y. The total uncertainty in k is based on an assumed 1 m s −1 uncertainty in the wind speed data and from the 0.43 • C uncertainty of the AVHRR satellite-based SST (Sect. 3.2), which impact the Schmidt number calculation and ranged between ±13 and ±41 m month −1 . The uncertainty in K 0 is obtained by summing a 0.5 % uncertainty in the computation of K 0 (McGillis and Wanninkhof, 2006)  pCO 2 results from the quadratic sum of the total uncertainty in pCO 2, MLR (15.9 µatm in permanently well-mixed provinces and 17.1 µatm in seasonally stratified provinces; Sect. 4.1) and an assumed ±1 µatm uncertainty in the atmospheric pCO 2 . The covariances were obtained by holding all parameters at their mean values, except SST, and by using the changes in K 0 , k and pCO 2 due to SST to determine these covariances. Depending on the provinces, Cov K0k varied between −12.96 and −27.61, Cov K0 pCO 2 between 3.15 × 10 −5 and 1.05 × 10 −4 and Cov k pCO 2 between −3.94 × 10 −4 and −1.05 × 10 −4 . For further details on the calculation of the terms used in Eq. (7) see Omar et al. (2007). The random error associated with the mean annual flux estimates was calculated as the quadric sum of the uncertainties in the monthly mean fluxes. Thus, air-sea CO 2 fluxes in the paper and figures are given with their computed uncertainties.

Performance of MLR
We performed MLRs to estimate surface pCO 2 based on SST, Chl a, the time variable TI, WND and PAR in the sWEC, and by including MLDr (MLD / depth ratio) and without WND in the nWEC. Table 2 Fig. 6a and b). It is worth noting that when excluding Chl a in the MLR, R 2 decreased by 0.03 and 0.05 and RMSE increased by 1.4 and 2.1 µatm in sWEC and nWEC, respectively. Overall, the RMSE accounted for less than 10 % of the amplitude of the pCO 2 signal (approximately 200 µatm). For each variable and each MLR, we calculated the p values, which were all smaller than 0.001 (not shown in Table 2), meaning that all of the variables were statistically significant in the MLR.
From the normalized coefficients, we calculated the percentages of variability explained by each variable. The seasonal pCO 2 signal followed an average dynamic close to a sinusoidal signal, with maximal values in fall and minimal values in spring, with transitional values in winter and summer. Therefore, the time variable TI, which is a sinusoidal function, contributed to more than half of the variability of the pCO 2 signal, highlighting the strong seasonality observed on this signal (Fig. 5a, Table 2). Beside TI, the most significant variables in terms of relative contribution were SST and PAR, with 22 and 15 % in sWEC and both with 15 % in nWEC, respectively. Chl a contributed 7 and 6 % in the sWEC and nWEC, respectively, a relatively low value considering that, as reported by Marrec et al. (2013), biological processes are the main driver of seasonal pCO 2 variability in the WEC. The addition of MLDr improved the performance of the MLR despite its relatively small contribution compared to the other normalized coefficients. Due to the complexity of the algorithms, a quantitative interpretation of non-normalized coefficients is difficult. For example, according to our model, pCO 2 decreases by 14.4 µatm when SST increases by 1 • C (Table 2). This value contradicts the expected thermodynamic relationship between SST and pCO 2 from Takahashi et al. (1993). The goal of this study was to develop suitable algorithms to predict pCO 2 variability in Biogeosciences, 12, 5371-5391, 2015 www.biogeosciences.net/12/5371/2015/ continental shelf seas by maximizing the performance of the MLR and not to define empirical relationships between the variables and pCO 2 . Figures 5a, b and c show the monthly gridded (0.05 • of latitude) pCO 2 , pCO 2 predicted from MLR coefficients (pCO 2, MLR ) and associated residuals (pCO 2, obs -pCO 2, MLR ) from January 2011 to January 2014 between Roscoff and Plymouth. In the sWEC, pCO 2 values lower than 350 µatm were observed during spring 2011, whereas in the same period in 2012 and 2013, pCO 2 remained close to the atmospheric equilibrium, between 350 and 400 µatm. Thus, the dynamics of pCO 2 in both provinces presented important inter-annual variability. As the pCO 2 simulation by the MLR is mainly driven by a seasonal cycle (TI), which is the same every year, these inter-annual discrepancies can yield bias in the MLR simulation. For the sWEC the MLR model overestimated pCO 2 during spring and summer 2011 (residuals up to 30 µatm, Fig. 5b and c) and underestimated pCO 2 in spring 2012 (residuals down to −50 µatm, Fig. 5b and c) by simulating an average decrease of pCO 2 both years. In Fig. 6c and d, residuals are plotted vs. observed pCO 2 in the sWEC and nWEC, and in Fig. 6e and f monthly mean residuals over each province are plotted vs. months from January 2011 to December 2013. In the sWEC, when observed pCO 2 (pCO 2, Obs ) values were below 350 µatm, as in spring 2011, pCO 2, MLR values were much higher than pCO 2, Obs and residuals were highly negative. In the sWEC, residuals as a function of the observed pCO 2 were not homogeneously distributed, especially during year 2011 (Fig. 6), with high negative residuals when pCO 2 was below 350 µatm and high positive residuals when pCO 2 was over 450 µatm. In the nWEC, the distribution of residuals was more homogeneous (Fig. 6c and d); the less pronounced inter-annual variability was responsible for the better performance of the algorithms (R 2 ) in this part of the WEC. Shadwick et al. (2010) and Signorini et al. (2013) undertook similar studies on the Scotian shelf and the northeastern American continental shelf, respectively. They estimated pCO 2 as a function of SST, Chl a and k, and from SST, salinity, Chl a and a time variable (TI), respectively, using MLR. Based on 14 monthly mean values from a high-frequency data set at a moored buoy, the algorithm developed by Shadwick et al. (2010) attained a R 2 of 0.81 with an associated standard error of 13 µatm. They extrapolated this algorithm over the entire Scotian shelf region to investigate pCO 2 and air-sea CO 2 fluxes from remotely sensed data from 1999 to 2008. Signorini et al. (2013) reported R 2 and associated RMSE ranging from 0.42 to 0.87 and from 22.4 to 36.9 µatm, respectively. They divided the northeastern American continental shelf into five distinct regions according to their physwww.biogeosciences.net/12/5371/2015/ Biogeosciences, 12, 5371-5391, 2015  ical and biogeochemical attributes. Their studies were based on SOCAT surface ocean pCO 2 , and the environmental variables used to predict pCO 2 came from remotely sensed and modeled data. The performances of our MLRs are within the same range as those in these previous studies. We developed our algorithms based on a 3-year data set obtained during highly contrasting years, which contributed to the robustness of our model to predict a representative seasonal cycle of pCO 2 as seen in the nWEC. However, the WEC is a highly dynamic continental shelf ecosystem characterized by strong inter-annual variations. Very exceptional events, inherent to continental shelf areas, remain difficult to simulate with our method, which explain the lower performances of our MLR for the sWEC, as seen from the high residuals observed during 2012. We compared air-sea CO 2 fluxes (Eq. 5) calculated from observed pCO 2 and from pCO 2 simulation ( Fig. 7 and Table 3). Figure 7 shows the air-sea CO 2 flux variation in the sWEC and the nWEC based on pCO 2, obs and pCO 2, MLR from January 2011 to January 2014. Fluxes were computed from the mean monthly pCO 2 of each province, and the standard deviation on MLR fluxes corresponds to MLR fluxes computed plus and minus the RMSE obtained in the respective provinces (Table 2). Seasonal air-sea CO 2 flux cycles were well described by the algorithm-defined pCO 2 , particularly for the nWEC, with both provinces acting as a sink of atmospheric CO 2 during spring and summer and as source of CO 2 to the atmosphere during fall and winter. The interannual variability of pCO 2 observed in the sWEC during spring and summer was also reflected in the flux computations, the fluxes based on MLR overestimating the CO 2 sink in spring 2012. Table 3 reports the annual flux estimates in both provinces based on in situ pCO 2 observations and pCO 2, MLR . On an annual scale, the seasonally stratified nWEC waters acted as a sink of atmospheric CO 2 at a rate of 0.0 to 0.5 mol C m −2 yr −1 based on in situ pCO 2 measurements. Fluxes computed from pCO 2, MLR also indicated that Biogeosciences, 12, 5371-5391, 2015 www.biogeosciences.net/12/5371/2015/ Table 3. Air-sea CO 2 fluxes (in mol C m −2 yr −1 ) calculated from observed pCO 2 and from pCO 2 obtained by MLR along the ferry track in nWEC and sWEC using Nightingale et al. (2000) k parameterization. In the main text and figures, all these fluxes are given with their respective uncertainties as computed in Sect. 3.5. Here we provide in brackets values computed using Wanninkhof et al. (1992) and Wanninkhof and McGillis (1999) k parameterizations to give a range of computed air-sea CO 2 fluxes. the nWEC acts as a sink of atmospheric CO 2 , but we observed some discrepancies between the magnitude of in situ and MLR-based fluxes. The permanently well-mixed sWEC waters acted as a source of CO 2 to the atmosphere from 2011 to 2013 ranging between 0.5 and 0.8 mol C m −2 yr −1 , and annual CO 2 fluxes computed from observed and modeled pCO 2 were in good agreement. The performances of our algorithms to estimate monthly surface pCO 2 allowed us to compute suitable air-sea CO 2 fluxes in WEC provinces during three contrasted years.

Spatial and temporal extrapolation of the algorithms
We applied the previous algorithms (Sect. 4.1) over our study area (Fig. 2) from monthly mean remotely sensed and modeled parameters in the permanently well-mixed sWEC, IS and LE and in the seasonally stratified nWEC, sCS and nCS. The pCO 2 values computed from these variables were averaged by province from January 2003 to December 2013 (Fig. 8). The available SOCAT and LDEO observed pCO 2 data ( Fig. 4 and Table 1) were binned into 0.05 • × 0.05 • grid cells and averaged over the provinces. Observed pCO 2 monthly mean data were superimposed on the algorithm pCO 2 time series (Fig. 8).
For the nWEC (Fig. 8d), the modeled data followed the main features of the seasonal cycle described by the observed data and were in relatively good quantitative agreement (see Fig. S2 for more details). Spring pCO 2 minima were in the same range, despite the discrepancy of timescales. During fall and winter, maximum values were not always reached, suggesting a relatively small overestimation of modeled pCO 2 values at this time. In the sCS (Fig. 8c), where observed data covered most months from 2003 to 2013, predicted data fitted reasonably well the observed pCO 2 during spring and summer. During fall and winter, only few observed pCO 2 data were above equilibrium values, but the model predicted surface water pCO 2 oversaturation compared to atmospheric equilibrium. Despite the limited number of observed data available for the nCS, IS and LE, the model predictions were in good agreement with the observed pCO 2 (Fig. 8 and S2). For the sWEC (Fig. 8e), the predicted pCO 2 values were higher than the observed data, our algorithm thus mainly overestimating the pCO 2 . During spring, the minimal observed pCO 2 values were not always reached by the model. Figure 4 shows that the observed pCO 2 data available for the sWEC were acquired along the Ushant front (Pingree et al., 1975;Morin, 1984;Sournia et al., 1990) at the border of the province (sWEC, nWEC and sCS) delimited based on summer SST (Fig. 2). This border is a frontal zone between well-mixed and stratified systems with enhanced biological activity due to the constant supply of nutrients from the deep layer of stratified systems, especially in summer when the winter nutrient stock is totally depleted (Holligan, 1981;Morin, 1984, Le Fèvre, 1986, Le Boyer et al., 2009. This enhanced productivity might induce biological consumption of CO 2 , which would explain the overestimation of modeled pCO 2 in the frontal zone. The SOCAT and LDEO data were not representative of a homogeneous system, hindering a direct comparison.
www.biogeosciences.net/12/5371/2015/ Biogeosciences, 12, 5371-5391, 2015 Directly comparing monthly mean pCO 2 values obtained from algorithms and the SOCAT and LDEO pCO 2 data could generate an important bias because of the timescale difference between these data sets. Monthly gridded SO-CAT and LDEO data were mainly based on measurements performed at daily scales. Computed pCO 2 values were representative of the average monthly pCO 2 variability, which tends to smooth extreme values obtained at shorter timescales and prevent any observation of short-term processes. Despite this timescale discrepancy the mean differences between observed and predicted pCO 2 were −1 ± 27 in the sCS, −3 ± 25 in the nWEC, 3 ± 26 in the nCS, −7 ± 26 in the IS and −3 ± 29 µatm in LE on an annual scale. Considering the uncertainties relative to the MLR of 17 µatm (Sect. 4.1), these results were very promising and allowed us to validate the extrapolation of our method over our study area. The results obtained in the sWEC, with an annual mean difference of −14 ± 26 µatm, were less promising as explained above and in Sect. 4.1. The comparison of SOCAT and LDEO data with the model predictions provided indications on the MLR performance on a wider spatial scale. For the first time, we thus computed the seasonal and long-term dynamics of pCO 2 and associated air-sea CO 2 fluxes over a decade for this part of the northwestern European continental shelf (Sect. 4.3) despite the relative uncertainties inherent to the method.

Seasonal and biogeochemical controls of pCO 2 in stratified ecosystems
Figures 9 to 11 show the monthly values of Chl a, computed pCO 2 and associated air-sea CO 2 fluxes in the stratified and homogeneous regions of our study area defined in Fig. 2. Based on in situ MLD data at fixed station E1 (Western Channel Observatory of Plymouth, Fig. 1), on MLD from Armor-3D L4 Analysis observation products and on modeled MLD (Sect. 3.2, Fig. S1), we generally observed an onset of stratification in the nWEC and CS from April to October. Modeled MLD data indicated that water column stratification generally started one month earlier and ended one month later in the CS than in the nWEC. The formation of shallow surface layers (≈ 30 m in the CS and 15 m in the nWEC) triggers the initiation of spring phytoplankton blooms in the CS and nWEC (Pingree, 1980). The earlier onset of stratification in the CS than in the nWEC, due to less intense tidal streams (Pingree, 1980), is consistent with the preliminary signs of the spring bloom observed firstly in the CS (Fig. 9) and resulted in pCO 2 values below 350 µatm first observed in CS in April and one month later in the nWEC (Fig. 9). In the CS, pCO 2 values remained below the atmospheric equilibrium until October despite the apparent lack of biological activity in surface waters (Chl a concentrations < 1 µg L −1 ). Subsurface phytoplankton blooms can occur within the thermocline  at the interface with the deep cold water pool, which is not depleted in nutrients (Pemberton et al., 2004;Southward et al., 2005;Smyth et al., 2009;Hickman et al. 2012), and maintained low pCO 2 values. In the nWEC, spring Chl a values remained between 1 and 2 µg L −1 , with particularly elevated values in July, resulting in CO 2 undersaturated waters with respect to the atmosphere until September. The breakdown of stratification and associated remineralization of organic matter started one month earlier in the nWEC than in the CS and resulted in surface waters oversaturated in CO 2 in the two provinces from September and October, respectively. From end of fall to the end of the winter (March of the following year), surface waters remained undersaturated in CO 2 with respect to the atmosphere ( Fig. 8 and 10) in seasonally stratified systems due to dominant thermodynamical control (Marrec et al., 2013).

Seasonal and biogeochemical controls of pCO 2 in permanently well-mixed ecosystems
The interpretation of the seasonal dynamics of pCO 2 linked to Chl a based on satellite observations in the all-year wellmixed sWEC, LE and IS is more complex than in adjacent seasonally stratified systems. In the IS, we obtained abnormally high Chl a satellite estimates based on the OC3 algorithm (Sect. 3.2) most of the year caused by elevated suspended particles and colored dissolved organic matter concentrations (McKee and Cunningham, 2006). Moreover, the areas defined as sWEC and LE are not only representative of homogeneous systems, but they also include tidal mixing frontal zones. These frontal regions host higher biological production than well-mixed systems (Pingree et al., 1975), which enhance the Chl a signal. The latter had a minor contribution (7 %, Table 2) in the computation of pCO 2 in homogeneous systems and did not have a large effect on pCO 2 prediction. Instead, as reported by previous studies (Boalch et al., 1978;L'Helguen et al., 1996;Wafar et al., 1983), the main factor controlling phytoplankton production in homo-geneous ecosystems is the light availability, represented in the MLR by the PAR. The PAR contributed to 16 % of the variability of computed pCO 2 and was concomitant with the observed CO 2 undersaturation from April to July. In June, when day length was the longest and meteorological conditions were generally favorable, we observed peaks in Chl a values associated with the lowest pCO 2 values ( Fig. 8 and  10) of 320 µatm as biology exerted a dominant control on the pCO 2 signal (Marrec et al., 2013). The productive period is shorter in all-year well-mixed systems than in seasonally stratified areas (Marrec et al., 2013(Marrec et al., , 2014. Surface pCO 2 values were below the atmospheric equilibrium from March to July in the sWEC, LE and IS, whereas these patterns were observed from March to September/October in the CS and the nWEC (Fig. 8 and 10). This is also due to the fact that at the end of summer the organic matter remineralization processes started earlier in homogeneous systems and surface seawaters were oversaturated in CO 2 compared to the atmosphere in August in these ecosystems. The pCO 2 values reached maximum values around 450 µatm during fall in sWEC, LE and IS. Similarly to stratified systems, the wellmixed surface waters remained undersaturated in CO 2 with respect to the atmosphere ( Fig. 8 and 10) during the following winter months due to dominant thermodynamical control (Marrec et al., 2013).
The alternate control of the pCO 2 seasonality by biological processes during spring, summer and fall and by thermodynamic in winter-summer transitions is representative of temperate coastal ecosystems in Europe (Borges et al., 2006;Bozec et al., 2005;Bozec et al., 2006). The separation in fixed representative provinces based on physical properties provides insight on the different phases and amplitudes of the pCO 2 seasonality in these poorly explored ecosystems. The sharp boundaries of pCO 2 (Fig. 10) between permanently well-mixed and seasonally stratified systems can appear as surprising, especially between August and October. However, previous studies in the WEC (Marrec et al., 2013(Marrec et al., , 2014 showed important pCO 2 gradient from 50 µatm to 100 µatm between sWEC and nWEC surface waters. In August and September 2014 we performed crossings from a newly exploited VOS, the Pont-Aven (Brittany Ferries), between Roscoff and Cork (Ireland) and had access to pCO 2 data from DIC/TA measurements in the CS (Fig. S3). By comparing these in situ pCO 2 data with mean pCO 2 data along the ferry tracks calculated from our MLR from 2003 to 2013, we clearly observed the presence of these sharp boundaries between sWEC and nWEC provinces and between nWEC, LE and CS waters. We choose to fix the delimitation of these provinces based on the averaged July and August SST from 2003 to 2013 (Sect. 2, Fig. 2). These regional boundaries represent the shifting area of thermal fronts. As we could not estimate pCO 2 using our algorithms in frontal zones, such an approach appeared suitable. Table 4. Annual air-sea CO 2 fluxes (in mol C m −2 yr −1 ) in the seasonally stratified (nCS, sCS and nWEC) and permanently mixed provinces (sWEC, LE and IS) of our study area between 2003 and 2013 and the mean annual fluxes over the decade using Nightingale et al. (2000) k parameterization. Scaled annual fluxes over province areas (Table 2) in Tg C yr −1 were calculated from the mean annual fluxes over the decade. Fluxes in brackets were calculated using Wanninkhof et al. (1992) and Wanninkhof and McGillis (1999)

Variability of air-sea CO 2 fluxes over the shelf and the decade
On an annual scale, the permanently well-mixed sWEC, IS and LE acted as net sources of CO 2 to the atmosphere at a mean rate (from 2003 to 2013) of 0.2 ± 0.2, 0.3 ± 0.2 and 0.3 ± 0.3 mol C m −2 yr −1 , respectively, (Table 4) whereas the seasonally stratified systems acted as net sinks of atmospheric CO 2 , with mean values over 11 years of −0.6 ± 0.3, −0.9 ± 0.3 and −0.5 ± 0.3 mol C m −2 yr −1 for the nCS, sCS and nWEC, respectively (Table 4). Air-sea CO 2 fluxes computed from predicted pCO 2 corroborate the hypothesis of Borges et al. (2005), with permanently well-mixed systems acting as sources of CO 2 to the atmosphere and seasonally stratified systems acting as a sink of atmospheric CO 2 . The only previous flux estimate for the CS was based on a study by Frankignoulle and Borges (2001), reported in Borges et al. (2006), which indicated that the CS acts as sink of CO 2 of −0.8 mol C m −2 yr −1 . In the sCS we obtained an averaged flux value of −0.9 ± 0.3 mol C m −2 yr −1 , which is in agreement with this previous study. Further, we report what is, to the best of our knowledge, the first estimate of air-sea CO 2 flux in the IS of 0.3 ± 0.2 mol C m −2 yr −1 . These values are all on the same order as the mean annual air-sea CO 2 flux value of −1.9 mol C m −2 yr −1 for European coastal waters reported by Borges et al. (2006). The good agreement between the compilations of annually integrated fluxes computed from field measurements by Borges et al. (2006) and the results found in this study further support the robustness of our MLR.
Our study provides a first assessment of the seasonality of pCO 2 and air-sea CO 2 fluxes over 11 years as well as of the inter-annual and multi-annual variability. Monthly surface ocean pCO 2 derived from algorithms ( Fig. 8) showed important inter-annual variability for the seasonal cycle of CO 2 in each province. Monthly air-sea CO 2 fluxes (Fig. 12) followed the same trend as pCO 2 , resulting in significant interannual differences in the intensity and/or direction of annual fluxes (Table 4 and Fig. 12). The IS and LE remained overall annual sources of CO 2 to the atmosphere from 2003 to 2013, except in 2007 when they acted as sinks of atmospheric CO 2 . In the sWEC, the annual 11-year average flux value of 0.2 ± 0.2 mol C m −2 yr −1 corresponds to annual values ranging from −0.7 ± 0.2 to 0.8 ± 0.3 mol C m −2 yr −1 . The sWEC acted as a sink of atmospheric CO 2 in 2006 and 2007, and as a source of CO 2 to the atmosphere or neutral for the other years. In 2007, in permanently well-mixed systems, a particularly intense spring phytoplankton bloom (data not shown) occurred, which resulted in important CO 2 undersaturation and a CO 2 sink. The CO 2 outgassing during fall 2007 was one of the lowest observed over the decade (Fig. 12), due to relatively weak wind speeds at this time (data not shown) and resulting low k values. The association of these two features explained the annual CO 2 sink obtained in 2007. Seasonally stratified systems showed variability in the intensity of annual air-sea CO 2 fluxes but remained sinks of atmospheric pCO 2 over the decade, except the nWEC in 2009. In addition to the changes of ocean-atmosphere pCO 2 gradient, the wind-dependent gas transfer velocity has a strong influence on air-sea CO 2 fluxes. For example, during fall 2009, monthly pCO 2 values were in the same range as the other years ( Fig. 8) but we observed peaks of CO 2 outgassing in response to more intense monthly wind speeds (> 10 m s −1 ). As mentioned above in Sect. 4.1, our method precluded establishment of empirical relationships between the variables and pCO 2 , and it is therefore difficult to quantitatively and  directly interpret the influence of each variable in the pCO 2 simulation. We scaled the mean annual fluxes over province areas (Tables 1 and 4) and obtained air-sea CO 2 fluxes of −0.45 ± 0.21, −0.69 ± 0.23 and −0.07 ± 0.04 Tg C yr −1 in the nCS, sCS and nWEC, and of 0.02 ± 0.03, 0.02 ± 0.02 and 0.06 ± 0.04 Tg C yr −1 in the sWEC, LE and IS, respectively. These fluxes correspond to absorption of −1.11 ± 0.32 Tg C yr −1 over our study area (the total uncertainty was calculated as the quadratic sum of uncertainties in each province). Borges et al. (2006) estimated the CO 2 sink over the European continental shelves at −68.1 Tg C yr −1 , while Chen and Borges (2009) reported air-sea CO 2 flux of −16.1 Tg C yr −1 in the northeastern Atlantic continental shelf region. The contribution of our study area, which represents 5 % of the European continental shelf reported by Borges et al. (2006), appears rather small because of the large extent of well-mixed ecosystems. However, considering the lack of investigation of air-sea CO 2 fluxes in the CS, IS and to a lesser extent in the WEC, our study allowed for the first time the estimation of the spatiotemporal dynamics of air-sea CO 2 fluxes in these provinces of the northwestern European shelf.

Concluding remarks and perspectives
Based on a 3-year data set of pCO 2 measurements acquired on a VOS in the WEC, we estimated surface ocean pCO 2 and air-sea CO 2 fluxes in the northwestern European continental shelf waters using MLRs from remotely sensed SST, Chl a, PAR and wind speed (in the sWEC), from modeled MLD (in the nWEC) and from a time variable TI. For the first time, seasonal and long-term dynamics of pCO 2 and air-sea CO 2 fluxes over this part of the northwestern European continental shelf were evaluated over a decade, despite the relatively high uncertainties inherent to such method. We thus provide the first estimate of air-sea CO 2 flux in the poorly documented CS, IS and WEC. As mentioned above, very few data are currently available in these coastal seas. However, the number of surface in situ pCO 2 data grows exponentially and these data are now easily available on data portals such as SOCAT and LDEO to develop and validate such algorithms. For example, in situ data of the CO 2 system are currently acquired during seasonal cruises within the CANDYFLOSS project (NERC, collaboration National Oceanographic Center/Biological Station of Roscoff) in the CS and the IS, and with a FerryBox system operating between Roscoff (France) and Cork (Ireland). In the future, these data will improve and allow further developments of our algorithms with an adequate division of the shelf area in representative biogeochemical provinces and by developing specific algorithms in each province.
The reconstructed decadal data sets highlighted the importance of multi-annual study of air-sea CO 2 fluxes in continental shelf seas. As mentioned by Keller et al. (2014), it can be difficult to detect relevant trends in the seawater pCO 2 signal, particularly in coastal areas with high inter-and intraannual variability. Beaugrand et al. (2000) and Tréguer et al. (2014) demonstrated that coastal marine systems of western Europe are connected to large-scale North Atlantic atmospheric circulation, the North Atlantic Oscillation (NAO), and there is a consensus that these coastal systems are highly sensitive to natural and anthropogenic climate change (Goberville et al., 2010. Thomas et al. (2008) investigated the influence of the NAO on air-sea CO 2 fluxes in the North Atlantic and suggested that multi-annual variability of the ocean CO 2 system was linked to the NAO phasing. Salt et al. (2013) demonstrated the connection between NAO forcings and pH and CO 2 variability in the North Sea, another shelf sea of the northwestern European continental shelf. We did not attempt an evaluation of the long-term trend of our CO 2 signal as we believe our algorithm needs to be further improved with more in situ data as mentioned above. In the future, a similar approach could be applied on our data set to investigate the possible links between large-scale climatic indices and the multi-annual variability of pCO 2 and airsea CO 2 on this part of the northwestern European continental shelf, which is closely connected to North Atlantic open ocean waters.
The Supplement related to this article is available online at doi:10.5194/bg-12-5371-2015-supplement.