CMEMS-LSCE: a global, 0.25°, monthly reconstruction of the surface ocean carbonate system

Observation-based data reconstructions of global surface ocean carbonate system variables play an essential role in monitoring the recent status of ocean carbon uptake and ocean acidification, as well as their impacts on marine organisms and ecosystems. So far, ongoing efforts are directed towards exploring new approaches to describe the complete marine carbonate system and to better recover its fine-scale features. In this respect, our research activities within the Copernicus Marine Environment Monitoring Service (CMEMS) aim to develop a sustainable production chain of observation-derived global ocean carbonate system datasets at high space–time resolutions. As the start of the long-term objective, this study introduces a new global 0.25^∘ monthly reconstruction, namely CMEMS-LSCE (Laboratoire des Sciences du Climat et de l'Environnement) for the period 1985–2021. The CMEMS-LSCE reconstruction derives datasets of six carbonate system variables, including surface ocean partial pressure of CO₂ (pCO₂), total alkalinity (A_T), total dissolved inorganic carbon (C_T), surface ocean pH, and saturation states with respect to aragonite (Ω_ar) and calcite (Ω_ca). Reconstructing pCO₂ relies on an ensemble of neural network models mapping gridded observation-based data provided by the Surface Ocean CO₂ ATlas (SOCAT). Surface ocean A_T is estimated with a multiple-linear-regression approach, and the remaining carbonate variables are resolved by CO₂ system speciation given the reconstructed pCO₂ and A_T; 1σ uncertainty associated with these estimates is also provided. Here, σ stands for either the ensemble standard deviation of pCO₂ estimates or the total uncertainty for each of the five other variables propagated through the processing chain with input data uncertainty. We demonstrate that the 0.25^∘ resolution pCO₂ product outperforms a coarser spatial resolution (1^∘) thanks to higher data coverage nearshore and a better description of horizontal and temporal variations in pCO₂ across diverse ocean basins, particularly in the coastal–open-ocean continuum. Product qualification with observation-based data confirms reliable reconstructions with root-mean-square deviation from observations of less than 8 %, 4 %, and 1 % relative to the global mean of pCO₂, A_T (C_T), and pH. The global average 1σ uncertainty is below 5 % and 8 % for pCO₂ and Ω_ar (Ω_ca), 2 % for A_T and C_T, and 0.4 % for pH relative to their global mean values. Both model–observation misfit and model uncertainty indicate that coastal data reproduction still needs further improvement, wherein high temporal and horizontal gradients of carbonate variables and representative uncertainty from data sampling would be taken into account as a priority. This study also presents a potential use case of the CMEMS-LSCE carbonate data product in tracking the recent state of ocean acidification. The data associated with this study are available at https://doi.org/10.14768/a2f0891b-763a-49e9-af1b-78ed78b16982 (Chau et al., 2023).

1 Introduction

Between 1750 and 2019, the ocean took up an estimated 25 % (or 170 ± 20 Pg C) of total cumulated anthropogenic CO₂ (685 ± 75 Pg C) emitted to the atmosphere (IPCC AR6 – the Sixth Assessment Report of the United Nations Intergovernmental Panel on Climate Change, Canadell et al., 2021). While the uptake of anthropogenic CO₂ mitigates global warming, it also profoundly modifies seawater chemistry in a suite of well-understood reactions (Orr et al., 2005), leading to an increase in hydrogen ion concentration ([H⁺]), as well as a decrease in carbonate ion concentration ([${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$]) and in the saturation state of seawater (Ω) with respect to calcium carbonate minerals (CaCO₃). The increase in hydrogen ion concentration ([H⁺]) is commonly reported as a decrease in pH (pH = −log [H⁺]) and referred to as ocean acidification.

Changes in carbonate chemistry impact calcifying plankton and benthos as a direct result of decreasing seawater saturation state with respect to CaCO₃ (Fabry et al., 2008; Thomsen et al., 2015). Ocean acidification also modifies the production of marine trace gases exchanged at the air–sea interface (Hopkins et al., 2020), the availability of nutrients fuelling primary production (Doney et al., 2009), and the speciation of pollutants (Millero et al., 2009; Hoffmann et al., 2012). These chemical changes interact with warming and ocean deoxygenation to drive major changes in marine ecosystems (Doney et al., 2020) and to alter global biogeochemical cycles with the potential for feeding back to radiative forcing (Gehlen et al., 2011; Hopkins et al., 2020). The likelihood of major disruptive impacts of ocean acidification on marine ecosystems, if future CO₂ emissions were to go unabated, is reflected by the Sustainable Development Goal 14.3 (SDG 14.3) – namely to “Reduce Ocean Acidification: minimize and address impacts of ocean acidification” (https://www.globalgoals.org/14-life-below-water, last access: 20 March 2023). While not specifically mentioned, moving towards SDG 14.3 implies the understanding of historical and contemporary carbonate chemistry and its mean state, trends, and variability.

In situ time series have played an important role in monitoring ocean acidification over the last decades (Bates et al., 2014; Lauvset et al., 2015; Sutton et al., 2019; Pérez et al., 2021; Leseurre et al., 2022; Skjelvan et al., 2022). At these sites, seawater pH (Ω) has been either directly measured or calculated from measurements of other carbonate system variables. These variables include surface ocean partial pressure of CO₂ (pCO₂), total alkalinity (A_T), and dissolved inorganic carbon (C_T). While changes in the time series of carbonate system variables reflect well the impacts of enhanced anthropogenic CO₂ uptake on ocean chemistry at a local scale (Steinberg et al., 2001; González-Dávila and Santana-Casiano, 2009; Dore et al., 2009; Bates et al., 2014; Pérez et al., 2021), the reliable upscaling to large ocean regions or entire basins requires a significant extension of the existing observing network (Lauvset et al., 2015; Bakker et al., 2016; Sutton et al., 2019; Lauvset et al., 2022a).

Time series data are completed by bottle data from international cruises. These data are synthesized by the Global Ocean Data Analysis Project v2.2022 (GLODAPv2.2022) and include about 1.4 million measurements of surface-to-interior ocean pH, A_T, C_T, and other parameters (Lauvset et al., 2022b,  https://www.glodap.info/, last access: 30 September 2022). Likewise, underway measurements of near-surface CO₂ fugacity, i.e. pCO₂ corrected for non-ideal gas behaviour, have been compiled in the Surface Ocean CO₂ ATlas (SOCAT) since its first release in 2011 (Pfeil et al., 2013). The latest version, SOCATv2022, yields approximately 33.7 million high-quality controlled data (Bakker et al., 2022,  http://www.socat.info/, last access: 17 June 2022). Despite millions of observations being available, data coverage is still modest; e.g. CO₂ fugacity samples over the global ocean cover less than 2 % of its surface for each month in the last 3 decades (Bakker et al., 2016; Hauck et al., 2020). Mapping methods have thus become an essential tool in ocean carbon cycle research, allowing us to interpolate or extrapolate these sparse measurements into space–time-varying fields of carbonate system variables.

Recent years have seen the rapid development of machine learning approaches to map global surface ocean pCO₂ (see Rödenbeck et al., 2013; Landschützer et al., 2016; Denvil-Sommer et al., 2019; Gregor et al., 2019; Chau et al., 2022b,  for instance). Thanks to these efforts, the carbon cycle community can now draw on an ensemble of reconstructions for the observation-based assessment of the ocean carbon sink (Friedlingstein et al., 2022). However, only a few global observation-based reconstructions are available for pH, A_T, C_T, and Ω with respect to calcite and aragonite (see  Gregor and Gruber, 2021, for a review). The reconstruction of global distributions of these variables is hampered by an insufficient amount of direct measurements (Bakker et al., 2016; Lauvset et al., 2022a). Alternatively, the complete carbonate system can be obtained by speciation given the information of any coupling of pCO₂, pH, A_T or C_T together with chemical (e.g. phosphate, silicate, nitrate) and physical variables (e.g. temperature, salinity), as well as corresponding dissociation constants (Park, 1969; Lewis and Wallace, 1998; Dickson et al., 2007).

Regardless of the developments in different observation-based estimation methods, Takahashi et al. (2014), Iida et al. (2021), and Gregor and Gruber (2021) propose global climatologies or monthly varying fields of all variables of the carbonate system, i.e. pCO₂, pH, A_T, C_T, and Ω. These data products have a spatial resolution of 1^∘ (∼ 100 km × 100 km) or even coarser. Nevertheless, the variations in carbonate system variables over the coastal regions where their instantaneous gradients are driven by smaller-scale features like ocean upwelling, wind turbulence, eddies, water runoff, and sharp biological productivity (Jones et al., 2012; Bakker et al., 2016; Laruelle et al., 2017) are poorly described at such spatial resolutions. Here, we improve on existing studies by providing a global, 0.25^∘, monthly observation-based product consisting of datasets of six core surface ocean carbonate system variables of the marine carbonate system and their associated 1σ uncertainty. This high-resolution data product covers the years from 1985 to 2021. Laboratoire des Sciences du Climat et de l'Environnement (LSCE) is in charge of the product within the European Copernicus Marine Environment Monitoring Service (CMEMS). Our product is referred to as CMEMS-LSCE hereafter.

The reconstruction of surface ocean carbonate system variables starts with the reconstruction of surface ocean pCO₂ and A_T in each regular grid of 1 month × 0.25^∘ × 0.25^∘ in the period 1985–2021 (444 months in total). Next, the variables of pH (reported on total scale), C_T, and Ω are derived by speciation. The advantages of the combination of pCO₂ and A_T over others for the speciation of the carbonate system are as follows: (1) pCO₂ is the most extensively measured parameter; (2) A_T can be accurately predicted from salinity, temperature, and nutrient concentrations; and (3) the combination of these two prior variables results in the slightest uncertainty in pH estimates (Zeebe and Wolf-Gladrow, 2001; Lauvset and Gruber, 2014; Takahashi et al., 2014; Orr et al., 2018). The three main successive modules used in the CMEMS-LSCE production chain are summarized as follows:

• i.

  Reconstruction of pCO₂ (Sect. 3.1). The CMEMS-LSCE-FFNN ensemble-based approach (Chau et al., 2022b) is modified to map gridded datasets of SOCATv2022 CO₂ fugacity and predictors in order to reconstruct pCO₂ at a finer-spatial-scale resolution of 0.25^∘ for every month in the period 1985–2021 (444 months in total). The spatial resolution of new feed-forward neural networks (FFNNs) is 16-fold higher than the original. By design, the ensemble of model outputs allows us to yield the best model estimate (i.e. ensemble mean) and model uncertainty (i.e. ensemble standard deviation) at each 0.25^∘ grid cell and for each month. Global monthly reconstructions of pCO₂ proposed by this study complement the previous climatological product by Landschützer et al. (2020), i.e. a combination of the two existing datasets covering, respectively the open ocean at 1^∘ (Landschützer et al., 2016) and the coastal sector at 0.25^∘ (Laruelle et al., 2017).

• ii.

  Reconstruction of A_T (Sect. 3.2). Locally interpolated alkalinity regression (LIAR;  Carter et al., 2016, 2018) is chosen to estimate monthly total alkalinity over the global surface ocean. Various reconstruction methods for A_T exist (see Carter et al., 2016; Broullón et al., 2019; Gregor and Gruber, 2021,  for a review), but we choose LIAR due to its global applicability, simplicity in setting, and accuracy compared to other published approaches (Carter et al., 2018; Gregor and Gruber, 2021). Importantly, LIAR allows us to determine reconstruction uncertainty propagated from multiple sources of input uncertainties at desired model resolutions.

• iii.

  Reconstruction of pH, C_T, and saturation states with respect to aragonite (Ω_ar) and calcite (Ω_ca) (Sect. 3.3). CO2SYS (Lewis and Wallace, 1998; Van Heuven et al., 2011) is a standard software used for the speciation of carbonate parameters in the marine CO₂ system (see  Olsen et al., 2016; Bresnahan et al., 2021; Gregor and Gruber, 2021; Woosley, 2021, for a few of said parameters). A complement to the CO2SYS software (CO2SYS.v2) developed by Orr et al. (2018) is used to quantify the uncertainty associated with these carbonate system variables. All the input data uncertainties are propagated through the CO2SYS.v2 processing chain.

The global, monthly, 0.25^∘ resolution datasets of surface carbonate variables are intensively evaluated against different observation-based products independent from our model fitting at a global scale to in situ locations. In Sect. 4, multiple metrics are proposed for product analyses and assessments. Results are presented in Sect. 5 with an emphasis on the evaluation of the best reconstruction and associated model uncertainty for each variable (Sect. 5). This section also highlights the advantages obtained with an increase in spatial resolution and presents an application of the CMEMS-LSCE product in tracking ocean acidification over the last 3 decades. Section 7 summarizes key results, discusses the potential for future model upgrades, and introduces possible product use cases. The high-resolution data product described in this paper (netCDF format) can be accessed via repository under the following data DOI: https://doi.org/10.14768/a2f0891b-763a-49e9-af1b-78ed78b16982 (Chau et al., 2023).

Bakker et al. (2022)Good (2022); Good and Worsfold (2022)Nardelli et al. (2016); Droghei et al. (2018)Menemenlis et al. (2008)Maritorena et al. (2010)Chevallier et al. (2005, 2010); Chevallier (2013)Takahashi et al. (2009)Garcia et al. (2019)

Table 1 Input data used in the reconstructions of CMEMS-LSCE carbonate system variables over the global ocean in 1985–2021.

Note: Last access was on 15 April 2022 for all input databases, except for SOCATv2022 data (17 June 2022) and WOA18 data (30 July 2022). Data products 1–8 are used in the pCO₂ reconstruction. Products 2–3 and 9–11 are used to compute A_T, C_T, pH, Ω_ar, and Ω_ca. Global values of product uncertainty (data^† or analysis errors expressed as σ) have been reported for $f{{\mathrm{CO}}_{\mathrm{2}}}^{†}$ (< 5 µatm), SST (0.15 ^∘C), SSS (0.2 PSU), SSH (0.02 m), MLD (–), Chl a (0.03 mg m⁻³), xCO₂ (–), $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{clim}}$ (10 µatm), NO₃ (1.8 µmol kg⁻¹), SiO₂ (3.6 µmol kg⁻¹), and PO₄ (0.12 µmol kg⁻¹).

2 Data used and reprocessing

2.1 Input data products for surface ocean carbonate system reconstructions

Many observation-based products are used as predictors of our target carbonate system variables (Table 1). Global ocean maps of sea surface temperature (SST), sea surface salinity (SSS), sea surface height above geoid (SSH), and chlorophyll a (Chl a) come from the Copernicus Marine Environment Monitoring Service (CMEMS; Good et al., 2020; Nardelli et al., 2016; Droghei et al., 2018; Maritorena et al., 2010). Mixed-layer depth (MLD) fields belong to Estimating the Circulation and Climate of the Ocean project Phase II (ECCO2; Menemenlis et al., 2008). CO₂ mole fractions (xCO₂) are derived from the CO₂ atmospheric inversion of the Copernicus Atmosphere Monitoring Service (CAMS; Chevallier et al., 2005, 2010; Chevallier, 2013). Surface ocean concentrations of nitrate (NO₃), silicate (SiO₂), and phosphate (PO₄) are extracted from the World Ocean Atlas 2018 (WOA18; Garcia et al., 2019). The climatological pCO₂ ($p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{clim}}$) product is provided by Lamont Doherty Earth Observatory (LDEO; Takahashi et al., 2009). Details of these products, including resource access, data coverage, and resolutions, are presented in Table 1.

With the exception of xCO₂, nutrient concentrations, and $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{clim}}$, these input data products have original resolutions equivalent to or even finer than a spatial resolution of 0.25^∘ and a temporal resolution of monthly. When mismatches in data resolutions appear, input data products are interpolated to fit the pre-defined model resolutions. The monthly average is applied for daily (hourly) data products, and the Climate Data Operators (CDO) remap operator is used to regrid the input data fields with spatial resolutions different from 0.25^∘. The datasets of SST and xCO₂ – the two key variables driving global pCO₂ changes (Bates et al., 2014; Gruber et al., 2019; Landschützer et al., 2019; Chau et al., 2022b; Friedlingstein et al., 2022) – cover the full learning period and the whole globe, as expected. The other predictor data are not available before the 1990s, when new types of satellite measurements started, and one of them (i.e. Chl a) does not cover the high latitudes of the winter hemisphere. We therefore gap-fill the time series in an ad hoc manner, as in previous studies (Landschützer et al., 2016; Gregor et al., 2019; Chau et al., 2022b). Monthly climatologies of SSS, Chl a, and MLD computed based on the available data are used for each missing year. Likewise, climatologies plus linear trends of SSH following global warming effects serve for the pre-1993 period. Missing Chl a data in the high latitudes of the winter hemisphere are replaced by the minimum concentration of Chl a over the available data for the same grid cell (∼ 0.01 mg m⁻³). WOA18 nutrients and LDEO $p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{clim}}$ are already climatologies per se, and we apply them for all the analysis years (1985–2021). The sine function is applied to convert latitude, while both the sine and cosine are used to transform longitude to conserve their periodical behaviours.

CO₂ fugacity from SOCATv2022 (Bakker et al., 2022) is used as the target data in our monthly pCO₂ reconstructions. The SOCAT project collects and qualifies underway observations via international vessels, moorings, or autonomous platforms. It grids the observations at spatial resolutions of 1^∘ or 0.25^∘, resulting in the two major SOCAT gridded data products. The temporal resolution of these two products is monthly. While the 1^∘ data product (SOCATv2022r100) covers the global ocean, the 0.25^∘ product covers solely the coastal regions. The SOCAT coastal areas are within 400 km from the shoreline (Sabine et al., 2013; Bakker et al., 2016); see Fig. A1a in the Appendix for an illustration. To merge the two resolutions, we first duplicate the 1^∘ open-ocean SOCATv2022 data (∼ 2 × 10⁵ data points) over its 16 0.25^∘ sub-cells. These 0.25^∘ open-ocean data are then combined with the 0.25^∘ coastal–ocean SOCATv2022 data (∼ 4 × 10⁵ data points) to generate a global monthly 0.25^∘ ocean data product fed to our reconstruction model of pCO₂ (Sect. 3.1). The merged SOCATv2022 product at monthly and 0.25^∘ resolutions is referred to as SOCATv2022r025 hereafter. The assumption of open-ocean data homogeneity of pCO₂ within 1^∘ grid boxes (∼ 100 km × 100 km) does not degrade the reconstruction skill over the global open ocean (see Sect. 5 for results), where pCO₂ observations are spatially auto-correlated within a global median distance of 400 ± 250 km (Jones et al., 2012). The data distribution of SOCATv2022 CO₂ fugacity before and after combining is shown in Fig. A2 and Table 3.

2.2 Product qualification and comparison

The monthly, 0.25^∘ resolution reconstructions of carbonate system variables are qualified with gridded observation-based datasets and in situ time series which are not used in our model fitting (Table 2).

• The SOCAT data in each reconstruction month are excluded from the model fitting, which avoids overfitting and ensures fairness in the model evaluation (Chau et al., 2022b). The global monthly, 0.25^∘ resolution CMEMS-LSCE-FFNN pCO₂ fields can therefore be evaluated against the pCO₂ data converted from SOCATv2022 CO₂ fugacity (Körtzinger, 1999) at the same resolution. Doing this, CMEMS-LSCE-FFNN pCO₂ is assessed with more than 32 × 10⁵ open-ocean data and 4 × 10⁵ coastal–ocean data (Table 3). The SOCATv2022 data have low random uncertainty (2–5 µatm), but the spatio-temporal sampling bias from the month and grid centres is significant (Bakker et al., 2016). The 0.25^∘ data reconstruction is also compared to its previous version with a spatial resolution of 1^∘ (Chau et al., 2022a, b).

• The monthly, 0.25^∘ reconstructions of A_T, C_T, and pH are qualified based on GLODAPv2.2022 data (Lauvset et al., 2022b). GLODAP provides non-gridded datasets of ocean carbon variables which have been compiled and bias-corrected from water samples taken at various depths. The measurement uncertainty is 4 µmol kg⁻¹ in A_T and C_T and between 0.01–0.02 in pH (Lauvset et al., 2022a). Only direct measurements at depths shallower than 10 m and with a flag of 2 (best quality control) are selected for this evaluation. Measurements in each box of 1 month × 10 m × 0.25^∘ × 0.25^∘ are averaged to obtain representative data of surface A_T, C_T, and pH at 0.25^∘ grid cells for months in the period 1985–2021. This results in roughly 16 × 10³ data points for A_T and C_T over the global ocean (Table 4). Only half of that amount stems from direct pH measurements. Another half, referred to as indirect measurements (i.e. pH calculated with A_T and C_T), is excluded from this data evaluation. Over 30 % of GLODAP data are distributed along the coasts.

• In situ time series of direct measurements of carbonate system variables (pCO₂, A_T, C_T, and pH) are used to qualify our product at the local scale (Table 2).

  • a.

    Sutton et al. (2019) present data over multiple sites equipped with autonomous moorings that have been measuring surface ocean pCO₂ and pH since 2004. These stations cover a wide range of oceanic conditions from the tropics to high latitudes (Fig. A1b). More than half of 42 stations are distributed over the continental shelves, and many of them observe pCO₂ and pH in the regime of coral reefs (Tables A2 and A3). Measurement uncertainty is up to 2 µatm for pCO₂ and 0.02 for pH.

  • b.

    For A_T and C_T, we consider eight accessible time series to provide insights into changes in the surface ocean carbonate system over the recent decades (Bates et al., 2014; Coppola et al., 2020; Gregor and Gruber, 2021; Pérez et al., 2021; Leseurre et al., 2022). Two of them belong to the biogeochemical observing systems located in the subtropical Atlantic: the Bermuda Atlantic Time Series (BATS, Michaels and Knap, 1996; Steinberg et al., 2001) and the European Station for Time Series in the Ocean Canary islands (ESTOC, González-Dávila and Santana-Casiano, 2009). The other two provide direct measurements in the same ocean basin but in the subpolar region: the Irminger Sea and Iceland Sea (IRMINGER and ICELAND; Olafsson et al., 2010). Further mentions are made of time series distributed in specific conditions, including a high-Arctic fjord in Svalbard Fischer et al. (AWIPEV, 2017); Gattuso et al. (AWIPEV, 2023a) and the Mediterranean basin (DYFAMED, Coppola et al., 2021). Hawaii Ocean Time series in the subtropical North Pacific (HOT, Dore et al., 2009) and OISO time series in the subpolar Southern Ocean (KERFIX, Metzl and Lo Monaco, 1998) allow us to complete the reconstruction qualification in different ocean basins. Measurement uncertainty (for instance, from replication experiments) reported at these sites is below 3 µmol kg⁻¹. Except for HOT and ESTOC (surface ocean measurements), data over all the stations are extracted from seawater depths shallower than 10 m. A monthly average is applied for all the mentioned time series in order to be compatible with output from the CMEMS-LSCE chain of models.

Bakker et al. (2022)Chau et al. (2022a, b)Lauvset et al. (2022a, b)Sutton et al. (2019)Michaels and Knap (1996)Steinberg et al. (2001)Coppola et al. (2020, 2021)González-Dávila and Santana-Casiano (2009)Dore et al. (2009)Fischer et al. (2017)Gattuso et al. (2023b)Olafsson et al. (2010)Metzl and Lo Monaco (1998)

Table 2 Data sources used in product evaluation and comparison. Values in brackets for each variable present measurement-based data uncertainty or analysis^‡ uncertainty.

3 Reconstruction methods

3.1 Ensemble pCO₂ mapping feed-forward neural networks

The CMEMS-LSCE-FFNN (Chau et al., 2022b) is based on an ensemble of 100 feed-forward neural network (FFNN) models mapping SOCAT CO₂ fugacity (fCO₂) and predictor variables (Eq. 1).

$$\begin{array}{}\text{(1)}& \begin{array}{rl}f{\mathrm{CO}}_{\mathrm{2}}=& \mathrm{FFNN}(\mathrm{SST},\mathrm{SSS},\mathrm{SSH},\mathrm{Chl}-a,\mathrm{MLD},\\ & x{\mathrm{CO}}_{\mathrm{2}},f{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{clim}},\mathrm{latitude},\mathrm{longitude})\end{array}\end{array}$$

The predictors of fCO₂ include sea surface temperature (SST), salinity (SSS), surface height (SSH), chlorophyll a (Chl a), mixed-layer depth (MLD), CO₂ mole fraction (xCO₂), fCO₂ climatologies (f${\mathrm{CO}}_{\mathrm{2}}^{\mathrm{clim}}$), and the geographical coordinates (latitude and longitude). The datasets of SOCAT fCO₂ and predictors are first reprocessed to match model-fitting requirements (Sect. 2.1). After excluding the data in the reconstruction month, the data within the 3-month window are randomly separated into FFNN training and validation subsets with a ratio of 2:1. The subsampling process is repeated for every 100 FFNN runs, which results in 100 different datasets for model fitting. The excluded SOCATv2022 datasets are used in model evaluation. The CMEMS-LSCE-FFNN approach was originally developed for pCO₂ reconstructions at monthly, 1^∘ resolutions where pCO₂ is converted from fCO₂ following the formulation by Körtzinger (1999). The model best estimate and its uncertainty are defined as the ensemble mean (μ) and ensemble spread (σ) of 100 model outputs of pCO₂.

This study slightly modifies the CMEMS-LSCE-FFNN ensemble approach by Chau et al. (2022b) to achieve pCO₂ reconstructions at monthly, 0.25^∘ resolutions. Some of the input datasets presented here (Table 1) are different from those presented in Chau et al. (2022b) (Table S1 in the Supplement). The up-to-date input datasets have higher resolutions and a better coverage over the coastal ocean, as well as in the high latitudes. Furthermore, the new CMEMS data resources offer space–time-varying uncertainty fields which are important in quantifying carbonate system variable uncertainties.

For comparable evaluations in this study, we execute 100-member ensembles of FFNN models at spatial resolutions of both 1^∘ (FFNNr100) and 0.25^∘ (FFNNr025) using the same lot of input data resources (Table 1). Be reminded that the training data of fCO₂ are extracted from the SOCATv2022r100 product for FFNNr100, while they come from the SOCATv2022r025 product (i.e. the merged product of the 1^∘ open-ocean dataset and the 0.25^∘ coastal–ocean dataset) for FFNNr025. All input datasets are reprocessed with respect to each model resolution (Sect. 2.1). Section 3.1 compares these two CMEMS-LSCE-FFNN versions and highlights the skill of the finer-resolution data product.

3.2 Locally interpolated alkalinity regression

Locally interpolated alkalinity regression (LIAR;  Carter et al., 2016, 2018) is an ensemble-based regression method developed for the global reconstruction of total alkalinity (A_T). Regression coefficients were learnt based on GLODAPv2 data (Olsen et al., 2016) binned within regular windows of 5^∘ × 5^∘. For prediction, the LIAR software interpolates between the regression coefficients to arbitrary resolutions specified by the users. This study employs eight LIAR models (Carter et al., 2018, Table 2) for calculating A_T at monthly, 0.25^∘ resolutions. Each model represents a combination of predictor variables (see the full presentation in Eq. 2):

$$\begin{array}{}\text{(2)}& A_{\mathrm{T}}=\mathrm{LIAR}(\mathrm{SSS},\mathrm{SST},{\mathrm{NO}}_{\mathrm{3}},{\mathrm{SiO}}_{\mathrm{2}}).\end{array}$$

The eight regression models include salinity (SSS) – the predominant predictor of A_T – and some combinations of temperature (SST), nitrate (NO₃), and silicate (SiO₂). The model which has the smallest propagation uncertainty is chosen to provide the best estimate of A_T.

Global monthly total alkalinity and 1σ uncertainty are estimated with given input data from the monthly CMEMS SSS and SST fields and from the WOA18 datasets of nutrient concentrations (Table 1). Uncertainty in terms of the A_T field is estimated systematically through input uncertainty propagation along the processing chain (Carter et al., 2018). Here we define the input uncertainty of predictors in terms of standard deviations (1σ). Input uncertainty fields associated with the monthly CMEMS SSS and SST are products' analysis errors (see e.g. Fig. A7), while uncertainties in terms of the WOA18 NO₃ and SiO₂ climatologies are set to 15 % of data values per cell. The 15 % quantity refers to the median percentage of standard analysis errors against climatological means of nutrient concentrations (see product standard errors in Table 6, Garcia et al., 2019). The WOA18 standard analysis errors are defined as misfits between their interpolated data and GLODAPv2 bottle data (Olsen et al., 2016). The spatial distribution of the error percentage of the WOA18 nutrient concentrations at the ocean surface is illustrated in Fig. A8.

3.3 Carbonate system speciation

The CO2SYS speciation software was first developed by Lewis and Wallace (1998) to determine carbonate system parameters in the marine CO₂ system based on a set of equilibrium equations (Dickson et al., 2007). Here we use the speciation program written by Van Heuven et al. (2011) and its extension (CO2SYS.v2) with uncertainty propagation proposed by Orr et al. (2018). To obtain a complete description of the ocean carbonate system, the CO2SYS.v2 initialization requires the following input conditions:

• i.

  values of any coupling of the parameters pCO₂, A_T, C_T, and pH

• ii.

  temperature and pressure

• iii.

  total concentrations of all the non-CO₂ acid-base systems

• iv.

  equilibrium constants used to describe seawater acid-base chemistry.

The third condition involves total concentrations of both conservative and non-conservative constituents in the non-CO₂ acid-base systems. The amount of conservative constituents such as borate, fluoride, and sulfate in surface seawater is estimated with salinity. The total concentration of non-conservative constituents (nutrients) is computed approximately with silicate (SiO₂) and phosphate (PO₄). Further information on the carbonate system speciation can be found in Dickson et al. (2007) and Dickson (2010).

With the reconstructions of pCO₂ and A_T (Sects. 3.1 and 3.2), the CO2SYS.v2 speciation software is used to derive pH (on total scale), C_T, Ω_ar, and Ω_ca and to determine their uncertainty over the ocean surface at a resolution of 0.25^∘. Equation (3) expresses all input–output variables of CO2SYS.v2 for this study. Note that the estimates for other carbonate system variables such as hydrogen ion (H⁺) concentration and Revelle factor (RF) – a measure of the carbonate buffer capacity – are also available (Figs. A4 and A6) but are beyond the scope of our data evaluation.

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{pH},C_{\mathrm{T}},{\mathrm{\Omega }}_{\mathrm{ar}},{\mathrm{\Omega }}_{\mathrm{ca}}=& \text{CO2SYS.v2}(p{\mathrm{CO}}_{\mathrm{2}},A_{\mathrm{T}},\\ & \mathrm{SST},\mathrm{SSS},P,{\mathrm{SiO}}_{\mathrm{2}},{\mathrm{PO}}_{\mathrm{4}},\text{constants})\end{array}\end{array}$$

The FFNN best estimate (ensemble mean) of pCO₂ reconstructions (Sect. 3.1) and the LIAR outputs of A_T (Sect. 3.2) are used as the prior inputs of the CO2SYS.v2 at each grid cell for every month in the period 1985–2021. We take the same data products of SST, SSS, and nutrient concentrations as for the previous reconstructions (Table 1). Pressure (P) is assumed to be 0 dbar at the ocean surface. For equilibrium constants, we choose the best empirical values recommended by Dickson et al. (2007) and Dickson (2010). These settings include the dissociation constants K₁ and K₂ from Lueker et al. (2000) and $K_{{\mathrm{HSO}}_{\mathrm{4}}}$ from Dickson (1990) in combination with the total boron ratio salinity formulation by Uppstrom (1974).

The uncertainty of the CO2SYS.v2 variables is estimated by error propagation (Orr et al., 2018). Inputs for the CO2SYS.v2 error propagation include the reconstruction uncertainty of pCO₂ (FFNN ensemble standard deviation) and of A_T (LIAR error propagation). The uncertainties of SST, SSS, and nutrient concentrations are set to the same values as in the previous section (Sect. 3.2). Equilibrium constants' standard errors are default values (see Table 1, Orr et al., 2018). As for FFNN and LIAR, uncertainty values of each carbonate system variable are computed for each month in 1985–2021 and at each 0.25^∘ grid box over the global surface ocean.

4 Evaluation metrics

4.1 Model best estimate and uncertainty quantification

The 100 FFNN models result in an ensemble of 100 estimates of global monthly, 0.25^∘ surface ocean pCO₂ fields (Sect. 3.1). Specifying any $t=\mathrm{1}:\mathrm{444}$ (month), $i=\mathrm{1}:\mathrm{180}$ (latitude), and $j=\mathrm{1}:\mathrm{360}$ (longitude), the best estimate (μ_tij) and uncertainty (σ_tij) at time t and grid cell ij are deduced from 100 FFNN pCO₂ estimates $\left(X\right(t,i,j,m){)}_{m=\mathrm{1}}^{m=\mathrm{100}}$ as follows:

$$\begin{array}{}\text{(4a)}& {\displaystyle }{\mathit{\mu }}_{tij}={\displaystyle \frac{{\sum }_{m=\mathrm{1}}^{m=\mathrm{100}}X(t,i,j,m)}{\mathrm{100}}},\text{(4b)}& {\displaystyle }{\mathit{\sigma }}_{tij}=\sqrt{{\displaystyle \frac{{\sum }_{m=\mathrm{1}}^{m=\mathrm{100}}{\left[X(t,i,j,m)-{\mathit{\mu }}_{tij}\right]}^{\mathrm{2}}}{\mathrm{100}}}}.\end{array}$$

For pH, A_T, C_T, Ω_ar, and Ω_ca, the best estimates and associated uncertainties (μ_tij and σ_tij) are obtained directly from the LIAR and CO2SYS.v2 speciation tools and their error propagation (Sects. 3.2 and 3.3).

To assign representatives of μ and σ estimates for carbonate system variables at a specific space–time window, we define statistics with respect to each of the three following cases:

• i.

  a representative over a period of time (T months)

  $$\begin{array}{}\text{(5a)}& {\displaystyle }{\mathit{\mu }}_{ij}& {\displaystyle }={\displaystyle \frac{{\sum }_t{\mathit{\mu }}_{tij}}{T}}\text{(5b)}& {\displaystyle }{\mathit{\sigma }}_{ij}& {\displaystyle }=\sqrt{{\displaystyle \frac{{\sum }_t{\mathit{\sigma }}_{tij}^{\mathrm{2}}}{T}}}\end{array}$$
• ii.

  a representative over a region (e.g. ocean basins and sub-basins, the global ocean)

  $$\begin{array}{}\text{(6a)}& {\displaystyle }{\mathit{\mu }}_t& {\displaystyle }={\displaystyle \frac{{\sum }_{ij}{\mathit{\mu }}_{tij}\times A_{ij}}{{\sum }_{ij}{A}_{ij}}}\text{(6b)}& {\displaystyle }{\mathit{\sigma }}_t& {\displaystyle }=\sqrt{{\displaystyle \frac{{\sum }_{ij}{\mathit{\sigma }}_{tij}^{\mathrm{2}}\times A_{ij}}{{\sum }_{ij}{A}_{ij}}}}\end{array}$$
• iii.

  a representative over a period of time and a region.

  $$\begin{array}{}\text{(7a)}& {\displaystyle }\mathit{\mu }& {\displaystyle }={\displaystyle \frac{{\sum }_{t,ij}{\mathit{\mu }}_{tij}\times A_{ij}}{T\times {\sum }_{ij}{A}_{ij}}}\text{(7b)}& {\displaystyle }\mathit{\sigma }& {\displaystyle }=\sqrt{{\displaystyle \frac{{\sum }_{t,ij}{\mathit{\sigma }}_{tij}^{\mathrm{2}}\times A_{ij}}{T\times {\sum }_{ij}{A}_{ij}}}}\end{array}$$

In the above, t is in a time period with length T, and A_ij is the area of each grid cell in a desired region. It is noteworthy that the statistics in Eqs. (5b)–(7b) are not the standard deviation associated with the mean quantities in Eqs. (5a)–(7a), but they stand for the best representative of uncertainty estimates over an ocean basin and/or time period. These statistics also offer support for the comparison with model–observation deviations (e.g. Eq. 10), which has typically been used in the calculation of standard uncertainty, as proposed in the previous studies (Jiang et al., 2019; Iida et al., 2021; Gregor and Gruber, 2021). Subscripts in the notations of the best model estimates (μ) and model uncertainties (σ) in Eqs. (4)–(7) are dropped for general situations.

The best model estimate (μ) is assessed against model uncertainty (σ) through the σ-to-μ ratio (%):

$$\begin{array}{}\text{(8)}& R(\mathit{\mu },\mathit{\sigma })=\mathrm{100}\mathit{\%}{\displaystyle \frac{\mathit{\sigma }}{\left|\mathit{\mu }\right|}}.\end{array}$$

The σ-to-μ ratio allows an evaluation of the significance of the model estimate. Model estimates of carbonate variables are reliable with R(σ,μ) values of less than 20 % (Rose, 2013).

4.2 Model performance in comparison with evaluation data

Assume that μ_tij and O_tij are the best model estimate and an observation (or its gridded data) available at time t, and grid cells ij, μ, and O are their means over the total number of evaluation data (N). Model skills are assessed against observation data (Table 2) with the following metrics:

• mean model–observation differences (bias)

  $$\begin{array}{}\text{(9)}& \mathrm{Bias}={\displaystyle \frac{{\sum }_{t,ij}({\mathit{\mu }}_{tij}-O_{tij})}{N}}\end{array}$$

• root-mean-square deviation (RMSD)

  $$\begin{array}{}\text{(10)}& \mathrm{RMSD}=\sqrt{{\displaystyle \frac{{\sum }_{t,ij}({\mathit{\mu }}_{tij}-O_{tij}{)}^{\mathrm{2}}}{N}}}\end{array}$$

• coefficient of determination (r²).

  $$\begin{array}{}\text{(11)}& r^{\mathrm{2}}={\displaystyle \frac{\left[{\sum }_{t,ij}\right({\mathit{\mu }}_{tij}-\mathit{\mu })\times (O_{tij}-O){]}^{\mathrm{2}}}{{\sum }_{t,ij}({\mathit{\mu }}_{tij}-\mathit{\mu }{)}^{\mathrm{2}}\times {\sum }_{t,ij}(O_{tij}-O{)}^{\mathrm{2}}}}\end{array}$$

5 Results

5.1 Surface ocean pCO₂

This section presents the reconstruction of surface ocean pCO₂ at monthly and 0.25^∘ resolutions. The reconstruction skill is evaluated against SOCATv2022 data from global to ocean basin scales and at the level of grid cells (Sect. 2.2). We compare the novel reconstruction at a higher spatial resolution to the one at a coarser spatial resolution (Chau et al., 2022b). Emphasis is put on the skill to reproduce spatial and temporal variations of pCO₂ across a variety of coastal regions and time series stations.

Figure 1 presents global maps at a 0.25^∘ resolution of long-term averages of pCO₂ and corresponding uncertainty estimates. Reconstructed pCO₂ distributions reveal well-documented large-scale structures. Values are high over upwelling regions (e.g. Equatorial Pacific, California Boundary Current, western Arabian Sea). Low pCO₂ is associated with increased CO₂ solubility in cold high-latitude seawater (e.g. Arctic), strong biological production (e.g. China Sea), or the combination of both (e.g. subpolar North Atlantic, Southern Ocean between 35–50^∘ S). Spatial structures appear to be coherent from small to large spatial scales, both along the coast and moving towards the open ocean (see also Figs. 2–4). The combination of a down-scaled version of open-ocean and higher-resolution coastal SOCATv2022 data (Sect. 2.1) yields pCO₂ distributions without discontinuities. The uncertainty map (Fig 1b) represents the confidence level in surface ocean pCO₂ estimates (Fig 1a). Predominantly low uncertainty estimates (σ<5 µatm) indicate the global stability of the ensemble reconstruction. Exceptions are found in many coastal regions, open-ocean areas with sparse data coverage (e.g. Southern Pacific, Indian Ocean), or regions with substantially high or low surface ocean pCO₂ (e.g. Arctic, eastern equatorial Pacific). However, pCO₂ is reconstructed with a high degree of confidence over most of the global ocean with a σ-to-μ ratio (Eq. 8) below 10 % (Fig. A9a).

Figure 1 CMEMS-LSCE-FFNN pCO₂ over the global ocean at a spatial resolution of 0.25^∘. Temporal means of the model best estimate and 1σ uncertainty per grid cell over 1985–2021 are calculated by using Eq. (5).

Figure 2 Comparison of CMEMS-LSCE-FFNN mapping pCO₂ at 1^∘ (r100) and 0.25^∘ (r025) resolutions over four permanent upwelling regions associated with the eastern boundary currents (California, Humboldt, Canary, and Benguela; see Fig. A1 (A, B, G, and H) for geographical locations). For each region, spatial distributions of pCO₂ (μ) and uncertainty (σ) estimates and the coastal–ocean RMSD of pCO₂ averaged over 1985–2021 (Eqs. 5 and 10) are shown. Metrics presented in the legend for each of the third rows include the number of coastal–ocean SOCATv2022 data (N), regional RMSD (Eq. 10), and r² (Eq. 11).

Figure 3 Comparison of CMEMS-LSCE-FFNN mapping pCO₂ at 1^∘ (r100) and 0.25^∘ (r025) resolutions over four regions characterized by low pCO₂ values (Labrador, western South Atlantic, northern Europe, and Japan; see Fig. A1 (C, E, F, and J) for geographical locations). For each region, spatial distributions of pCO₂ (μ) and uncertainty (σ) estimates and the coastal–ocean RMSD of pCO₂ averaged over 1985–2021 (Eqs. 5 and 10) are shown. Metrics present in the legend for each of the third rows include the number of coastal–ocean SOCATv2022 data (N), regional RMSD (Eq. 10), and r² (Eq. 11).

Figure 4 Comparison of CMEMS-LSCE-FFNN mapping pCO₂ at 1^∘ (r100) and 0.25^∘ (r025) resolutions over the mouth of the river Amazon and the Western Arabian Sea (see Fig. A1 (D and I) for geographical locations). For each region, spatial distributions of pCO₂ (μ) and uncertainty (σ) estimates, and coastal–ocean RMSD of pCO₂ averaged over 1985–2021 (Eqs. 5 and 10) are showed. Metrics present in the legend for each of the 3rd row include the number of coastal–ocean SOCATv2022 data (N), regional RMSD (Eq. 10) and r² (Eq. 11).

Table 3 Evaluation statistics for monthly CMEMS-LSCE-FFNN reconstructions of pCO₂ at 1^∘ (r100) and 0.25^∘ (r025) spatial resolutions computed over the period 1985–2021. r100 → 025 and r025 → 100 refer to the versions downscaled or upscaled from the original CMEMS-LSCE-FFNN pCO₂ at 1^∘ and 0.25^∘ resolutions. SOCATv2022 gridded data independent from CMEMS-LSCE-FFNN training are used as benchmarks for model evaluation (see text for details). Statistics including total numbers of data, RMSD (Eq. 10), and r² (Eq. 11) are reported for both the open ocean (O) and coastal region (C). ^∗ Marks results with respect to the primary product proposed in this study.

Skill scores of the monthly, 0.25^∘ resolution reconstruction are presented in Table 3 (columns marked by an asterisk). The global RMSD (Eq. 10) between the best reconstruction and SOCATv2022r025 pCO₂ over the entire period is 14.3 µatm for the open ocean and 27.6 µatm for the coastal ocean. These two model errors are lower than 4 % and 8 % of the global mean pCO₂ (Table 6). Moreover, variability present in observation-based data is reproduced by the CMEMS-LSCE-FFNN with high values of r² (open ocean: 0.83, coast: 0.74). The reconstruction quality is similar among major ocean basins. Spatial distributions of SOCATv2022 data, bias, and RMSD are shown in Figs. A2b and d and A3b, d, f, and h. Estimation skills are low in the ocean basins, with sparse data coverage and significant space–time variability of pCO₂ (e.g. Arctic, eastern equatorial Pacific, land–ocean continuum).

Table 3 also presents statistics for the monthly FFNN products of surface ocean pCO₂ at spatial resolutions of 0.25^∘ (r025) and 1^∘ (r100) together with their variants (r100 → 025 and r025 → 100). The latter are, respectively extrapolation and interpolation versions of the original r100 and r025 datasets. We used the Climate Data Operators (CDO) remap operator to regrid FFNN model outputs to a finer or coarser spatial resolution. For compatibility, we compare statistics between

• i.

  FFNN(r025) and FFNN(r100 → 025) by using SOCATv2022r025 as evaluation data

• ii.

  FFNN(r025 → 100) and FFNN(r100) by using SOCATv2022r100 as evaluation data.

The FFNN(r025) central to this study yields a lower RMSD and a higher correlation to the SOCAT data than the FFNN(r100 → 025). As expected, the improvement in reconstruction skill with higher model resolution is larger over coastal regions than in the open ocean. The FFNN(r025) product after interpolating to a coarser resolution, i.e. FFNN(r025 → 100), agrees with the original 1^∘ resolution data product over the whole of the ocean.

The motivation to increase the spatial resolution of the reconstruction is to improve the representation of horizontal gradients of pCO₂ at fine scales. Figures 2–4 exemplify spatial distributions for the two reconstructions (r025 and r100) over the coastal–open-ocean continuum. A total of 10 distinct oceanic regions are considered (see Fig. A1a and Table A1 for the 10 locations), and these can be classified into three groups:

• permanent eastern-boundary-current upwelling systems with relatively high pCO₂ (California Current System – CCS, Humboldt Current System – HCS, Canary Current System – CnCS, and Benguela Current System – BCS)

• regions characterized by low pCO₂ values driven by cold water temperatures and strong biological production (Labrador Sea, western South Atlantic, northern Europe, and Sea of Japan)

• other regions either under the influence of strong river runoff (Amazon mouth) or monsoon-driven upwelling (western Arabian Sea).

The legend of Figs. 2–4 includes regional RMSD and r² computed between the best estimates of two models and coastal–ocean SOCATv2022r025 data. The coarser spatial resolution product is co-located at the same 0.25^∘ grid cells for this analysis. These figures illustrate important discrepancies in pCO₂ data density between coastal regions with poorly monitored regions (e.g. HCS, BCS, Amazon mouth) in contrast with areas with higher data coverage (e.g. northern Europe, Sea of Japan).

Over 7 out of the 10 analysed regions with reconstruction at monthly, 0.25^∘ resolutions yield RMSDs below 10 % of the global mean of coastal–ocean pCO₂ estimates (Table 6) and r² values higher than 0.3 – e.g. northern Europe (RMSD = 32.4 µatm, r²=0.81), Sea of Japan (RMSD = 20.8 µatm, r²=0.71), and CnCS (RMSD = 28.8 µatm, r²=0.4). The CMEMS-LSCE-FFNN model projections of pCO₂ lack skill over the HCS (RMSD = 50.0 µatm, r²=0.33), the region under influence of the Amazon River (RMSD=45.7 µatm, r²=0.38) and the western Arabian Sea (RMSD = 44.8 µatm, r²=0.49). In nearshore sectors of these coastal areas, pCO₂ estimates are also subject to a substantial amount of uncertainty (σ > 20 µatm). The lack of model skill reflects the combination of low data density and strong pCO₂ gradients driven by multiple underlying physical and biogeochemical processes. The HCS, for instance, is characterized by the highest pCO₂ levels (Fig. 2) among the four eastern boundary current systems, with interannual variability being amplified with the El Niño–Southern Oscillation (ENSO) events (Feely et al., 1999; Landschützer et al., 2016). Similarly, high pCO₂ levels with substantial seasonal variability are observed over the western Arabian Sea (Fig. 4b), the key driver being monsoonal upwelling (Sabine et al., 2002; Sarma et al., 2013, 2023). In contrast to the two aforementioned coastal regions, high CO₂ undersaturation and strong pCO₂ gradients (Fig. 4a) are found in the area under the influence of Amazon River discharge (Olivier et al., 2022). Extreme values and large variability of pCO₂ challenge any approach to estimating pCO₂ data over these regions (Ibánhez et al., 2015; Bakker et al., 2016; Landschützer et al., 2020).

The two FFNN reconstructions (r025 and r100) share similarities in terms of their overall structures of pCO₂ over the coastal–open-ocean continuum (Figs. 2–4). However, the higher spatial resolution outperforms its lower-resolution counterpart in reproducing fine-scale features of pCO₂ in the transition from nearshore regions to the adjacent open ocean. The increase in model spatial resolution translates into a greater spatial coverage of the continental shelves such as the Labrador Sea, northern Europe, and the Sea of Japan (Fig. 3) and thus an increase in the number of data over the coastal domain. The increase in spatial resolution allows a gain in the prediction probability of pCO₂ variations on the order of roughly 2 % over the eastern boundary currents to 7 % over the western South Atlantic (Figs. 2 and 3b).

Figure 5 Time series of surface ocean pCO₂ (µatm) at coastal observing stations (Table A2 and Fig. A1b): model best estimate (curve), 1σ uncertainty (envelope), and monthly average of in situ observations (point). The reconstructed data at 1^∘ (r100) and 0.25^∘ (r025) resolutions are co-located to in situ observations provided by Sutton et al. (2019). Means of the best estimate and 1σ uncertainty (μ±σ) calculated over the observing time span are shown in brackets. Statistics include number of months with observations (N), bias, RMSD, and r² computed for the two reconstructions. ${\mathit{\sigma }}_{p{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{t}}$ stands for temporal standard deviation from monthly averages of pCO₂ observations.

At the local scale, the reconstruction of in situ pCO₂ (Sutton et al., 2019) over the open ocean is on the high order of confidence (Table A3). Low RMSD (between 7.8 and 23.5 µatm) and sustainably high r² (from 0.45 to 0.98) dominate evaluation statistics over the 18 open-ocean stations. Obviously, CMEMS-LSCE-FFNN has less skill in the coastal sector, and model–observation deviation varies depending on a wide range of pCO₂ conditions. However, coastal–ocean RMSD can be smaller than 10 % of the station climatology (e.g. KILONALU, KANEOHE, ALAWAI), and the reproduction availability of temporal variations of pCO₂ possibly exceeds 70 % (e.g. SEAK, KODIAK, DABOB). Through Fig. 5, we further assess seasonal to inter-annual variability reproduced at the eight coastal sites (see Fig. A1b and Table A2 for station location) where measurements are available for both pCO₂ and pH (analysed in Sect. 5.3) and where they are poorly constrained by the 1^∘ reconstruction (Chau et al., 2022b). The temporal variability of pCO₂ reported for these time series sites reflects a combination of processes (Sutton et al., 2019), e.g. the California Current System (CAPEARAGO and CCE2), the western coastal upwelling (CAPEELIZABETH), eutrophication enhancing respiration of CO₂ (FIRSTLANDING), and multiple stressors on coral reef environments (CHEECAROCKS, GREYREEF). As shown in Fig. 5 (scattered points for observations), time series of coastal pCO₂ are still short. The longest time series covers 127 months of pCO₂ monitoring since 2010 (CCE2), while the shortest one contributes 17 months with observations (FIRSTLANDING).

Analysing the station time series, we have found that data have been sampled within a few days, with an average offset of about 1 week from the month centre. At these coastal sites, the temporal standard deviation from monthly averages of pCO₂ (${\mathit{\sigma }}_{p{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{t}}$) exceeds analytical errors (2 µatm, Sutton et al., 2019). ${\mathit{\sigma }}_{p{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{t}}$ ranges from 20.1 µatm at GREYREFF to values as large as 65.6 µatm at CAPEARAGO or 70 µatm at FIRSTLANDING. The monthly average of pCO₂ might not be adequately represented by discreet samples at sites with a large temporal standard deviation of pCO₂. The misfit between the monthly reconstruction and discreet observations is exacerbated in dynamical coastal environments and might explain, in part, the large RMSD of reconstructions of monthly coastal pCO₂ (e.g. GREYREEF: 38.3 µatm, CAPEARAGO: 79.9 µatm, FIRSTLANDING: 77.3 µatm) for the r025 reconstruction. The RMSD is mostly lower for the FFNN reconstruction at 0.25^∘ resolution compared to the FFNN at 1^∘ resolution by 2.2 µatm (CCE2) to 23.2 µatm (COASTALMS). Similarly, r² increases between 7 %–23 % at higher resolutions. Overall, seasonal to interannual variations of coastal–ocean pCO₂ are better reproduced in the reconstruction at 0.25^∘ resolution (Fig. 5).

Figure 6 CMEMS-LSCE A_T and C_T over the global ocean at a spatial resolution of 0.25^∘. Temporal means of the model best estimate and 1σ uncertainty per grid cell over 1985–2021 are calculated by using (Eq. 5).

5.2 Total alkalinity and dissolved inorganic carbon

This section presents and analyses global ocean surface reconstructions of total alkalinity (A_T) and dissolved inorganic carbon (C_T) at monthly, 0.25^∘ resolutions over 1985–2021. GLODAPv2.2022 bottle data (Sect. 2.2) serve as reference data for model evaluation. Model reconstruction skill is further assessed at the eight Eulerian time series sites: AWIPEV, BATS, DYFAMED, ESTOC, HOT, ICELAND, IRMINGER, and KERFIX (Table 2).

Figure 6 shows spatial distributions of the climatological mean and uncertainty (Eq. 5) for A_T and C_T. Despite being, in part, influenced by common biological and physical processes, both properties have contrasting distributions due to the strong correlation between surface ocean A_T and salinity (Lee et al., 2006; Broullón et al., 2019), as well as the contribution of air–sea gas exchange and biological productivity on surface ocean C_T levels (Feely et al., 2001; Takahashi et al., 2014). Over subtropical Atlantic gyres and the Mediterranean Sea, oceanic areas with net evaporation, A_T exceeds 2400 µmol kg⁻¹. Total alkalinity falls below 2150 µmol kg⁻¹ in regions where precipitation, river freshwater runoff, or seasonal sea ice melting dilute surface water salinity (e.g. subpolar North Pacific, Arctic, and equatorial river outflows). The distribution of C_T is relatively uniform between the Atlantic, Pacific, and Indian Ocean basins but shows pronounced latitudinal gradients. High concentrations of C_T are found throughout the Southern Ocean (C_T >2100 µmol kg⁻¹) where strong upwelling brings up subsurface water enriched in CO₂ and nutrients. The inefficient utilization of nutrients in this high-nutrient, low-chlorophyll region limits the biological drawdown of C_T, allowing the massive C_T input to be spread horizontally by westerlies (Key et al., 2004; Menviel et al., 2018). Levels of C_T below 1900 µmol kg⁻¹ are reconstructed over the equatorial Pacific, the equatorial eastern Atlantic, the eastern Indian Ocean, and coastal areas on the Arctic Ocean. While low C_T levels associated with equatorial upwelling reflect gas exchanges across the air–sea interface and enhanced biological production, the interaction between physical and biogeochemical processes at work in the Indian Ocean are less well understood (Takahashi et al., 2014). Low C_T levels found close to river mouths reflect outgassing of CO₂ across the salinity gradient, as well as enhanced biological uptake fuelled by river nutrient inputs. Representation uncertainty (Fig. 6c and d) associated with monthly alkalinity and C_T reconstructions is lower than 20 µmol kg⁻¹ throughout the open ocean. The open-ocean σ-to-μ ratio (Eq. 8) ranges between 0.5–1.5 %, which is relatively small (Fig. A9c and d). C_T uncertainty is computed through CO2SYS.v2 error propagation with reconstruction uncertainties of pCO₂ and A_T set as inputs. The largest values (σ > 30 µmol kg⁻¹) appear nearshore and at surrounding oceanic islands (Fig. 6d). A similar feature is found on the field of A_T (Fig. 6c), inherited from input uncertainty associated with the CMEMS salinity product (Fig. A7a).

Figure 7 Monthly time series of A_T and C_T at BATS, DYFAMED, ESTOC, and HOT stations (Table 2 and Fig. A1b): model best estimate (curve), 1σ uncertainty (envelope), and monthly average of surface (0–10 m) observations (point). Means of the best estimate and 1σ uncertainty (μ±σ) calculated over the observing time span are shown in brackets if accessible. Statistics include the number of months with observations (N), bias, RMSD, and r². ${\mathit{\sigma }}_{A_{\mathrm{T}}}^{\mathrm{t}}$ (${\mathit{\sigma }}_{C_{\mathrm{T}}}^{\mathrm{t}}$) stands for the temporal standard deviation from monthly averages of A_T (C_T) observations.

Figure 8 Spatial distribution of reconstruction skills for A_T, C_T, and pH over 1985–2021. Mean model–data difference (bias) and root-mean-square deviation (RMSD) between the reconstruction and GLODAPv2.2022 surface data (0–10 m) at a spatial resolution of 0.25^∘. The size of grid cells is scaled upon a better visualization.

Table 4 Skill scores computed between CMEMS-LSCE and GLODAPv2.v2022 in A_T, C_T, and pH over the period 1985–2021. Total numbers of data, RMSD (Eq. 10), and r² (Eq. 11) are reported for both the open ocean (O) and coastal region (C). Basin identification is shown in Fig. A1.

We qualify monthly, 0.25^∘ reconstructions of A_T and C_T with GLODAPv2.2022 data (Lauvset et al., 2022a) for the 37-year period (Table 4 and Fig. 8). The global open-ocean reconstruction scores an RMSD of 22.1 µmol kg⁻¹ and a r² of 0.9 in A_T. Similar numbers are found for C_T (RMSD = 22.7 µmol kg⁻¹ and r²=0.9). The model scores the good fit in the open Indian Ocean with RMSD smaller than 15.5 µmol kg⁻¹ and r² above 0.92 for both variables. The reconstruction deviates from GLODAP data in the western North Atlantic, subpolar North Pacific, tropics, and nearby major rivers (Fig. 8a, b, c, and d).

A_T and C_T are underestimated in the continental shelves of north Alaska and the northeastern Atlantic, the Mediterranean Sea, the South China Sea, and nearby river plumes (Fig. 8a and c). The Arctic yields the poorest estimations among all the ocean basins with a global RMSD over 100 µmol kg⁻¹ (Table 4). The prediction probability of variability in A_T (C_T)] is relatively large for the open ocean at 79 % (71 %) but is rather unsatisfying over the coastal ocean (46 % [40 %]). Extrapolating these carbonate variables towards the shore remains challenging, with much higher errors and uncertainty estimates obtained over the continental shelf compared to the open-ocean reconstruction (Table 4, Figs. 6c, d and 8a–d). The coastal–ocean errors are on the order of 10 % of the global mean values of A_T and C_T (Table 6).

The reconstruction of A_T distributions relies on LIAR coefficients fit with GLODAPv2 data (Olsen et al., 2016) covering the years before 2015. These data are also part of the latest version, GLODAPv2.2022 (Lauvset et al., 2022a). They therefore do not correspond to an independent dataset for the evaluation data of the CMEMS-LSCE reconstruction. To accomplish a cross-validation, reconstructions of A_T and C_T are compared to observations for eight time series stations: AWIPEV, BATS, DYFAMED, ESTOC, HOT, ICELAND, IRMINGER, and KERFIX (see Table 2 and Fig. A1b for data sources and station locations). Table A4 presents the evaluation statistics for all the stations, and Fig. 7 illustrates the comparison between monthly time series of A_T and C_T extracted from the CMEMS-LSCE datasets and measurements at the four sustained long-term monitoring sites. More than 270 (80) months in the years 1988–2021 (1995–2009 and 1998–2017) include measurements of A_T and C_T at BATS and HOT (ESTOC and DYFAMED)]. As shown in Table A4, the Arctic site (AWIPEV) provides 52-month data in 2015–2020, while the three other stations sparsely observed A_T and C_T at the surface layer, resulting in fewer than 30 monthly mean data in 2014–2021 (IRMINGER and ICELAND) and 1992–2018 (KERFIX). The reconstructed time series fit monthly averages of in situ measurements well. Mean estimates of A_T (C_T) over the observing period are about 2283 (1983) µmol kg⁻¹ at KERFIX (HOT) to 2420 (2219) µmol kg⁻¹ at DYFAMED (AWIPEV). At all the stations (DYFAMED and AWIPEV excepted) and for the two variables, the model–observation misfit is small (Bias <11 µmol kg⁻¹, RMSD < 14 µmol kg⁻¹) relative to the aforementioned mean estimates (Table A4). The highest offset between the CMEMS-LSCE estimation and observations for all the stations is found at DYFAMED (A_T: −145.1 µmol kg⁻¹, C_T: −124.7 µmol kg⁻¹). DYFAMED provides long-term time series of A_T and C_T measurements in the northwestern Mediterranean Sea (Fig. A1b). Salinity and alkalinity have substantial values due to the net evaporation (Coppola et al., 2020). The average of A_T in the Mediterranean Sea exceeds that for the global ocean by 10 % (Palmiéri et al., 2015). These characteristics set the Mediterranean Sea aside from the ocean basins. Although the bias between reanalysed SSS and observations (Fig. A10) is relatively small, LIAR (Carter et al., 2018) was trained on GLODAPv2 (Olsen et al., 2016), including only a few observations in this area. The distinct relationship between alkalinity and salinity prevailing in the Mediterranean Sea is likely not reproduced by LIAR, leading to an underestimation of A_T and a systematic bias in relation to C_T at DYFAMED (Fig. 7). ESTOC is located close to the North Atlantic east coast and is under the influence of the Canary Current System (CCS, Fig. A1). Spatial gradients and temporal variability are higher in the CCS (Fig. 2c) compared to BATS and HOT, which are both located in the centre of subtropical gyres. Despite showing good estimates of A_T and C_T in terms of RMSD at ESTOC, temporal variability of observations is reconstructed with the lowest r² (Table A4). Particularly, seasonality to multi-year variations in C_T are predicted at r²=0.47 for ESTOC compared to r² > 0.7 for AWIPEV, ICELAND, IRMINGER, BATS, and HOT. Over all the stations, the model underestimates temporal changes in A_T (Fig. 7a; BATS: r²=0.33, DYFAMED: r²=0.12, ESTOC: r²=0.03, HOT: r²=0.32), which can be attributed to the large discrepancy in variability between in situ measurements and the CMEMS time series of salinity (Fig. A10a; BATS: r²=0.33, DYFAMED: r²=0.19, ESTOC: r²=0.03, HOT: r²=0.35). Model uncertainty (1σ envelop) of monthly A_T and C_T estimates (Fig. 7a) is also inflated somewhat proportionally to the CMEMS salinity product uncertainty (Fig. A10a).

5.3 Surface ocean pH and saturation state with respect to carbonate minerals

Surface ocean pH and saturation states with respect to aragonite (Ω_ar) and calcite (Ω_ca) are critical indicators used to measure ocean acidification. This section first presents an overall evaluation of these variables. We then introduce estimates involved in the monitoring of ocean acidification in 1985–2021 as an essential application of the CMEMS-LSCE surface ocean carbon product.

Figure 9 CMEMS-LSCE pH and Ω_ar over the global ocean at a spatial resolution of 0.25^∘. Temporal means of the model best estimate and 1σ uncertainty per grid cell over 1985–2021 are calculated by using Eq. 5.

5.3.1 General analysis and evaluation

The spatial distribution of surface ocean pH reported on the scale of total hydrogen ion (H⁺) is shown in Fig. 9 (the corresponding figure for H⁺, Fig. A4, is included in the Supplement). Both temporal means of the best model estimate and 1σ uncertainty of pH share spatial patterns with pCO₂ (Fig. 1). The variables of pH and pCO₂ correlate closely through equilibrium relationships of dissolved CO₂ in seawater: an increase in pCO₂ generally corresponds to a decrease in pH. The distribution of the climatological mean of pH displays a gradient with latitudes between 8.03 and 8.11 pH units across most of the basins (Fig. 9a). Values of pH below 8 are associated with the upwelling of CO₂-rich waters (e.g. eastern equatorial Pacific, western Arabian Sea). pH exceeds 8.15 in sub- and polar cold surface water and in the regions with high biological productivity (e.g. Labrador Sea, Nordic seas, Southern Ocean between 35–50^∘ S).

The saturation state of surface ocean waters with respect to calcium carbonate minerals of aragonite and calcite is defined as the ratio of the product of the concentrations of calcium ions (Ca²⁺) and carbonate ions (${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$) to the solubility of the respective calcium carbonate mineral (CaCO₃) in surface seawater. With aragonite being the more soluble polymorph, its degree of saturation (Ω_ar) is smaller than that of calcite (Ω_ca) (Mucci, 1983). With the exception of this offset, the spatial distributions of their climatological means share common spatial patterns over the global ocean (Figs. 9b and A5a). Surface seawater is generally supersaturated, i.e. Ω_ar and Ω_ca greater than 1. The magnitude of surface ocean calcium carbonate saturation state varies with latitude. Values as large as 3.7–4.5 (5–7) for aragonite (calcite) are reconstructed in subtropical and tropical regions. Ω_ar and Ω_ca decrease toward the poles. In the Southern Ocean, surface seawater enriched in CO₂ from vertical mixing has Ω_ar (Ω_ca) values in the range of 1.5–2.1 (2–3.4). Low saturation states are also computed in the Arctic and for waters of upwelling regimes (Fig. 9b). Locally, Ω_ar drops below 1.3 and even falls under the CaCO₃ dissolution threshold of 1 (Gattuso and Hansson, 2011) in the Arctic water runoff and Baltic Sea.

The uncertainty (1σ) of pH, Ω_ar, and Ω_ca that propagated the speciation of the CO₂ system takes into account the ensemble spread of pCO₂ estimates and analysis errors of other variables (Sect. 3.3). Monthly pH uncertainty estimates fall in the 95 % confidence interval of [0.008,0.036], with a global mean value of 0.011. These estimates are in close agreement with the global uncertainty between 0.01–0.022 pH units calculated by Jiang et al. (2019), Iida et al. (2021), and Gregor and Gruber (2021). pH uncertainty is typically larger than 0.03 in the Arctic and in coastal regions (Fig. 9c). In contrast, the reconstructions of Ω_ar and Ω_ca are subject to high uncertainty (σ>0.175) between 30^∘ S–30^∘ N (Fig. 9d and A5b). Regarding the σ-to-μ ratio, mean uncertainty estimates per cell for the saturation states in the (sub-) tropical band are relatively small compared to the mean of the best monthly estimates (Fig. A9e and f). The Arctic and the coastal oceans remain the regions with the largest reconstruction uncertainties for Ω_ar and Ω_ca, as well as for pCO₂ and pH (Fig. A9a and b). Excluding these regions, R(σ,μ) (Eq. 8) is less than 0.3 % for pH and 8 % for Ω_ar and Ω_ca.

The monthly CMEMS-LSCE reconstruction at 0.25^∘ resolution is assessed against pH measurements from GLODAPv2.2022 bottle data (Table 2). For the period 1985–2021, the global RMSD amounts to 0.022 (0.060) pH units and r² scores at 0.70 (0.45) over the open (coastal) ocean (Table 4). Model bias lies within $[-\mathrm{0.01},\mathrm{0.01}$] pH units, and RMSD is below 0.02 pH units over the open ocean, except for high latitudes over 60^∘ (Fig. 8e and f). At the local scale, the time series from Sutton et al. (2019) are used for further evaluation (Tables A2 and A3). There exist many fewer evaluation data for pH than for pCO₂, e.g. only 2 months of monitoring pH at COASTALLA and FIRSTLANDING or devoid of pH measurements at equatorial observing systems. CMEMS-LSCE reconstructs pH over the open ocean with rather high scores, e.g. at BOBOA (RMSD = 0.011 and r²=0.71) and KEO (RMSD = 0.014 and r²=0.86). Referring to the eight coastal sites evaluated for pCO₂ in Sect. 5.1, RMSD can be as small as 0.035 and 0.04 pH units at CCE2 and GRAYSREEF, while it is over 0.05 pH units at the other stations (e.g. COASTALLA: 0.068, CAPEARAGO: 0.069). Similarly to the pCO₂ time series (Fig. 5, Sect. 5.1), pH has been monitored with low sampling frequency (roughly a few days in the tracking month), and the temporal sampling deviation of instantaneous observations from monthly averages (${\mathit{\sigma }}_{\text{pH}}^{\mathrm{t}}$) is significant. This temporal sampling uncertainty of pH contributes to the mismatch between model estimates and observations. For example, ${\mathit{\sigma }}_{\text{pH}}^{\mathrm{t}}$ amounts to 0.048 pH units at CCE2 and 0.020 pH units at GRAYSREEF and reaches the highest values of 0.078 pH units at COASTALLA and 0.086 pH units at CAPEARAGO. Although model–observation misfit and model uncertainty remain high over the coastal sector (see also Figs. 8e and f and 9c), their estimates do not surpass 1 % of the global mean pH (8.082). The reconstructed pH time series reproduce measurement variability with relatively high correlation, r² in [0.21,0.94], which reinforces the reliability of CMEMS-LSCE pH data.

5.3.2 Ocean acidification: key features from global to local scales

The monthly, 0.25^∘ CMEMS-LSCE datasets of pH, Ω_ar, and Ω_ca are at the basis of two CMEMS ocean indicators monitoring surface ocean acidification from 1985 to 2021: (1) annual global means and (2) global trend maps.

In Fig. 10, we present annual global means of surface ocean pH and saturation states for aragonite (Ω_ar). An illustration of calcite (Ω_ca) is provided in the Appendix (Fig. A13a). For each variable, the calculation of annual global area-weighted means of best estimates (line) and 1σ uncertainties (envelope) follows Eq. 7. The trends reported in the legend result from linear least-squares regression on annual global means of 100 ensembles of the carbonate system variables. These ensembles are generated with a Gaussian distribution having the mean and variance as the best model estimate μ and squared uncertainty (σ²) at monthly time steps and 0.25^∘ grid cells, respectively. pH decreases from 8.110 ± 0.017 in 1985 to 8.049 ± 0.014 in 2021 with a descent rate of −0.017 ± 0.004 per decade. Similar trends are found for the surface ocean saturation states with respect to calcium carbonate minerals. The global mean estimates of Ω_ar (Ω_ca) amount to 3.141 ± 0.198 (4.807 ± 0.302) and 2.862 ± 0.174 (4.372 ± 0.266) for the open and coastal oceans. The saturation state declines at a rate of −0.080 ± 0.029 per decade with respect to aragonite, while the reduction is steeper for calcite (−0.114 ± 0.045 per decade).

Figure 10 Yearly global area-weighted mean of surface seawater pH reported on total scale (a) and surface ocean saturation states with respect to aragonite (b). Global means of the best estimate (μ, plain line) and of uncertainty (σ, envelop) are computed with Eq. (7a). Trend and uncertainty in the legend are computed with linear regressions on the 100-member ensemble of yearly global means for each variable.

Figure 11 Global trend maps of surface seawater pH reported on total scale (a) and surface ocean saturation states with respect to aragonite (b). Linear trend of CMEMS-LSCE pH and Ω_ar is estimated per 0.25^∘ grid cell over 1985–2021. Cross-hatching covers the regions with uncertainty over 10 % (20 %) of pH (Ω_ar) trend estimates.

Global trend maps of surface ocean pH, Ω_ar, and Ω_ca over the entire period are illustrated in Figs. 11 and A13b. Linear least-squares regression is used to estimate secular trends at every 0.25^∘ grid cell. The linear fits of each variable against time rely on the 100-member ensemble generated with the best estimates and propagated uncertainties of pH, Ω_ar, and Ω_ca (see Fig. A14 for examples). Regression slope and residual standard deviation estimates are defined as linear trend and uncertainty of pH, Ω_ar, and Ω_ca. Hatched area represents pH (Ω_ar and Ω_ca) trend estimates (μ) with the highest uncertainties (σ), i.e. σ-to-μ ratio (Eq. 8) above 10 % (20 %). These regions include a portion of the Arctic, Antarctic, equatorial Pacific, and coastal ocean (Figs. 11, A11, and A12); 95 % of pH trend estimates over the global ocean are in the range of [$-\mathrm{0.022},-\mathrm{0.012}$] per decade (Fig. 11a). In the broad open ocean of the tropics and subtropics, pH has been declining around −0.018 per decade to −0.012 per decade. Faster decrease rates are found in the Indian Ocean and Southern Ocean, with values between −0.022 and −0.018 per decade. The fastest reductions are computed for the eastern equatorial Pacific and the Arctic with rates exceeding −0.025 per decade. A similar magnitude of pH trends over these regions is also found in (Lauvset et al., 2015; Leseurre et al., 2022; Ma et al., 2023). The spatial distribution of saturation states with respect to calcium carbonate minerals generally shows the opposite latitudinal pattern (Figs. 11b and A13b). The magnitude of Ω_ar (Ω_ca) trends over the 30^∘ S–30^∘ N band can be as large as −0.086 per decade (−0.134 per decade) to the greatest extent of −0.186 per decade (−0.275 per decade) (e.g. eastern equatorial Pacific). Trends of Ω_ar and Ω_ca computed in polar and subpolar Northern Hemisphere regions are not significant.

Bates et al. (2014)Bates et al. (2014)

Table 5 Secular trend estimates of pH and Ω_ar at seven time series stations (Bates et al., 2014). Trend and uncertainty estimates are reported as μ±σ. Monthly time series in the CMEMS-LSCE datasets are extracted at the grid box nearest to each station location (Fig. A1b). For the first three stations, this study calculates linear trends starting in the year 1985. Brackets show values computed over the full period 1985–2021.

^∗ Stations with direct observations of pH.

Trend estimates derived from reconstructions of pH and Ω_ar are evaluated at seven time series stations (Bates et al., 2014) in Table 5. Time series locations are shown in Fig. A1b. With the exception of CARIACO and HOT, for which pH measurements are available, long-term trend estimates by Bates et al. (2014) rely on time series of pH and Ω_ar calculated via speciation from measurements of A_T and C_T. The 100-member ensembles of monthly time series of pH and Ω_ar are extracted from the 0.25^∘ grid box nearest to each monitoring station. Linear least-squares regression is then used to infer estimates of their secular trends and associated uncertainties (see Fig. A14 for illustration). Trend estimates derived from CMEMS-LSCE reconstructions at HOT, BATS, ESTOC, and Munida are in line with previous studies for both pH and Ω_ar (Dore et al., 2009; González-Dávila and Santana-Casiano, 2009; Bates et al., 2014). The magnitude of the trend estimate at the Irminger Sea for 1985–2012 (pH: −0.014 ± 0.001 per decade, Ω_ar: −0.006 ± 0.011 per decade) is smaller than that determined by Bates et al. (2014). However, the CMEMS-LSCE pH trend is consistent with the estimate by Pérez et al. (2021) (−0.017 ± 0.002 per decade). Moreover, 1σ uncertainty reported for both pH and Ω_ar trend estimates by Bates et al. (2014) is large at this station (pH: −0.025 ± 0.006 per decade, Ω_ar: −0.080 ± 0.040 per decade), highlighting the associated uncertainty. Long-term trends of pH and Ω_ar are also underestimated at the Iceland Sea monitoring site, but the bias is not as large as at the Irminger Sea (Table 5). Low data-sampling frequency at these two stations (Table 1, Bates et al., 2014) could be on account of trend estimate deviation. At CARIACO, the CMEMS-LSCE time series yields a decrease in Ω_ar of −0.059 ± 0.053 per decade, relatively close to Bates et al. (2014) (−0.066 ± 0.028 per decade). The decrease in pH derived from CMEMS-LSCE is, however, larger than in Bates et al. (2014).

Table 6 Summary of global evaluation statistics for CMEMS-LSCE surface ocean carbonate system datasets at monthly, 0.25^∘ resolutions over the period 1985–2021. μ and σ stand for the global area-weighted means of monthly best estimates and 1σ uncertainties for each variable (Eq. 7). RMSD (Eq. 10) and r² (Eq. 11) are computed with SOCATv2022 for pCO₂ and GLODAPv2.2022 for pH, A_T, and C_T. The division between the coastal (C) and open (O) oceans is at 400 km distance from the shore line (Fig. A1a).

6 Data availability

The CMEMS-LSCE datasets (netCDF format) of six carbonate system variables have been delivered to the European Copernicus Marine Environment Monitoring Service (CMEMS; product ID: MULTIOBS_GLO_BIO_CARBON_SURFACE_REP_015_008, DOI: https://doi.org/10.48670/moi-00047). The high-resolution data product described in this paper can be accessed via the repository under the following data DOI: https://doi.org/10.14768/a2f0891b-763a-49e9-af1b-78ed78b16982 (Chau et al., 2023).

7 Summary

This study presents the CMEMS-LSCE product, a dataset of six carbonate system variables (Table 6) covering the global surface ocean at a spatial resolution of 0.25^∘ for every month in the period 1985–2021 (444 months). Datasets of individual carbonate system variables are built on the combination of the three methods. First, we adapt an ensemble of 100 feed-forward neural network models (CMEMS-LSCE-FFNN;  Chau et al., 2022b) to estimate surface ocean partial pressure of CO₂ (pCO₂) at the pre-defined data resolution. Second, the high-resolution total alkalinity (A_T) reconstruction is obtained by using locally interpolated alkalinity regression (LIAR;  Carter et al., 2016, 2018). Finally, surface ocean pH, total dissolved inorganic carbon (C_T), and saturation states with respect to aragonite (Ω_ar) and calcite (Ω_ca) are calculated with the carbonate system speciation software (CO2SYS.v2; Lewis and Wallace, 1998; Van Heuven et al., 2011; Orr et al., 2018), given the global monthly reconstructions of pCO₂ and A_T and other environmental input data (Sect. 3). Results are 2D fields of the best estimate and associated uncertainty (1σ) of carbonate system variables available at each grid box of 1 month × 0.25^∘ × 0.25^∘; 1σ uncertainty is referred to as the ensemble standard deviation of 100 FFNN outputs for pCO₂ while it is propagated through the processing chain of LIAR and CO2SYS.v2, taking into account different uncertainty sources of input parameters for other variables.

Multiple observation-based datasets, which are not used for the CMEMS-LSCE reconstructions at monthly and 0.25^∘ resolutions, serve as benchmarks in the assessments of product quality from global to local scales (e.g. Tables 3, 4, and A3; Figs. 2–5, 7, and 8). A summary of the primary statistics for all the six carbonate variables is presented in Table 6. Over the full period of 1985–2021, CMEMS-LSCE yields global RMSDs of 14.3 µatm and 27.6 µatm in comparison with SOCATv2022 pCO₂ for the open and coastal oceans, respectively. Temporal variability of observation-based data is well reproduced, with r² of 0.83 for the open ocean and 0.74 for the coastal domain. In comparison to CMEMS-LSCE at monthly and 1^∘ resolutions (Chau et al., 2022b), the reconstructions over coastal areas are improved at higher resolution. The monthly, 0.25^∘ reconstruction outperforms its 1^∘ counterpart in reproducing horizontal and temporal gradients of pCO₂ over a variety of oceanic regions, as well as at nearshore time series stations (Figs. 2–5). Evaluations with GLODAPv2022 bottle data and time series stations results in good reconstruction skills for A_T, C_T, and pH at monthly and 0.25^∘ resolutions (Tables 4 and A3, Figs. 7 and 8). At the global scale, the open-ocean reconstruction scores an RMSD smaller than 23 µmol kg⁻¹ and a r² of 0.9 in A_T and C_T. The model–observation deviation is higher in the coastal zone. However, it does not exceed 5 % of the global mean values, and r² is above 0.6 for both coastal A_T and C_T. Regarding pH, the CMEMS-LSCE reconstruction provides estimates with RMSD = 0.022 (0.060) and r²=0.7 (0.45) over the open (coastal) ocean. From the statistics in Tables 3 and 4, the Indian Ocean and the Southern Ocean have poor data density (Fig. A2) but generally show the best global reconstruction among the ocean basins. Thus, model evaluation with different numbers of observation data might not reflect a fair comparison of skill scores (e.g. RMSD and r²) between regions. Data density is much higher in the Arctic, Atlantic, and Pacific than in the Indian and Southern oceans. The increased data density reveals stronger spatio-temporal variability, for instance, related to coastal dynamics or upwelling than what is resolved in the two latter basins. RMSD and r² computed on the lower data variability result in better model scores.

The spatial distribution of long-term mean 1σ uncertainty estimates (Figs. 1b, 6c, and d and 9c and d) indicates higher confidence levels for open-ocean estimates than over the coastal sector. The evaluation of temporal mean 1σ uncertainty estimates relative to climatological mean values μ (Figs. 1a, 6a and b, and 9a and b) results in σ-to-μ ratios (Eq. 8) below 5 % and 8 % for pCO₂ and Ω_ar, 2 % for A_T and C_T, and 0.4 % for pH over the open ocean (Fig. A9). The σ-to-μ ratio reaches values as high as 10 % to 20 % for pCO₂ and Ω_ar in the coastal domain. The global means of open-ocean 1σ uncertainty estimates (Eq. 7a) for CMEMS-LSCE pCO₂ (8.5 µatm), A_T (16.7 µmol kg⁻¹), C_T (15.8 µmol kg⁻¹), pH (0.011), and Ω_ar (0.180) are in line with those reported by previous studies despite being derived from different statistics. For instance, Iida et al. (2021) calculated 1σ uncertainty based on the median absolute deviation of regression model fits from open-ocean observations. Their approach yielded global σ averages of 17.8 µatm, 11.5 µmol kg⁻¹, 0.018, and 0.110 for pCO₂, normalized C_T, pH, and Ω_ar, respectively. In Gregor and Gruber (2021), the authors propagated the sum squared errors (global RMSD and measurement uncertainties) of pCO₂ (15 µatm) and A_T (22 µmol kg⁻¹), obtaining global uncertainty estimates of 19 µmol kg⁻¹ in C_T and 0.022 in pH. Mean uncertainty estimates over the coastal region are on the order of 2-fold, computed for the open ocean for these four variables (Table 6), corroborating results by Gregor and Gruber (2021) (Fig. 7).

Our high-resolution carbon data product opens the door to various analyses of the marine carbonate system from the global to local scale. This study exemplifies an application of the data for monitoring ocean acidification over recent years. The monitoring indicators derived from the monthly, 0.25^∘ surface ocean CMEMS-LSCE product consist of (1) yearly global means of surface ocean pH and saturation states with respect to aragonite Ω_ar and calcite Ω_ca and (2) global maps of multi-annual trends of surface ocean pH, Ω_ar, and Ω_ca (Figs. 10, 11, and A13). In 1985, the global mean surface ocean pH was 8.110 ± 0.017. It was 8.049 ± 0.014 in 2021 (Fig. 10a). Over the same 37-year time period, Ω_ar decreased from 3.141 ± 0.198 to 2.862 ± 0.174 (Fig. 10b). The rate of decline of surface ocean pH and Ω_ar was, respectively, −0.017 ± 0.004 per decade and −0.080 ± 0.029 per decade from 1985 (see also results for Ω_ca in Sect. 5.3.2). Estimates of pH trend lie between [$-\mathrm{0.022},-\mathrm{0.012}$] per decade across most of the open ocean (Fig. 11a). In general, surface ocean pH decreased more rapidly in the Indian Ocean and Southern Ocean than in the tropics and subtropics. These findings are in close agreement with the suggestions by Lauvset et al. (2015) and Ma et al. (2023), but future studies would need to include analyses of underlying drivers to provide insight into regional differences in pH changes. By contrast, the greatest reduction in surface ocean saturation states (Fig. 11b) was computed for the two latter regions. The global trend maps of pH and Ω_ar highlight the eastern equatorial Pacific as one of the vulnerable regions with respect to ocean acidification. In this area, the decline rate of pH exceeds −0.025 per decade and −0.186 per decade for Ω_ar. The comparison of multi-annual trends of pH and Ω_ar at time series stations (Table 5 and Fig. A14) highlighted the consistency between CMEMS-LSCE estimates and previous studies (Dore et al., 2009; González-Dávila and Santana-Casiano, 2009; Bates et al., 2014; Pérez et al., 2021). For most of these sites, the trends evaluated for 1985–2021 are greater than those relative to the sub-period before the year 2012. The faster rate of ocean acidification over the full period compared to pre-2012 probably reflects a steeper acceleration in ocean uptake of anthropogenic CO₂ in the last decade. The global maps of CMEMS-LSCE pH, Ω, and their trend estimates would be potential indicators for ocean acidification along with the SDG 14.3.1 – namely “average marine acidity (pH) measured at an agreed suite of representative sampling stations” (https://sdgs.un.org/goals/goal14, last access: 31 July 2023).

The production chain of CMEMS-LSCE carbonate system variables will be maintained, and further improvements with the aim of reducing model–observation misfit and improving the quantification of model uncertainty are on the way forward. Being at the core of the chain, model upgrades of CMEMS-LSCE-FFNN will be tackled first. At this time, SOCAT does not provide open-ocean data of CO₂ fugacity gridded at monthly, 0.25^∘ resolutions. Our ensemble-based approach thus draws on two SOCATv2022 data sources: a downscaled version of the 1^∘ open-ocean data and the 0.25^∘ coastal–ocean data (see Sect. 2.1). Open-ocean SOCAT datasets gridded at finer regular resolutions (if accessible) will be updated to gain more accuracy in our model fitting. Selections of data products for predictors needed for model input are equally important. For instance, the CMEMS SSS product used here results in a globally good reconstruction of total alkalinity (Table 4). However, the temporal variability in CMEMS SSS data does not match that in observations (Fig. A10), and this feature is retained in time series of total alkalinity (Fig. 7). Despite our best efforts in determining the overall product uncertainty in estimates of carbonate system variables, part of the input uncertainty is still not taken into account or only partially quantified due to a lack of time–space-varying uncertainty fields associated with predictor variables (e.g. SSH, Chl a, MLD, nutrient concentrations). Moreover, temporal sampling bias in pCO₂ and pH is likely to contribute to deviations between observations and model output (Fig. 5 and Table A3). The total analytical error uncertainty should be considered with great care during reconstruction and model output evaluation.

The CMEMS-LSCE approach leads as the first series of long-term reconstructions of pCO₂, pH, A_T, C_T, Ω_ca, and Ω_ar extending seamlessly from the global open ocean to coastal regions at monthly, 0.25^∘ resolutions. Future-use cases recommended for this high-resolution product include (1) estimation of monthly to interannual variations, long-term trends of carbonate system variables, and air–sea CO₂ exchanges at the surface layer from the local scale to large ocean basins; (2) analyses of interactions between these variables and the effects of other physical and biogeochemical factors on ocean acidification and changes in the marine carbonate system; (3) assessments of horizontal and temporal gradients of carbonate system variables in the coastal–open-ocean continuum; (4) evaluation or combination with other model- or observation-based products (e.g. Biogeochemistry Argo, Southern Ocean Carbon and Climate Observations and Modeling); and (5) improvements in coastal reconstructions based on observation system simulation experiments (e.g. with finer spatio-temporal model resolutions). The CMEMS-FFNN surface ocean carbon product at monthly, 0.25^∘ resolutions will be accessible through the CMEMS data portal (see Sect. 6).

Appendix A

Figure A1 (a) Ocean basins (https://github.com/RECCAP2-ocean/RECCAP2-shared-resources/tree/master/data/regions, last access: 11 July 2022): coastal mask (grey, approximately 400 km from the shoreline) and feature regions analysed in this study (cyan box, Table A1). (b) Location of time series stations recording in situ observations used in data evaluation (Table 2): blue stars for ocean acidification (Bates et al., 2014), black stars for A_T and C_T (Metzl and Lo Monaco, 1998; Coppola et al., 2021; Gattuso et al., 2023a), and other coloured scattered objects for pCO₂ and pH (Sutton et al., 2019). Asterisk (∗) marks the two stations with additional A_T and C_T observations (Olafsson et al., 2010) available for assessments.

Figure A2 Spatial distribution of total months in 1985–2021 containing SOCATv2022 gridded data: (a, c) 1^∘ data product (r100), (b, d) 0.25^∘ data product (r025). Open-ocean data (O) in each 0.25^∘ grid box are created by setting conservatively the open-ocean SOCATv2022 data at the 1^∘ grid box containing it. The coastal–ocean SOCATv2022 data (C) are assigned within 400 km from the shoreline (Fig. A1a).

Figure A3 Global maps of mean model–data difference (bias, a–d) and root-mean-square deviation (RMSD; e–h) between the reconstruction and SOCATv2022 pCO₂ (µ atm) over 1985–2021. (a, c, e, g) CMEM-LSCE-FFNN with a resolution of 1^∘ (r100); (b, d, f, h) CMEM-LSCE-FFNN with a resolution of 0.25^∘ (r025). Open-ocean data (O) in each 0.25^∘ grid box used for evaluation are created by setting conservatively the open-ocean SOCATv2022 data value at the 1^∘ grid box containing it. Coastal–ocean data (C) are extracted from each of the two SOCATv2022 gridded data products.

Figure A4 CMEMS-LSCE H⁺ over the global ocean at a spatial resolution of 0.25^∘. Temporal means of the model best estimate and 1σ uncertainty per grid cell over 1985–2021 are calculated by using Eq. (5).

Figure A5 CMEMS-LSCE Ω_ca over the global ocean at a spatial resolution of 0.25^∘. Temporal means of the model best estimate and 1σ uncertainty per grid cell over 1985–2021 are calculated by using Eq. (5).

Figure A6 CMEMS-LSCE Revelle factor (RF) over the global ocean at a spatial resolution of 0.25^∘. Temporal means of the model best estimate per grid cell over 1985–2021 are calculated by using Eq. (5).

Figure A7 Spatial distribution of CMEMS SSS and SST product uncertainty over the global ocean at a spatial resolution of 0.25^∘; 1σ uncertainty is computed per grid cell by using Eq. (5) over 1985–2021.

Figure A8 Median percentage of analysis error uncertainty against climatological mean of surface WOA18 nutrient data: phosphate (PO₄), nitrate (NO₃), and silicate (SiO₂).

Figure A9 Spatial distribution of R(σ,μ) [%] (Eq 8), i.e. the ratio of model uncertainty (σ) against model best estimate (μ).

Figure A10 Monthly time series of SSS and SST at BATS, DYFAMED, ESTOC, and HOT stations (Table 2 and Fig. A1b): CMEMS data (curve), 1σ uncertainty (envelope), and monthly average of observations (point). Means of CMEMS data and 1σ uncertainty (μ±σ) calculated over the observation time span are shown in brackets if accessible. Statistics include the number of months with observations (N), bias, RMSD, and r²; ${\mathit{\sigma }}_{\mathrm{SSS}}^{\mathrm{t}}$ (${\mathit{\sigma }}_{\mathrm{SST}}^{\mathrm{t}}$) stands for temporal standard deviation from monthly averages of SSS and SST observations. Temporal variations in SSS observations are poorly described in the CMEMS SSS time series (Table 1) used in CMEMS-LSCE reconstructions.

Figure A11 (a) Global surface seawater pH trend over the period 1985–2021, (b) 1σ uncertainties associated with trend estimates, (c) σ-to-μ ratio R(σ,μ) [%] (Eq. 8) between uncertainty estimates (b) and the best trend estimates (a), and (d) mask applied over the regions where $R(\mathit{\mu },\mathit{\sigma })>\mathrm{10}\mathit{\%}$.

Figure A12 (a) Global surface seawater Ω_ar trend over the period 1985–2021, (b) 1σ uncertainties associated with trend estimates, (c) σ-to-μ ratio R(σ,μ) [%] (Eq. 8) between uncertainty estimates (b) and the best trend estimates (a), and (d) mask applied over the regions where $R(\mathit{\mu },\mathit{\sigma })>\mathrm{20}\mathit{\%}$.

Figure A13 (a) Yearly global area-weighted mean of surface ocean saturation states with respect to calcite (Ω_ca): global means of the best estimate (μ, plain line) and uncertainty (σ, envelop) are computed with Eq. (7a). (b) Global trend maps of Ω_ca over 1985–2021: cross-hatching covers the regions with uncertainty of a trend estimate over 20 % of the trend value.

Figure A14 Linear trend estimates (per decade) learnt on 100-member ensembles (grey points) of yearly mean time series of pH and Ω_ar at different stations (Bates et al., 2014). μ±σ present linear slope and residual standard deviation. Black or blue lines stand for linear fits over the full or sub-period in 1985–2021 (see Table 5 for comparison).

Table A1 Information of feature regions analysed in this study (Fig. A1a – cyan boxes).

Table A2 Information of moored time series of surface–ocean pCO₂ and pH observations (Sutton et al., 2019).

Table A3 Statistics computed between CMEMS-LSCE datasets (0.25^∘) and time series of pCO₂ and pH measurements (Sutton et al., 2019) at open-ocean (O) and coastal stations (C): total numbers of monthly mean observations (N), temporal standard deviation of observations from their monthly averages (σ^t), RMSD (Eq. 10), and r² (Eq. 11). See Table A2 and Fig. A1b for stations' information and locations.

Table A4 Statistics computed between CMEMS-LSCE datasets (0.25^∘) and time series of A_T and C_T measurements (0–10 m depth): total numbers of monthly mean observations (N), temporal standard deviation of observations from their monthly averages (σ^t), RMSD (Eq. 10), and r² (Eq. 11). See Table 2 and Fig. A1b for stations' information and locations.