A novel sea surface pCO₂-product for the global coastal ocean resolving trends over 1982–2020

In recent years, advancements in machine learning based interpolation methods have enabled the production of high-resolution maps of sea surface partial pressure of CO₂ (pCO₂) derived from observations extracted from databases such as the Surface Ocean CO₂ Atlas (SOCAT). These pCO₂-products now allow quantifying the oceanic air–sea CO₂ exchange based on observations. However, most of them do not yet explicitly include the coastal ocean. Instead, they simply extend the open ocean values onto the nearshore shallow waters, or their spatial resolution is simply so coarse that they do not accurately capture the highly heterogeneous spatiotemporal pCO₂ dynamics of coastal zones. Until today, only one global pCO₂-product has been specifically designed for the coastal ocean (Laruelle et al., 2017). This product, however, has shortcomings because it only provides a climatology covering a relatively short period (1998–2015), thus hindering its application to the evaluation of the interannual variability, decadal changes and the long-term trends of the coastal air–sea CO₂ exchange, a temporal evolution that is still poorly understood and highly debated. Here we aim at closing this knowledge gap and update the coastal product of Laruelle et al. (2017) to investigate the longest global monthly time series available for the coastal ocean from 1982 to 2020. The method remains based on a two-step Self-Organizing Maps and Feed-Forward Network method adapted for coastal regions, but we include additional environmental predictors and use a larger pool of training and validation data with ∼18 million direct observations extracted from the latest release of the SOCAT database. Our study reveals that the coastal ocean has been acting as an atmospheric CO₂ sink of −0.40 Pg C yr⁻¹ (−0.18 Pg C yr⁻¹ with a narrower coastal domain) on average since 1982, and the intensity of this sink has increased at a rate of 0.06 Pg C yr⁻¹ decade⁻¹ (0.02 Pg C yr⁻¹ decade⁻¹ with a narrower coastal domain) over time. Our results also show that the temporal changes in the air–sea pCO₂ gradient plays a significant role in the long-term evolution of the coastal CO₂ sink, along with wind speed and sea-ice coverage changes that can also play an important role in some regions, particularly at high latitudes. This new reconstructed coastal pCO₂-product (https://doi.org/10.25921/4sde-p068; Roobaert et al., 2023) allows us to establish regional carbon budgets requiring high-resolution coastal flux estimates and provides new constraints for closing the global carbon cycle.

1 Introduction

The exchange of carbon dioxide (CO₂) between the atmosphere and the ocean mainly depends on the gradient between the partial pressure of CO₂ (pCO₂) at the surface of the ocean and that of the overlying air on the global average. Over the past decade, the number of high-quality measurements of sea surface pCO₂ collected by research field programs and ships of opportunities has considerably increased. Moreover, large-scale community efforts have led to the compilation of tens of millions of sea surface pCO₂ measurements into uniform quality-controlled databases such as SOCAT (Surface Ocean CO₂ Atlas; Bakker et al., 2014), allowing for the quantification of the global oceanic CO₂ sink. However, in spite of this tremendous increase in data coverage, once gridded monthly at a typical spatial resolution of 1^∘ for the open ocean and 0.25^∘ for the coastal ocean, pCO₂ measurements remain largely discontinuous in time and space. The remaining regions and periods of time devoid of data thus prevent one from fully quantifying the air–sea CO₂ exchange and its full spatiotemporal variability based on measurements alone.

Therefore, in parallel to the ongoing measurement synthesis efforts, another research branch aiming at developing robust interpolation techniques to circumvent the spatial and temporal gaps in the data products has emerged. These techniques allow creating maps of pCO₂ that are continuous in space and time, typically at the monthly resolution (e.g., Chau et al., 2022; Gloege et al., 2022; Gregor and Gruber, 2021; Landschützer et al., 2014; Rödenbeck et al., 2014, 2015). The resulting observation-based continuous products (called hereafter “pCO₂-products”), however, differ in their spatial resolutions (e.g., from 0.25^∘ × 0.25^∘ in Chau et al., 2023, over 1^∘ × 1^∘ in Landschützer et al., 2014, to 4^∘ × 5^∘ in Majkut et al., 2014), their temporal coverage and their method of interpolation. Several studies have relied on direct interpolations of available pCO₂ measurements (e.g., Jones et al., 2015; Rödenbeck et al., 2014; Shutler et al., 2016) while others have first established linear (e.g., Iida et al., 2015; Park et al., 2010; Schuster et al., 2013) or nonlinear (e.g., Landschützer et al., 2014; Nakaoka et al., 2013; Zeng et al., 2014) predictive regression equations between a set of environment parameters (available everywhere and at every time within the domain of interest) and observed pCO₂ to perform the spatiotemporal extrapolation. These complementary pCO₂-products provide a better quantification of the spatial and temporal variability of the global oceanic CO₂ sink and its associated uncertainty on different time scales, going from seasonal fluctuations to decadal trends through interannual variability, while providing much improved observation-based benchmarks against which outputs from global model results can be evaluated (e.g., Hauck et al., 2020).

While significant efforts have been invested by the community to develop pCO₂-products for the global ocean, leading to a growing number of assessments of the CO₂ sink, most of these pCO₂-products ignore the coastal ocean (e.g., Landschützer et al., 2014) or resolve it by simply combining the coast with the open ocean (Chau et al., 2023). Indeed, the spatiotemporal investigations are performed for the entire ocean using the full set of observed pCO₂ data (coast and open ocean) in such a way that the specific conditions characteristic of coastal settings are not accurately accounted for in these products (e.g., Chau et al., 2022; Rödenbeck et al., 2013). In response to this shortcoming, other continuous pCO₂-products have been developed at the regional scale for several well monitored coastal seas (e.g., Bai et al., 2015; Hales et al., 2012; Jamet et al., 2007; Ono et al., 2004; Sarma et al., 2006) such as the California Current system (Sharp et al., 2022), European shelves (Becker et al., 2021) or the West Florida shelf (Chen et al., 2016). At the global scale, a significant step forward was made by Laruelle et al. (2017) when the first global coastal pCO₂-product at high spatial resolution (0.25^∘) was released for the entire coastal domain. This product, which is, to date, the only one available specifically developed for the global coastal ocean, is based on gridded coastal pCO₂ observations and nonlinear predictive regression equations between a set of environmental variables (drivers) and observed pCO₂ to perform the spatiotemporal extrapolation (the Self-Organizing Maps and Feed-Forward Network coastal pCO₂-product, ULB–SOM–FFN–coastalv1; Laruelle et al., 2017). This global coastal pCO₂-product provided a climatological mean (period 1998–2015) which allowed unprecedented investigation of the spatial distribution of the CO₂ sources and sinks in the global coastal ocean, especially for regions lacking data or regional assessments. It also allowed resolution of the seasonal variability of the air–sea CO₂ exchange in the coastal domain (Roobaert et al., 2019). Moreover, it was recently merged with an open ocean product to obtain a global reconstruction of the ocean CO₂ sink (Landschützer et al., 2020) and has been subsequently used to reduce the spread in global reconstructions (Fay et al., 2021). However, the ULB–SOM–FFN–coastalv1 pCO₂-product remains limited in its applications because it only provides a climatology covering a relatively short period (1998–2015) and is thus not suitable to evaluate the interannual variability or the long-term trends of the coastal air–sea CO₂ exchange. Such questions currently are at the forefront of the coastal research community's preoccupations (Bauer et al., 2013; Lacroix et al., 2021a; Laruelle et al., 2018; Regnier et al., 2013; Resplandy et al., 2023; Wang et al., 2017) but, because of the lack of adequate product, our confidence in the extent to which humans have perturbed the coastal air–sea CO₂ exchange since pre-industrial times remains low (Regnier et al., 2022). Moreover, the limitations of the ULB–SOM–FFN–coastalv1 do not yet allow us to produce robust trends in coastal pCO₂ fields against which global model outputs can be evaluated (e.g., Resplandy et al., 2023).

To address these limitations, this study expands and improves upon the version of the global coastal pCO₂-product of Laruelle et al. (2017) by extending its temporal coverage to four decades (1982–2020) and updating the methodology to resolve longer-terms changes in pCO₂, as described in the following section. The evaluation of this new product (ULB–SOM–FFN–coastalv2; Roobaert et al., 2023; https://www.ncei.noaa.gov/archive/accession/0279118; last access: November 2023) is done both spatially and for each decade individually, which represents an improvement compared with Laruelle et al. (2017) where the ULB–SOM–FFN–coastalv1 evaluation was limited to spatial and climatological seasonal cycles only. Using ULB–SOM–FFN–coastalv2 that relies on ∼18 million coastal direct observations from the SOCATv2022 database, we recalculate the coastal air–sea CO₂ exchange (FCO₂) for the 1982–2020 period and briefly describe the long-term trend of the global coastal CO₂ sink over this time frame. The long-term trend or “multidecadal trend” in this study is defined as a linear trend that spans a period exceeding 10 years that is in our case a trend that encompasses the years 1982–2020, resulting in a total of 39 years of observations. This study does not discuss the decadal change (period of 10 years) and interannual variabilities (year-to-year fluctuations) of the global coastal sink. In the future, these updated pCO₂- and FCO₂-products can be used as benchmarks for global oceanic models resolving trends in the coastal CO₂ dynamics, fulfilling a key knowledge gap identified in the latest Regional Carbon Cycle Assessment and Processes coastal synthesis (RECCAP2, Resplandy et al., 2023).

2 Methods

This section first describes the 2-step interpolation method used to generate the new version of the coastal pCO₂-product (Sect. 2.1) and the different datasets involved in this two steps procedure (Sect. 2.2). We then describe how the coastal air–sea CO₂ exchange is calculated (Sect. 2.3) and finally explain the approach used to quantify the uncertainties associated with our new pCO₂- and FCO₂-products (Sect. 2.4).

2.1 Self-Organizing Maps and Feed-Forward Network

We build upon the method described in Laruelle et al. (2017) to construct an updated observation-based continuous monthly pCO₂-product for the coastal ocean (ULB–SOM–FFN–coastalv2) at a 0.25^∘ spatial resolution over the 1982–2020 period. The method is based on the application of two artificial neural networks (the Self-Organizing Maps (SOM) network and the Feed-Forward Network (FFN)). The SOM first clusters the global coastal ocean into provinces characterized by similar environmental properties. In each province, the FFN then establishes nonlinear relationships between the observed pCO₂ and a set of environmental drivers of the coastal pCO₂ dynamics (which may be different from those used by the SOM). These relationships are then used to perform the spatiotemporal extrapolation of pCO₂ in each region defined by the SOM. This method was originally developed for the open ocean and is extensively described in Landschützer et al. (2013, 2014). It was later adapted for the global coastal ocean by Laruelle et al. (2017). We thus provide only a brief description of the methodology and focus here on the modifications introduced in this study.

In the first step, the global coastal ocean is divided into 10 biogeochemical provinces using the SOM clustering algorithm. Each resulting province is characterized by similar spatiotemporal patterns of a set of environmental variables, or drivers. In this study, we use the same drivers as in Laruelle et al. (2017), which consist of the wind speed calculated at 10 m above the sea surface (U₁₀), the sea surface temperature (SST), the sea surface salinity (SSS), the bathymetry, and the rate of change in the sea-ice coverage (see Sect. 2.2 for a description of the datasets). The SOM uses a neural network to detect similarities within multivariate datasets and uses an iterative procedure to distribute them into a predefined number of clusters. For each environmental driver, continuous monthly maps at the spatial resolution of 0.25^∘ are used as inputs for the neural network and each 0.25^∘ cell is allocated to one of the 10 provinces defined by the SOM. This procedure aims at minimizing the Euclidean distance between all points within each neuron of the network (see Landschützer et al., 2013, for more details). The spatial extension of these provinces varies from one month to the other because of the seasonal variations of the environmental drivers in such a way that a fixed grid cell in space may be assigned to several provinces over the course of a year. The choice of 10 provinces in the SOM stems from a sensitivity analysis that minimizes the average deviation between the observed pCO₂ and those simulated by the FFN algorithm (see second step below) while ensuring the presence of a minimum number of grid cells (>100) that can be used for the validation in each province (Laruelle et al., 2017). While their spatial extent varies seasonally, each province remains associated with specific regions over the course of the entire 1982–2020 period and the province occurring most often in each grid cell is shown in Fig. 1. Broadly, these provinces represent: province 1 (P1), the Antarctic shelf; P2 and P3, two subpolar/temperate coastal provinces of the Southern Hemisphere; P4 and P6, the large tropical coastal provinces; P5, a temperate province of the Northern Hemisphere which includes the Mediterranean Sea and the Norwegian Sea; P7, P8 and P10, high latitudes of the Northern Hemisphere provinces that are seasonally partly covered by sea ice (with the Baltic Sea and the Hudson Bay in P8); and P9, a permanent and cold polar province.

Figure 1 Coastal biogeochemical provinces generated by the Self-Organizing Maps (SOM) clustering algorithm. The spatial extension of these provinces can vary from one month to another due to seasonal variations of the environmental drivers used during SOM. Here we present their modal spatial distribution (see Sect. 2.1 for further details).

Table 1 Environmental drivers used for the Self-Organizing Maps (SOM) and the Feed-Forward Network (FFN) artificial neural networks to reconstruct the coastal pCO₂. These datasets are also used to calculate the coastal air–sea CO₂ exchange.

In a second step, within each biogeochemical province identified in step 1 (SOM), a FFN algorithm establishes nonlinear relationships between the observed sea surface pCO₂ and independent variables, or drivers, that are known to control its spatial and temporal variability. For each province, the FFN algorithm calculates relationships between the observed target variable (here pCO₂ using pCO₂ observations from the SOCAT_a dataset; see below) and inputs (environmental drivers; see below and Table 1) by adjusting weighting factors of a sigmoid activation function (one sigmoid function per neuron in the hidden layer) following an iterative procedure, i.e., a Levenberg–Marquardt backpropagation algorithm. At the first iteration, the weights of neurons are randomly assigned and the reconstructed pCO₂ is compared with the actual pCO₂ observations. Based on the resulting mismatch, the network weights are iteratively updated in a way that the error function – in our case the mean squared error between network output and actual observations – gets minimized. For each iteration, the FFN algorithm uses a fraction of the pCO₂ observations for the actual training of the network (i.e., the adjustment of the neuron weights), while another randomly selected fraction of the dataset is used to independently evaluate the performance of the algorithm. The final coefficients are obtained when the reconstructed pCO₂ simulated from the validation data does not significantly improve relative to the pCO₂ observations, to prevent overfitting. The final neuron weights and thus the resulting input–output relationships are used to reconstruct pCO₂ in each cell and for each month during the 1982–2020 period.

The predictors used for the FFN are U₁₀, SST, SSS, the atmospheric pCO₂, the rate of change in sea-ice coverage (except in regions not covered by sea ice, i.e., in P2, P3, P4 and P6), the bathymetry and the chlorophyll a concentration (Chl a). The Chl a is expressed as log₁₀(Chl a) to minimize the influence of its skewed distribution (Wrobel-Niedzwiecka et al., 2022). In P1, P8 and P9, we do not use Chl a as a driver because of the poor data coverage resulting from recurring cloud and/or sea-ice coverage in those provinces (see Sect. 2.2). This incomplete data coverage for Chl a is incidentally the reason why this predictor is not used at the SOM stage, because it requires complete global datasets. Atmospheric pCO₂, which was not included in the ULB–SOM–FFN–coastalv1, is also used as a driver of multidecadal changes induced by the increasing atmospheric pCO₂ concentration. Finally, we smooth spatially the monthly-resolved coastal pCO₂ field generated by the FFN using a moving 3 by 3 pixel window to remove abrupt pCO₂ transitions sometimes occurring at the boundaries between provinces. This smoothing procedure is described by Landschützer et al. (2014) and was also used in the ULB–SOM–FFN–coastalv1.

The surface pCO₂ data are extracted from the SOCATv2022 database (Bakker et al., 2022) that originally contains ∼40 million pCO₂ measurements for the entire global ocean (open and coastal seas combined). We randomly divide this dataset into two independent datasets: a group of data used for the FFN algorithm (SOCAT_a; see below) and a group of data that we use to validate our reconstructed pCO₂ (SOCAT_b). To do so, from the SOCATv2022 database, we follow the recommendation of the SOCAT community and use their accuracy criteria to only retain the data with the highest accuracy. To do so, we first select sea surface measurements expressed in fugacity of CO₂ (fCO₂) with a quality flag ranging from A to D (which corresponds to an estimated accuracy better than 5 µatm) and a World Ocean Circulation Experiment (WOCE) flag of 2 (good dataset following SOCAT) for the 1982–2020 period. Following Laruelle et al. (2017), we also remove fCO₂ values <30 and >1000 µatm that are likely derived from estuarine or fresh water systems that are not included in our coastal domain. We then randomly divide this dataset rich of ∼32 million fCO₂ measurements into a group of data used for the FFN algorithm (“a”, 80 % of the original dataset) and a group of data that we use to validate our reconstructed pCO₂ (“b”, 20 % of the original dataset). The two sets of data (SOCAT_a and SOCAT_b) are then gridded for each month at 0.25^∘ using the average of all fCO₂ values in each cell. Values are then converted from fCO₂ to pCO₂ using the equation of Takahashi et al. (2019, p. 7) and a coastal mask is applied on both gridded pCO₂-products. In this study, the coastal domain (“wide coastal ocean” with a total surface area of 76×10⁶ km², Laruelle et al., 2017) excludes the Black Sea, estuaries as well as inland water bodies, and its outer limit is defined as whichever point is furthest from the shoreline between the 1000 m isobath and a fixed 300 km distance (roughly the outer edge of territorial waters), following the coarse SOCAT definition of the coastal oceanic domain. At the end of this entire procedure, a total of ∼14 million and ∼4 million discrete coastal data have been allocated to SOCAT_a and SOCAT_b, respectively. A more common delineation of the coastal ocean is also used in this study when discussing the air–sea CO₂ exchange (Sect. 3.1) using the shelf break as the outer limit of the coastal domain (“narrow coastal ocean”, 28×10⁶ km²). The depth of the shelf break is calculated using a high-resolution global bathymetric database and estimated by calculating the slope of the sea floor. The isobath for which the increase in slope is the maximum over the 0–1000 m interval, yet still inferior to 2 %, defines the outer limit of the shelf break (Laruelle et al., 2013).

2.2 Environmental variables

The observational SST and SSS fields used as inputs for the SOM–FFN algorithm are calculated as the monthly means of the daily NOAA OI SST V2 (Reynolds et al., 2007) and of the daily Hadley center EN4 SSS (Good et al., 2013), respectively (Table 1). For U₁₀, we use the monthly mean of the 0.25^∘ resolution product of the European Center for Medium-Range Weather Forecasts (ECMWF) ERA5 wind product (Hersbach et al., 2023), which has a native temporal resolution of 6 h. The monthly mean of the daily 0.25^∘ dataset of Reynolds et al. (2007) is used for the sea-ice coverage. The rate of change in the sea-ice coverage for a given month x is then calculated as the difference between the sea-ice coverages of months x+1 and x−1. The atmospheric pCO₂ is from the SeaFlux product (Fay et al., 2021) which is calculated from the dry air mixing ratio of CO₂ (xCO₂) provided by the ESRL surface marine boundary layer CO₂ product (Dlugokencky et al., 2019; https://www.esrl.noaa.gov/gmd/ccgg/mbl/data.php; last access: October 2023) with a vapor correction according to Dickson et al. (2007), using the ERA5 sea level pressure (Hersbach et al., 2023) and applying the method of Landschützer et al. (2014) and Fay et al. (2021). It should be noted that, due to the proximity to the continent, the coastal ocean might be more exposed to anthropogenic sources of CO₂ and thus might be exposed to higher atmospheric pCO₂ compared with the global oceanic average. The use of spatially resolved dry air mixing ratio of CO₂ datasets, such as the one from the NASA's Orbiting Carbon Observatory 2 Goddard Earth Observing System (OCO-2 GEOS; Eldering et al., 2017), instead of the product used in this study might be more appropriate to include this effect. However, OCO-2 GEOS only covers 2015–2022, which is too short for the purpose of our study. It is also expected that the choice of the atmospheric pCO₂ does not considerably influence our FCO₂ calculations as the air–sea pCO₂ difference is mainly controlled by the oceanic pCO₂ (see Sect. 2.4). We use the bathymetry from the 2 min global ETOPO2 database (US Department of Commerce, 2006) and the Chl a field derived from the monthly 4 km merged GlobColour product from the EU Copernicus Marine Service information (Copernicus Marine Service, 2023; https://doi.org/10.48670/moi-00281, last access: October 2023), which is the product with the longest Chl a temporal coverage (1998–2020). However, because of recurrent cloud coverage everywhere and sea-ice coverage at high latitudes, the monthly averaged Chl a field at a 0.25^∘ resolution is discontinuous with grid cells devoid of data. We fill these cells (9 % excluding high latitude coastal regions) using a cascade of interpolation methods, the order of which depends on data availability in time and space in the surrounding cells: for an empty cell of the month x, the interpolation is performed by computing in the order of rank (1) the mean of the next neighboring cells of the month x, (2) the mean value of the month x+1 and month x−1 of the same cell, (3) the monthly mean value x of the cell for the entire 1998–2020 period and (4) the annual mean value of the cell for the entire 1998–2020 period. At high latitudes, where none of these options are feasible because of large bands without any data, we assign the modal value of Chl a as the default value in order to ensure the continuity of the maps that is required for the FFN algorithm. From this continuous gridded Chl a product over the recorded period, we then calculate for each cell a monthly seasonal climatology and attribute the climatological values to the unrecorded period (1982–1997). This means that, in our calculations, Chl a does not contribute to long-term changes in pCO₂ before 1997. All observational fields are converted from their original spatiotemporal resolution to monthly 0.25^∘ gridded resolution for the 1982–2020 period (except for the bathymetry which is constant over time) to match the observational pCO₂-product (SOCAT_a) resolution (Table 1).

2.3 Air–sea CO₂ exchange

The pCO₂ field generated by the SOM–FFN algorithm (ULB–SOM–FFN–coastalv2) is used to calculate the air–sea CO₂ exchange for each grid cell at the monthly timescale over the 1982–2020 period following Eq. (1):

$$\begin{array}{}\text{(1)}& {\text{FCO}}_{\mathrm{2}}=k{K}_{\mathrm{0}}\mathrm{\Delta }{\mathrm{pCO}}_{\mathrm{2}}(\mathrm{1}-\text{ice}),\end{array}$$

where FCO₂ (mol C m⁻² yr⁻¹) represents the coastal air–sea CO₂ exchange, ΔpCO₂ (atm) is the gradient between the oceanic pCO₂ and the atmospheric pCO₂ and K₀ (mol m⁻³ atm⁻¹) is the CO₂ solubility in seawater which is a function of SST and SSS following the equation of Weiss (1974). k (m yr⁻¹) represents the gas exchange transfer velocity which is a function of the second moment of the wind speed at 10 m above sea level and is calculated using the equation of Ho et al. (2011) and using the Schmidt number based on the equation of Wanninkhof (2014). The sea-ice coverage for each grid cell is represented by the term ice and ranges from 0 (no ice cover) to 1 (100 % ice cover). By convention a positive FCO₂ value corresponds to a source of CO₂ for the atmosphere. We use the same U₁₀, SST, SSS, sea-ice coverage and atmospheric pCO₂ datasets (see Sect. 2.2) to calculate FCO₂ and perform the FFN pCO₂-reconstruction. Our reconstructed coastal FCO₂ is also compared with coastal FCO₂ estimates derived from a synthesis of 214 regional FCO₂ estimations (Dai et al., 2022) and from the FCO₂-product derived from the original ULB–SOM–FFN–coastalv1 pCO₂-product (Roobaert et al., 2019).

2.4 Uncertainties in the reconstructed coastal data products

The uncertainties associated with our reconstructed pCO₂ and FCO₂ coastal products are estimated using the method proposed by Landschützer et al. (2014, 2018) and used by subsequent authors (e.g., Roobaert et al., 2019; Sharp et al., 2022). The FCO₂ uncertainty results from four sources of uncertainties, which are considered independent and thus summed quadratically:

$$\begin{array}{}\text{(2)}& {\mathit{\sigma }}_{{\mathrm{FCO}}_{\mathrm{2}}}=\sqrt{{\mathit{\sigma }}_{\mathrm{\Delta }{\mathrm{pCO}}_{\mathrm{2}}}^{\mathrm{2}}+{\mathit{\sigma }}_k^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{wind}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{ice}}^{\mathrm{2}}},\end{array}$$

where ${\mathit{\sigma }}_{{\mathrm{FCO}}_{\mathrm{2}}}$ represents the total FCO₂ uncertainty (Pg C yr⁻¹), ${\mathit{\sigma }}_{\mathrm{\Delta }{\mathrm{pCO}}_{\mathrm{2}}}$ is the uncertainty of the air–sea pCO₂ gradient, σ_k is the uncertainty associated with the choice of the k formulation in Eq. (1) and σ_wind is the uncertainty associated with the choice of the wind speed product (see Roobaert et al., 2018). We also include the effect of the choice of sea-ice product on the FCO₂ uncertainty (σ_ice) which was not included in the original calculations of Landschützer et al. (2014, 2018) but has been identified as a potential source of uncertainty in global coastal reconstructions (e.g., Resplandy et al., 2023). All of these four sources of uncertainty are expressed in petagrams of carbon per year. We do not include the uncertainty associated with the solubility term (K₀) in our uncertainty assessment since this contribution is minimal (0.2 %) as suggested by Weiss (1974). σ_wind is calculated following the strategy described in Roobaert et al. (2018) which consists of using the standard deviation of global FCO₂ fields calculated with three different wind products: the ERA5 (Hersbach et al., 2023), the Cross-Calibrated Multi-Platform Ocean Wind Vector 3.0 (Atlas et al., 2011) and the NCEP/NCAR reanalysis 1 (Kalnay et al., 1996). Since these wind products cover different time periods, σ_wind is calculated for the overlap period (1991–2011) between products. σ_k is estimated as the standard deviation of global FCO₂ fields calculated with four different global scale k parametrizations with the same wind speed (ERA5). We use the formulations of Ho et al. (2011), Sweeney et al. (2007), Takahashi et al. (2009) and Wanninkhof (2014), all suited for global scale applications (e.g., see Roobaert et al., 2018). σ_ice is calculated as the standard deviation of global FCO₂ fields calculated with two different sea-ice products: the NOAA dataset of Reynolds et al. (2007) and the sea-ice dataset of Rayner et al. (2003). ${\mathit{\sigma }}_{\mathrm{\Delta }{\mathrm{pCO}}_{\mathrm{2}}}$ mainly results from the oceanic pCO₂ uncertainty since the atmospheric pCO₂ uncertainty is significantly lower (Landschützer et al., 2018). For instance, Roobaert et al. (2019) quantified that uncertainties in atmospheric pCO₂ only contribute to 6 % in the overall FCO₂ uncertainty.

Following Landschützer et al. (2014), the uncertainty over the oceanic pCO₂ can be obtained from the quadratic sum of three sources of uncertainties:

$$\begin{array}{}\text{(3)}& {\mathit{\theta }}_{{\mathrm{pCO}}_{\mathrm{2}}}=\sqrt{{\left({\displaystyle \frac{{\mathit{\theta }}_{\mathrm{obs}}}{\sqrt{N}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\theta }}_{\mathrm{grid}}}{\sqrt{N}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\theta }}_{\mathrm{map}}}{\sqrt{N_{\mathrm{eff}}}}}\right)}^{\mathrm{2}}},\end{array}$$

where ${\mathit{\theta }}_{{\mathrm{pCO}}_{\mathrm{2}}}$ represents the total uncertainty of the oceanic pCO₂ (µatm), θ_obs is the experimental uncertainty associated with the sampling in the field of the observations from the SOCAT database (µatm), θ_grid is the uncertainty associated with the gridding of the observations from SOCAT into 0.25^∘ monthly meshed maps (µatm) and θ_map is the uncertainty derived from the comparison between the reconstructed pCO₂ and the observed gridded pCO₂ from the SOCAT database (µatm). Following Sharp et al. (2022), we use accuracies that are attributed to each fCO₂ measurement by the SOCAT community (flags A–D) to calculate θ_obs. Flags “A” and “B” represent an estimated accuracy of 2 µatm while the accuracy of flags “C” and “D” is 5 µatm. We first calculate the mean of all fCO₂ flags in each grid cell for each month. We then calculate the global average gridded flags uncertainty of all cells for the year x (or the entire period). For θ_grid, we first calculate in each grid cell of the month x the standard deviations of all fCO₂ values from the SOCAT database used for the gridding. We then calculate the average of these standard deviations of all grid cells for the year x (or the entire period). θ_map is calculated as the root mean squared deviation between the reconstructed pCO₂ and the gridded pCO₂ observation from the training dataset (SOCAT_a). We also divide each source of uncertainty (i.e., θ_obs, θ_grid and θ_map) by the square root of the number of pCO₂ samples (N; for details see Landschützer et al., 2018; Roobaert et al., 2019). For θ_map, the value of N is corrected to account for the fact that all individual errors are not spatially independent. To this end, we calculate the effective sample size (N_eff; see Landschutzer et al., 2018) by randomly selecting 1000 samples (40 % of the samples if the total number of samples is <1000) that cover our study period and calculating a lag 1 autocorrelation coefficient following Eqs. (18) and (19) of Landschützer et al. (2018). As we only use a subset of 1000 samples, we perform Monte Carlo simulations where this procedure is repeated 10 times and our final N_eff is calculated as the median of the 10 iterations. Finally, the total uncertainty on FCO₂ associated with the reconstructed pCO₂ (${\mathit{\sigma }}_{\mathrm{\Delta }{\mathrm{pCO}}_{\mathrm{2}}}$, Pg C yr⁻¹) is obtained by applying ${\mathit{\theta }}_{{\mathrm{pCO}}_{\mathrm{2}}}$ (µatm) in Eq. (1). All these procedures are performed globally for each year and for the entire period of our study.

Figure 2 Global maps of the climatological (1982–2020) averaged (a) reconstructed coastal pCO₂-product and (b) gridded pCO₂ from the SOCATv2022 database and used for the FFN algorithm (SOCAT_a). Panel (c) shows the temporal mean of the residuals between the coastal pCO₂-product and SOCAT_a and (d) shows standard deviation. Values in all panels are expressed in µatm.

Figure 3 Histograms of the pCO₂ residuals (difference between the reconstructed coastal pCO₂-product with (a) SOCAT_a and (b) SOCAT_b) for four decades expressed in µatm.

3 Results and discussion

3.1 pCO₂-product evaluation

Globally, our reconstructed coastal pCO₂-product compares well with the observed pCO₂ used to train the FFN algorithm (SOCAT_a) and reproduces all the well-known global spatial pCO₂ patterns with generally low pCO₂ (<360 µatm) at temperate as well as high latitudes and high pCO₂ (> 360 µatm) at low latitudes (contrast Fig. 2a with Fig. 2b). The spatial distribution of the temporal mean residuals (i.e., difference between the coastal pCO₂-product and SOCAT_a in each grid cell for every month where observations are available) reveals that, in some regions, underestimations (negative mean residual; blue colors in Fig. 2c) or overestimations (red colors) of the pCO₂ can be generated by the coastal pCO₂-product. However, most of the calculated residuals fall within the −20 to 20 µatm range, accounting for 69 % of the grid cells, while 45 % of the grid cells have absolute residuals <10 µatm (Figs. 2c and 3a).

A global mean of the residuals (bias) value of 0 µatm and a coefficient of determination (r²) of 0.7 are calculated, as expected, since the algorithm minimizes the root mean square error (RMSE) between the reconstructed pCO₂ and target pCO₂ observations. The global RMSE is, however, substantially larger (29 µatm) yet still comparable to those calculated, at the regional scale, in previous coastal pCO₂ studies based on statistical interpolations (e.g., see Chen et al., 2016) and slightly lower than the 32 µatm global RMSE calculated by Laruelle et al. (2017). Large differences can be observed between our product and SOCAT_a locally in regions that are known to present large spatiotemporal variabilities in pCO₂ and/or in regions lacking data to train our FFN algorithm. For instance, residuals and/or standard deviations >20 µatm are encountered in the Baltic Sea (which is henceforth treated as an independent biogeochemical province for some calculations), in upwelling regions (e.g., along the Peruvian upwelling), in coastal seas under the influence of seasonal changes in sea-ice coverage (e.g., along the Antarctic shelf) as well as along the very nearshore coastal domain (e.g., along the California current coast; Fig. 2c and d).

The overall consistency between our coastal pCO₂-product and SOCAT_a is diagnosed over the entire time span of our study, as illustrated by the histograms of residuals calculated for each of the four decades of our calculation period (Fig. 3a) or when the calculations are performed for each individual year (Table S1 in the Supplement). This is an important test as a study by Gloege et al. (2022) suggests that decadal trends in pCO₂-products may be obscured by changing residual distributions over time. In spite of the highly heterogeneous distribution of the number of pCO₂ observations available through time (<500 gridded cells in the 1980s vs. >10 000 in the 2010s), the shape and spread of the four histograms of the residuals are closely similar between decades with a distribution centered on a global mean bias close to 0 µatm and most of the residuals falling in the −20 and 20 µatm range. This demonstrates the accuracy of the method over time despite the skewed distribution of the calibration data. The analyses performed for each individual year reveal that global biases do not exceed 5 µatm. Exceptions are observed in the 1980s where biases (e.g., absolute bias of 6 µatm in 1987) and RMSE (e.g., 42 µatm in 1989) can be larger and partly attributed to an exceptionally low pCO₂ observational coverage during these periods (e.g., see Bakker et al., 2014). In the first version of the coastal product (ULB–SOM–FFN–coastalv1), the evaluation of the pCO₂-product of Laruelle et al. (2017) was restricted to spatial and climatological seasonal cycles. In this study, we successfully extended this analysis to the entire time period and evaluated each year and decade individually.

Table 2 Statistical analyses (bias, RMSE and r²) of the reconstructed coastal pCO₂-product compared with pCO₂ observations from SOCAT_a and SOCAT_b for the different biogeochemical provinces. For each province, bias, RMSE and r² are calculated using all of the monthly cells of the province for the period 1982–2020.

* Numbers in parentheses for P8 represent statistics when the Baltic Sea is removed from P8 and defined as an isolated province.

At the regional scale, 7 out of the 10 biogeochemical provinces yield RMSEs against SOCAT_a close to 20 µatm or lower with the best fit in P2 (RMSE = 13 µatm; Table 2). This is a significant improvement over the ULB–SOM–FFN–coastalv1 product which only had three provinces with RMSE <20 µatm (Laruelle et al., 2017) and none <15 µatm. We attribute this improvement to our advanced setup of the method such as the inclusion of the atmospheric pCO₂ as a driver as well as an increased number of available observations to train our FFN algorithms. In three provinces (P1, P5 and P7), however, the RMSE exceeds 35 µatm. Such values can partly be explained by the complex dynamics of the sea ice in the Antarctic shelf (P1) and by the limited number of observational data combined with the inclusion of coastal regions that present large spatiotemporal variabilities and cover two disconnected temperate basins (P5 and P7) of the Northern Hemisphere. This discrepancy was also highlighted by Landschützer et al. (2020). High discrepancies are further observed in the Baltic Sea (Fig. 2c and d) which is analyzed as a separate province (RMSE value of 41 µatm; Table 2). Excluding the Baltic Sea from P8 considerably reduces RMSE from 41 to 23 µatm. The inclusion of the Baltic Sea in P8 can also explain the large RMSE calculated by Laruelle et al. (2017) for the corresponding province of the ULB–SOM–FFN–coastalv1 coastal pCO₂-product (∼47 µatm).

3.2 Validation against independent data

Our reconstructed coastal pCO₂-product is also validated against an independent dataset that is derived from pCO₂ observations from the SOCATv2022 that were not used for the training of the FFN algorithm (see Sect. 2.1). This dataset consists of a pool of 404 206 gridded cells that are uniformly distributed between both hemispheres (SOCAT_b; Fig. S1 in the Supplement), an essential criterion for training the network. Globally, a good match is observed between the coastal pCO₂-product and SOCAT_b with a global bias and an RMSE of 0 and 29 µatm, respectively. These values are similar to those derived from the statistical analysis performed against SOCAT_a. At the biogeochemical provinces scale, RMSEs generally do not exceed 23 µatm (maximum value in P6; Table 2), except where important RMSEs (34 µatm for P1, P5 and P7) had already been calculated during the comparison with SOCAT_a (i.e., in regions under the sea-ice coverage dynamics and poor data coverage and provinces which encompass regions with high spatiotemporal pCO₂ dynamics).

As with SOCAT_a, the analysis against SOCAT_b demonstrates a good performance of our reconstructed coastal pCO₂-product over time with the histograms of the residuals calculated for each of the four decades presenting the same shape and spread in spite of the marked decrease in grid cell numbers over time (1054 grid cells in 1980 vs. 248 626 grid cells in 2010; Fig. 3b). Each of these histograms shows a distribution centered on a value of 0 µatm with ∼50 % of the grid cell residuals falling between −10 and 10 µatm. This is also true at the scale of the biogeochemical provinces with the four histograms of the residuals revealing global mean biases of 0 µatm and ∼50 % of the residual falling in the −10 and 10 µatm range. Only P8 stands as an exception (where 50 % of the residuals are between −40 and 40 µatm; Fig. S2), mainly due to the presence of the Baltic Sea in this province.

Figure 4 Panel (a) shows the temporal coverage (in years) where pCO₂ measurements extend over x years in SOCAT_b. The location of the eight coastal sites for which we present pCO₂ times series (black boxes) is also shown. Panels (b–i) show pCO₂ times series (in µatm) from the reconstructed pCO₂-product (in black), from SOCAT_b (in red), and from the atmospheric pCO₂ (in blue). For each region, we only select grid cells that extend over 30 years of observations in SOCAT_b. Medians are represented by circles and the vertical bars represent the monthly pCO₂ intra-spatial variability in the region. For each region, we report the bias (µatm), RMSE (µatm) and number of cells for the calculation between the reconstructed pCO₂-product and SOCAT_b as well as their respective long-term pCO₂ trend (in µatm decade⁻¹, which is calculated first as the slope of a linear trend using the monthly median values of all the deseasonalized data). The Cascadian shelf has no value for the SOCAT_b trend since no significant trend is detected based on a Mann–Kendall statistical test.

We also present pCO₂ time series derived from our reconstructed coastal pCO₂-product and compare them to data extracted from SOCAT_b for eight coastal sites (Fig. 4). The choice of these coastal regions is motivated by their data coverage extending over 30 years (Fig. 4a) and the fact that it is possible in these grid cells to reconstruct a spatially complete seasonal climatological cycle (i.e., data are available for all 12 months). For each region, we only extracted cells for which observations extended >30 years and reconstructed their times series from the coastal pCO₂-product and SOCAT_b, respectively (Fig. 4b–i). For most of the regions, the pCO₂-product properly captures the temporal dynamics of pCO₂ derived from the observations, bearing slight underestimations or overestimations in specific areas such as along the Cascadian shelf and the east coast of Australia (Fig. 4d and g). For the eight coastal sites, absolute biases are all <10 µatm with a minimum absolute bias of 1 µatm in the Irminger Sea and a maximum absolute bias of 9 µatm along the New Zealand coast. Except along the Cascadian shelf, all coastal sites present RMSEs lower than ∼30 µatm (five of the eight regions show RMSEs $<\sim \mathrm{20}$ µatm) which falls in the range of our global and regional RMSE values. The largest RMSE is calculated along the Cascadian shelf (62 µatm) and can partly be explained by the large spatial pCO₂ variability in SOCAT_b (as shown by the vertical bars in Fig. 4d) because of the riverine influence in the region.

Finally, the reconstructed pCO₂ times series are compared with three buoys with measurements longer than 10 years located in Cape Elizabeth (NDBC Buoy 46041), in Gray's Reef (NDBC Buoy 41008) and in the Gulf of Maine (Coastal Western Gulf of Maine Mooring; Sutton et al., 2019). The three pCO₂ time series for each buoy location mentioned above are presented in Fig. S3. Although smaller amplitude variabilities are generally observed, results show that the reconstructed pCO₂ times series follow those of the observational data with values that are mainly between the buoys errors (Fig. S3a–c). We speculate that the smaller amplitude stems from the coarser 0.25^∘ grid resolution of our method compared with the point nature of the buoy data. Landschützer et al. (2016) drew a similar conclusion when they compared their open ocean pCO₂ data with open ocean at HOT and BATS and is further corroborated by the much smaller variabilities obtained when raw SOCAT data are averaged at the grid cell level (Fig. 4). The exception is the Gulf of Maine where a general underestimation of pCO₂ is observed compared with the buoy observations. The reconstructed pCO₂-product also reproduces the observed climatological seasonal cycles (including a relatively good timing of the seasonal maxima and minima) for the three buoys as shown in Fig. S3d–f. Absolute average bias values of 14, 4 and 45 µatm, and RMSE values of 50, 53 and 61 µatm, are calculated between the pCO₂-product and the observations for Cape Elizabeth, Gray's Reef and in the Gulf of Maine, respectively. These statistical error values are larger than those calculated when the comparison with SOCAT is performed, whether on a global or regional scale. This is not surprising since the reconstructed pCO₂-product is for a global application and is quite challenging to compare with specific coastal buoys that present high temporal and spatial variability such as shown by the violin in Fig. S3.

The good evaluation of the reconstructed pCO₂ compared with SOCAT_a, SOCAT_b and buoys data gives us confidence to identify the linear trends of both the pCO₂ and FCO₂ over different temporal scales. For example, our results show that for all the eight studied regions represented in Fig. 4, an increase in pCO₂ over time comprised between 12 and 20 µatm decade⁻¹ is calculated for the long-term trend with our pCO₂-product, a range in good agreement with the 12–18 µatm decade⁻¹ obtained with SOCAT_b. Although New Zealand shows the largest bias between SOCAT_b and the pCO₂-product, they both show that this region displays the fastest trend in terms of pCO₂ rise (18 and 20 µatm decade⁻¹ for SOCAT_ b and the pCO₂-product, respectively).

3.3 Air–sea CO₂ exchange

This section describes the coastal air–sea CO₂ exchange patterns that are calculated using our new reconstructed pCO₂-product over the 1982–2020 period. The spatial and seasonal FCO₂ patterns are only briefly discussed (Sect. 3.3.1) since those have been extensively discussed in previous studies (e.g., see Dai et al., 2022; Resplandy et al., 2023; Roobaert et al., 2019). We thus focus on the long-term FCO₂ trends (Sect. 3.3.2), which are still poorly understood and highly debated (e.g., Lacroix et al., 2021a; Laruelle et al., 2018; Resplandy et al., 2023).

Figure 5 Panel (a) shows the spatial distribution of the annual average air–sea CO₂ exchange (FCO₂, in mol C m⁻² yr⁻¹) and (b) seasonal FCO₂ variability (expressed as the root mean square (RMS, in mol C m⁻² yr⁻¹)) calculated with the reconstructed coastal pCO₂-product (ULB–SOM–FFN–coastalv2, 1982–2020 climatology). The latitudinal mean FCO₂ distribution (red line) and its associated longitudinal variability (red shading) is presented in panel (c). This latter is compared with the FCO₂ calculated with the ULB–SOM–FFN–coastalv1 pCO₂-product (in green; Roobaert et al., 2019) and against a synthesis of 214 regional FCO₂ estimations (blue dots; Dai et al., 2022). For consistency in the comparison in (c) we applied the same coastal delimitation as in Dai et al. (2022) and Roobaert et al. (2019) to the FCO₂ ULB–SOM–FFN–coastalv2 product, i.e., we used the shelf break as the outer limit of the coastal domain (narrow coastal ocean). Panel (c) is also reconstructed based on an overlap period between the three products (1998–2020; except FCO₂ ULB–SOM–FFN–coastalv1 which is limited to the 1998–2015 period).

3.3.1 Spatial and seasonal variations

The spatial distribution of the climatological mean coastal FCO₂ shows that coastal regions in temperate areas (between 40 and 60^∘ in both hemispheres) and at high latitude (beyond 60^∘ in both hemisphere) mainly act as CO₂ sinks while CO₂ sources are mainly encountered in the subtropical band (Fig. 5a) which is consistent with the global latitudinal pattern established by previous studies (e.g., Borges, 2005; Borges et al., 2005; Cai, 2011; Cao et al., 2020; Chen et al., 2013; Dai et al., 2022; Laruelle et al., 2010, 2014; Roobaert et al., 2019). Globally, with the coastal delineation used in this study (“wide coastal ocean”; 76×10⁶ km²), the coastal ocean absorbs on average 0.40 Pg C per year (with an uncertainty of ±0.03 Pg C yr⁻¹; see Sect. 3.3.3) over the 1982–2020 period. Using the shelf break as the outer limit of the coastal domain (“narrow coastal ocean”; 28×10⁶ km²) which is a more common delineation of the coastal ocean, the globally integrated coastal sink amounts to $-\mathrm{0.18}\pm \mathrm{0.01}$ Pg C yr⁻¹ which is consistent with the latest estimates (e.g., −0.2 Pg C yr⁻¹ in Roobaert et al., 2019, and −0.25 Pg C yr⁻¹ in Dai et al., 2022). It should be noted that these comparisons are not straightforward because they do not cover the same time periods (i.e., the 1998–2015 period in Roobaert et al. (2019), 1998 to the present in Dai et al. (2022) and 1982–2020 in this study) and older assessments often do not report an explicit calculation period (Regnier et al., 2022).

Most of the intense CO₂ sinks (absolute FCO₂ value >0.5 mol C m⁻² yr⁻¹) are encountered at high latitudes of the Northern Hemisphere and in the temperate regions of the Southern Hemisphere, while CO₂ sources in the tropical bands are moderate except along upwelling areas such as in the Arabian Sea (Fig. 5a and c). A large fraction (44 % and 53 % for the wide and narrow domain, respectively) of the global CO₂ uptake is taking place north of 60^∘ N, which was already suggested in Laruelle et al. (2010) and further confirmed in subsequent studies (e.g., Cai, 2011; Dai et al., 2022; Laruelle et al., 2014; Roobaert et al., 2019). The spatial distribution of coastal CO₂ sources and sinks also closely follows the latitudinal FCO₂ profile calculated by Roobaert et al. (2019) which is based on ULB–SOM–FFN–coastalv1 (red and green lines in Fig. 5c). These global pCO₂-products, however, predict less variability in flux density than a compilation of regional estimations as shown in Fig. 5c when comparing our climatological FCO₂ latitudinal profile with the synthesis of 214 regional FCO₂ estimates which was already pointed by Dai et al. (2022) when comparing their data synthesis with the latitudinal FCO₂ profile of Roobaert et al. (2019), suggesting strong FCO₂ heterogeneities for a same latitudinal band. Finally, the seasonal coastal FCO₂ variability (expressed as the root mean square (RMS) of the seasonal amplitude) agrees with the few studies performed at global scale (e.g., see Dai et al., 2022; Roobaert et al., 2019) with high seasonal FCO₂ amplitudes (RMS values >1.5 mol C m⁻² yr⁻¹) at temperate and high latitudes and a low amplitude over the subtropical band (Fig. 5b).

Figure 6 Long-term trend in (a) the coastal FCO₂ (in mol C m⁻² yr⁻¹ decade⁻¹), (b) the air–sea pCO₂ gradient (ΔpCO₂, in µatm decade⁻¹), (c) the wind speed at 10 m above the sea surface (m s⁻¹ decade⁻¹) and (d) the sea-ice coverage (decade⁻¹) from 1982 to 2020. For each panel, the long-term trend is calculated as the slope of a linear regression on the monthly median values of all the deseasonalized data from 1982 to 2020. We only present grid cells where a significant trend is detected based on a Mann–Kendall statistical test.

3.3.2 Long-term trends in the coastal CO₂ sink

The rate of change in coastal FCO₂ and the various parameters involved in the FCO₂ calculation (i.e., ΔpCO₂, wind speed and sea-ice coverage) from 1982 to 2020 are presented in Fig. 6. Our results reveal significant spatial heterogeneities between the long-term temporal FCO₂ trends (linear trends that span over 39 years) observed within different coastal regions, a finding consistent with the range of varying slopes (including changes in sign of the slopes) already reported in local regional and discontinuous global studies (e.g., Becker et al., 2021; Laruelle et al., 2018; Wang et al., 2017). Our results also show that the rates of changes in ΔpCO₂ and FCO₂ follow each other (compare Fig. 6a with 6b). Coastal regions with negative (positive) ΔpCO₂ slopes present negative (positive) FCO₂ slopes, which translate into a stronger sink/weaker source (weaker sink/stronger source). Most coastal regions (∼60 % of the grid cells that present a significant trend using a Mann–Kendall statistical test with a significance threshold of 95 %) exhibit negative ΔpCO₂ and FCO₂ slopes (i.e., stronger sinks or weaker sources; blue colors in Fig. 6a and b) in agreement with past studies (e.g., Laruelle et al., 2018; Resplandy et al., 2023; Wang et al., 2017). Positive ΔpCO₂ and FCO₂ slopes (weaker sink or stronger sources; red colors) can also be observed such as along the Mediterranean Sea or Southeast Asia. Stronger FCO₂ rates of change (absolute value >0.6 mol m⁻² yr⁻¹ decade⁻¹) are mainly observed in mid to high latitude coastal regions and along upwelling regions (e.g., Moroccan upwelling current), while low latitude coastal regions show weaker slopes.

Although our results suggest that the long-term change in FCO₂ intensity mainly results from that of the ΔpCO₂ (compare Fig. 6a with 6b), the rate of change in FCO₂ can be amplified or dampened in some regions by changes in wind speed patterns and/or sea-ice coverage (through their effect on Eq. 1), in agreement with recent findings by Resplandy et al. (2023). For most of the coastal ocean, an increase in wind speed has been observed over the study period (positive slope, red colors in Fig. 6c) with a median value for the rate of change value of 0.04 m s⁻¹ decade⁻¹. This increase in wind speed promotes the FCO₂ exchange through its effect on the gas exchange transfer velocity (stronger sinks/sources). Rate changes in sea-ice coverage reveal a general retreat of sea-ice in the Northern Hemisphere (negative slope, blue colors in Fig. 6d) and a gain along the Antarctic shelf (positive slope, red colors) in agreement with, for example, Serreze and Meier (2019). A decrease in sea-ice coverage favors air–sea CO₂ exchange over a larger coastal surface area and during longer periods of the year, both of which strengthen, for instance, the CO₂ sink in coastal regions at high latitudes of the Northern Hemisphere.

Figure 7 Evolution of the global coastal CO₂ sink (FCO₂, in Pg C yr⁻¹) over time using the reconstructed coastal pCO₂-product (solid black line) with its associated uncertainties (dash-dotted black line; see Sects. 2.4 and 3.4 for further details). We use a 300 km distance from the coast as the outer limit of the coastal domain (wide coastal ocean) in (a) and the shelf break as the outer limit of the coastal domain (narrow coastal ocean) in (b).

Globally integrated, our results indicate that today's coastal ocean has been acting as a CO₂ sink since the beginning of our study period (1982) both in the wide coastal ocean (Fig. 7a) and in the narrow domain (Fig. 7b). For both domains this CO₂ sink, however, increases over time. In the wide coastal ocean, the global CO₂ uptake amounted to 0.28 Pg C yr⁻¹ in the 1980s (median value over the 1982–1992 period) and reached 0.54 Pg C yr⁻¹ in the 2010s (mean value over the 2010–2020 period) with small interannual fluctuations (∼ 0.01 Pg C yr⁻¹) of the CO₂ sink intensity diagnosed by our algorithm. The overall intensification of the coastal sink that we observed in this study (0.06 Pg C yr⁻¹ decade⁻¹ (± 0.0009 Pg C yr⁻¹ decade⁻¹ with a p value <0.05) and 0.02 Pg C yr⁻¹ decade⁻¹ (± 0.0005 Pg C yr⁻¹ decade⁻¹ with a p value <0.05) for the wide and narrow coastal domain, respectively) supports the only two available observational coastal studies performed at the global scale (i.e., Laruelle et al., 2018; Wang et al., 2017) which were, however, significantly limited by the small fraction of the coastal ocean domain investigated (e.g., 6 % in Laruelle et al., 2018), and both predict an increase in efficiency of the global coastal CO₂ sink over the past three decades. Our results are also in agreement with the conceptual approach of Bauer et al. (2013) as well as modeling studies, either using box models (e.g., Mackenzie et al., 2004, 2012; Rabouille et al., 2001; Ver et al., 1999) or, more recently, global ocean biogeochemical models (Bourgeois et al., 2016; Lacroix et al., 2021a, b) that all predict an increase in efficiency of the global coastal CO₂ sink at the century scale.

The significant strengthening of this global coastal sink that we observed in this study has approximately doubled between 1982 and 2020 (wide coastal domain) and results from a general tendency towards an increase in the coastal CO₂ sink intensities (e.g., in the high latitude of the Northern Hemisphere; Fig. 6a) combined with decreases in intensity of several CO₂ sources such as along upwelling currents (e.g., in the Arabian Sea). However, since a large fraction of the global CO₂ uptake results from coastal regions >40^∘ of the Northern Hemisphere, and since these CO₂ sink regions present strong negative rates of change in FCO₂ (Fig. 6a), our results suggest that the primary driver of this 2-fold increase in the global coastal CO₂ sink is to be found in the high latitudes of the Northern Hemisphere which contribute disproportionately to the global scale coastal FCO₂ trend. Further studies should, however, be carried out to support this conclusion, given the paucity of observational pCO₂ data in those high latitude regions that translate into high uncertainties in our pCO₂-product, as for example in the Arctic Ocean. Taking also the large heterogeneity in the long-term FCO₂ trends, a quantitative analysis of the respective contributions of different coastal systems to the global strengthening of the coastal CO₂ sink should also be performed in the future, using a regionalized approach. Moreover, changes in wind speed and sea-ice coverage have likely not been constant over time, and further analysis of their influence on the rate change of FCO₂ should be analyzed for each decade individually to better understand the interplay between these different drivers. Overall, our results highlight the complex nature of the coastal FCO₂ dynamics and emphasize the need for further investigation and understanding of the specific factors influencing the FCO₂ trends in different coastal regions.

Table 3 Global FCO₂ and pCO₂ uncertainties calculated for the reconstructed data products using a wide delimitation and a narrow delimitation (given in parentheses) of the coastal domain.

* The numbers correspond to uncertainties calculated using a wide coastal delimitation, while those enclosed in brackets represent uncertainties calculated using a narrow coastal delimitation.

3.4 Uncertainties associated with the data products

The global coastal CO₂ sink of −0.40 Pg C yr⁻¹ (wide coastal domain) that we calculate in this study using the ULB–SOM–FFN–coastalv2 pCO₂-product is associated with a relative uncertainty that amounts to ∼10 % (value of 0.03 Pg C yr⁻¹ for ${\mathit{\sigma }}_{{\mathrm{FCO}}_{\mathrm{2}}}$; see Eq. 2 and Table 3). This global uncertainty mainly results from the uncertainty associated with the oceanic pCO₂ (${\mathit{\sigma }}_{\mathrm{\Delta }{\mathrm{pCO}}_{\mathrm{2}}}$, uncertainty of 0.02 Pg C yr⁻¹). The choice of the gas exchange transfer velocity formulation yields a 7 % difference (σ_k, uncertainty of 0.01 Pg C yr⁻¹) on the global FCO₂ calculation, whereas we calculate a ∼4 % difference on FCO₂ depending on the wind speed product choice (σ_wind, uncertainty of 0.01 Pg C yr⁻¹) or on the sea-ice product choice (σ_ice, uncertainty of 0.01 Pg C yr⁻¹) on the global FCO₂ calculation. It is noteworthy, though, that the uncertainty in the mean is substantially lower than that calculated for individual months (see Fig. 7) or regions (e.g., see discussion above regarding the Baltic Sea region) due to compensating errors as was also identified in a study by Gloege et al. (2022).

The total uncertainty of the oceanic pCO₂ (${\mathit{\theta }}_{{\mathrm{pCO}}_{\mathrm{2}}}$, value of 0.63 µatm; see Eq. 3) mainly results from the SOM–FFN mapping method to reconstruct the coastal pCO₂-product (θ_map=28.82 µatm; Table 3). This uncertainty falls within the range of values reported in the literature from different statistical interpolation methods to generate coastal pCO₂ data products (RMSE values generally between 10 and 35 µatm; see Chen et al., 2016) which are calculated from regional studies and would be expected to be smaller than those calculated for global scale analysis (or even the performance of our algorithm at the scale of its provinces, which generally cover a much larger surface area than most regional studies). The θ_map uncertainty calculated in this study is, however, higher than reported for the open ocean (typical RMSE values <20 µatm; e.g., Landschützer et al., 2014), mainly because of the complex biogeochemical dynamics and larger variability observed in the coastal seas compared with the open ocean. We calculate a global value of 3.08 µatm for θ_obs, the uncertainty on the sampling in the field of the observations from the SOCAT database, which is slightly higher than the value reported by Pfeil et al. (2013; value of 2 µatm). For θ_grid, the uncertainty associated with the meshing of the observations from SOCAT to gridded 0.25^∘ monthly maps, we calculate a global mean value of 5.70 µatm, which is close to the value reported by Sabine et al. (2013; 5 µatm) for the open ocean. It should be noted that all these uncertainties are calculated globally and can be larger at the regional scale (e.g., see Roobaert et al., 2019) as exemplified by the uncertainty associated with the choice of wind speed product on the FCO₂ calculation (see Roobaert et al., 2018). Moreover, due to the temporal heterogeneity of the data coverage in the SOCAT database, our FCO₂ uncertainties can also vary temporally. As shown in Fig. 7, the global FCO₂ uncertainties that we report for each year (dash-dotted black lines) are largest in the 1980s (e.g., global ${\mathit{\sigma }}_{{\mathrm{FCO}}_{\mathrm{2}}}$ value of 0.11 Pg C yr⁻¹ in 1987) because of the scarcity of pCO₂ measurements before 1990 in the SOCAT database and decrease over time. Our global uncertainties are also slightly larger along the nearshore domain of the coastal ocean. Using the narrow definition of the coastal domain (i.e., the shelf break as the outer limit), we calculate a global value of 0.01 Pg C yr⁻¹ for ${\mathit{\sigma }}_{{\mathrm{FCO}}_{\mathrm{2}}}$ (7 % uncertainty on the global FCO₂, which is consistent with the global FCO₂ uncertainty calculated by Roobaert et al. (2019; 10 %)). A 7 %, 2 % and 8 % FCO₂ differences are calculated depend on the k formulation used (σ_k value of 0.01), the wind product choice (σ_wind=0.002 Pg C yr⁻¹) and the sea-ice choice (σ_ice=0.01 Pg C yr⁻¹), respectively with the narrow coastal domain. For θ_map, θ_grid and θ_obs, we calculate values of 35 µatm, 8 µatm and 3 µatm, respectively.

4 Data availability

The ULB–SOM–FFN–coastalv2 pCO₂- and FCO₂-products can be found at https://doi.org/10.25921/4sde-p068 (Roobaert et al., 2023). The bathymetry is derived from the 2 min global ETOPO2 database (Table 1; US Department of Commerce, 2006), the Chl a from the monthly 4 km merged GlobColour product for the 1998–2020 period from the EU Copernicus Marine Service information (https://doi.org/10.48670/moi-00281, Copernicus Marine Service, 2023), the SST and SSS from the daily NOAA OI SST V2 (Reynolds et al., 2007) and from the daily Hadley center EN4 SSS (Good et al., 2013), respectively. We use the atmospheric pCO₂-product from the SeaFlux product (Fay et al., 2021) which is calculated from the dry air mixing ratio of CO₂ (xCO₂) provided by the ESRL surface marine boundary layer CO₂ product (Dlugokencky et al., 2019; https://www.esrl.noaa.gov/gmd/ccgg/mbl/data.php) with a vapor correction according to Dickson et al. (2007) and using the ERA5 sea level pressure (Hersbach et al., 2023; https://doi.org/10.24381/cds.f17050d7). The sea-ice coverage is derived from the monthly mean of the daily 0.25^∘ dataset of Reynolds et al. (2007) and the wind speed from the 6 h first moment of the 0.25^∘ resolution product of the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 (Hersbach et al., 2023). The pCO₂ observations are derived from the SOCAT database v2022 (Bakker et al., 2022). The wind speed is calculated from the monthly mean of the 0.25^∘ resolution product of the ECMWF ERA5 wind product (https://doi.org/10.24381/cds.f17050d7; Hersbach et al., 2023), which has a native temporal resolution of 6 h.

5 Conclusions and directions for future research

The release of the global coastal pCO₂-product in 2017 by Laruelle et al. (2017) was a significant step forward for the investigation of the spatial distribution of CO₂ sources and sinks as well as their seasonal variabilities in the shallow portion of the ocean. It was also instrumental to the completion and harmonization of global ocean air–sea CO₂ fluxes (Fay et al., 2021), hence supporting global carbon budget analyses (Friedlingstein et al., 2022). However, this product was not designed or evaluated regarding its ability to resolve the interannual and decadal variabilities and the long-term evolution of the coastal air–sea CO₂ exchange, which are still poorly understood (e.g., Bauer et al., 2013; Lacroix et al., 2021a; Laruelle et al., 2018; Regnier et al., 2013, 2022; Resplandy et al., 2023; Wang et al., 2017). In this study, we presented a new coastal pCO₂-product for the 1982–2020 period using ∼18 million direct coastal observations from the SOCATv2022 database (Bakker et al., 2022) combined with an updated version of the coastal two-step SOM and FFN method used by Laruelle et al. (2017). We also provided a new coastal air–sea CO₂ exchange product for the same period and examined the long-term trends, that is, the temporal evolution of the global coastal CO₂ sink over the past four decades. This analysis reveals that the long-term trend of the air–sea pCO₂ gradient drives most of the long-term evolution of the coastal CO₂ sink, wind speed and sea-ice coverage playing a significant role regionally. Trend analysis of the coastal FCO₂ has also been attempted using global ocean pCO₂-products that cover the coastal domain (Resplandy et al., 2023). However, these investigations have been inconclusive, likely because global ocean pCO₂-products cannot yet sufficiently well capture the specific and changing conditions occurring along the coastal domain (e.g., Chau et al, 2022; Rödenbeck et al., 2013; see also Resplandy et al., 2023). Our updated coastal pCO₂-product circumvents these limitations and provides a first robust assessment against which outputs from global oceanic modes results can be evaluated (e.g., Resplandy et al., 2023). It will thus help better constrain the anthropogenic perturbation of the global ocean carbon cycle. In the future, our machine learning approach could also be used to diagnose the main drivers of change in the global coastal ocean sink and, more specifically, changes in the long-term trend evolution of the coastal pCO₂ field. This approach, in conjunction with process-based simulations, is critically needed to evaluate and mitigate the impact of multiple anthropogenic perturbations (e.g., atmospheric pCO₂ increase, physical climate, eutrophication and hypoxia) on the global coastal carbon cycle and associated biodiversity loss and other marine stressors.