2.1 GRaCER algorithm

The GRaCER algorithm builds conceptually on a series of cluster-regression algorithms that have been successfully used for the interpolation and extrapolation of surface ocean pCO₂ (Sasse et al., 2013; Landschützer et al., 2013, 2016; Iida et al., 2015; Gregor et al., 2019). The main advantage of such a two-step approach is that the first clustering step organizes the variability regionally and temporally. This greatly enhances then the fidelity of the second step, i.e., the regression, as the size of the regression problem is reduced from the global domain to smaller, more homogeneous regions. A second advantage is that this clustering brings together regions with similar seasonality and similar co-variability with potential predictors, irrespective of the number of observations. The regression step explains the variability within each region over time and space dimensions, including interannual variability. Further, the clustering permits the regression to transfer information from spatially distant but geochemically similar regions, making the interpolation and extrapolation more robust in data-poor regions. The main innovation of the GRaCER algorithm relative to the previously used two-step approaches is its use of ensembles of cluster regressions, i.e., the generation of a whole series of clusters and corresponding regressions, which overcomes the boundary problems that are inherent in all two-step approaches.

For the clustering step, we use monthly climatological data of pCO₂ and TA and related parameters (Fig. 2a–c) to determine the main patterns of variability of the target variable and its co-variability with potential predictor variables. We opted for a clustering on climatological data rather than on the actual monthly data in order to clearly focus on the clustering step's role to isolate primarily regions with the same seasonal cycle. The alternative approach, i.e., to cluster on the monthly data, is also more prone to errors since the climatological distributions are better known than their month-to-month variations. Finally, clustering on monthly data would also take away signals from the regression step, which is actually better suited to capture the smaller level variations associated with the interannual variability and trends. The mini-batch K-means implementation in the Python scikit-learn package is used to perform the clustering due to its computational efficiency and scalability with large data sets (Pedregosa et al., 2011). A user-defined number of cluster centers are initiated, where cluster centers represent the mean of the points in a cluster. The K-means$++$ algorithm is used to initiate the cluster centers; it randomly selects the location of the first cluster center and then iteratively selects a best-guess location for the remaining cluster centers. Thereafter, the algorithm minimizes the distance between cluster centers and data points in the variable space. Once the clusters have been defined for the climatological domain, the co-located training data are assigned to the monthly clusters.

The regression is then performed individually for each of the clusters (Fig. 2d–f). The GRaCER method does not use a prescribed regression method – rather the appropriate algorithm for the particular use case is implemented. Importantly, the algorithm must be able to scale appropriately to the size of the problem. For example, the training data set for TA is $\mathrm{1}/\mathrm{20}$th of the size of the pCO₂ training data set; thus a more computationally expensive method can be used to predict TA.

The ensemble members are created by performing the cluster-regression step multiple times. Creating an ensemble is possible due to the fact that each clustering instance is slightly different (Fig. 2g–i). In practice, the spatial distribution of the clusters is similar; i.e., there is consistency in the typology of the clusters, particularly in regions where clusters are well defined, such as in the subtropical gyres and in the tropical eastern Pacific. However, there are regions that belong to different clusters; i.e., there is slight variance in the typology between ensemble members. The differences are due to the random initialization of the first cluster center in the K-means clustering step and the fact that clustering variables for some regions have weak gradients in spatial auto-correlation resulting in weak association with a cluster. In practice, this means that the locations of cluster boundaries vary between ensemble members; thus the ensemble mean does not have discrete boundaries (Fig. 2j).

2.2 Algorithm implementation

2.3 Data

Data are used to develop the two-step GRaCER model, i.e., clustering and regression, and to evaluate the estimates. Table 1 provides an overview of all data employed and the purposes for which they are used, and Table 2 shows the corresponding source of the data. We describe each data set by parameter and use.

Table 1Variables used as the clustering features and predictor variables for regression. Details about these data are given in the text. Note that clustering features are all resampled to monthly climatologies. ^a TA^map is the exception where a quasi-annual mean is used. The first column for regression shows the target variable. All machine learning models use the same variables to train and predict the final estimates with the exception of ^b SSS, where ungridded GLODAP data are used to train and SODA salinity is used to predict. All other references to SSS refer to the gridded SODA product.

Download Print Version | Download XLSX

Table 2Data sources used in this study.

Download Print Version | Download XLSX

2.4 Computation of DIC, pH, and Ω

The remaining parameters of the marine carbonate system, i.e., DIC, pH, and Ω, are computed using the Python version of CO2SYS (Humphreys et al., 2020) originally developed by Lewis et al. (1998). In addition to pCO₂ and TA, CO2SYS requires the input of sea-surface temperature, sea-surface salinity, pressure (assumed 0 dbar at the surface), ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$, and H₄SiO₄. We use the same data sources described in Sect. 2.3.2 and Table 1. Climatologies of H₄SiO₄ and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ are repeated for each year, thus assuming no interannual variability. The impact of this assumption is minimal (≪ 2 µmol kg⁻¹), as can be shown by varying these nutrients over seasonal cycle. The dissociation constants by Dickson et al. (1990) for $K_{{\mathrm{HSO}}_{\mathrm{4}}}$ and the total boron–salinity relationship by Uppström (1974) were used, as recommended by Orr et al. (2015) and Raimondi et al. (2019). For further details on the calculation and the full description of the marine carbonate system, see Dickson et al. (2007).

An important consideration in these calculations is the internal consistency of the marine carbonate system, i.e., the error due to uncertainties in the equations and coefficients that describe the marine carbonate system. Raimondi et al. (2019) pointed out that the pCO₂–TA pair has the lowest error in the calculation of pH (0.003 ± 0.008 pH units) using the dissociation constants by Mehrbach et al. (1973) as refitted by Dickson and Millero (1987). However, using the same pair and the same dissociation constants resulted in an estimate of Ω with respect to aragonite that is very different from that computed using the DIC–TA pair. But since Raimondi et al. (2019) lacked direct measurements of Ω, it remains unclear which pair is actually better for Ω. We cannot resolve this here but need to acknowledge that this inconsistency adds some additional uncertainty to our computed Ω values.

3.1 Sources of errors for TA and pCO₂

We identify three sources of errors that contribute to the total uncertainty for pCO₂ and TA, namely the uncertainty stemming from the measurement (M), representation (R), and prediction (P) errors. Assuming independence of the three error sources, the total uncertainty (E) for the TA and pCO₂ estimates can thus be expressed as the root of the squared sum of the uncertainties from the three error sources:

The measurement error reflects the combination of potential biases (systematic errors) from sampling and measurement as well as random errors associated with sampling and the imprecise nature of the measurement system. Since both TA and pCO₂ are being measured against certified reference materials and have undergone extensive secondary quality control, we assume that they have no systematic error, i.e., that their bias is zero. We also assume the sampling error to be small, so that the uncertainty M associated with the measurement error can be well approximated by the precision of the employed measurement methodology (Dickson et al., 2007).

The representation error, R, is a result of the fact that we develop our statistical model (GRaCER) on a grid that is in many places coarser in time and space than the typical scales of variability of TA and pCO₂. As a result, any given observation may not be representative for the $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ monthly grid cell used as a basis for our regression, leading to a bias in the estimated mean relative to the true spatial and temporal mean. This problem is particularly severe if the number of samples within any grid cell is low and the spatiotemporal variability is high. This is often the case. For example, more than 90 % of the data in the monthly gridded pCO₂ product are based on a single day of sampling within the month. The situation is more dire for TA, for which there are 10-fold fewer observations than pCO₂. Since we are lacking full knowledge of the spatial and temporal variability of TA and pCO₂, we cannot fully quantify the representation error. Instead, we approximate it using the few regions where we have sufficient observations or then using closely related parameters for which we have more observations. For simplicity, we make the assumption that the uncertainty associated with this error is, on a global average, normally distributed with a bias of 0.

The uncertainty associated with the prediction error, P, is determined by the test scores from the evaluation of the statistical model vis-à-vis the independent test data. The test scores describe the error incurred in the prediction of the subset of data that are not used in the training step. This error includes also the propagated uncertainty associated with the predictor variables.

We summarize these uncertainties with mean biases ($\frac{{\sum }_{i=\mathrm{1}}^N({\widehat{y}}_i-y_i)}{N}$) and root mean squared error (RMSE, $\sqrt{\frac{{\sum }_{i=\mathrm{1}}^N({\widehat{y}}_i-y_i{)}^{\mathrm{2}}}{N}}$), where y is the target value, $\widehat{y}$ is the predicted value, and N is the number of observations. We separate the coastal and open-ocean regions using the COastal Segmentation and related CATchments (COSCATs) mask (Laruelle et al., 2013) in order to reflect their very different levels of spatiotemporal variability.

Table 3Summary of the uncertainties of total alkalinity and pCO₂ from the different error sources (see Table 2) separately evaluated for the open ocean and for coastal regions (defined by the COSCATs regions, Laruelle et al., 2013). See text for details on how the different sources were quantified.

Download Print Version | Download XLSX

3.1.1 Uncertainty for total alkalinity

We adopt an uncertainty M associated with the measurement error of TA of ± 5 µmol kg⁻¹ based on the laboratory intercomparison by Bockmon and Dickson (2015). This is only half the accuracy of ± 10 µmol kg⁻¹ or 0.5 % reported by GLODAPv2 (Olsen et al., 2019). We consider this to be an overly conservative estimate, since Bockmon and Dickson (2015) pointed out that the majority of the laboratories involved in the round-robin exercise achieved an accuracy of better than ± 5 µmol kg⁻¹. We thus opted for this lower value that is more representative of the majority of the data.

Owing to the sparseness of the TA observations, we cannot estimate the uncertainty R associated with the representation error directly. Instead, we use the high correlation between TA and salinity. This permits us to determine the representation error for TA indirectly from an estimate of the representation error of sea-surface salinity. Concretely, we compare the test RMSE of TA predicted with GLODAPv2's in situ salinity with the RMSE of TA predicted with the satellite-based SODA salinity (see Table 1). Since the latter salinity is supposed to reflect the true time–space average over each grid cell, the difference between these two salinities is a direct estimate of the uncertainty associated with the representation error for salinity. Consequently, the difference in TA from these two estimates is an estimate of the uncertainty associated with the representation error for TA. The resulting estimates for the open and coastal ocean are summarized in Table 3.

The uncertainty P associated with the prediction error is based on the model's RMSE score calculated from test data and is listed in Table 3 and Fig. 3a. The global mean prediction error for the open ocean amounts to 13 µmol kg⁻¹, with some regional differences. The prediction error is more than twice this number in the coastal regions, i.e., 28 µmol kg⁻¹ (coastal regions are defined by the COSCATs regions, Laruelle et al., 2013). We find especially high prediction errors, for example, in the highly dynamic Amazon outflow region or the Gulf of Maine in the northwestern Atlantic. However, in such regions, one can expect that part of the high prediction error is actually stemming from a representation error, as we are not using directly co-measured variables when we train our regression model. In Fig. A2, we show the spatially and climatologically mapped test errors for TA in Fig. A2 using the GRaCER approach.

While the global bias of the TA product of OceanSODA-ETHZ is near zero (0.5 µmol kg⁻¹), confirming our assumption about the unbiased nature of our prediction error, this is not the case regionally. For example, OceanSODA-ETHZ tends to consistently overestimate TA in the southeastern Atlantic and underestimate TA in the southern Indian Ocean (Fig. 3b). A seasonal breakdown of the biases into DJF and JJA reveals that the winter period of each hemisphere has biases in the high latitudes, though the paucity of data means that we can place less weight on this finding.

Figure 3Test metrics for total alkalinity (a, b) and pCO₂ (c, d). The left-hand column (a, d) shows the root mean squared error (RMSE) compared with the target data. Similarly, the middle column (b, e) shows bias compared to the respective training data sets (GLODAP v2.2019 and SOCAT v2019). The right-hand column shows the zonally averaged biases for June, July, and August (JJA, blue) and December, January, and February (DJF, orange). A 2D spatial convolution was first applied to the $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ pixels to make regional patterns in the biases and RMSE clearer, and data were then aggregated into $\mathrm{4}{}^{\circ }\times \mathrm{4}{}^{\circ }$ pixels for clearer visualization.

A good check on the model prediction error is provided by comparing the estimated TA against independent observations. To this end, we use data from the Hawaii Ocean Time-series (HOT), the Bermuda Atlantic Time-series Study (BATS), and the Irminger station shown in Fig. 4a, b, e, f, i, and j and Table 4. For the period 1990–2018, the bias for BATS is 3 µmol kg⁻¹ and for HOT −2 µmol kg⁻¹, indicating that the method captures the overall structure and variability of TA well at these subtropical stations. Further, the mean seasonal cycle is relatively well represented at HOT and BATS, being within 1 standard deviation of the interannual variability when averaged as a climatology (Fig. 4b, f). However, the results are not as good for the Irminger station in the Atlantic high latitudes (∼ 65^∘ N), where OceanSODA-ETHZ has a large negative bias (−10 µmol kg⁻¹) when compared to TA computed from the observed pCO₂ and DIC. OceanSODA-ETHZ also overestimates the weak seasonal cycle of TA at the Irminger station, contributing to the large bias that is particularly strong from December to May. The RMSE at Irminger station is 15 µmol kg⁻¹, less than 5 µmol kg⁻¹ larger than the RMSE for HOT and BATS stations (10 µmol kg⁻¹ respectively) owing to the small interannual and seasonal amplitude at Irminger. The RMSEs are thus smaller than the mean prediction error (13 µmol kg⁻¹) at the subtropical stations yet exceed this mean estimate at the high-latitude station.

Figure 4A comparison of a subset of measurements from long-term observation stations (gray) with predicted total alkalinity (TA) (a, b, e, f, i, j) and partial pressure of CO₂ (pCO₂) (c, d, g ,h, k, l). The top row (a–d) shows data for the Bermuda Ocean Time-series Study (BATS), the middle row (e–h) shows data for the Hawaii Ocean Time-series (HOT), and the bottom row (i, l) shows data for the Irminger station. The narrow panels show the average of the seasonal climatology for the time series. The gray shading shows the standard deviation of the observations for the period 1990 to 2018, while the orange/blue lines show the average estimate. TA for the Irminger station is calculated from pCO₂ and DIC, and pCO₂ is calculated for BATS and HOT using DIC and TA, as described in Sect. 2.4.

Download

Table 4Comparison of training and independent data sources with various methods for the open-ocean region using the COSCATs coastal mask by Laruelle et al. (2013). GLODAP refers to the GLODAP v2 2019 data (Olsen et al., 2019), HOT refers to the Hawaii Ocean Time-series (Dore et al., 2009), BATS refers to Bermuda Atlantic Time-series Study (Bates and Peters, 2007), SOCAT is the 2019 version of the Surface Ocean Carbon Atlas (Bakker et al., 2016), LDEO is the Lamont-Doherty Earth Observatory pCO₂ data set for points not present in SOCAT (Takahashi et al., 2020), and SOCCOM is the pH measured by autonomous floats from the Southern Ocean Carbon and Climate Observations and Modeling project (Johnson et al., 2016). Statistical outliers were excluded in the calculation of LDEO RMSE. OS-ETHZ is the OceanSODA-ETHZ data from this study, NNGv2 is from Broullón et al. (2018), LIARv2 from Carter et al. (2018), CMEMS-FFNNv2 from Denvil-Sommer et al. (2019), and SOMFFN from Landschützer et al. (2016). NNGv2 and LIARv2 predictions are made with SODA salinity and OSTIA sea-surface temperature, resulting in different estimates to the original publications (Broullón et al., 2018; Carter et al., 2018). Note that the full data set is used for OceanSODA-ETHZ, unlike Table 3, which presents the errors for test years.

Download Print Version | Download XLSX

3.1.2 Uncertainty for pCO₂

For the uncertainty M associated with the measurement error of pCO₂, we adopt a value of ±2 µatm. This reflects the fact that 80 % of the data we have used from SOCAT (flags A and B) have a precision better than that number and an accuracy of similar magnitude. The remaining data we used (SOCAT flags C and D) have a measurement precision and accuracy of less than ±5 µatm.

We estimate the uncertainty R associated with the representation error of pCO₂ on the basis of a spatiotemporal gradient analysis. To this end, we compare the pCO₂ in our regular grid that has a resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ by 1 month, with the pCO₂ binned to a grid with twice this resolution, i.e., $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ by 15 d. In regions with high spatiotemporal coverage, the difference in the average of adjacent grid cells represents the potential change that can occur within the coarser $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ by a 1-month grid cell. The spatial and temporal gradients are calculated separately, and we take the average of these two elements. Using this analysis, we estimate a representation error of our pCO₂ estimates of 7 and 17 µatm for the open and coastal ocean respectively (Table 3).

From the RMSE of our test data, we estimate an uncertainty P associated with the prediction error of pCO₂ of 12 µatm for the open ocean and 28 µatm for the coastal ocean (Table 3). Within the open ocean (Fig. 3e), the eastern tropical Pacific has the highest RMSE, but this is also the region with the highest variance in the observations. The strong horizontal gradients in the region increase the errors, particularly at cluster boundaries. The high RMSE for the coastal region stems primarily from coastal Antarctica as well as some coastal regions in the higher latitudes of the Northern Hemisphere. The former is due to large uncertainties in pCO₂ during the summer months when retreating ice and ensuing rapid net primary production result in large gradients (Bakker et al., 2008). The climatologically mapped errors for pCO₂ are shown in Fig. A2b and d, created by mapping the test errors to the clusters and averaged over the ensemble.

The comparison between the regression estimated and observed pCO₂ reveals regional biases (Fig. 3e), despite the global bias being close to zero (−0.37 µatm). Some of the highest biases are found, again, in the eastern tropical Pacific, where strong horizontal gradients drive the observed juxtaposed biases. The large negative biases in winter in the Southern Ocean are likely driven by the paucity of data in this region (Gregor et al., 2019; Gray et al., 2018; Bushinsky et al., 2019).

The time series comparisons show that the seasonal cycle is well represented at BATS and HOT with r² scores of 0.89 and 0.82 respectively (Fig. 4d, h). Low biases (<2 µatm in absolute terms) further demonstrate that pCO₂ estimates are reliable in the subtropics. The seasonal cycle is also well captured at the Irminger station in the high latitudes, but a lower r² score and larger bias (−8.0 µmol kg⁻¹) allude to the dampened amplitude of the seasonal cycle particularly in the winter months (Fig. 4).

Figure 5Root mean squared error (RMSE – a, d) and biases (b, c, e, f) for dissolved inorganic carbon (DIC, a, b, c) and pH (d, e, f) compared with in situ GLODAP v2.2019 data. The two subplots in the rightmost column compare the zonally averaged bias for JJA and DJF. Data were processed for plotting as described in Sect. 4.

Figure 6A comparison of observations from long-term observation stations (gray) with predicted dissolved inorganic carbon (DIC). The top row (a, b) shows data for the Bermuda Ocean Time-series Study (BATS), the middle row (c, d) shows data for the Hawaii Ocean Time-series (HOT), and the bottom row (e, f) shows data for the Irminger station. The narrow panels on the right show the average of the seasonal climatology for the time series, where the gray shading shows the standard deviation of the observations for the period 1990 to 2018, while the green line shows the average estimate.

Download

3.2 Uncertainties of the calculated parameters

We determine the uncertainties of the calculated parameters in two ways. First, we propagate the uncertainties of pCO₂ and TA through pyCO2SYS (Orr et al., 2018; Humphreys et al., 2020) onto the computed parameters DIC, pH, and Ω. This yields an expected uncertainty that we refer to as a bottom-up total uncertainty estimate. Second, we obtain a top-down uncertainty estimate by comparing the calculated DIC and pH with independent measurements (see Table 4). We first describe the top-down estimates and then the bottom-up estimates. For the top-down estimate, we use the GLODAP DIC and pH data as our independent test of the method's performance. This is because these data are not used in any way in our estimation.

In the global mean, the computed DIC in OceanSODA-ETHZ has a very low bias compared with in situ GLODAP measurements (0.5 µmol kg⁻¹) (Table 4). Spatially, the DIC biases reveal a more nuanced picture, with large positive biases in the western equatorial Pacific and negative biases in the western equatorial Atlantic. The bias in the western equatorial Atlantic matches the negative bias in the pCO₂ in the same region; however, the source of the DIC bias in the Pacific is not clear from the pCO₂ and TA test data biases. The global mean of the top-down uncertainty estimate for DIC (16.3 µmol⁻¹ for the open ocean) is then perhaps a better reflection of the uncertainty as positive and negative values do not cancel each other out.

It is interesting to point out that the computed DIC in OceanSODA-ETHZ compares very favorably to directly estimated DIC products, such as that provided by NNGv2. Our uncertainty associated with the prediction error for DIC of 16.3 µmol⁻¹ is substantially better than theirs (23.1 µmol kg⁻¹) across all independent data.

The comparison of the DIC time series data (BATS, HOT, and Irminger stations) supports the findings of the global top-down estimates (Fig. 6). The biases are relatively low for HOT and BATS (−1.0 and 0.4 µmol kg⁻¹ respectively), but the bias is much larger at Irminger station (−6.9 µmol kg⁻¹), which is at a much higher latitude (64^∘ N) compared to HOT and BATS in the subtropics (<35^∘ N). However, this bias is not reflected in the zonal average of the seasonal biases (Fig. 5c). Similarly, the RMSE is also larger at Irminger station (14 µmol kg⁻¹) compared to HOT and BATS (∼ 9 µmol kg⁻¹) in the subtropics. Despite these differences in the top-down error, the same amount of variability is represented by the OceanSODA DIC for all three stations (∼ 0.76) owing to the larger seasonal cycle at Irminger station.

The pH comparison with the GLODAP pH measurements shows that OceanSODA-ETHZ has a negligible bias (0.001). As with DIC, regional biases in pH are larger than the global average, with the coastal and high-latitude oceans contributing significantly to the regional biases. The RMSE of pH with respect to GLODAP is also low (0.024) but is slightly outperformed by the RMSE of pH calculated with LIARv2 TA and FFNNv2 pCO₂ (made available in the FFNNv2 data set, 0.023). But given the uncertainty of GLODAP pH (0.01), the difference is negligible.

For the bottom-up estimate, we propagate the total uncertainty of pCO₂ and only the uncertainties associated with the measurement and prediction errors for TA. We do this to avoid including the representation error twice in the bottom-up estimate, as we hypothesize that the representation error of TA is largely accounted for by the representation error of pCO₂ and vice versa. We choose to use the representation error of pCO₂ rather than TA as the larger number of samples gives us greater confidence in the estimate. Moreover, we feel that the assumptions we make in the estimate of the TA representation error are larger than those for pCO₂, thus further justifying our choice in using only the representation error of pCO₂.

Figure 7A comparison of propagated uncertainties with independent uncertainties as an assertion of the validity of uncertainty estimates. The map (a) shows the separation between coastal and open ocean, panel (b) shows the error contributions in the open ocean, and panel (c) shows those in the coastal ocean. The total uncertainty has been broken into the three different components. Note that the values represented by the bar plots are not equivalent to values in Table 3 as the latter shows pCO₂ and TA total uncertainties for test data only, while the bar charts show total uncertainties for all data; further, the breakdown of the error contributions is proportional to the contribution of the sum of the squares (see Eq. 2).

The comparison of top-down vs. bottom-up uncertainty estimates for open and coastal oceans is shown in Fig. 7, with values of the total uncertainties also shown. The top-down and bottom-up uncertainty estimates for DIC are relatively accurate, with the estimates being within 5 % of each other in the open ocean and 9 % in the coastal ocean. However, the estimates are not as coherent for pH where the bottom-up error is 23 % smaller than the top-down uncertainty in the open ocean. The difference is even bigger in the coastal ocean where the bottom-up uncertainty is 31 % smaller than the top-down uncertainty.

4.1 Comparison with other climatologies

A spatial comparison between OceanSODA-ETHZ and existing products might reveal potential biases in our product if the bias is present in all comparisons. We compare TA against LIARv2 and NNGv2 (Carter et al., 2018; Broullón et al., 2018) by first taking the difference for each month and then calculating the average for these differences. The same approach is used to compare pCO₂ with SOMFFN and FFNNv2 (Landschützer et al., 2016; Denvil-Sommer et al., 2019).

The differences in TA between OceanSODA-ETHZ and NNGv2 and LIARv2 are on the same order of magnitude in the open ocean as the prediction error (13 µmol kg⁻¹) but are slightly larger on average for NNGv2 than for LIARv2 (Fig. 8a, b). There is also some agreement in the spatial pattern of the distribution, particularly in data-sparse parts of the Pacific Ocean and Indian Ocean. The larger differences may stem from the fact that both LIARv2 and NNGv2 are more constrained by spatial coordinates than GRaCER. NNGv2 uses latitude and longitude as predictors, while the LIARv2 approach interpolates the linear regression coefficients for every 5^∘ grid cell (Broullón et al., 2018; Carter et al., 2018). The divergence in data-sparse regions is thus not surprising. Although, it must be emphasized that this comparison serves more as a sanity check than as a ground-truthing exercise.

Figure 8A comparison of the mean differences between TA (a, b) and pCO₂ (c, d) for OceanSODA-ETHZ and other published methods: (a) LIARv2, (b) NNGv2, (c) CMEMS-FFNNv2, and (d) MPI-SOMFFN. The markers in the North Pacific show the locations used in the comparison of the climatology in Fig. 11.

The differences between OceanSODA-ETHZ pCO₂ and FFNNv2 and SOMFFN are smaller and not as spatially coherent compared to those for TA. The differences are marginally larger for the OceanSODA-FFNNv2 than for OceanSODA-SOMFFN (Fig. 8c, d). This is consistent with the data in Table 4, where the metrics for the SOMFFN are very similar to OceanSODA-ETHZ. The smaller difference between the latter should not come as a surprise as the GRaCER approach is built on the two-step cluster-regression approach of the SOMFFN, while the FFNNv2 approach includes spatial coordinates (Landschützer et al., 2016; Denvil-Sommer et al., 2019). In general, the dissimilarity between the differences is encouraging as it indicates that OceanSODA-ETHZ is not consistently biased relative to SOMFFN and FFNNv2.

Figure 9A basin–mean comparison of OceanSODA-ETHZ pCO₂ with four gap-filling methods: MPI-SOMFFN, JENA-MLS, CMEMS-FFNNv2, and CSIR-ML6. The thin lines show the differences compared to the individual methods, while the thick line shows the mean difference across the four methods. We do not show the Arctic Ocean as OceanSODA-ETHZ covers only 23 % of the region. The vertical dashed line in each figure marks the year 1990, where estimates prior to this period show biases in panels (b) and (c).

Download

We also show the temporal evolution of the basin–mean differences between OceanSODA-ETHZ pCO₂ and other gap-filling methods (Fig. 9). In the Atlantic (Fig. 9a), OceanSODA-ETHZ pCO₂ is <2 µatm lower than the mean of the other gap-filling methods for the period 1990 to 2008. Thereafter, the difference is <1 µatm. In the Indian Ocean, our pCO₂ estimates have a persistent negative difference of ∼2 µatm (Fig. 9c). The comparison in the Pacific (Fig. 9b) is the most consistent with the other methods, with a slight positive difference in the beginning of the period (pre-1990). The OceanSODA-ETHZ estimates of pCO₂ in the Southern Ocean (Fig. 9d) have a large positive difference prior to 1990 – up to 6 µatm for one of the ensemble members. This difference quickly diminishes and is near zero by 1990. There is also a negative difference later in the period (2004 to 2015); however, the ensemble spread over this period is large.

The comparison with other methods illustrates that while gap-filling methods are converging on a global scale, there are regional differences. Further, large differences in pCO₂ between methods prior to 1990 indicate high uncertainty for this period.

4.2 Seasonal climatologies

The climatological mean spatial distribution of TA, pCO₂, DIC, pH, and Ω, obtained by averaging the estimates from 1985 through 2018 reveals a very rich and diverse pattern of variability with commonalities and differences (Fig. 10). The climatological maps are accompanied by Hovmöller diagrams that show the zonal average of the seasonal cycle for each of the variables (Fig. 11). We also show climatological time series for each of the variables at high-latitude (55^∘ N, 170^∘ W), midlatitude (30^∘ N, 170^∘ W), and low-latitude (10^∘ N, 170^∘ W) locations.

Figure 10Mean maps of the GRaCER-based estimates for the period 1985–2018 for (a) total alkalinity and (c) pCO₂, as well as those of the computed variables, (b) dissolved inorganic carbon, (d) pH, (e) Ω_calc (saturation state with regard to calcite), and (f) Ω_arag (saturation state with regard to aragonite). The three black markers in each plot show the locations chosen for the seasonal analysis in Fig. 11f–j.

Total alkalinity shows the largest differences between basins, with the mean alkalinity being much higher in the saltier Atlantic than in the Pacific and Indian basins (Boutin et al., 2018) (Fig. 10a). The spatial variability of TA on a global scale exceeds the variability of the seasonal cycle (Fig. 11a, e). Seasonal variability of TA is on the order of 20 µmol kg⁻¹ at the chosen mid- and high-latitude locations, while the latitudinal gradient is as large as 150 µmol kg⁻¹. However, much of the TA seasonality is driven by seasonal changes in salinity due to precipitation and ice melt in the respective regions.

Dissolved inorganic carbon is more homogeneous across the basins but has a much larger meridional gradient than TA, amounting to more than 150 µmol kg⁻¹ (Fig. 10b). This meridional gradient is seasonally substantially more modified than is the case for TA, particularly in the high latitudes of the Northern Hemisphere where the seasonal cycle is as large as 100 µmol kg⁻¹ (Fig. 11b, g). The larger seasonal cycle in DIC is due to the carbon uptake by springtime phytoplankton blooms and stratification during the warmer seasons (Siegel et al., 2002). The magnitude of the springtime blooms is dampened by iron limitation in the Southern Ocean, visible by a smaller seasonal cycle amplitude (∼ 40 µmol kg⁻¹) (Watson et al., 2000; Tagliabue et al., 2017). However, the background DIC concentration in the Southern Ocean is much larger due to upwelling of DIC-rich circumpolar deep waters driven by the persistent westerlies south of 50^∘ S (Marshall and Speer, 2012).

The spatial distribution and seasonal cycles of pCO₂ and pH (Figs. 10c, d and 11c, d, h, i) are strongly negatively correlated due to the inverse stoichiometric relationship between dissolved aqueous CO₂ and [H⁺] (Dickson et al., 2007). The reduction in DIC is concomitant with the reduction in pCO₂ in the high latitudes due to the biological uptake of [CO₂] (Fig. 11g, h). However, this relationship does not hold true in the midlatitudes, where the slight decrease in DIC is contrasted by a relatively strong increase in pCO₂. This is due to the positive temperature dependence of pCO₂ (the opposite is true for pH) (Takahashi et al., 1993), which will be elaborated on in the discussion.

Figure 11Hovmöller plots (a–e) showing the zonally averaged seasonal climatologies for (a) total alkalinity, (b) dissolved inorganic carbon, (c) pCO₂, (d) pH, and (e) aragonite saturation state (Ω_Ar). The second row of panels (f)–(j) shows the corresponding variables for a high-latitude (55^∘ N, 180^∘ E, blue), midlatitude (30^∘ N, 180^∘ E, orange), and low-latitude (10^∘ N, 180^∘ E, green) locations. The units for panels (f)–(j) correspond with the units in panels (a)–(e).

Download

The spatial distribution of Ω (Figs. 10e, f and 11e, j) strongly reflects the concentration of the carbonate ion, which can be well approximated by the difference between TA and DIC (Sarmiento and Gruber, 2006). Given that the seasonal cycle of TA is much weaker than DIC, the latter dominates the seasonal cycle of Ω_Ar (Fig. 11e, j). This would also be true for Ω_Ca, which only differs from Ω_Ar in magnitude and not in distribution or seasonality (the latter is not shown).

4.3 Global- and basin-scale long-term trends

The OceanSODA-ETHZ data set can provide important novel constraints on the long-term trends in ocean acidification. We determine the long-term trends by a linear regression approach, restricting the period to 1990 through 2018, thus leaving out the 1980s, where the estimates are much more uncertain.

The global- and basin-scale trends for pCO₂ are remarkably similar (varying by only about 0.5 µatm per decade around the global mean of 16.5 µatm per decade) (Table 5. The ocean trends are also very consistent with the trends in atmospheric CO₂ (∼ 18.6 µatm per decade), reflecting the fact that the CO₂ system in the surface ocean is following the increase in the atmosphere very closely globally as well as in all ocean basins. The atmospheric trend is slightly steeper than those in the ocean. This is as expected since the atmospheric CO₂ increase forces the ocean, leading to a slight delay, a reflection of an increase in the air–sea disequilibrium over time (see, e.g., discussion in Matsumoto and Gruber, 2005). This growth in the air–sea disequilibrium is the driving force behind the increase in the oceanic sink strength for anthropogenic CO₂ over time (Gruber et al., 1996).

The basin-scale consistency holds true for pH as well (−0.016 units per decade) and Ω_ar (−0.07 units per decade), where global values are very similar to those averaged across the different ocean basins. In contrast, there is more spatial variability in the DIC trends, but the overall magnitude of about 9 µmol kg⁻¹ per decade largely represents the uptake of anthropogenic CO₂ from the atmosphere. The spatial differences in the DIC trends are also mirrored in the spatial variations in the TA trends, with regions with higher DIC trends having higher TA trends. This makes sense in terms of positive trends in TA increasing the buffering capacity, hence permitting DIC to grow faster for the same increase in seawater pCO₂. The trends in TA themselves are driven almost entirely by salinity with a basin-scale correlation of 0.99 (see Cheng et al., 2020).

Table 5Linear trends and their standard errors for OceanSODA-ETHZ variables for the period 1990 to 2018. All columns show increases per decade (decade⁻¹). All trends in the table are significant (P<0.05). We exclude the Arctic as the OceanSODA-ETHZ product only covers 23 % of this region and may thus give spurious trends. The ocean basins are defined by the map shown in Fig. A4.

Download Print Version | Download XLSX

5.1 Choosing the appropriate machine learning configuration

Here we consider two notable decisions that have a large impact on the final estimates: (1) the use of the ensemble approach and (2) the choice of regression algorithm. For details on the minor choices, see Sect. A3 in the Appendix.

As previously motivated, we opt for the cluster-regression approach that is able to generalize estimates in sparse regions due to information sharing within a cluster. However, cluster boundaries are often semi-discrete, resulting in artifactual boundaries. This makes the output of cluster-regression approaches less suitable for studies where gradients over short time periods or distances are assessed, e.g., for the detection of extreme events. Our approach removes these boundaries and improves the robustness of the estimates by eliminating, to a large extent, the sensitivity of the regression to the clustering algorithm.

The second major consideration is the choice of regression algorithm. Our choice of different algorithms for TA and pCO₂ may seem peculiar. However, this decision was informed by the nature of the problems. While regressing TA and pCO₂ are conceptually similar, the size of each data set and the distribution of training samples set them apart. When gridded to a monthly by $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ resolution, the number of data are ∼300 000 for SOCAT pCO₂ and ∼16 000 for GLODAP TA – i.e., nearly 20 times more data for the former. However, the strong linear correlation between TA and salinity (r=0.96) compensates for the poor sampling distribution – that is, as long as the regression method is able to extrapolate – a criterion that the support vector regression (SVR) method meets.

pCO₂ is poorly correlated to any proxy, suggesting that the regression problem is more complex, thus requiring a method that is appropriately non-linear. Gradient boosted decision trees (GBDTs) and feed-forward neural networks (FFNNs) meet this criterion. But why not just use one of these approaches? Work by Gregor et al. (2019) found that an ensemble of methods (SVR, GBDTs, and FFNN) outperformed each individual member. And while SVR performed well in Gregor et al. (2019), the method does not scale to larger problems, which GBDTs and FFNNv2 are capable of, leading to our choice of the latter two. One critique of GBDTs in this application may be that, being a tree-based method, it is not able to extrapolate. However, we feel that the cluster-regression approach combined with the large number of training data for pCO₂ compensates for this shortcoming. GBDTs provide also provide useful diagnostics, such as feature importances, that, when combined with the GRaCER approach, provide useful information about the spatial and seasonal importance of the proxies.

5.2 Why are pH uncertainties less well constrained?

One of the novel contributions of this study is that we are able to assert the validity of our results by comparing the bottom-up (propagated) with the top-down uncertainties (in situ comparisons). Using this approach, we show that the uncertainty estimates of DIC are remarkably well constrained, with the top-down estimate being within 5 % of the bottom-up uncertainty estimate for the open ocean (Fig. 7b). The same statistics for pH yields a 23 % difference between the two budgets. The question is thus, why are pH bottom-up and top-down uncertainty estimates much less consistent?

To assess this problem from the bottom-up perspective, we need to consider the uncertainties of pCO₂ and TA. But given that the DIC uncertainty budgets are well constrained, we can, with some certainty, rule out the bottom-up uncertainty estimate as the source for the larger mismatch in pH.

The source of the mismatch must thus be driven primarily by uncertainties in the top-down perspective, where it may be that the representation error of pH is larger than for DIC. We rule out the measurement error as a contributor to the mismatch, as the bias of the measurements (provided accurate calibration to reference samples) should be normally distributed around zero. Thus, the representation error is the most likely candidate, due to the temperature- and pressure-sensitive nature of pH compared to the conservative nature of DIC with respect to the same variables (Dickson et al., 2007). This is important given that our surface pH data from GLODAP can be as deep as 30 m, while pCO₂ is typically measured at the intake depth of most ships, which is typically at 5 m (Bakker et al., 2016). A basic sensitivity study of the variability of pH in the surface layer shows that the median standard deviation of surface pH at a single station is 0.004 units (even when normalized to a standard temperature and depth). The same approach for DIC yields a standard deviation of 2 µmol kg⁻¹. This suggests that the vertical mismatch error of GLODAP alone already explains a good deal of the discrepancy in the uncertainty budgeting. One might reduce this uncertainty by placing tighter constraints on the definition of surface pH, limiting the depths between 3 and 15 m, for example, but at the loss of valuable test samples. The representation error (described in Sect. 3.1) could further explain the remaining disparity in the pH error budget.

5.3 Can we reduce the total uncertainty?

The last two decades have seen major improvements in the accuracy and precision of the TA and pCO₂ measurements, leading to substantial reductions in the measurement errors. The introduction of certified reference materials for TA and a standardized approach for measuring pCO₂ with reference gases means that the uncertainties associated with the measurements are low.

In contrast, the prediction uncertainty is the largest contributor to the total uncertainty for both DIC and TA, suggesting that this could be a fruitful avenue to pursue. However, current literature suggests that this is unlikely. Gregor et al. (2019) showed that, within a selection of six gap-filling methods, all achieved similar accuracy scores when compared with independent data, upon which the authors suggested that we have hit a wall.

This leaves the representation error, which contributes a moderate fraction to the total pCO₂ and TA uncertainties in the open ocean and even less in the coastal ocean. A back-of-the-envelope calculation shows that increasing the resolution of pCO₂ 4-fold (from monthly by 1^∘ × 1^∘ to 8-daily by 0.25^∘ × 0.25^∘) could decrease the representation error by 2.5 µatm (35 %) for the open ocean and 3.2 µatm for the coastal ocean (20 %). This is perhaps not as much as expected, but these small gains are larger than those that are currently being made with regard to the prediction errors (Gregor et al., 2019). This is not applicable for TA, where ungridded values are already used to train GRaCER and where the selection of the predictor variables is likely more important.

Why are these gains smaller than hoped? It may be that our gradient approach for calculating the representation error breaks down as the resolution increases due to the decreasing number of adjacent grid points. This is hardly surprising considering that 78 % of grid cells in the SOCAT v2019 monthly gridded product are represented by sampling on a single day that falls within that period (Bakker et al., 2016). Another possibility is that decreasing both the spatial and temporal resolution exposes the sharp mesoscale gradients that are otherwise averaged over at larger resolutions (Resplandy et al., 2009; Monteiro et al., 2015). With the currently available data, it seems as if the reduction to be made in the total uncertainty by increasing the resolution of the prediction data will be small.

Figure 12Attribution of DIC, TA, and temperature to the seasonal cycle of pCO₂. The top two figures show the seasonal cycle of pCO₂ and the drivers thereof for (a) high latitudes, (b) midlatitudes, and (c) low latitudes. These locations are shown with the markers in the Pacific Ocean in panel (c). The map (d) shows the dominant driver of the seasonal cycle for each region calculated as the value with the maximum seasonal amplitude.

5.4 Regional sensitivity of pCO₂ to driver variables

Here we demonstrate one of the possible ways in which the OceanSODA-ETHZ data can be used to gain further insight into the marine carbonate system.

We decompose and attribute the mean seasonal cycle variability of pCO₂ to its drivers, namely TA, DIC, temperature, and salinity. Past studies using observation based products have been limited to a simpler thermal/non-thermal decomposition of pCO₂ due to the lack of DIC and TA (Takahashi et al., 2002; Landschützer et al., 2015). This is thus the first time that a full decomposition of pCO₂ has been applied to observation-based data for a global domain. The decomposition is performed with pyCO2SYS by keeping all but one of the drivers constant (to the average) and assessing the influence on pCO₂ (Humphreys et al., 2020). This is similar to previous studies that applied an analogous decomposition to ocean simulation output (Lovenduski et al., 2007; Fassbender et al., 2018; Gallego et al., 2018). An important note to make, considering that we intend that OceanSODA-ETHZ is used for ocean acidification studies, is that the decomposition of pH would result in virtually the same contribution of the drivers (Dickson et al., 2007).

The seasonal amplitude of pCO₂ is driven predominantly by changes in DIC in the high and equatorial latitudes and by temperature in the midlatitudes (Fig. 12d). At first glance, these results may seem similar to those that the simpler thermal decomposition might result in, requiring only temperature. However, in regions where the seasonal amplitude of pCO₂ is smaller, the importance of TA becomes more apparent. For example, at the midlatitude station (Fig. 12b), TA and DIC synergistically act to dampen the impact of temperature on pCO₂. Conversely, at the equatorial station (Fig. 12c) the effect of TA on pCO₂ opposes that of the more dominant DIC. Further, there are regions in the tropics where TA is the dominant driver due to the weak seasonal cycle of both temperature and DIC.

5.5 Recommendations for use

In order to use the OceanSODA-ETHZ product in an optimal manner, it is important to be aware of its strengths and weaknesses.

The primary use of the OceanSODA-ETHZ data set is to determine and assess the seasonality, the interannual variations, and trends of ocean acidification thanks to its containing all relevant parameters of the marine carbonate system (Landschützer et al., 2015, 2016; Gregor et al., 2018; Keppler and Landschützer, 2019). However, users of the OceanSODA-ETHZ product should be aware of the fact that data prior to the 1990s should be treated with care due to the paucity of SOCAT pCO₂ training data during this period (Rödenbeck et al., 2015; Watson et al., 2020). This was recently demonstrated by Watson et al. (2020), who used an ensemble of various regression approaches to show that the spread of pCO₂ estimates prior to the 1990s is large due to the paucity of data. Similarly, Gregor et al. (2019) showed that pCO₂ estimates prior to 1990 tend to have a slightly positive bias.

The product is also very well suited for assessing models. Thanks to the spatially resolved estimates of uncertainty for TA and pCO₂ (Fig. A2), one cannot only assess the model–observation mismatches, but also weigh them with the appropriate uncertainties.

A strength of the OceanSODA-ETHZ product is that it extends further into the coastal margin than most previous studies (Iida et al., 2015; Landschützer et al., 2016; Denvil-Sommer et al., 2019). This is achieved (i) by including coastal observations during the training and (ii) by using a larger number of clusters compared to other clustering approaches (Landschützer et al., 2016; Watson et al., 2020). This permits us to better separate open-ocean and coastal variability through the inclusion of suitable variables in the clustering step (e.g., Chl a for pCO₂, and see Fig. A5 to see a representation of cluster boundaries). This gives us confidence in the coastal estimates of the climatological seasonal cycle. Our product is therefore comparable to that of Landschützer et al. (2020), who blended separate coastal and open-ocean pCO₂ products into a single climatological product with a monthly resolution (Landschützer et al., 2016; Laruelle et al., 2017).

The total uncertainties of our estimates in the coastal ocean are considerably larger compared to the open-ocean estimates (Fig. 7). This reflects the much higher spatiotemporal variability of the physical and chemical environment in the coastal ocean, leading to much higher variations in the marine carbonate system (Laruelle et al., 2017). Since our predictor variables are only partially reflecting this variability, a large portion of the high total uncertainty is due to a high representation error (Table 3). Increasing the resolution of the products may improve the total uncertainty of coastal estimates as done by Laruelle et al. (2017). Until we arrive at this point, the OceanSODA-ETHZ data should be used with care in the coastal ocean. Further, we recommend that researchers interested in the investigation of interannual variability and trends in the coastal ocean using the OceanSODA-ETHZ product should also look a the underlying in situ data to gain a better understanding of the variability, trends, and uncertainties for the coastal region of interest.