The northern European shelf as increasing net sink for CO 2

. We developed a simple method to reﬁne existing open ocean maps and extending them towards different coastal seas. Using a multi linear regression we produced monthly maps of surface ocean f CO 2 in the northern European coastal seas (North Sea, Baltic Sea, Norwegian Coast and in the Barents Sea) covering a time period from 1998 to 2016. A comparison with gridded 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 5 f CO 2 , pH and air-sea CO 2 ﬂux. The surface ocean f CO 2 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 f CO 2 increase by 2 µ atm yr − 1 , which is similar to the atmospheric trend. The Baltic Sea does not show a signiﬁcant trend. Here, the variability was much larger than the expected trends. Consistently, the pH trends were smaller than expected for an increase of f CO 2 in pace with the rise of atmospheric CO 2 levels. The calculated air-sea CO 2 ﬂuxes revealed that most regions were net sinks for CO 2 . Only 10 the southern North Sea and the Baltic Sea emitted CO 2 to the atmosphere. Especially in the northern regions the sink strength increased during the studied period.

. Overview of trends in surface ocean CO2 reported in the literature.

Reference
Time range dpCO2/ dt /µatm yr −1 North Sea Thomas et al. (2007Thomas et al. ( ) 2001Thomas et al. ( -2005 normalized to 16 • North Sea Salt et al. (2013Salt et al. ( ) 2001Salt et al. ( -2005 Yasunaka et al. (2018Yasunaka et al. ( ) 1997Yasunaka et al. ( -2013 Barents Sea Yasunaka et al. (2018Yasunaka et al. ( ) 1997Yasunaka et al. ( -2013 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 f CO 2 (CO 2 fugacity) maps based on already existing open ocean pCO 2 maps, in four example regions: North Sea, Baltic Sea, Norwegian Coast and the Barents Sea. A multi linear and the maps were then used to investigate trends in coastal f CO 2 and pH in the entire region from 1998 to 2016. Finally, we used the f CO 2 maps to determine the air-sea CO 2 exchange and its temporal and spatial patterns.

Study area
This work focuses on the northern European continental shelf and marginal seas. As we want to show the performance of the 5 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 (Figure 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 on 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 10 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.

Data handling
The CO 2 data used in this study were extracted from SOCAT version 5 (Bakker et al., 2016). Their coverage is shown in Figure   15 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 2 values from the closest 20 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 2 (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 2 values were extracted from two different products (Rödenbeck et al. (2014) (version oc_v1.5) and Landschützer et al. (2017Landschützer et al. ( , 2016   For preparing the input data for the MLR, observations closest to each SOCAT f CO 2 data point in time and space were extracted from the 3D fields with the driver data. This produces, for each of the driver data, a vector as long as the SOCAT f CO 2 observations. After this, the f CO 2 data as well as all extracted driver data were binned on a monthly 0.125 • x0.125 • 5 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 f CO 2 data. If a grid box for example only contains f CO 2 observations from the first week of the month and the northwestern corner, we make sure, that also the gridded driver data only contains 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 f CO 2 maps produced in this work. These 10 driver data were used for determining the MLRs.
For producing the final maps, a second set of the driver data was prepared, in the following called field data. Here the driver data were directly regridded to a monthly 0.125 • x0.125 • grid, providing full spatial and temporal coverage and a homogeneous average in each grid box. The field data were used to produce the f CO 2 maps based on the equation derived MLR equations.

15
The multi linear regression models were constructed by forward and backward stepwise regression using the driver data as predictor variables to model the f CO 2 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 1 ...X n being the driver data and a 1 ...a nn the regression coefficients. y = a 0 + a 1 · X 1 + ... + a n · X n + a 12 · X 1 X 2 + ... + a mn · X m X n + a 11 · X 2 1 + ... + a nn · X 2 n (1) The pCO 2 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 5 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., 2016;Lind et al., 2018).

10
The three linear fits were compared to each other in each region by taking into account the R 2 and the root mean square error (RMSE) of the fit, and 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 (I) would fit the observations: 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 5 were also compared to the gridded SOCAT v5 data, resulting in an average offset and standard deviation. In order to compare the maps against data that were not used to produce the maps, we predicted the f CO 2 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 f CO 2 observations in SOCAT v2019, gridded on a monthly 0.125 • x0.125 • grid. We also compare the maps directly with observations from repeated sampling locations in the North Sea and the Baltic Sea. 10

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 after 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 15 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), KSO − 4 dissociation constants after Dickson (1990) and the boron-salinity relation after 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 of AT over the last 25 20 years (Müller et al., 2016).

Calculation of trends
For calculating trends of f CO 2 and ocean acidification, the data in every grid box were deseasonalised by subtracting the long-term averages of the respective months. Then a linear fit was applied to the deseasonalised 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 25 resulted in significant trends, we decided to limit the trend analysis to starting years 1998 through 2006 and ending years 2008 through 2016.

Flux calculation
The air-sea disequilibrium was calculated as the difference between our mapped f CO 2 values and atmospheric f CO 2 in each grid cell and time step. The atmospheric f CO 2 was determined by converting the xCO 2 from the NOAA Marine Boundary 30 Layer Reference product from the NOAA GMD Carbon Cycle Group into f CO 2 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 2 flux (F) according to Equation 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) (Equation 4). The Schmidt number, Sc, was calculated according to Wanninkhof (2014) and the solubility coefficient for CO 2 , K 0 , after Weiss (1974).
In our calculations, we used 6-hourly winds of the NCEP-DOE Reanalysis 2 product. The coefficient a q in Equation 4 is strongly dependent on the used wind product (Roobaert et al., 2018). We determined it to be a q = 0.16 cm h −1 for the 6-10 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%. 15 3 Results

Maps of f CO 2
The skill assessment metrics for MLR 1, MLR 2 and MLR 3 are presented in Table 5. It shows the the R 2 and RMSE of the fit, the ME, as well as the average offset and standard deviation to the gridded SOCAT data.  Figure 3 shows, from left to right, the spatial distribution of the average difference between the predicted f CO 2 by MLR 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).
to slightly overestimate the f CO 2 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 f CO 2 well. However, in the Gulf of Finland it usually predicts too low f CO 2 values during May/June while it slightly underestimates events of very high f CO 2 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 Figure 3 (MLR 3 only)), 5 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 shows the monthly averages of f CO 2 produced by MLR 1 for February, May, August and November. In all regions, the highest f CO 2 values occur in the winter, while the lowest f CO 2 occur in summer. The largest seasonal cycle could be observed in the Baltic Sea, where f CO 2 reached well below 200 µatm in mid summer and over 500 µatm during the winter. 10 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 2 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 f CO 2 distribution in the North Sea, their specific shape here, are also a consequence of the influence of the Rödenbeck maps on our 15  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 2 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 5 also visible in subsequent pH and trend maps.

Maps of pH
The monthly average of pH calculated from MLR 1 f CO 2 ranges from about 8 during winter to 8.15 during summer in the North Sea and at the Norwegian coast ( Figure 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  , 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.

Performance of the pCO 2 maps
The performance of the MLR and the maps are evaluated in different ways: (1) the R 2 and the RMSE of the fit, (2) the average 5 deviation and its standard deviation, as well as the ME between the produced f CO 2 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 f CO 2 maps to observations from two time series stations. (2) -(4) will be shown for the time period covered by the driver data (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) and for the prediction of the f CO 2 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 10 dataset. Please note that the comparability of the model performance between the different regions is limited. All used statistical parameters 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) 5 mapping coastal f CO 2 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 influence of these spatial correlations is relatively small in regions with a high data density (as the european :::::::: European shelf) 20 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 Figure 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 f CO 2 in winter and early spring, while it overestimates the f CO 2 in summer. In all other regions, 25 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.
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 (Figure 7). In the other regions, the agreement is somewhat reduced 30 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 f CO 2 values (see Figure 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.  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 f CO 2 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 10 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 Figure 2). The MLR method can also lead to unrealistic All regions with questionable f CO 2 are also questionable in their pH data. There is a number of very high pH in the Barents Sea ( Figure 5), that are associated with also very low f CO 2 (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.

Trends in f CO 2 and pH
Trends in surface ocean f CO 2 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 5 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 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.
The monthly maps of f CO 2 from 1998 to 2016 enable us now to estimate the trend in surface ocean f CO 2 for the entire 10 region, equally distributed over the seasons (Figure 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 Figure 9. . It has to be noted that there was almost no measured f CO 2 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 −1 ), 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 −1 ). In the North Sea this is consistent with a recent study of Omar et al. (2019), 25 which is directly based on observations. These low trends will result in an increase in the strength of the ocean carbon sink with time.
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, of which many were covering periods shorter than 10 years (Table 1). In order to compare the general patterns of f CO 2 trends determined from our 30 maps with those directly determined from observations over a similar time range, we estimated the f CO 2 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   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 −1 in the Barents Sea, (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 as the Baltic Sea, experience a high anthropogenic pressure. Anthropogenic impacts other than rising atmospheric CO 2 concentrations influence the ocean carbon 5 system, for example the nutrient load of rivers can effect coastal ecosystems and primary production through eutrophication.
This will result in lower f CO 2 in summer and higher f CO 2 in winter (Borges and Gypens, 2010;Cai et al., 2011). Another important process that influences the carbon system in the Baltic Sea are inflow events from the North Sea. In between such events, CO 2 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 of surface f CO 2 because of the input 10 of DIC rich waters from below. Another reason is the observed change in alkalinity with time. This affects the f CO 2 though changes in the buffer capacity of the inorganic carbon system (Müller et al., 2016).
In most other regions, the sea surface f CO 2 trends were typically smaller than the rise in the atmospheric CO 2 concentration.
A possible explanation is an earlier onset of the spring bloom. For example, in the North Sea a significant drawdown in f CO 2 has been observed as early as February in some years, but there is also a large variability . The bloom timing 15 and onset in the North Sea after the 1990s has 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 was found to be significantly earlier in the 2010s compared to the previous decades (Desmit et al., 2020). Even if the trend in winter f CO 2 was following the atmospheric xCO 2 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 f CO 2 in each month in the four regions considered. Particularly in the North Sea and Baltic Sea, very low f CO 2 trends are observed in February -May, suggesting that changing timing of the spring bloom might be important here. Investigating the seasonal f CO 2 in more detail (Figure 10b), revealed an earlier and deeper f CO 2 drawdown in the second decade of our analysis (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) than in the first (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(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 f CO 2 growth rates in this region, which is among the ones with the smallest f CO 2 trends in the North Sea (about 1 µatm yr −1 , Fig. 9). In the other 5 regions, no such changes could be established with confidence. Future investigations should aim at generating f CO 2 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 2 trend. In contrast, the Norwegian coast and the Barents Sea experience a robust trend much lower than the 10 atmospheric trend (Norwegian Coast: 1 − 1.5 µatm yr −1 , Barents Sea: around 1 µatm yr −1 ). Here we can also see a stable pattern of warming over time scales of 10 to 15 years. The warming in itself would result in an increase of f CO 2 with time, in addition to the atmospheric forcing. As we are observing a trend smaller than the atmospheric trend, temperature effects can't 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 leads to a longer productive season, and this will result in more 15 months with strong undersaturation in CO 2 (Oziel et al., 2016). In the Baltic Sea the pattens 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 −1 ). Although slightly smaller our results broadly agree with trend estimates based on measurements of 4.6 -6.1 µatm yr −1 over 2008-2015  (Schneider and Müller, 2018). Finally, it also needs to be noted that the uncertainty of the f CO 2 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.
For pH, the trend in most regions is around -0.002 yr −1 (Figure 11) The air-sea CO 2 fluxes and their trends ( Figure 13) show that most regions are a net and increasing sink for CO 2 . 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 2 , have been described in the literature 10 before (Thomas et al., 2004), but the gradient between them as represented here, may be a too sharp (Section 3.1). However, there is a large interannual variability in the f CO 2 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). 5 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 10 dependence on the air-sea disequilibrium, however (not shown). This indicates that the seasonal and interannual variability in our f CO 2 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 mmol m −2 d −1 (Barents Sea) and -4 to -8 mmol m −2 d −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 15 winds. This is also why the Baltic Sea is a net source region. Although Parard et al. (2017) (Table 5). A number of studies addresses 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 mmol C d −1 m −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 25 has largest influence in summer when also the disequilibrium is large. In contrast, the uncertainty in ∆f CO 2 will cause larger errors in winter, when the wind speeds are high.
There is an ongoing discussion, 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 2 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 30 disequilibrium along the Norwegian Coast in 2010, we could not find any significant correlation between the CO 2 fluxes and the NAO index.

Conclusions
The MLR approach presented in this work is a relatively easy and straight forward method to produce monthly f CO 2 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 present in observations in the different regions for both f CO 2 and pH. The surface seawater f CO 2 trends were mostly lower than the atmospheric trends and also lower than the trends 5 found in the open North Atlantic. Results show that the northern European shelf is an increasing net sink for CO 2 . 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.
Longterm 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 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 insitu ::::: 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 CDOM, can also lead to an improved performance of the regressions, especially in regions with a high load of terrestrial dissolved organic 5 carbon.
While MLR derived sea surface f CO 2 maps provide a coherent picture of the entire region, they have clear limitations and should be interpreted with caution in regions with few or none observations. Both, for producing high quality maps, as well for their validation a large number of observations is essential.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 10 addition to that, our work neglects the areas closest to land due to unavailability of CO 2 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 f CO 2 observations Couldrey, M. P., Oliver, K. I. C., Yool, A., Halloran, P. R., and Achterberg, E. P.: On which timescales do gas transfer velocities control North Desmit, X., Nohe, A., Borges, A. V., Prins, T., Cauwer, K. D., Lagring, R., Zande, D. V. d., and Sabbe, K.: Changes in chlorophyll concentration and phenology in the North Sea in relation to de-eutrophication and sea surface warming, Limnology