A global monthly climatology of oceanic total dissolved inorganic carbon: a neural network approach

Anthropogenic emissions of CO₂ to the atmosphere have modified the carbon cycle for more than 2 centuries. As the ocean stores most of the carbon on our planet, there is an important task in unraveling the natural and anthropogenic processes that drive the carbon cycle at different spatial and temporal scales. We contribute to this by designing a global monthly climatology of total dissolved inorganic carbon (TCO₂), which offers a robust basis in carbon cycle modeling but also for other studies related to this cycle. A feedforward neural network (dubbed NNGv2LDEO) was configured to extract from the Global Ocean Data Analysis Project version 2.2019 (GLODAPv2.2019) and the Lamont–Doherty Earth Observatory (LDEO) datasets the relations between TCO₂ and a set of variables related to the former's variability. The global root mean square error (RMSE) of mapping TCO₂ is relatively low for the two datasets (GLODAPv2.2019: 7.2 µmol kg⁻¹; LDEO: 11.4 µmol kg⁻¹) and also for independent data, suggesting that the network does not overfit possible errors in data. The ability of NNGv2LDEO to capture the monthly variability of TCO₂ was testified through the good reproduction of the seasonal cycle in 10 time series stations spread over different regions of the ocean (RMSE: 3.6 to 13.2 µmol kg⁻¹). The climatology was obtained by passing through NNGv2LDEO the monthly climatological fields of temperature, salinity, and oxygen from the World Ocean Atlas 2013 and phosphate, nitrate, and silicate computed from a neural network fed with the previous fields. The resolution is $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ in the horizontal, 102 depth levels (0–5500 m), and monthly (0–1500 m) to annual (1550–5500 m) temporal resolution, and it is centered around the year 1995. The uncertainty of the climatology is low when compared with climatological values derived from measured TCO₂ in the largest time series stations. Furthermore, a computed climatology of partial pressure of CO₂ (pCO₂) from a previous climatology of total alkalinity and the present one of TCO₂ supports the robustness of this product through the good correlation with a widely used pCO₂ climatology (Landschützer et al., 2017). Our TCO₂ climatology is distributed through the data repository of the Spanish National Research Council (CSIC; https://doi.org/10.20350/digitalCSIC/10551, Broullón et al., 2020).

1 Introduction

The ocean is the major carbon reservoir of the Earth. Most of this carbon occurs as dissolved inorganic carbon (TCO₂, also known as DIC or C_T; Ciais et al., 2013; Tanhua et al., 2013). Three species make up TCO₂: dissolved CO₂ (generally considered to be the sum of the dissolved CO₂ itself, CO₂(aq), and carbonic acid, H₂CO₃), bicarbonate ion (${\mathrm{HCO}}_{\mathrm{3}}^{-}$), and carbonate ion (${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$). The relative concentrations of these species with respect to each other determine the seawater pH (Zeebe and Wolf-Gladrow, 2001). The seawater CO₂ chemistry system can be represented as a set of chemical equilibria reactions that describe the speciation of the various ions of TCO₂ as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{g}\right)\rightleftharpoons {\mathrm{CO}}_{\mathrm{2}}\left(\text{aq}\right)\\ {\displaystyle }& {\displaystyle }{\mathrm{CO}}_{\mathrm{2}}\left(\text{aq}\right)+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\rightleftharpoons {\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}\\ {\displaystyle }& {\displaystyle }{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}\rightleftharpoons {\mathrm{H}}^{+}+{\mathrm{HCO}}_{\mathrm{3}}^{-}\\ {\displaystyle }& {\displaystyle }{\mathrm{HCO}}_{\mathrm{3}}^{-}\rightleftharpoons {\mathrm{H}}^{+}+{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}.\end{array}$$

Since the Industrial Revolution, the concentration of TCO₂ in the global ocean has increased, generally to a certain depth level (depending on the particular processes in each ocean area) due to the entry of CO₂ into the seawater from the atmosphere (Sarmiento and Gruber, 2002; Doney et al., 2009; Vázquez-Rodríguez et al., 2009; Bates et al., 2012; Sallée et al., 2012; Khatiwala et al., 2013). The uptake is driven by the increasing partial pressure of CO₂ (pCO₂) in the atmosphere relative to the ocean, generated by the anthropogenic emissions of CO₂ that cause an annual net flux of this gas into the ocean (Le Quéré et al., 2018). Accompanying the change in TCO₂, the pH and carbonate ion concentration have been declining because of the anthropogenic process previously mentioned, these changes being reflected in the proportions of the chemical species of TCO₂ (Kleypas and Langdon, 2000; Orr et al., 2005). These changes in seawater chemistry framed in the ocean acidification process can negatively influence various processes involving marine organisms such as calcification, growth, and survival (Orr et al., 2005; Fabry et al., 2008; Hendriks et al., 2010; Hoegh-Guldberg and Bruno, 2010; Kroeker et al., 2013).

In addition to the secular trends driven by the uptake of anthropogenic CO₂, ocean TCO₂ varies both temporally and spatially as a consequence of several natural processes. This variability may reach values of 15 % of the mean TCO₂ value in the ocean (Lee et al., 2000). The processes that increase TCO₂ are net flux of CO₂ from the atmosphere to the ocean, organic matter remineralization, and the dissolution of calcium carbonate (CaCO₃). The processes that reduce TCO₂ are net flux of CO₂ from the ocean to the atmosphere, primary production, and calcification. Advection and mixing also influence the variability of TCO₂ in these two ways (Sabine et al., 2002). In the surface ocean, the main variables influencing the variability of TCO₂ are temperature and salinity (Weiss et al., 1982; Lee et al., 2000; Wu et al., 2019), through the modification of the solubility of CO₂, affecting the seawater pCO₂ (which is almost instantaneous) and thus the air–sea CO₂ flux, which eventually drives the change in TCO₂ over time. Nutrients and oxygen can also reflect the processes that modify the concentration of TCO₂ through their consumption and release, like during the cycling of organic matter (Körtzinger et al., 2001; Bauer et al., 2013). From products generated with measured data (Key et al., 2004; Takahashi et al., 2014; Lauvset et al., 2016) and in modeling studies (e.g., Doi et al., 2015), it is known that the global surface distribution of TCO₂ follows a zonal gradient: there is a reduction in its concentration from the poles to the Equator, reflecting the processes that control its variability. Key et al. (2004) emphasize that this distribution is associated with the distribution pattern of nutrients. Recently, Wu et al. (2019) found that the distribution of surface-salinity-normalized TCO₂ (nDIC) has two main drivers: temperature and upwelling. At depth, the variation shown in almost any measured profile of TCO₂ mainly reflects the remineralization of organic matter and, to a lesser extent, the dissolution of CaCO₃ (Millero, 2007), resulting in an increase in TCO₂ from the surface to intermediate depths.

Understanding the distribution and variability of TCO₂ in the ocean and its secular trends driven by anthropogenic carbon uptake is needed to assess the magnitude and possible impacts of ocean acidification. It is also necessary for the evaluation of numerical models that include the carbon cycle and their estimates of past, current, and future ocean carbon cycle behavior (e.g., Yool et al., 2013; Aumont et al., 2015; Butenschön et al., 2016; Le Quéré et al., 2016;) Goris et al., 2018). Seasonality of TCO₂ and the horizontal and vertical variability underscore the necessity to design a climatology with both monthly and spatial resolutions according to the processes that influence this variable on a global scale. The existing climatologies of TCO₂ do not include all these characteristics collected together. Key et al. (2004) and Lauvset et al. (2016) built an annual climatology at 33 depth levels using interpolation techniques with data from the Global Ocean Data Analysis Project version 1 (GLODAPv1; Key et al., 2004) and version 2 (GLODAPv2; Key et al., 2015; Olsen et al., 2016), respectively. Takahashi et al. (2014) published a monthly climatology for the surface ocean computed from climatologies of pCO₂ and total alkalinity (A_T). Other studies used the covariability between TCO₂ and other more commonly measured variables discussed above for mapping and gap-filling via empirical regressions and neural networks. Lee et al. (2000) used temperature and nitrate to compute surface nDIC with an area-weighted error of ±7 µmol kg⁻¹. Sauzède et al. (2017) and Bittig et al. (2018) trained neural networks with GLODAPv2 data to compute TCO₂ over the depth range 0–8000 m with an accuracy of ±9 and ±7.1 µmol kg⁻¹, respectively. The input variables used in those studies were location, pressure, temperature, salinity, dissolved oxygen, and time.

In the present study, we introduce the use of neural networks for going one step further in the design of a climatology. We have generated a climatology of TCO₂ with a resolution consistent with that of the climatology of A_T of Broullón et al. (2019): horizontal resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$, 102 depth levels between 0 and 5500 m, and a monthly (0–1500 m) and annual (1550–5500 m) temporal resolution. The availability of global databases containing variables of the seawater CO₂ system with more and more data – e.g., GLODAPv2.2019; the Lamont–Doherty Earth Observatory (LDEO) database (Takahashi et al., 2017); the Surface Ocean CO₂ Atlas (SOCAT; Bakker et al., 2016) – and the great ability of the neural networks to interpolate as shown in other climatological studies about CO₂ system variables (Landschützer et al., 2014; Broullón et al., 2019) show the appropriateness of this approach for generating a global monthly climatology covering more than the surface ocean.

2 Methodology

2.1 Neural network design

A feedforward neural network was configured to compute TCO₂ in the global ocean and to create a global climatology based on the good results previously obtained with this method in similar studies (e.g., Broullón et al., 2019). Briefly, a neural network of this type (Fig. S1a in the Supplement) is used to extract relationships between a set of input variables and a target one through a training process. At this stage, the inputs are passed through different parallel layers composed of a tunable number of neurons to reach values as closest as possible to the target ones (Fig. S1a). Initially, all inputs enter each neuron of the first layer, where they are multiplied by different weights depending on the neuron they go to. Inside the neurons (Fig. S1b), the results of the previous operation are summed, and a bias is added. The obtained value inside each neuron is passed through an activation function, which yields an output. The outputs of each neuron in each layer go to the following layer undergoing the same process described to this point. In the last layer, which is composed of one neuron, a unique value for the target variable is calculated for each pair consisting of inputs and a target. This value is compared to the desired one, and the difference between both values is backpropagated through the entire network in order to adjust the weights and biases and to start the processes again and reach an accurate output value after multiple iterations. A complete description of the most common algorithms used to backpropagate and minimize the errors can be found in Rumelhart et al. (1986), Levenberg (1944), and Marquardt (1963).

The method used here is equivalent to that fully described by Broullón et al. (2019) for A_T. In addition to the target variable (TCO₂ instead of A_T), the main changes in the present study compared to that of Broullón et al. (2019) are the inclusion of the input variable “year”, accounting for the anthropogenic increase in the TCO₂ pool, and the use of the pCO₂ database from the LDEO (Takahashi et al., 2017) in addition to the extended GLODAPv2.2019 (Olsen et al., 2019) to enable more robust TCO₂ estimates in the surface ocean. Similar to Broullón et al. (2019), the neural networks were trained using the Levenberg–Marquardt method (Levenberg, 1944; Marquardt, 1963) through the trainlm function (detailed in Beale et al., 2018) in MATLAB. The splitting of the database used in the present study (see Sect. 2.2) into the sets needed for training and testing the network is depicted in Fig. 1. The data were randomly associated with each dataset to capture (training) and evaluate (test) all possible variability. The input variables are temperature, salinity, phosphate, nitrate, silicate, oxygen, sample position, and year (Fig. S1a). The number of neurons tested in the unique hidden layer to find the best neural network was 16, 32, 64, 128, and 256. Ten networks were trained for each number of neurons. The criteria for selecting the final number of neurons are based on a trade-off between the root mean square error (RMSE; between the measured TCO₂ and that estimated by the neural network) on the one hand and the generalization of the network (to prevent overfitting, maintaining a similar error in the training and in the test sets) on the other hand. Furthermore, an additional criterion based on the influence of each input variable on the TCO₂ extracted with the connection weight approach (Olden and Jackson, 2002) was followed to ensure that biogeochemical input variables have a larger influence on the TCO₂ estimates than the input variables related to sample position for selecting a proper network. The influence of each input variable on the computed TCO₂ was obtained from Eq. (1):

$$\begin{array}{}\text{(1)}& C_i=\sum _{k=\mathrm{1}}^H{w}_{ik}\cdot w_k,\end{array}$$

where C_i is the relative importance of the input variable i, H the number of neurons in the hidden layer, w_ik the weight of the connection between the variable i and the neuron k of the hidden layer, and w_k the weight of the connection between the neuron k of the hidden layer and output layer.

Figure 1 Division of the complete database in the datasets needed to train the neural network. The percentages in each level are relative to the number of data in the previous one. Data in the datasets of the first level are always the same for each network. Data in the sets of the second level are randomly associated with each set for each network to find the best network weights because of the different starting points in the error weight space of the training process (see also Broullón et al., 2019).

2.2 Data

We included the LDEO database version 2016 (Takahashi et al., 2017; https://www.nodc.noaa.gov/ocads/oceans/LDEO_Underway_Database, last access: 13 November 2017) because it contains significantly more data in the surface layer than GLODAPv2.2019. Since the higher variability in the surface layer may lead to high errors in modeling variables of the seawater CO₂ system (e.g., Carter et al., 2018; Bittig et al., 2018; Broullón et al., 2019), including the LDEO database should force the network to reach a more robust fit. The idea is that these additional data probably have more different relationships between input variables and TCO₂ to help the neural network to adequately capture spatiotemporal variability. The pCO₂, temperature, and salinity data from the LDEO were monthly averaged for each year in a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid. The points where the standard deviation (SD) of the averaged pCO₂, temperature, and salinity were greater than ±20 µatm, 1.5 ^∘C, and 0.5, respectively, were discarded since the objective is to capture the monthly variability, and therefore an extremely high submonthly variability could lead to errors. To obtain TCO₂ values from the LDEO data, an additional variable of the CO₂ system is necessary, for which we take A_T computed using the neural network NNGv2 of Broullón et al. (2019). The input variables required by NNGv2 were obtained from (1) temperature and salinity from the LDEO; (2) filtered oxygen from the World Ocean Atlas version 2013 (WOA13; see Broullón et al., 2019); and (3) phosphate, nitrate, and silicate computed with CANYON-B (Bittig et al., 2018) using the previous variables as inputs. Finally, TCO₂ was calculated from this A_T and the averaged pCO₂ using the MATLAB version of the CO_2SYS program (van Heuven et al., 2011); we used the dissociation constants of Mehrbach et al. (1973; as refit by Dickson and Millero, 1987) and the borate dissociation constant of Dickson (1990). Note that we used this software and set of constants for all seawater CO₂ chemistry calculations in the present study. The TCO₂ calculated this way and the associated input variables were used as a part of the training and testing data for the neural networks created here. The final number of data points derived from the LDEO was 54572.

To represent interior ocean conditions, the GLODAPv2.2019 database (Olsen et al., 2019) was added to the LDEO dataset for training and testing the neural network. Only samples which had data for all input variables and TCO₂ were used. This database was included in two ways: (1) only samples where all variables passed the second quality control (n=287 953; Olsen et al., 2016; Olsen et al., 2019; hereafter abbreviated Gv2QC) and (2) all samples (n=321 647; hereafter abbreviated Gv2). Therefore, two neural network options were trained and tested: NNGv2QCLDEO and NNGv2LDEO, respectively.

2.3 Comparison of methods

We compared our method with CANYON-B of Bittig et al. (2018), where TCO₂ values were also computed from multiple input variables. Both methods are based on neural networks but with certain differences as summarized in Table 1.

Table 1 Differences between the methods used in the present study and in CANYON-B (Bittig et al., 2018).

An error analysis was carried out in the same areas for which this was done by Broullón et al. (2019) for A_T and in several depth ranges (0–50, 50–200, 200–500, and 500–1000 m and 1000 m–bottom) for the two methods (our method and CANYON-B) and for the two datasets (Gv2QC and LDEO). The Gv2QC database was analyzed in this section instead of Gv2 because in the designing of CANYON-B only quality-controlled data were included. The analysis of CANYON-B using the LDEO dataset is useful to evaluate the validity of the approach followed by converting pCO₂ to TCO₂ since CANYON-B has not been trained with this dataset.

Computed pCO₂ from A_T and TCO₂ derived from neural networks was also evaluated in the LDEO dataset to assess the adequacy of including this dataset in our approach and to assess the ability of NNGv2 of Broullón et al. (2019) and the present TCO₂ neural network to compute other variables of the seawater CO₂ system. Furthermore, we compared the magnitude of the errors with the ones obtained by Landschützer et al. (2014), in which pCO₂ is computed directly with a neural network, to evaluate the accuracy of our computed pCO₂.

2.4 Validation

In addition to the ability to compute TCO₂ using the Gv2 and LDEO test sets, the neural network has been tested using independent data from 10 ocean time series, located in different regions of the world ocean (data were obtained from https://www.nodc.noaa.gov/ocads/oceans/time_series_moorings.html, last access: 4 June 2019): Hawaii Ocean Time-series (HOT ALOHA and HOT ALOHA SURFACE; Dore et al., 2009), Bermuda Atlantic Time-series Study (BATS; Bates et al., 2012), European Station for Time-series in the Ocean at the Canary Islands (ESTOC; González-Dávila et al., 2010), Iceland Sea time series (ICELAND; Olafsson et al., 2010), Irminger Sea time series (IRMINGER; Olafsson et al., 2010), Kyodo North Pacific Ocean time series (KNOT; Wakita et al., 2010), K2 (Wakita et al., 2010), Ocean Weather Station Mike (OWS; Gislefoss et al., 1998), and Kerguelen Islands in the Indian sector of the Southern Ocean (KERFIX; Jeandel et al., 1998). CANYON-B was also used to compute TCO₂ in the time series to show the differences between that method and ours. The TCO₂ values were obtained by feeding the neural networks with the measured values of the input variables at each time series. The data from these time series allow us to test the ability of the neural network to reconstruct not only the seasonal variability of TCO₂ at the various locations and depths sampled but also its long-term trends. For the trend analyses, the measured and estimated TCO₂ values were deseasonalized following Bates et al. (2014).

As an additional test, the measured pCO₂ or the pCO₂ calculated from measured TCO₂ and A_T at the time series stations was compared with pCO₂ calculated from the neural-network-generated values of A_T and TCO₂. This provides insight into the combined performance of the NNGv2 of Broullón et al. (2019) and the neural network designed in the present study. Furthermore, we compared the magnitude of the errors to that obtained by Landschützer et al. (2014) for some of the time series.

2.5 Climatology of TCO₂

We used the selected network, based on the results of the analyses described above, to construct a climatology of TCO₂. Climatologies of the input variables were passed through the network to obtain the climatological fields of TCO₂. The spatiotemporal resolution of the product is determined by that of the climatologies used as inputs: $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ horizontal resolution, 102 upper depth levels of the WOA13, and monthly (for 0–1500 m depth) to annual (for 1550–5500 m depth) temporal resolution. Temperature and salinity climatologies were obtained from objectively analyzed WOA13 fields (Locarnini et al., 2013; Zweng et al., 2013; https://www.nodc.noaa.gov/OC5/woa13/woa13data.html, last access: 6 February 2017). Oxygen, phosphate, nitrate, and silicate climatologies were taken from Broullón et al. (2019; https://doi.org/10.20350/digitalCSIC/8644, last access: 1 August 2019). These climatologies of nutrients were created using the objectively analyzed climatologies of temperature, salinity, and oxygen (Garcia et al., 2014; oxygen climatology from WOA13 has been filtered by applying a fifth-order one-dimensional median filter through the depth dimension; see Broullón et al., 2019) from the WOA13 in CANYON-B (Bittig et al., 2018). As a year input is needed, we decided to center the TCO₂ climatology around 1995 based on the time distribution of the data used to create the WOA13 climatologies: the World Ocean Database 2013 (Boyer et al., 2013).

The computed climatological values were compared with those from measured data to assess the uncertainty of the climatology since the WOA13 does not offer an uncertainty field with the objectively analyzed climatologies. Unfortunately, only two locations have enough measured data to calculate a pure climatological value of TCO₂ for each month: HOT ALOHA and BATS. The measured values were monthly averaged at several depth levels, and the anthropogenic carbon as calculated by Lauvset et al. (2016) was added or subtracted to correct the data to the reference year of the climatology according to

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\text{TCO}}_{\mathrm{2}}^{{\text{year}}_{\mathrm{2}}}={\text{TCO}}_{\mathrm{2}}^{{\text{year}}_{\mathrm{1}}}-{\mathrm{C}}_{{\text{ant}}_{\mathrm{2002}}}\left[\right(\mathrm{1}+\mathrm{0.0191}{)}^{({\text{year}}_{\mathrm{1}}-\mathrm{2002})}\\ \text{(2)}& {\displaystyle }& {\displaystyle }-(\mathrm{1}+\mathrm{0.0191}{)}^{{\text{year}}_{\mathrm{2}}-\mathrm{2002}}],\end{array}$$

where ${\text{TCO}}_{\mathrm{2}}^{{\text{year}}_{\mathrm{2}}}$ is the TCO₂ corrected for year₂, which is the reference year of the climatology; ${\text{TCO}}_{\mathrm{2}}^{{\text{year}}_{\mathrm{1}}}$ is the TCO₂ measured in year₁, C${}_{{\text{ant}}_{\mathrm{2002}}}$ is the anthropogenic carbon for 2002; and 0.0191 is the annual increase rate derived from the scaling factor determined by Gruber et al. (2019) for the global ocean between 1994 and 2007.

We compared our climatology with previously published climatologies of TCO₂. The monthly surface climatology created by Takahashi et al. (2014) was used to assess the spatiotemporal differences in the surface layer. The annual climatology of Lauvset et al. (2016) was used to evaluate the spatial differences in the deeper parts of the ocean. For the comparisons, the climatologies of Takahashi et al. (2014) and Lauvset et al. (2016) were adjusted to the year 1995, subtracting the anthropogenic carbon (C_ant) of Lauvset et al. (2016) as in Eq. (2).

Finally, a surface climatology of pCO₂ was computed from the TCO₂ climatology of the present study and the A_T climatology of Broullón et al. (2019) to assess the potential of computing climatologies of other variables of the seawater CO₂ system. For comparison, the updated monthly pCO₂ climatology from Landschützer et al. (2016, 2017) was used. The values between 1981 and 2010 were averaged to obtain the climatological year 1995. The variable selected from Landschützer et al. (2017) was that labeled as spco2_raw (sea surface pCO₂) in the netCDF file.

It should be noted that the RMSE and the bias were obtained for all the comparisons, the last statistic being computed as the difference between the measured (or computed by the method to compare) TCO₂ and the one obtained with the neural network of the present study.

3 Results

3.1 Neural network analysis

Following the established criteria to obtain the optimal number of neurons, the configuration with 128 neurons in the hidden layer was selected. From the 10 networks trained with this number of neurons for each approach (NNGv2LDEO and NNGv2QCLDEO), the ones with the lowest influence of the position input variables were selected. These two networks present a similar RMSE in both training and test datasets, showing there is no overfitting. Because in Gv2QC, both NNGv2LDEO and NNGv2QCLDEO produce the same global RMSE (6.1 µmol kg⁻¹), it is likely that the Gv2 dataset contains high-quality measurements, and the possible errors in the non-QC data of this dataset are clearly avoided by the network; otherwise NNGv2LDEO should have a higher RMSE in the test dataset than NNGv2QCLDEO because of an overfitting of the errors in the Gv2 dataset. The same holds for the LDEO dataset. The network properly fitted TCO₂ derived from the LDEO dataset since it does not significatively increase the global RMSE relative to a network only trained with Gv2. Therefore, we decided to continue with NNGv2LDEO only since it has fitted more relationships between variables (e.g., Gv2 has more data points than Gv2QC in the Mediterranean Sea), providing a more robust fitting. For this network, the influence of each input variable on the computed TCO₂ is depicted in Fig. S2. The position variables together (latitude, clongitude, slongitude, and depth) have no more than 30 % influence, allowing biogeochemical variables to be the main ones responsible for the variability of TCO₂. Furthermore, the input variable year has an influence lower than 5 %. This is probably responsible for capturing the positive interannual trend due to the TCO₂ increase derived from anthropogenic emissions of CO₂ to the atmosphere (see Sect. 3.2).

Figure 2 Regression of TCO₂ computed using NNGv2LDEO and TCO₂ in Gv2 and the LDEO dataset. The graph is divided into pixels. The color of each pixel is determined by the number of points inside it. Note the logarithmic scale of the pixels accounting for the large number of data.

The global RMSE is quite low for the Gv2 dataset and for the LDEO dataset (Fig. 2). The measured and the computed data are highly correlated (Fig. 2), and the bias is negligible in both datasets. The higher RMSE in the LDEO dataset likely results from the higher variability of TCO₂ in the surface layer and from uncertainties in its calculation from pCO₂.

Table 2 RMSE (bias) by area and depth for TCO₂ and pCO₂, computed with CANYON-B and NNGv2LDEO in Gv2QC and LDEO datasets. For each depth range, the RMSE (bias) in each area was weighted by the contribution of its data to the total. Units are micromoles per kilogram (µmol kg⁻¹) for TCO₂ and microatmospheres (µatm) for pCO₂.

The RMSE by area and depth for NNGv2LDEO and CANYON-B in Gv2QC is shown in Table 2. The highest errors for the two methods are in the 0–50 m layer for the Gv2QC dataset and the LDEO dataset. These errors get smaller with increasing depth for all areas, and the depth-weighted RMSE of the two methods is not significantly different below 50 m. In the LDEO dataset, NNGv2LDEO produces a lower error than CANYON-B, except for two areas: East GIN (Greenland, Iceland, and Norwegian) Seas and the Bengal Basin (Table 2), although there are only 9 and 13 data points, respectively, in each area. Interestingly, CANYON-B is able to reproduce the TCO₂ data derived from the complete LDEO dataset with a lower error than the one it obtains for the complete Gv2QC dataset in the surface ocean (RMSE LDEO: 16.4 µmol kg⁻¹; RMSE Gv2QC, 0–5 m: 17.8 µmol kg⁻¹), supporting the approach of computing reliable TCO₂ values from the pCO₂ of the LDEO dataset and the A_T computed with NNGv2 (Broullón et al., 2019) since CANYON-B was not trained with the LDEO database. A similar result was obtained for NNGv2LDEO but with a higher difference between the two errors (RMSE LDEO: 11.4 µmol kg⁻¹; RMSE Gv2QC, 0–5 m: 17.1 µmol kg⁻¹). Finally, the surface RMSE towards LDEO data of NNGv2LDEO is clearly lower than that of CANYON-B. This shows the value of including pCO₂-derived surface TCO₂ among the training data, through which there are more fitted relations in our new method.

For data from Gv2 where no QC was performed for at least one of the variables used in the present study (Gv2noQC), the RMSE also decreases with increasing depth ( <50 m: 22.5 µmol kg⁻¹; 50–200 m: 9.8 µmol kg⁻¹; 200–500 m: 7 µmol kg⁻¹; 500–1000 m: 5.4 µmol kg⁻¹; > 1000 m: 5.4 µmol kg⁻¹). Thus, the error in Gv2noQC is similar to that in the areas with the highest error in Gv2QC (Table 2; except in Beaufort Sea, where the error is considerably higher). However, the higher error in Gv2noQC is mainly caused by the samples located in the Arctic Ocean since cruises in the Atlantic and Pacific oceans are modeled with a very low error. Therefore, using Gv2noQC does not imply the introduction of low-quality data in our study; otherwise the network would not compute TCO₂ with low errors in Gv2QC because of an overfitting of the possible low-quality data that Gv2noQC could contain.

Figure 3 Differences between (a) Gv2 TCO₂ and NNGv2LDEO TCO₂ (0–30 m) and (b) LDEO TCO₂ and NNGv2LDEO TCO₂ (0 m). This figure was made with Ocean Data View (Schlitzer, 2016).

Figure 4 Histograms and box plots of differences between measured and neural-network-computed TCO₂ in (a) Gv2 and (b) LDEO. ^* TCO₂ computed from measured pCO₂ and neural-network-derived A_T.

In general, the highest differences between measured and estimated TCO₂ occurs in the high-latitude surface oceans (Figs. 3 and 4). In Gv2, 40 % of the samples with differences beyond ±3RMSE (3 times RMSE; threshold selected to refer to samples with large residuals) are at latitudes greater than 70^∘ N. In the LDEO dataset, 39 % of the samples with differences beyond ±3RMSE are from latitudes south of 70^∘ S. These samples where RMSE is high are 7.5 % of the total north of 70^∘ N in Gv2 and 42 % of the total south of 70^∘ S in the LDEO dataset. The samples with low salinities have the highest errors (Fig. 4). A total of 41.5 % of the samples in Gv2 and 43 % in the LDEO dataset with differences beyond ±3RMSE have salinities below 33. Furthermore, in the LDEO dataset, the number of samples with residuals beyond ±3RMSE increases with increasing SD of both pCO₂ and salinity in the monthly averaging in each pixel in the LDEO subset (Fig. S3). This result shows the difficulty of modeling areas with a high submonthly variability in pCO₂ and salinity and supports the exclusion of the averaged LDEO data with a high SD since it could cause the network to interpret the submonthly variability as monthly variability (note that the purpose of this study is to capture the monthly variability).

Like for modeling A_T (Takahashi et al., 2014; Broullón et al., 2019), the Arctic Ocean is one of the regions with the highest RMSE of neural-network-estimated TCO₂. The major Arctic rivers contribute TCO₂ concentrations ranging between 400 and 3600 µmol kg⁻¹ (estimated by Tank et al., 2012), derived mainly from carbonated rocks in the watersheds. Other areas like the Okhotsk Sea also show a high RMSE (Table 2 and Fig. 3), probably because of the high riverine input of TCO₂ (Watanabe et al., 2009). An input variable accounting for the contribution of the rivers to the TCO₂ pool would improve the neural network performance in areas like these, but it is not available.

The errors of the pCO₂ computed in the LDEO dataset with TCO₂ from NNGv2LDEO and A_T from NNGv2 (Broullón et al., 2019) are similar to the errors obtained by Landschützer et al. (2014) for the SOCAT database in some of the areas (10–16 µatm; Table 2). This result shows the potential of computing pCO₂ values with neural networks trained for other variables of the seawater CO₂ system, at least in some ocean regions. The global error of the pCO₂ in the LDEO dataset is clearly higher than that obtained by Landschützer et al. (2014) for the SOCAT dataset (22 vs. 12 µatm, respectively), although the critical areas are mainly the same (Fig. S4): equatorial Pacific upwelling system, Arctic and subarctic waters around the Alaska Peninsula, the Southern Ocean, the Gulf Stream, and the North Atlantic Current. At this point, the following should be considered: (1) the pCO₂ computed in the present study is derived from A_T and TCO₂ and not from specific modeling for pCO₂, and therefore it contains errors associated with this computation (∼6 µatm; Millero, 1995) and the neural network estimates of A_T and TCO₂; (2) the present study includes the Arctic region where the highest errors occur (Table 2; Beaufort Sea and High Arctic areas); and (3) there is a longer temporal range in the present study (1973–2016). The analysis of Landschützer et al. (2014) in the LDEO dataset for data that differ from SOCAT shows a global error higher than the one obtained in the present study for all LDEO data between 1998 and 2011 (25.9 vs. 21.3 µatm, respectively). The error between 40^∘ S and 40^∘ N is similar in the two studies (Landschützer et al., 2014: 16.5 µatm; NNGv2LDEO: 16.4 µatm). Although it is not the main objective of this work, these last two results show how NNGv2LDEO and NNGv2 (Broullón et al., 2019) have the potential to compute pCO₂ values between 40^∘ S and 40^∘ N with similar errors as the method with the lower error in the pCO₂ modeling to obtain a climatology and with lower errors in high latitudes for the LDEO dataset, even taking into account the inclusion of the critical area of the Arctic in the computation of the error of the pCO₂ from the present study (it is not included in Landschützer et al., 2014) and the higher number of data from high latitudes in the present study (15 479 vs. 3799).

3.2 Time series validation

The good generalization of the network in the test dataset containing data from Gv2 and the LDEO dataset with a similar RMSE as the one reached in the training set is also evidenced through independent time series data (Table 3). Except for KERFIX, where the number of data points is very low and Olsen et al. (2019) suggested an adjustment to the original data of −39 µmol kg⁻¹, TCO₂ computed using NNGv2LDEO and CANYON-B at the time series locations is characterized by low errors and biases (Table 3). NNGv2LDEO computes TCO₂ with a lower RMSE and bias than CANYON-B for most of the time series stations (Table 3). CANYON-B reaches a lower RMSE in HOT ALOHA SURFACE and ESTOC than NNGv2LDEO, but the bias is considerably higher in these time series for CANYON-B.

Table 3 RMSE and bias between measured and computed TCO₂ concentrations in several time series. The comparison was done using only water samples where all the input variables for NNGv2LDEO and the TCO₂ were measured in the same water sample.

The seasonal variability is well captured by NNGv2LDEO, showing its great potential to design a monthly climatology. In the surface layer, where the seasonal variability is the highest, the computed values are strongly correlated with the measured TCO₂ in all the time series (Fig. 5). In addition, the high correlation holds for all depths (Table S1). The TCO₂ computation with a low error in these time series located in different oceanic regimes as well as in the areas of Table 2 shows the good performance of NNGv2LDEO in almost any region of the ocean.

Figure 5 Measured (blue line) and computed (orange line) TCO₂ with NNGv2LDEO for the depth range 0–15 m (0–30 m in panel b) for several time series. (a) BATS, (b) HOT ALOHA SURFACE, (c) ESTOC, (d) ICELAND, (e) IRMINGER, (f) KNOT, (g) K2 and (h) OWS.

Table 4 RMSE and bias between measured pCO₂ (and in some cases, computed from measured A_T and TCO₂ in time series where pCO₂ was not measured) and computed pCO₂ with A_T from NNGv2 (Broullón et al., 2019) and TCO₂ from NNGv2LDEO in several time series. The time period for pCO₂ from this study is the same as in Table 3. Consult Table 2 in Landschützer et al. (2014) for its time period. The depth range is 0–15 m. Only time series with more than 30 data points are included. RMSE and bias for computed A_T with NNGv2 (Broullón et al., 2019) and TCO₂ with NNGv2LDEO are included to show the errors in the variables used to compute TCO₂.

Assessing the potential of neural networks to obtain values of other variables of the seawater CO₂ system in the time series, pCO₂ calculated with A_T from NNGv2 (Broullón et al., 2019) and TCO₂ from NNGv2LDEO compared quite well with pCO₂ as measured or calculated from A_T and TCO₂ at the time series stations (Table 4). Except for BATS, the pCO₂ obtained in the present study has a lower error than that reported by Landschützer et al. (2014; Table 4). In contrast, the bias in the present study is higher, except for ESTOC. Considering the error involved in the calculation of pCO₂ from A_T and TCO₂ (∼6 µatm); Millero, 1995), and the error in the computed A_T and TCO₂ with the neural networks (Table 4), our results demonstrate again the ability of NNGv2 and NNGv2LDEO to calculate other variables of the seawater CO₂ system with a relatively low error.

Table 5 Long-term trends (seasonally detrended) of the measured and computed TCO₂ and pCO₂ from neural networks at time series locations in the depth range 0–15 m.

^* Computed from measured A_T and TCO₂ in time series where pCO₂ was not measured.

Using NNGv2LDEO, it is also possible to reproduce the secular trends in TCO₂. Using seasonal detrending to enhance the multiannual changes, similar trends in the longer time series are found for the measured TCO₂ and the neural-network-computed TCO₂ (Table 5). The same holds for pCO₂ (Table 5), although at the IRMINGER site the trend obtained from the neural-network-generated data is significantly lower than that from measured data. The neural networks seem to capture the anthropogenic influence in the seawater CO₂ system and thus the ocean acidification process (Fig. 6). Furthermore, using NNGv2LDEO increases the number of TCO₂ data where the various inputs were measured but not TCO₂ itself. This allows for the evaluation of high-frequency changes (Fig. 6) and for the calculation of interannual trends with a low error (as temporal sampling biases are reduced).

Figure 6 Time series of TCO₂ using NNGv2LDEO at (a) BATS and (b) HOT ALOHA locations. The water column shows a higher concentration of TCO₂ year by year. This figure was made with Ocean Data View (Schlitzer, 2016).

3.3 Climatology

Using NNGv2LDEO we have demonstrated its ability to compute TCO₂ values with low errors and especially to capture the monthly variability of this variable. In addition, the climatologies of the input variables used to create the climatology of TCO₂ have been satisfactorily evaluated previously for the construction of an A_T climatology (Broullón et al., 2019). Considering these results, a monthly climatology of TCO₂ is obtained by passing the input climatologies through NNGv2LDEO.

The spatial distribution of the surface annual mean climatology of TCO₂ (Fig. 7a) is similar to two recent climatologies: those of Takahashi et al. (2014) and Lauvset et al. (2016). The largest surface TCO₂ concentrations occur in the Southern Ocean, subpolar North Atlantic, Nordic Seas, and Mediterranean Sea (note that the latter is not included in these other climatologies). In general, surface TCO₂ decreases from high to low latitudes. The Indian and the Pacific oceans are characterized by lower concentrations of TCO₂ at higher latitudes than the Atlantic, the latter being the ocean with the highest surface TCO₂ by area. TCO₂ increases with depth in all oceans, in particular in the upwelling regions, where this increase is expanded eastwards with depth (Fig. 7b and video at https://doi.org/10.20350/digitalCSIC/10551, Broullón et al., 2020). Depending on the area, the values reach a maximum at certain intermediate depths, and below it the concentration gradually decreases or remains almost constant (Fig. S5).

Figure 7 Annual mean climatology of TCO₂ at (a) 0 m, (b) 100 m, and (c) 1000 m. This figure was made with Ocean Data View (Schlitzer, 2016).

Figure 8 Seasonal amplitude of TCO₂ at (a) 0 m, (b) 100 m, and (c) 1500 m. The contour lines of 25, 50, 75, and 100 µmol kg⁻¹ are shown. This figure was made with Ocean Data View (Schlitzer, 2016).

The largest seasonal variability occurs at the surface at high latitudes, in the Pacific upwelling region, the equatorial African coasts, and in the area under influence of the Amazon River (Fig. 8a). At depth, the seasonal variability decreases, except for the Pacific upwelling region, where it increases and moves progressively northward between 30 and 150 m (Fig. 8b). This increase is correlated with the high seasonal variability of the climatologies of nutrients, oxygen, and temperature at these depths. Czeschel et al. (2012) also showed an increase in the subsurface variability of oxygen from measured profiles. Similar increases also occur in the Indian Ocean north of 20^∘ S between 50 and 100 m and in the equatorial Atlantic Ocean in the same depth range. At the 1500 m level, the seasonal variability is below 10 µmol kg⁻¹ in most of the ocean (Fig. 8c). This last result shows that an annual climatology below 1500 m is sufficient.

Table 6 Comparison of four annual mean surface climatologies of TCO₂. Numbers in the lower left corner represent RMSE. Numbers in the upper right corner represent r².

^* The domain analyzed is the same as in Takahashi et al. (2014) for coherency reasons.

Although the surface patterns of the annual mean of the TCO₂ climatology are very similar to those of the other recent climatologies (Takahashi et al., 2014; Lauvset et al., 2016), differences do occur. The annual mean climatology of the present study is closest to that of Takahashi et al. (2014; Table 6). The largest differences between these two climatologies are located in the Arctic, North Pacific, Peru upwelling area, western South Pacific, and the area of influence of the Antarctic Circumpolar Current (Fig. S6a). The Atlantic and the Indian oceans do not show significant differences. Our climatology shows more deviations from that of Lauvset et al. (2016), compared in the grid of Takahashi et al. (2014; Table 6). The highest differences are found in the North Pacific, around Antarctica, the Nordic Seas, the South and North Atlantic, and in several less localized areas around the oceans (Fig. S6b). When the climatology of Takahashi et al. (2014) is compared to that of Lauvset et al. (2016), the differences are even higher (Table 6), and the critical areas are the same of those of the previous comparison. Although it is clear that discrepancies between the three climatologies are derived from the different methods used, the higher similarity between ours and the one of Takahashi et al. (2014) is probably due to the influence of the same source used to create them, the World Ocean Atlas.

The comparison of our climatology with that of Lauvset et al. (2016) at the 33 depth levels of Lauvset et al. (2016) shows a reduction in the RMSE with depth. Between 0 and 1000 m, the RMSE is reduced from ∼32 to 7 µmol kg⁻¹ (Table S2; note the higher RMSE at surface compared to the one obtained for the grid of Takahashi et al., 2014, because of the inclusion of areas which are not included in the latter's grid, and the difficulty of modeling TCO₂ in some areas, like the Arctic and the Mediterranean Sea). This reduction with depth is probably due to the reduction in the variability in most of the ocean below the surface. The surface values in Lauvset et al. (2016) are likely characteristic of months in which most of the sampling was carried out. Because of the lower variability of TCO₂ at depth, the values are closer to the annual mean, and therefore the two compared climatologies are more similar at depth than at surface depth levels. Below 1000 m, the differences between the two climatologies are not significant, with an RMSE around 5 µmol kg⁻¹ and a bias around 0.5 µmol kg⁻¹ at each depth level (Table S2).

Our monthly climatology shows a high correspondence with that of Takahashi et al. (2014), although the RMSE values show that there are also large differences in certain areas (Table 7). These areas are mainly the same as those in the comparison of the annual mean climatologies, but some other small regions with high differences appear for each month all through the ocean (Fig. S7).

Table 7 Comparison of the monthly TCO₂ climatology of Takahashi et al. (2014) and the one of the present study.

Figure 9 Comparison of the monthly climatological profiles of TCO₂ computed from measured data (red profile; shadow area is the SD of the averaged values at each depth level) and those from the TCO₂ climatology at the (a) BATS and (b) HOT ALOHA locations. The units on the x axis are micromoles per kilogram (µmol kg⁻¹).

Unfortunately, the uncertainty of the TCO₂ climatology cannot be assessed globally and robustly. As Broullón et al. (2019) stated, the unavailability of an uncertainty field associated with the objectively analyzed WOA13 climatologies does not allow us to perform a proper global uncertainty assessment. Therefore, the analysis is relegated to the areas where repeated sampling of TCO₂ has been carried out monthly over a long period, that is, the HOT ALOHA and BATS time series stations. The climatology of TCO₂ from NNGv2LDEO is consistent with the monthly climatological values at these two places (Fig. 9). In general, the profiles from the TCO₂ climatology are within the variability range (shadow area in Fig. 9) of the monthly averaged measured data for each depth level. In the upper 30 m of the water column, the climatology of TCO₂ differs from the measured BATS data from May to August. This difference is mainly explained by the surface error of the network shown for this time series in Fig. 5a, where the computed TCO₂ decreases from maximum to minimum sooner than the measured TCO₂. For HOT ALOHA, the RMSE of the profiles of the TCO₂ climatology oscillates between 3.6 and 9.2 µmol kg⁻¹ with a mean value of 6.3 µmol kg⁻¹ (bias range: −3.8 to 1.2 µmol kg⁻¹; mean bias: −1.4 µmol kg⁻¹). At BATS, the RMSE is lower than for HOT ALOHA: 1.1 to 8.5 µmol kg⁻¹ with a mean value of 4.4 µmol kg⁻¹ (bias range: −1.4 to 7.6 µmol kg⁻¹; mean bias: −2.5 µmol kg⁻¹). Furthermore, the seasonal variability of the TCO₂ climatology is quite similar to that of the measured data at BATS and HOT ALOHA (Fig. S8). Although in other time series there are not enough measured data to obtain climatological values, these pseudoclimatological values also correlate very well with the TCO₂ climatology (data not shown). These results suggest that the climatology is robust in different oceanographic regimes and adequately captures the seasonal cycle of TCO₂.

It has been demonstrated in this study that pCO₂ and possibly other variables of the seawater CO₂ system can be computed from A_T and TCO₂ derived from neural networks with a relatively low error in different datasets (LDEO in Sect. 3.1 and time series in Sect. 3.2). The pCO₂ climatology (Fig. S9) computed from the TCO₂ climatology of the present study and the A_T climatology of Broullón et al. (2019) is very similar to that of Landschützer et al. (2017). The differences between the annual mean climatology of the two studies are below 15 µatm in most of the ocean (RMSE: 8.3 µatm; bias: 2.9 µatm; r²: 0.82). The differences above this threshold are mainly located in the Pacific equatorial upwelling system, the eastern part of the South Pacific Gyre, Nordic Seas, Labrador Sea, Atlantic section of the Southern Ocean, Bay of Bengal, and the waters surrounding the eastern margin of Asia (Fig. S10). In most of these areas, both methods have the greatest errors (Figs. 2 and 4 in Landschützer et al., 2014, and Fig. S4 of the present study).

On a monthly basis, the RMSE between the two pCO₂ climatologies is between 13.6 and 15.6 µatm, and the correlation is lower than for the annual mean comparison (r²: 0.55–0.72 vs. 0.82). The areas with the higher differences are the same as in the annual comparison, but other small regions appear along the ocean month by month (Fig. S11). Furthermore, the seasonal variability in the two climatologies matches in a great extension of the ocean, although there are areas with notable differences (Fig. S12). In general, the pCO₂ climatology is quite similar to that of Landschützer et al. (2017), and this result helps show that both the TCO₂ climatology of the present study and the A_T climatology of Broullón et al. (2019) are mostly robust and suggest that climatologies of other seawater CO₂ system variables can be confidently computed.

4 Data availability

The climatologies of TCO₂ and pCO₂ as well as NNGv2LDEO designed in this study are available at the data repository of the Spanish National Research Council (CSIC; https://doi.org/10.20350/digitalCSIC/10551; Broullón et al., 2020).

5 Conclusions

We presented a tool for computing TCO₂ in the global ocean. Compared to previous methods, the uncertainties in such computations have been reduced. Including two updated datasets containing thousands of measurements of inorganic carbon variables across the ocean in the training of the neural network, we were able to capture a wide range of variability of TCO₂. The low errors obtained in independent subsets at time series stations are further evidence of the potential of the network in computing TCO₂.

Our global monthly climatology created with a neural network is the first that covers the oceans from the surface to the abyss at such temporal resolution. In addition to the accuracy of the network, the low uncertainty of the climatology in different regions and its usefulness in creating climatologies of other seawater CO₂ chemistry variables (i.e., pCO₂) show its robustness. Therefore, we present the global climatology of TCO₂ to the scientific community to complement the recently designed climatology of A_T by Broullón et al. (2019) for its use in the initialization and evaluation of models or any other analysis related to the carbon cycle.