The northern European shelf as an increasing net sink for CO₂

We developed a simple method to refine existing open-ocean maps and extend them towards different coastal seas. Using a multi-linear regression we produced monthly maps of surface ocean fCO₂ in the northern European coastal seas (the North Sea, the Baltic Sea, the Norwegian Coast and the Barents Sea) covering a time period from 1998 to 2016. A comparison with gridded Surface Ocean CO₂ Atlas (SOCAT) v5 data revealed mean biases and standard deviations of 0 ± 26 µatm in the North Sea, 0 ± 16 µatm along the Norwegian Coast, 0 ± 19 µatm in the Barents Sea and 2 ± 42 µatm in the Baltic Sea. We used these maps to investigate trends in fCO₂, pH and air–sea CO₂ flux. The surface ocean fCO₂ trends are smaller than the atmospheric trend in most of the studied regions. The only exception to this is the western part of the North Sea, where sea surface fCO₂ increases by 2 µatm yr⁻¹, which is similar to the atmospheric trend. The Baltic Sea does not show a significant trend. Here, the variability was much larger than the expected trends. Consistently, the pH trends were smaller than expected for an increase in fCO₂ in pace with the rise of atmospheric CO₂ levels. The calculated air–sea CO₂ fluxes revealed that most regions were net sinks for CO₂. Only the southern North Sea and the Baltic Sea emitted CO₂ to the atmosphere. Especially in the northern regions the sink strength increased during the studied period.

1 Introduction

For facing global challenges, such as predicting and tracking climate change, it is important to improve our understanding of the ocean carbon sink and its variability. Open oceans, especially those of the Northern Hemisphere, are relatively well understood and described in their large-scale variability (Gruber et al., 2019; Landschützer et al., 2018, 2019; Fay and McKinley, 2017). Reliable autonomous systems for measuring carbon dioxide partial pressure from commercial vessels were developed in the early 2000s (Pierrot et al., 2009) and have since been deployed on a large number of vessels (e.g., Bakker et al., 2016). This has resulted in sufficient data to develop methods to interpolate the data and to describe large-scale air–sea CO₂ exchange and its variability (Landschützer et al., 2014, 2013; Rödenbeck et al., 2013; Jones et al., 2015). These methods apply a wide variety of approaches, such as linear interpolation, machine learning and model-based estimates. By comparing the different results, it is possible to achieve a good estimate of the uncertainty associated with the respective methods (Rödenbeck et al., 2015).

Thomas et al. (2007)Salt et al. (2013)Salt et al. (2013)Fröb et al. (2019)Omar et al. (2019)Omar et al. (2019)Laruelle et al. (2018)Laruelle et al. (2018)Wesslander et al. (2010)Schneider and Müller (2018)Laruelle et al. (2018)Yasunaka et al. (2018)Laruelle et al. (2018)

Table 1 Overview of trends in surface ocean CO₂ reported in the literature.

Meyer et al. (2018)Kitidis et al. (2019)Parard et al. (2017)Yasunaka et al. (2018)Yasunaka et al. (2018)

Table 2 Overview of air–sea CO₂ fluxes reported in the literature. Negative sign denotes flux from atmosphere to ocean.

Despite coastal seas covering 7 %–10 % of the world's oceans (Bourgeois et al., 2016), their contribution to the oceanic carbon sink is not yet fully constrained. Whether coastal seas are a net sink or source for atmospheric CO₂ and how their role will change in a changing climate is still under debate (Bauer et al., 2013; Laruelle et al., 2010). Compared to the open ocean, longer time series and higher spatial and temporal resolution of the observations are needed in order to capture all relevant coastal processes. Small-scale circulation patterns governed by topographic features; thermal and haline stratification; or mixing through tidal cycles, upwelling or internal waves result in a need for more complex maps with a higher resolution (Bricheno et al., 2014; Lima and Wethey, 2012; Blanton, 1991). These physical drivers are not the only reasons for coastal seas being more complicated to understand. Generally, coastal regions are more productive than open-ocean regions due to better availability of nutrients (e.g., mixing at continental margins, river runoff). While deeper coastal regions are seasonally stratified, shallow regions are vertically mixed, allowing for exchange between the benthic and pelagic parts of the ecosystem (Griffiths et al., 2017; Wollast, 1998). Together with strong gradients of productivity this leads to spatial and temporal heterogeneity in surface CO₂ content.

The different maps developed for describing the open-ocean surface pCO₂ (CO₂ partial pressure) dynamics and air–sea CO₂ fluxes are not directly applicable in coastal regions. Many exclude data from continental shelves completely while all of them have too coarse a spatial resolution (typically between 1 and 5^∘) to properly resolve coastal seas with their small-scale variability. A few studies have described coastal carbon dynamics, but most of them have strong regional or temporal limitations. Table 1 shows an overview of studies with estimated pCO₂ trends over the northern European shelf, while Table 2 presents available flux estimates. Laruelle et al. (2017) used a neural network approach to produce a global pCO₂ climatology for coastal seas, describing more distinct seasonal variability in the Northern Hemisphere than in the southern Pacific and Atlantic. A global climatology covering both open-ocean and coastal regions was recently constructed by combining this product with the open-ocean product of Landschützer et al. (2016, 2020). Laruelle et al. (2018) published trend estimates in regions with high data coverage based on winter data spanning up to 35 years. Their findings is that the pCO₂ rise in coastal regions tends to lag the atmospheric rise in CO₂. However, few studies attempted to constrain coastal air–sea fluxes before. Kitidis et al. (2019) estimated fluxes between 0 and −15 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in the North Sea, depending on the season (more negative during summer than during winter) and the region (more negative fluxes in the northern North Sea compared to the south). For the Baltic Sea, Parard et al. (2016, 2017) used a neural network approach to produce surface ocean pCO₂ maps from 1998 to 2011 and estimated an air–sea flux of 1.2 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. Yasunaka et al. (2018) estimated a flux of 8–12 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in the Barents Sea and along the Norwegian Coast using a self-organizing map technique. Most of the other available studies on the trends in coastal pCO₂ are based on data from either summer or winter. Estimates based on summer-only data typically show large interannual variations (Thomas et al., 2007; Salt et al., 2013), which led to the conclusion that here the interannual variability masks the actual long-term trend. The approach to use winter-only data (Fröb et al., 2019; Omar et al., 2019), however, is based on the assumption that during this season the influence of biological processes is negligible and therefore winter data can be used to establish a baseline trend. However, also using winter-only data has its drawbacks. In particular the choice of which months to include can cause biases, and the optimal selection can differ from region to region.

In this study we present a new approach to develop monthly fCO₂ (CO₂ fugacity) maps based on already existing open-ocean pCO₂ maps in four example regions: the North Sea, the Baltic Sea, the Norwegian Coast and the Barents Sea. A multi-linear regression (MLR) was used to fit driver data against fCO₂ observations. Based on the resulting fCO₂ maps and a salinity–alkalinity correlation we also produced monthly maps of coastal pH. The performance of the produced maps was evaluated, and the maps were then used to investigate trends in coastal fCO₂ and pH in the entire region from 1998 to 2016. Finally, we used the fCO₂ maps to determine the air–sea CO₂ exchange and its temporal and spatial patterns.

2 Method

2.1 Study area

Figure 1 The study area and the location of the four different regions: North Sea (purple), Norwegian Coast (red), Barents Sea (green) and Baltic Sea (blue).

This work focuses on the northern European continental shelf and marginal seas. As we want to show the performance of the MLR method we picked a number of regions with very different characteristics: the North Sea, the Baltic Sea, the Norwegian Coast and the western Barents Sea (Fig. 1). We decided to concentrate on these regions because (1) the data coverage in these regions is fairly high and (2) the authors have strong knowledge of the specific regions. This is important in order to properly evaluate the maps and to assess whether or not the output is realistic. The four regions were defined based on the COastal Segmentation and related CATchments (COSCAT) segmentation scheme (Laruelle et al., 2013). The threshold for defining a region as coastal sea was set to a depth limit of 500 m. By using this definition, we produce an overlap to the open-ocean maps, allowing our maps to be merged with the open-ocean maps. Please note that this study concentrates on the continental shelf area. The near-coastal zones (e.g., intertidal zones) are not included due to the limited availability of driver data in these regions.

2.2 Data handling

Figure 2 The number of months with fCO₂ data from SOCAT v5 in each grid box. The data cover a range of 20 years (240 months).

Cavalieri et al. (1996)Center (2006)Rödenbeck et al. (2014)Landschützer et al. (2017)

Table 3 Products used as driver data in the MLR and the maps.

The CO₂ data used in this study were extracted from Surface Ocean CO₂ Atlas (SOCAT) version 5 (Bakker et al., 2016). Their coverage is shown in Fig. 2. A newer version of the SOCAT database (SOCATv2019) was used for validating the maps against independent data. An overview over the reanalysis products used as driver data is given in Table 3. We use as basic driver data sea surface temperature (SST), sea surface salinity (SSS), chlorophyll a concentration (Chl a), mixed layer depth (MLD), bathymetry (BAT), distance from shore (DIST), ice concentration (ICE) and the change in ice concentration from the month to month (prior to current). Chl a values during the dark winter season were set to 0. In addition to the reanalysis data, pCO₂ values from the closest coastal grid cell of the open-ocean map were used as a driver in our MLR. We neglect the approximately 1 µatm difference between partial pressure (reported in the mapped products) and fugacity of CO₂ (reported in SOCAT) as it is much smaller than the accuracy of the data extracted from SOCAT v5 (2 to 10 µatm) and the uncertainty associated with the open-ocean maps. The open-ocean pCO₂ values were extracted from two different products (Rödenbeck et al., 2014, version oc_v1.5; and Landschützer et al., 2017, 2016, version 2016). Rödenbeck et al. (2014) is based on a data-driven diagnostic model of mixed layer ocean biogeochemistry fitted against surface pCO₂ observations, while Landschützer et al. (2016) uses a two-step neural network (a feed-forward network coupled with self-organizing maps, FFN-SOM) trained with the pCO₂ observations. Please note that the Rödenbeck open-ocean map contains data in coastal grid boxes, while the Landschützer open-ocean map is restricted to the open-ocean regions. The MLR models based on these two are called MLR 1 (based on the coastal pCO₂ values from the Rödenbeck map) and MLR 2 (based on the nearest open-ocean pCO₂ values of the Landschützer map), respectively. To determine the extent to which the regressions benefit from the information in the open-ocean maps, a third MLR, MLR 3, was determined. Here, we do not use any of the open-ocean maps as a driver, but to account for the annual rise in CO₂, year is included in the set of driver data.

For preparing the input data for the MLR, observations closest to each SOCAT fCO₂ data point in time and space were extracted from the 3-D fields with the driver data. This produces, for each of the driver data, a vector as long as the SOCAT fCO₂ observations. After this, the fCO₂ data as well as all extracted driver data were binned on a monthly 0.125^∘ × 0.125^∘ grid covering 1997 to 2016. These procedures ensure that the driver data have the same bias in space and time within each grid box as the fCO₂ data. If a grid box for example only contains fCO₂ observations from the first week of the month and the northwestern corner, we make sure that also the gridded driver data only contain values from the first week and the northwestern corner of the grid box, and not an average over the entire month and grid box. This is mostly important for the chlorophyll driver data, which are available at a very high resolution compared to the fCO₂ maps produced in this work. These driver data were used for determining the MLRs.

For producing the final maps, a second set of the driver data was prepared, which is called field data in the following. Here the driver data were directly regridded to a monthly 0.125^∘ × 0.125^∘ grid, providing full spatial and temporal coverage and a homogeneous average in each grid box. The field data were used to produce the fCO₂ maps based on the MLR equations.

2.3 Multi-linear regression

Table 4 Driver used in the different regressions.

The multi-linear regression models were constructed by forward and backward stepwise regression using the driver data as predictor variables to model the fCO₂ observations. In each step of this regression procedure, the model's tolerance to addition or exclusion of a variable is tested. This decision on whether to add or remove a term is based on the p value of the F statistic with or without the term in question. The entrance tolerance was set to 0.05 and the exit tolerance to 0.1. The model includes constant, linear, and quadratic terms as well as products of linear terms. Equation (1) gives the basic equation, with X₁ … X_n being the driver data and a₁ … a_nn the regression coefficients.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}y=& a_{\mathrm{0}}+a_{\mathrm{1}}\cdot X_{\mathrm{1}}+\mathrm{\dots }+a_n\cdot X_n+a_{\mathrm{12}}\cdot X_{\mathrm{1}}{X}_{\mathrm{2}}\\ & +\mathrm{\dots }+a_{mn}\cdot X_m{X}_n+a_{\mathrm{11}}\cdot X_{\mathrm{1}}^{\mathrm{2}}+\mathrm{\dots }+a_{nn}\cdot X_n^{\mathrm{2}}\end{array}\end{array}$$

The pCO₂ value of the respective open-ocean maps was used for MLR 1 and MLR 2, while year was always used as a driver variable in MLR 3. Inclusion of stationary drivers (such as month, latitude and longitude) in the MLR increased the performance of MLR 2 and MLR 3. However, these were still not better than MLR 1, and we therefore decided to limit this analysis to dynamic parameters. Using dynamic drivers only assures a dynamic description of the conditions in the field and gives us the possibility to reproduce changes caused by a regime shifts, for example the ongoing Atlantification of the Barents Sea (Oziel et al., 2017; Lind et al., 2018).

2.4 Validation

The three linear fits were compared to each other in each region by taking into account the R² and the root mean square error (RMSE) of the fit, as well as the Nash–Sutcliffe method efficiency (ME) (Nondal et al., 2009). The method efficiency compares how well the model output (E_n) fits the observations (I_n) for every data point n to how well a simple monthly average ($\stackrel{\mathrm{‾}}{I}$) would fit the observations:

$$\begin{array}{}\text{(2)}& \text{ME}={\displaystyle \frac{{\sum }_n(I_n-E_n{)}^{\mathrm{2}}}{{\sum }_n(I_n-\stackrel{\mathrm{‾}}{I}{)}^{\mathrm{2}}}}.\end{array}$$

A method efficiency > 1 means that using just monthly averages of all data in the region would fit better to measured data than the respective model. Generally, a method efficiency > 0.8 is considered bad. Besides the statistics of the fit itself, the final maps were also compared to the gridded SOCAT v5 data, resulting in an average offset and standard deviation (SD). In order to compare the maps against data that were not used to produce the maps, we predicted the fCO₂ for the years 2017 and 2018 (i.e., we applied the trained multi-linear model to driver data from 2017 and 2018) and compared these maps to fCO₂ observations in SOCAT v2019, gridded on a monthly 0.125^∘ × 0.125^∘ grid. We also compare the maps directly with observations from repeated sampling locations in the North Sea and the Baltic Sea.

2.5 Ocean acidification

For calculating the pH, alkalinity (AT) was estimated in the North Sea, along the Norwegian Coast, and in the Barents Sea via a salinity–alkalinity correlation following Nondal et al. (2009). Alkalinity describes the capacity of the sea water to buffer changes in pH. As the concentration of most of the weak bases in seawater is strongly dependent on the salinity, alkalinity can in many regions be estimated from salinity. However, in regions with a high amount of organic bases in seawater, for example in strong blooms or at river mouths, deviations from the alkalinity–salinity relationship can occur. The carbonate system was calculated using the CO2SYS program (van Heuven et al., 2009) with carbonic acid dissociation constants of Mehrbach et al. (1973) as refitted by Dickson and Millero (1987), ${\mathrm{KSO}}_{\mathrm{4}}^{-}$ dissociation constants following Dickson (1990) and the boron–salinity relation following Uppström (1974). For the Baltic Sea, we did not calculate pH as the alkalinity–salinity relationship in this region is complex due to different AT–S relations in different sub-regions of the Baltic Sea and a non-negligible increase in AT over the last 25 years (Müller et al., 2016).

2.6 Calculation of trends

For calculating trends of fCO₂ and ocean acidification, the data in every grid box were deseasonalized by subtracting the long-term averages of the respective months. Then a linear fit was applied to the deseasonalized time series. For illustrating the influence of interannual variability we calculated the trend for different time ranges. As a time range less than 10 years barely resulted in significant trends, we decided to limit the trend analysis to starting years 1998 through 2006 and ending years 2008 through 2016.

2.7 Flux calculation

The air–sea disequilibrium was calculated as the difference between our mapped fCO₂ values and atmospheric fCO₂ in each grid cell and time step. The atmospheric fCO₂ was determined by converting the xCO₂ from the NOAA marine boundary layer reference product from the NOAA GMD Carbon Cycle Group into fCO₂ by using monthly SST and SSS data (Table 3) and monthly air pressure data from the NCEP-DOE Reanalysis 2 (Kanamitsu et al., 2002). We calculated the air–sea CO₂ flux (F) according to Eq. (3), such that negative fluxes are into the ocean. The gas transfer coefficient k was determined using the quadratic wind speed (u) dependency of Wanninkhof (2014) (Eq. 4). The Schmidt number, Sc, was calculated according to Wanninkhof (2014) and the solubility coefficient for CO₂, K₀, following Weiss (1974).

$$\begin{array}{}\text{(3)}& {\displaystyle }F& {\displaystyle }=k\cdot K_{\mathrm{0}}\cdot (f{\mathrm{CO}}_{\mathrm{2},\mathrm{sw}}-f{\mathrm{CO}}_{\mathrm{2},\mathrm{atm}})\text{(4)}& {\displaystyle }k& {\displaystyle }=a_q\cdot \langle u^{\mathrm{2}}\rangle \cdot {\left({\displaystyle \frac{\text{Sc}}{\mathrm{660}}}\right)}^{-\mathrm{0.5}}\end{array}$$

In our calculations, we used 6-hourly winds of the NCEP-DOE Reanalysis 2 product. The coefficient a_q in Eq. (4) is strongly dependent on the wind product used (Roobaert et al., 2018). We determined it to be a_q = 0.16 cm h⁻¹ for the 6-hourly NCEP 2 product following the recommendations of Naegler (2009) and by using the World Ocean Atlas sea surface temperatures (Locarnini et al., 2018). The barrier effect of sea ice on the flux was taken into account by relating the flux to the ice cover extent following Loose et al. (2009). As the gas exchange in areas that are considered 100 % ice covered from satellite images should not be completely neglected, we use a sea ice barrier effect for a 99 % sea ice cover in all grid cells where the sea ice coverage exceeded 99 %.

3 Results

3.1 Maps of fCO₂

Table 5 Statistical evaluation of the MLR 1, MLR 2 and MLR 3 in comparison to the open-ocean maps of Rödenbeck et al. (2015) and Landschützer et al. (2017) for each region. The data for the open-ocean map of Landschützer et al. (2017) are in parentheses since this is based on an extrapolation of the nearest open-ocean grid cell towards the coast. The number of grid cells containing data is given behind the region abbreviations.

The skill assessment metrics for MLR 1, MLR 2 and MLR 3 are presented in Table 5. It shows the R² and RMSE of the fit, the ME, and the average offset and SD to the gridded SOCAT data. The coefficients for MLR 1, MLR 2 and MLR 3 are provided in the Supplement. The MLRs substantially improve the predictions of the open-ocean maps in all studied regions, showing a better average offset and SD to SOCAT v5 and ME than the coarser-resolution open-ocean maps (for example the Rödenbeck map: North Sea, 0 ± 95 µatm; Norwegian Coast, 2 ± 17 µatm; Barents Sea, 22 ± 40 µatm; Baltic Sea, 4 ± 48 µatm; MLR1: North Sea, 0 ± 26 µatm; Norwegian Coast, 0 ± 16 µatm; Barents Sea, 0 ± 19 µatm; Baltic Sea, 2 ± 42 µatm). In all regions MLR 1 has the best model efficiency, the highest R² and the smallest RMSE of the fit, while these statistics are worse for MLR 2 and MLR 3. This can be explained by the fact that the Rödenbeck map contains information about the continental shelf and the Barents Sea, while for MLR 2 the closest open-ocean grid cell of Landschützer et al. (2017) was used. The fact that MLR 3 showed the weakest performance shows the value of using information from the open-ocean maps in the regression.

Figure 3 Average regional differences between MLR 1 and gridded SOCAT v5 data, the Rödenbeck map and gridded SOCAT v5 data, MLR 1 and the Rödenbeck map, and MLR 3 and the gridded SOCAT v5 data (from left to right).

Figure 3 shows, from left to right, the spatial distribution of the average difference between the predicted fCO₂ by MLR 1 and the gridded SOCAT v5 data, the Rödenbeck map and the gridded SOCAT v5 data, the difference between MLR 1 and the Rödenbeck map, and, for comparison, the difference between MLR 3 and the SOCAT v5 data. In the North Sea, MLR 1 seems to slightly overestimate the fCO₂ in the constantly mixed region at the entrance of the English Channel and the area off the Danish North Sea coast. In the Baltic, MLR 1 generally describes the spatial variability in fCO₂ well. However, in the Gulf of Finland it usually predicts fCO₂ values that are too low during May/June while it slightly underestimates events of very high fCO₂ in December/January. Regardless, the spatial biases vs. SOCAT are clearly smaller for MLR 1 than for the original Rödenbeck map. Further, as the predictions of MLR 2 and 3 are clearly inferior to those of MLR 1 (Table 5 and Fig. 3, MLR 3 only), we will use MLR 1 results for the further analyses. An extended validation of the MLR 1 maps can be found in the discussion section.

Figure 4 The average fCO₂ of MLR 1 (1998–2016) for 1 example month in each season (February, May, August and November).

Figure 4 shows the monthly averages of fCO₂ produced by MLR 1 for February, May, August and November. In all regions, the highest fCO₂ values occur in the winter, while the lowest fCO₂ values occur in summer. The largest seasonal cycle could be observed in the Baltic Sea, where fCO₂ reached well below 200 µatm in midsummer and over 500 µatm during the winter.

We notice that the gradients that exist between the grid cells in the Rödenbeck map are still visible in our maps in some regions, for example the sharp gradient in the southern North Sea in February or the east–west and north–south gradients in the entire North Sea in August. Such gradients are also evident in directly mapped pCO₂ data (Kitidis et al., 2019); however, here they are strongly meridional and latitudinal in their extent. As such, while these gradients do reflect actual features of the fCO₂ distribution in the North Sea, their specific shape here is also a consequence of the influence of the Rödenbeck maps on our estimates – from the use of these maps as a driver in the MLR and their importance in improving the statistical performance vs. the MLR that did not use these values as a driver (MLR 1 vs. MLR 3, Table 5). Also, they do reflect the uncertainty of – and our level of confidence in – the estimated pCO₂ values, being approximately similar to or slightly larger than the RMSE of MLR 1 (Table 5). Any smoothing would be completely artificial and, while being more visually pleasing, would not better reflect the truth in any meaningfully quantifiable extent. We have therefore chosen to leave them untouched. These gradients are therefore also visible in subsequent pH and trend maps.

3.2 Maps of pH

Figure 5 The average pH based on MLR 1 (1998–2016) for 1 example month in each season (February, May, August and November).

The monthly average of pH calculated from MLR 1 fCO₂ ranges from about 8 during winter to 8.15 during summer in the North Sea and at the Norwegian Coast (Fig. 5). Towards the Barents Sea the pH maximum during summer increases to 8.2. The pH of 8.00–8.15 in regions with a large influence from the Atlantic, such as the northern North Sea and the Norwegian Coast, is in good agreement with the range of pH determined for the open North Atlantic (Lauvset and Gruber, 2014; Lauvset et al., 2015). In the North Sea, the pH is in the same range as reported in Salt et al. (2013), and it also shows the same distribution in August/September, with higher pH in the northern North Sea and lower pH in the southern part.

4 Discussion

4.1 Performance of the pCO₂ maps

The performance of the MLR and the maps is evaluated in different ways: (1) using the R² and the RMSE of the fit; (2) using the average deviation and its SD, as well as the ME between the produced fCO₂ maps and the gridded observations as a regional average; (3) showing the median deviation between the MLR and the gridded observations on a monthly level; and (4) by comparing the data from the fCO₂ maps to observations from two time series stations. (2)–(4) will be shown for the time period covered by the driver data (1998–2016) and for the prediction of the fCO₂ values for 2017 and 2018. These predicted values are compared with data from the newest SOCAT release (SOCATv2019) and provide a comparison with an independent dataset. Please note that the comparability of the model performance between the different regions is limited. All statistical parameters used are influenced by characteristics that can vary substantially between the different regions, such as range of the data, their variability or the amount of grid cells with data. For example, in a diverse region with many measurements the amount of variability captured by these measurements is most likely larger and thus will lead to a weaker correlation.

Generally, the uncertainty of MLR 1 is in the same range as in other studies (Laruelle et al., 2017; Yasunaka et al., 2018) mapping coastal fCO₂ dynamics: 25 µatm in the North Sea, 16 µatm along the Norwegian Coast, 12 µatm in the Barents Sea and 39 µatm in the Baltic Sea (based on the RMSE in Table 5). In the Baltic Sea, which has a large variability in itself, Parard et al. (2016) obtained lower SDs through dividing the area in smaller sub-regions.

One clear drawback of the here presented MLR 1 is the clearly visible grid pattern of the open-ocean pCO₂ product that was used as input data with its grid size of 5^∘ × 4^∘, mentioned in Sect. 3.1. There are two ways how one could get rid of this artifact in a future release. A finer resolution of the open-ocean maps used will lead to a better representation of the actual gradients in our mapped product. Rödenbeck et al. just released a newer, finer resolution of their open-ocean product (2.5^∘ × 2^∘) that we intend to use in a future version of this data product. However, this will not be sufficient to eradicate the artifact completely. Another approach, running the MLR without an open-ocean pCO₂ product, can provide a coastal pCO₂ product without this artifact. While in principle it is preferential to have coastal maps that are independent of the open-ocean products, MLR 3, which is running without open-ocean pCO₂ as a driver, clearly did not reach the same accuracy as MLR 1 (Table 5). New and better input fields or a different regression method could help improve the independent coastal maps in the future. Another impact that the open-ocean pCO₂ product of Rödenbeck et al. can have on MLR 1 is the potential introduction of patterns from regions further away as the spatial correlations used in producing the Rödenbeck et al. pCO₂ just ignore land barriers. However, the influence of these spatial correlations is relatively small in regions with a high data density (as the European shelf) and the multi-linear regression used to produce MLR 1 corrects for this.

The seasonal differences between MLR 1-determined values and the SOCAT v5 data for each region are shown in Fig. 6. This comparison shows a very good agreement. For MLR 1, the seasonal variations of the median bias are small in the North Sea, along the Norwegian Coast and in the Baltic Sea. In the Barents Sea, however, the bias varies seasonally. Here, MLR 1 slightly underestimates the fCO₂ in winter and early spring, while it overestimates the fCO₂ in summer. In all other regions, the median seasonal bias is smaller than the uncertainty of the maps. The larger seasonal bias in the Barents Sea is most likely caused by the larger seasonal bias in the number of available observations. There are no data available in October, December and January.

Figure 6 Boxplots showing the median deviation of MLR 1 from the gridded SOCAT 5 data for each region (red line). The boxes show the upper and lower 75th percentiles. Ninety-nine percent of the data lie within the range of the purple whiskers. Extremes are shown as gray crosses.

Figure 7 Boxplots showing the median deviation between MLR 1 (based on observations until 2016) and measured fCO₂ values in 2017 and 2018. The boxes show the upper and lower 75th percentiles. Ninety-nine percent of the data lie within the range of the purple whiskers. Extremes are shown as gray crosses. The numbers of grid cells with data available were 5047 for the North Sea, 1543 for the Norwegian Coast, 2312 for the Barents Sea and 5414 for the Baltic Sea.

Figure 8 Time series of ship-of-opportunity (SOOP) data from Nuka Arctica (a, blue) and Finnmaid (b, blue) compared with MLR 1 at the same location (red). In light blue the predictive MLR output for the years 2017 and 2018 is shown.

When comparing all observations from the years 2017 and 2018 to the predictions by the MLR 1, we find a good agreement in the North Sea (2±20 µatm) and no seasonal bias (Fig. 7). In the other regions, the agreement is somewhat reduced compared to the years 1997–2016 (−9 ± 39 µatm (Norwegian Coast), −5 ± 29 µatm (Barents Sea) and 28 ± 58 µatm (Baltic Sea)). In these regions we also observe a seasonal bias in the years 2017 and 2018. At least for the Baltic Sea this could be a result of the extraordinary warm and dry summer in 2018 that led to very low fCO₂ values (see Fig. 8 and the data in SOCAT, Bakker et al., 2016). Please note that for this comparison the MLR was extrapolated in time. Only observations until December 2016 were used to produce the MLR. Therefore accuracy of the maps itself is reduced.

In a second test to investigate to which extent MLR 1 can reproduce observations we compared the MLR output with time series data from two voluntary observing ship lines in two very different regions with a good data coverage: M/V Nuka Arctica in the northern North Sea (0–2^∘ E, 58–60^∘ N) and M/V Finnmaid in the Baltic Sea (23–24^∘ E, 59–60^∘ N) (Fig. 8). The agreement between the MLR 1 and the observations is very good. MLR 1 reproduces the general seasonality and some of the interannual variability, also in the years 2017 and 2018, the observations of which were not used in the regression.

When performing interpolation exercises it is always important to be aware of the fact that the resulting maps might come with biases and do not represent all regions equally well. While the here presented maps give a good general overview about the surface ocean fCO₂ variability in regions with a relatively large amount of data, the reliability, however, is limited in regions where the data coverage is more scarce. This is especially the case when the region with scarce data coverage is showing different characteristics in, for example, temperature and salinity, compared to the rest of the region. One such example is the Gulf of Bothnia in the Baltic Sea region where almost all data used to derive the MLR is from south of 60^∘ N, i.e., not in the Gulf of Bothnia but in the Baltic Proper and western Baltic Sea (see Fig. 2). The MLR method can also lead to unrealistic extreme values and even negative fCO₂. Some such values occur in the northeastern Barents Sea as well as in parts of the Baltic Sea (about 0.01 % of the grid cells in each region). As pH cannot be calculated for negative fCO₂, we excluded all negative fCO₂ values for the calculation of pH. Excluding the negative values resulted in a change in the average fCO₂ of 0.05 µatm (Baltic Sea) and 0.3 µatm (Barents Sea) so they are of negligible importance for the flux estimates. While the negative values are easy to spot and discard, there are most likely other unrealistically low values in spring and summer in the very north and northeastern Barents Sea as well as some parts of the Baltic Sea. However, there are no data available in SOCAT v5 or available elsewhere to validate this.

All regions with questionable fCO₂ are also questionable in their pH data. There is a number of very high pH regions in the Barents Sea (Fig. 5) that are associated with also very low fCO₂ (Fig. 4) that might not be realistic. In addition, estimated pH values in low-salinity regions where the actual alkalinity–salinity deviates strongly from the Nondal et al. (2009) one used here (e.g., river mouths in the southern North Sea or the Skagerrak) should be interpreted with caution.

4.2 Trends in fCO₂ and pH

Trends in surface ocean fCO₂ in coastal regions are often difficult to assess because of the scarcity of data relative to the highly dynamical character of these regimes and their large interannual variability. For example, the start of the productive season can range from February to April even within a small area, such that even restricting the analysis to specific seasons (e.g., winter) can be challenging. Also, due to a lack of data, especially winter data, most observational studies are based on summer data. Further, the fact that these measurements typically do not take place every year adds even more uncertainty to the estimated trend, as interannual variability can mask the trend signal.

Figure 9 The trend in surface ocean fCO₂ estimated from deseasonalized fCO₂. Panel (a) shows the spatial distribution of the trend over the time period from 1998 to 2016. Grid boxes without a significant trend are denoted with a black dot. Panels (b)–(e) show the trends in different time periods in four regions, from the various years on the y axis to the various years on the x axis. Non-significant trends were left blank. Significant trends in sea surface temperature are indicated with crosses/circles. The color bar is centered on the approximate annual fCO₂ rise in the atmosphere (2 µatm yr⁻¹).

The monthly maps of fCO₂ from 1998 to 2016 enable us now to estimate the trend in surface ocean fCO₂ for the entire region, equally distributed over the seasons (Fig. 9, left). All trends were computed from deseasonalized data. The interannual variability of the trend estimates in each region is shown in the panels on the right hand side in Fig. 9. We exclude the northern Baltic Sea from the trend map because we do not expect to have a robust trend estimate in that region as there are only very few data from that region in the regression. Based on the linear regression the significant trends in fCO₂ have an average uncertainty of 0.5 µatm yr⁻¹ (North Sea), 0.4 µatm yr⁻¹ (Norwegian Coast), 0.4 µatm yr⁻¹ (Barents Sea) and 0.7 µatm yr⁻¹ (Baltic Sea), while the shorter time periods shown have a higher uncertainty; no time periods longer than 1998–2016 (for which the given uncertainties of the trend apply) are shown. For pH trends, the average uncertainties of the regressions over 1998–2016 are 5 × 10⁻⁴ (North Sea) and 7 × 10⁻⁴ (Norwegian Coast and Barents Sea).

In most of the regions addressed in this study, the trend in the surface ocean is lower than the trend in atmospheric xCO₂ (global average 2.02 ppm yr⁻¹ (Cooperative Global Atmospheric Data Integration Project, 2015)). Trends exceeding the atmospheric values in the period from 1998 to 2016 can only be observed at the entrance of the English Channel, in Storfjorden/Svalbard and the Gulf of Finland (2.5–3 µatm yr⁻¹). It has to be noted that there was almost no measured fCO₂ as MLR input from Storfjorden. Therefore, these trends are highly uncertain. The trend in the western North Sea is only slightly lower than the trend in the atmosphere (1.5–2 µatm yr⁻¹), while the trends in the eastern North Sea, along the Norwegian Coast and in the Barents Sea are lower (0.5–1.5 µatm yr⁻¹). In the North Sea this is consistent with a recent study of Omar et al. (2019), which is directly based on observations. These low trends will result in an increase in the strength of the ocean carbon sink with time.

Table 6 fCO₂ trend calculated from gridded, deseasonalized SOCAT v5 observations.

Generally, only few regressions over time ranges of less than 10 years turned out to be significant. This is an important finding when comparing the trends determined from our maps with the trends reported in literature, many of which were covering periods shorter than 10 years (Table 1). In order to compare the general patterns of fCO₂ trends determined from our maps with those directly determined from observations over a similar time range, we estimated the fCO₂ trends also from the SOCAT v5 observations that were used to produce the MLR (Table 6). We gridded and deseasonalized the SOCAT v5 data and divided the entire region into nine subregions. A figure showing the fits and the data coverage can be found in Appendix A. These observation-based trends show similar general patterns as those based on our maps (Fig. 9, 1998–2016): (1) largest trends in the southern North Sea, (2) decreasing towards the north with trends around the atmospheric trend in the northern North Sea and trends around 1 µatm yr⁻¹ in the Barents Sea, and (3) close to atmospheric trends in the Baltic Sea.

The observation that some subareas (the Baltic Sea or along the shore of the western North Sea) did not show a significant trend can be explained by the fact that coastal sea systems, especially enclosed areas such as the Baltic Sea, experience a high anthropogenic pressure. Anthropogenic impacts other than rising atmospheric CO₂ concentrations influence the ocean carbon system; for example the nutrient load of rivers can affect coastal ecosystems and primary production through eutrophication. This will result in lower fCO₂ in summer and higher fCO₂ in winter (Borges and Gypens, 2010; Cai et al., 2011). Another important process that influences the carbon system in the Baltic Sea is inflow events from the North Sea. In between such events, CO₂ accumulates in deeper water layers, causing an increasing gradient of dissolved inorganic carbon (DIC) across the halocline. Whenever deep winter mixing occurs, this will then lead to a large increase in surface fCO₂ because of the input of DIC-rich waters from below. Another reason is the observed change in alkalinity with time. This affects the fCO₂ though changes in the buffer capacity of the inorganic carbon system (Müller et al., 2016).

Figure 10 (a) The long-term trend (1998–2016) in surface ocean fCO₂ each month. (b) The average seasonality in fCO₂ for the periods 1998–2007 (green) and 2007–2016 (purple) in the northeastern North Sea (58–60^∘ N, 3–8^∘ E), normalized to December. The SD for each month is shown as the shaded area.

In most other regions, the sea surface fCO₂ trends were typically smaller than the rise in the atmospheric CO₂ concentration. A possible explanation is an earlier onset of the spring bloom. For example, in the North Sea a significant drawdown in fCO₂ has been observed as early as February in some years, but there is also a large variability (Omar et al., 2019). The bloom timing and onset in the North Sea after the 1990s have been shown to be mainly triggered by the spring–neap tidal cycle and the air temperature (Sharples et al., 2006). The bloom timing and onset were found to be significantly earlier in the 2010s compared to the previous decades (Desmit et al., 2020). Even if the trend in winter fCO₂ was following the atmospheric xCO₂ increase, such a change in bloom timing and onset would lead to a trend lower than in the atmosphere when averaging over the entire year. Figure 10a shows the annual trends in fCO₂ in each month in the four regions considered. Particularly in the North Sea and Baltic Sea, very low fCO₂ trends are observed in February–May, suggesting that changing timing of the spring bloom might be important here. Investigating the seasonal fCO₂ in more detail (Fig. 10b) revealed an earlier and deeper fCO₂ drawdown in the second decade of our analysis (2007–2016) than in the first (1998–2007) in the northeastern North Sea (58–60^∘ N, 3–8^∘ E). This strongly suggest that an earlier and stronger spring bloom is lowering the annual fCO₂ growth rates in this region, which is among the ones with the smallest fCO₂ trends in the North Sea (about 1 µatm yr⁻¹, Fig. 9). In the other regions, no such changes could be established with confidence. Future investigations should aim at generating fCO₂ maps with higher temporal resolution, as changes in the timing of the spring bloom might be a matter of days or weeks, which would not be fully resolved by the monthly maps presented here.

When looking at the interannual variability, it becomes obvious that the trend in the North Sea is slightly smaller than the atmospheric CO₂ trend. In contrast, the Norwegian Coast and the Barents Sea experience a robust trend much lower than the atmospheric trend (Norwegian Coast: 1–1.5 µatm yr⁻¹; Barents Sea: around 1 µatm yr⁻¹). Here we can also see a stable pattern of warming over timescales of 10 to 15 years. The warming in itself would result in an increase in fCO₂ with time, in addition to the atmospheric forcing. As we are observing a trend smaller than the atmospheric trend, temperature effects cannot be the driver here. The lower trend stems most likely from an earlier onset of spring bloom. It has been shown that the Atlantification and the reduced ice coverage of the Barents Sea lead to a longer productive season, and this will result in more months with strong undersaturation in CO₂ (Oziel et al., 2017). In the Baltic Sea the patterns are different. Here the variability is much larger, while most of the time periods show a trend larger than the atmospheric trend (3–3.5 µatm yr⁻¹). Although slightly smaller our results broadly agree with trend estimates based on measurements of 4.6–6.1 µatm yr⁻¹ over 2008–2015 (Schneider and Müller, 2018). Finally, it also needs to be noted that the uncertainty of the fCO₂ maps was highest in the Baltic Sea. This makes it also more difficult, if not impossible, to properly detect these small differences in the trends.

Figure 11 The trend in surface ocean pH estimated from deseasonalized pH. (a) The spatial distribution of the trend over the time period from 1998 to 2016 is shown. Grid boxes without a significant trend are denoted with a black dot. Panels (b)–(e) show the trends in different time periods in three regions, from the various years on the y axis to the various years on the x axis. Non-significant trends were left blank.

For pH, the trend in most regions is around −0.002 yr⁻¹ (Fig. 11). As expected, regions with the strongest trend in fCO₂ also show the highest trend in pH, such as the southern North Sea. The trend in the northern North Sea and along the Norwegian Coast is in good agreement with the pH trends found in studies focusing on the open Atlantic Ocean (−0.0022 yr⁻¹, Lauvset and Gruber, 2014) and the North Atlantic and Nordic Seas (−0.002 yr⁻¹, Lauvset et al., 2015).

4.3 CO₂ disequilibrium and flux

Figure 12 The average air–sea CO₂ disequilibrium over the period 1998–2016 (a, red colors indicate average undersaturation, while blue colors indicate average oversaturation). For every region average disequilibria are shown as seasonal averages (b) and time series of annual disequilibria (c). Blue line: North Sea; red line: Norwegian Coast; yellow line: Barents Sea; purple line: Baltic Sea

The average air–sea CO₂ disequilibrium (ΔfCO₂ = fCO_2,sea − fCO_2,atm) is shown in Fig. 12. The only region showing an average supersaturation is the southern North Sea. Towards the north, the surface ocean becomes more and more undersaturated, with the lowest values in the Barents Sea. The values in the Barents Sea (−60 to −80 µatm in the southern Barents Sea and less than −100 µatm around Svalbard) are in agreement with those estimated by Yasunaka et al. (2018). The seasonal cycle of ΔfCO₂ follows a biologically driven pattern with higher values in the winter and lower values from April to August. The seasonal cycle is largest in the Baltic and smallest in the Barents Sea.

Figure 13 The average air–sea CO₂ flux over the period 1998–2016 (a, red colors indicate sink regions, while blue colors indicate source regions). For every region average fluxes are show as seasonal averages (b) and time series of annual fluxes (c).

The air–sea CO₂ fluxes and their trends (Fig. 13) show that most regions are a net and increasing sink for CO₂. The only net source regions are the southern North Sea and the Baltic Sea. The two different regimes in the North Sea, with the southern, nonstratified part being a source and the northern temporarily stratified part a sink for CO₂, have been described in the literature before (Thomas et al., 2004), but the gradient between them as represented here may be a too sharp (Sect. 3.1). However, there is a large interannual variability in the fCO₂ disequilibrium (Omar et al., 2010), and studies based on different years find conflicting results regarding the direction of the flux (Kitidis et al., 2019; Schiettecatte et al., 2007; Thomas et al., 2004). This large interannual variability was also present in our maps. During some years, larger parts of the North Sea were a net source, while during other years also the southern North Sea acted as net sink (not shown).

The seasonal variations in the air–sea flux are driven by a combination of the changes in the disequilibrium, the wind strength and the ice cover. As there is less wind during summer, when the disequilibrium is large, but a smaller disequilibrium during winter, when the wind strength is high, the seasonal variability in the flux is often less clear than that in the disequilibrium. This can be seen in the Barents Sea and Norwegian Coast. Yasunaka et al. (2018) found the seasonal and interannual variation in the Barents Sea and the Norwegian Sea mostly corresponded to the wind speed and the sea ice concentration. We see the strongest dependence on the air–sea disequilibrium, however (not shown). This indicates that the seasonal and interannual variability in our fCO₂ maps is larger than in the maps generated by Yasunaka et al. (2018). Still, our average fluxes fit well with those reported by Yasunaka et al. (2018) of −8 to −12 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Barents Sea) and −4 to −8 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (Norwegian Coast). In the North Sea there is almost no net flux during winter, as the surface ocean is more or less in equilibrium with the atmosphere. In the Baltic Sea, there are high fluxes into the atmosphere during winter as here a large oversaturation coincides with strong winds. This is also why the Baltic Sea is a net source region. Although Parard et al. (2017) did find slightly smaller seasonal fluxes (+15 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ during winter and −8 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ during summer), the annual air–sea CO₂ fluxes are in good agreement (0 to +4 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ between 1998 and 2011).

The uncertainty in the calculated fluxes is a result of the uncertainties in the fCO₂ observations, ΔfCO₂ maps, the gas exchange parameterization and the wind product. The uncertainty of the ΔfCO₂ is mostly driven by the uncertainty of the MLR, resulting in an error between 12 and 39 µatm, according to the RMSE values of MLR1 for the different regions (Table 5). A number of studies address the uncertainty of gas exchange parameterizations and the wind products (Couldrey et al., 2016; Gregg et al., 2014; Ho and Wanninkhof, 2016). For this study, we apply an uncertainty of the gas transfer velocity of 20 % (Wanninkhof, 2014). This will result in an uncertainty of the air–sea flux of about 2 $\mathrm{mmol}\mathrm{C}{\mathrm{d}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{2}}$. It has to be kept in mind that the absolute uncertainty in k increases with increasing wind speed but that the uncertainty in the wind speed has the largest influence in summer when also the disequilibrium is large. In contrast, the uncertainty in ΔfCO₂ will cause larger errors in winter, when the wind speeds are high.

There is an ongoing discussion of how and to which extent the dominant climate mode in the North Atlantic, the North Atlantic Oscillation (NAO), is driving the variability in the CO₂ fluxes (Salt et al., 2013; Tjiputra et al., 2012; Watson et al., 2009). Even though some features in the time series seem to coincide with very extreme states of the NAO, such as a very large disequilibrium along the Norwegian Coast in 2010, we could not find any significant correlation between the CO₂ fluxes and the NAO index.

5 Conclusions

The MLR approach presented in this work is a relatively easy and straight forward method to produce monthly fCO₂ maps with a high spatial resolution in coastal seas, and the use of available open-ocean maps improved the coastal maps significantly. The maps reproduce nicely the main spatial and temporal patterns that are present in observations in the different regions for both fCO₂ and pH. The surface seawater fCO₂ trends were mostly lower than the atmospheric trends and also lower than the trends found in the open North Atlantic. Results show that the northern European shelf is an increasing net sink for CO₂. Only the Baltic Sea is a net source region. This method clearly has the potential to be extended to a larger region. However, it should be handled with care in regions with only a small number of observations as the MLR can lead to unrealistic values.

Long-term observations with a high temporal resolution are extremely important for developing maps such as presented here. While a decent spatial coverage exists for the open North Atlantic, most coastal regions are still undersampled, in particular relative to their high variability in time and space. To further understand and interpret the trends in fCO₂ and pH it is necessary to increase our knowledge and understanding of the interaction between primary production, respiration in the water column and the sediments, mixing and gas exchange, and their influence on the carbon cycle.

Besides an increased amount of in situ observations, also improved, higher-resolution driver data hold the potential to enable a better representation of spatial gradients. Including not only satellite-derived chlorophyll data but also colored dissolved organic matter (CDOM) can also lead to an improved performance of the regressions, especially in regions with a high load of terrestrial dissolved organic carbon.

While MLR-derived sea surface fCO₂ maps provide a coherent picture of the entire region, they have clear limitations and should be interpreted with caution in regions with few or no observations. A large number of observations is essential both for producing high quality maps and for their validation. Also, observations of a second parameter of the carbon system would be beneficial for deriving pH maps – to reduce and quantify the error introduced by estimating alkalinity from salinity. In addition to that, our work neglects the areas closest to land due to unavailability of CO₂ data and reanalysis products in those areas. For adding their contribution to the flux estimates, new platforms specialized on measurements directly at the land–ocean interface need to be developed.

Appendix A: Trend in surface ocean fCO₂ observations

Figure A1 Trend in surface ocean fCO₂ in deseasonalized, gridded observation data (SOCAT v5).