The Global Ocean Data Analysis Project version 2 (GLODAPv2) was an effort from the international community to consolidate all data from ocean bottle samples collected as part of many oceanic cruises (Olsen et al., 2016). The GLODAPv2 database (available at http://cdiac.ornl.gov/oceans/GLODAPv2/) provides a single high-quality internally consistent global data product that contains CO₂-relevant ocean interior measurements from ship-based surveys. The GLODAPv2 database includes samples of core variables such as salinity, oxygen, macronutrients, and seawater CO₂ chemistry from 724 oceanic cruises. In this study, we focused on seven variables representative of the macronutrients and of the seawater carbonate system: nitrate (NO₃⁻), phosphate (PO₄³⁻), silicate (Si(OH)₄), pH on the total scale (pH_T), total alkalinity (A_T), total dissolved inorganic carbon (C_T), and partial pressure of CO₂ (pCO₂). Note that we estimated this last variable from the A_T and C_T measurements available in GLODAPv2 (see details below).

Initially, GLODAPv2 was instigated to prepare a unified, bias-corrected interior ocean data product. Thus, a high quality control, QC, based on two steps (i.e., primary and secondary QC) was applied to each data (Olsen et al., 2016). The primary QC was carried out following routines outlined in Sabine et al. (2005) and Tanhua et al. (2010), essentially based on inspection of property-property plots. The secondary QC for salinity, oxygen, nutrients, C_T, and A_T was more complex, and carried out through crossover (i.e., comparing data where two different cruises crossed or came close to each over) and inversion analyses (i.e., calculation of corrections required to minimize all cruise-by-cruise offsets). This two step-based method was introduced by Gouretski and Jancke (2000) and Johnson et al. (2001). For the secondary QC applied to the GLODAPv2 database, the crossover offsets were calculated using the running-cluster crossover routine (Tanhua et al., 2010), with data from beneath 2,000 m to minimize effects of real variations. For pH_T, crossover analysis was not possible because data only exist for a small fraction of the cruises. To pass the secondary QC, pH_T measurements had to be concomitant with C_T and/or A_T for calculating offsets (see details in Olsen et al., 2016). For Mediterranean Sea data, the secondary QC always failed because none of the cruises inside the Mediterranean had an overlap with other cruises (e.g., outside the Mediterranean) thus preventing the crossover analysis. Hence, only “high-quality” GLODAPv2 data that passed the secondary QC were used to train and validate the CANYON methods except for the Mediterranean Sea where we used data that only passed the primary QC.

The subset of the GLODAPv2 database used for our study (i.e., the data that passed the secondary QC, except for the Mediterranean Sea data as explained above) contained 37,863 concurrent profiles of water-column (from the surface to a maximum sample depth of 8,000 m) hydrological properties together with nutrients concentration and/or parameters of the carbonate system (see Figure 1). These data were collected between 1972 and 2013 and were representative of the diversity of oceanic regions, i.e., 25 % were collected in the North Atlantic, 10 % in the South Atlantic, 22 % in the North Pacific, 12 % in the South Pacific, 10 % in the Indian Ocean, 13 % in the Southern Ocean, 7 % in the Arctic Ocean, and ~0.2 % in the Mediterranean Sea (geographic boundaries are provided in Figure S1). On the temporal scale, most of the data were acquired since the 1990's and more data were available for the spring and summer months (Figure S2). There was a sampling bias according to latitude as data from autumn and winter months (i.e., December to March for the Northern hemisphere, and May to August for the Southern hemisphere) were less represented at high than low latitudes (i.e., >45° North and South, respectively) in the GLODAPv2 database (Figure S2).

All the data used to train and validate CANYON were measurements recorded in the GLODAPv2 database, except the pCO₂ estimates that we calculated from A_T and C_T measurements using the R package “seacarb” (Gattuso et al., 2015, 2016). The carbonate system parameters were computed using the carbonic acid dissociation constants of Lueker et al. (2000), the hydrogen fluoride dissociation constant of Perez and Fraga (1987), the dissociation constant for bisulfate of Dickson (1990), and a ratio of total boron to salinity derived from Uppström (1974). In situ measurements of salinity, temperature, hydrostatic pressure as well as the concentrations of PO₄³⁻ and Si(OH)₄ were used to calculate pCO₂. When not available in the GLODAPv2 database, the concentrations of PO₄³⁻ and Si(OH)₄ were estimated using our CANYON algorithm (see associated accuracies in Section Overall CANYON Performance).

For the neural network development, vertical profiles of nutrients and carbonate system parameters from eight independent zones of the GLODAPv2 database (squares of 10° latitude × 10° longitude) were first removed from the general database to provide a more “independent data set” used for an independent validation of the algorithm developed in this study. These zones were chosen in several major oceanic basins and were representative of the Sub-Equatorial Pacific, the Sub-Equatorial Indian, the North Atlantic Subtropical Gyre, the North Atlantic Subpolar Gyre, the North Pacific, the South Atlantic, the South Indian, and the South Pacific (Figure 1). The remaining profiles were then split into two subsets with 80 % and 20 % of the data, the so-called training and validation datasets, respectively (see the number of data for each variable in Table 1).

Using the validation database (i.e., 20 % of the general database minus the eight independent zones), we evaluated the performance of the method by comparing the CANYON-retrieved nutrient concentrations and carbonate system parameters with the measurements in the GLODAPv2 database using the statistics from Section Statistical Evaluation of Method Performance. Scatterplots of CANYON-retrieved variables vs. GLODAPv2 measurements (Figure 3) show that the CANYON method predicts nutrient concentration and carbonate system parameters with good accuracy (i.e., of 0.93, 0.066, and 3.0 μmol kg⁻¹ for the concentrations of NO₃⁻, PO₄³⁻, and Si(OH)₄, respectively, and of 0.019, 7 μmol kg⁻¹, 10 μmol kg⁻¹ and 5.1 % or 20 μatm at 400 μatm for pH_T, A_T, C_T, and pCO₂, respectively). The determination coefficients of the seven linear models between the CANYON-retrieved and GLODAPv2-variables are comprised between 0.982 and 0.996 with slopes ranging from 0.986 to 0.999. In Figure 3 only very few data points diverge from the 1:1 line. A higher scatter is observed for low CANYON-retrieved NO₃⁻ (and PO₄³⁻), which mostly corresponds to low surface nutrient concentrations inside and at the edges of the subtropical gyres. Moreover, the higher scatter observed near the surface than deeper for most variables is probably due to the higher inherent variability in the surface.

To go further, the final CANYON-accuracies for the seven variables can be estimated using the merged accuracies of CANYON estimations and GLODAPv2 measurements (RMSE_final = √[RMSE_CANYON² + RMSE_{RMSEGLODAPv2²}]). The accuracies of GLODAPv2 measurements are 0.46, 0.033, and 1.1 μmol kg⁻¹, for NO₃⁻, PO₄³⁻, Si(OH)₄, respectively, and 0.005, 6 μmol kg⁻¹, and 4 μmol kg⁻¹, for pH_T, A_T, and C_T, respectively from Olsen et al. (2016) and 5.6 % (22 μatm at 400 μatm) for pCO₂ from uncertainty propagation of the carbonate system calculations using seacarb errors (Gattuso et al., 2015, 2016). Thus, ultimately, the final global accuracies of CANYON are 1.04, 0.074, and 3.2 μmol kg⁻¹ for NO₃⁻, PO₄³⁻, and Si(OH)₄ concentrations, respectively, and 0.020, 9 μmol kg⁻¹, 11 μmol kg⁻¹ and 7.6 % (30 μatm at 400 μatm) for pH_T, A_T, C_T, and pCO₂, respectively.

The training and validation datasets of the neural networks used to retrieve carbonate system parameters were smaller than the datasets for the retrieval of nutrient concentrations (Figure 3 and Table 1). It is thus possible that the carbonate system networks are less robust than the nutrient ones. In any case, all the MLPs could be updated in the future as more data become available; this seems especially important for the pH_T database, which is presently the least populated. In order to assess the importance of this potential weakness, we developed a special neural network using all pH_T data available, i.e., all the data that passed the primary quality control (see details in Section The GLODAPv2 Database). The results of this special CANYON algorithm, based on more but a priori less accurate data used for training, are not improved when compared to our initial results (i.e., RMSE of 0.030). Given this and in order to maintain consistency among CANYON algorithms and their retrieval performance, all neural networks were trained using data that had passed the secondary quality control (except for the Mediterranean Sea, see details in Section The GLODAPv2 Database).

To identify possible trends, errors were plotted against each input variable for each retrieved nutrient concentration (Figure 4) and carbonate system parameter (Figure 5). Some errors are larger because of the small numbers of data (the intensity of the shading in each box refers to the number of data). Nevertheless, some clear trends are observed. The range of errors for Si(OH)₄ retrieval seems to be higher at high latitudes in the Southern hemisphere (i.e., latitudes <−60°) and low temperatures, i.e., the ranges of box-plot whiskers increased with both decreasing latitude and temperature in Figure 4. This corresponds to regions of significant Si(OH)₄, mostly in the Southern Ocean, which suggests that the CANYON method is less accurate for Si(OH)₄ retrievals in the Southern Ocean than in other areas. The C_T estimates show a high error at high temperatures and unusually low (<34 psu) salinities. pCO₂ estimates exhibit an increased uncertainty at the extremes of the O₂ input (Figure 5), i.e., the range of errors increases at low O₂ concentrations corresponding to high pCO₂ (Figure S3), and at high O₂ concentrations corresponding to cold, low salinity polar surface waters. For CANYON-estimated NO₃⁻, PO₄³⁻, and C_T, the upper layer (i.e., ≤700 m) has the broadest error range. The Si(OH)₄ displays an opposite trend with error larger for deeper than surface waters (i.e., ≥700 m). The CANYON-estimated A_T, pH_T, and pCO₂ errors are not affected by depth inputs. Finally, Figure S3 shows that the CANYON retrieved variables are not biased against the range of in situ values to retrieve, except for a few extreme values where few data (i.e., lightly shaded boxes) were available in the training and validation databases.

The above results indicate that hydrographic data and information about season (doy) and geolocation can certainly predict some aspects of the dynamics of biogeochemical variables at the surface (e.g., nutrient and C_T drawdown during the spring bloom, seasonal reset to preformed nutrient concentrations with winter ventilation) with O₂ being the most important predictor in CANYON for production and remineralization, particularly in the ocean interior. The year is the only variable that accounts in CANYON for the increase in anthropogenic CO₂. As a consequence, when the year was not included, CANYON-estimated pH_T, C_T, and pCO₂ showed a clear remaining trend due to the missing information about the long-term changes in seawater CO₂-carbonate chemistry (data not shown).

The smoothed mean differences between CANYON-retrieved and in situ measurements were plotted as vertical profiles for each of the seven variables and the eight independent zones (Figure 6). In general, the accuracy (i.e., RMSE) for each variable is comparable to the accuracy determined on the validation data set (Table 3). Beyond this general agreement, there are a few discrepancies. The errors appear to be higher in the 0 – 200 m layer than below. This is maybe due to a larger variability in this upper layer, caused by not only biogeochemical processes but also air-sea exchange of O₂ that act to decouple O₂ from the CANYON outputs (see also Figure 3). The Sub-Equatorial Pacific, North Pacific and South Atlantic display higher RMSE than calculated from the 20 % validation data (Table 3). These results suggest that the GLODAPv2 data set for these specific zones could have been underrepresented in the training database with respect to the regional variability in nutrient concentrations and carbonate system parameters.

The above comparisons and the identification of spatio-temporal domains where CANYON performed less accurately could be used to identify more objectively the periods or regions that require more intense acquisition of discrete high-quality measurements. Indeed, the periods and regions where CANYON fails to reproduce well in situ observations are those where the variability (mostly seasonal) was not well-represented in the training dataset. From Table 3 and Figure 6, these regions are mainly the Sub-Equatorial Pacific, the North Pacific and the South Atlantic, and the periods are most likely late autumn and winter, when harsh sea conditions generally prevent ship-based collection of high-quality measurements. In fact, the GLODAPv2 database seems to be biased in this respect (Section The GLODAPv2 Database and Figure S2).

Here we present further results important for specific applications. Indeed, one of a potential application of the CANYON method is the calibration of NO₃⁻ and pH_T sensors mounted on BGC-Argo profiling floats (because the corresponding sensors may drift over long-term deployment). To overcome this problem, CANYON could be used to compute deep (e.g., ≥1,000 m) reference measurements each time the float profiles. With this aim in mind, specific MLPs were developed to retrieve NO₃⁻ concentration and pH_T only at depths between 950 and 2,050 m, i.e., different from the full water-column MLPs used above. Results show that specific, deep MLPs do not significantly improve the quality of the results compared to the MLPs developed for the entire water column applied to this layer (i.e., RMSE for NO₃⁻ of 0.50 and 0.54 μmol kg⁻¹, respectively, and for pH_T of 0.013 units for the two approaches). The CANYON method can therefore be used to retrieve NO₃⁻ concentration and pH_T for specific applications focused on deep layers with excellent accuracy.

For pCO₂ applications, most studies focus on the surface layer of the open ocean and on CO₂ air-sea exchange (e.g., Takahashi et al., 2009). The CANYON method could be used to address the regional and seasonal variability of air-sea CO₂ fluxes in view of comparing it with the results of previous studies based on neural networks (Landschützer et al., 2014, 2016). For this application, it is important to ascertain that the general MLP performs as adequately as a specific MLP developed to retrieve pCO₂ for the surface layer only (i.e., ≤100 m). Here again, the CANYON algorithm developed for the entire water column is as robust as a specific, surface-focused MLP (i.e., RMSE of 34 and 33 μatm, respectively).

The CANYON method was further applied outside the GLODAPv2 domain, on which it was trained, using an independent dataset that contained similar variables. The Hawaii Ocean Time Series (HOT) database contains monthly vertical profiles of hydrological properties, O₂, nutrient concentrations and carbonate system parameters from the deep-water station ALOHA since 1994 (Karl and Lukas, 1996; Dore et al., 2003). The HOT temperature, salinity, and O₂ measurements from HOT were used as input variables to estimate the concentrations of nutrients and the carbonate system parameters. The overall agreement between the CANYON-simulated variables and their measured in situ counterparts is satisfactory, as shown by the absolute differences between the two datasets and the corresponding comparison statistics (Figure 7 and Table 4). There are few systematic biases, such as the model underestimation of NO₃⁻ from 450 to 550 m (i.e., the depth horizon of the nitracline) and the overestimation of PO₄³⁻ and C_T in the upper 150 m (Figures 7b,d,l, respectively).

Interestingly, running CANYON on this dataset unraveled O₂ outliers in the database. A first run of CANYON produced abnormally high differences between the retrieved and in situ values of the seven CANYON-estimated variables for profiles from 25 February 2003 to 29 March 2003. This corresponded to five O₂ profiles with abnormally high concentrations of O₂ (see black triangles in Figure S4). These profiles were subsequently removed from the analysis to avoid contaminating the CANYON-retrieved data shown in Figure 7 and the corresponding statistics in Table 4.

CANYON also provided a way to fill a gap in the HOT dataset. Indeed, pH_T had not been acquired during the 1999–2003 period, but the input variables needed to run CANYON had been measured. pH_T could thus be estimated during that period (Figure 7g) with a mean accuracy of 0.033 units (i.e., RMSE in Table 4). However, it is obviously not possible to compare the values predicted by CANYON with (non-existent) corresponding in situ measurements from 1999 to 2003 (Figure 7h).

Using the HOT temperature, salinity and O₂ measurements from the last year, i.e., 2015, we estimated pH_T for the 15 years to come (by changing only the year input in CANYON), and a decline is found in pH_T of 0.024 units over this 15-year period (i.e., decrease of 0.0016 ± 0.0004 units year⁻¹). This value is consistent with the decrease in pH_T of 0.0019 ± 0.0002 year⁻¹ reported in the central North Pacific (Dore et al., 2009) and more generally of 0.0013 to 0.0026 units year⁻¹ units during the 20–30 last years (Bates et al., 2014). This suggests that the CANYON approach could perhaps also be used outside the temporal range of training for the carbonate system parameters for estimating near future changes in pH_T, C_T, and pCO₂, thanks to the use of the year among the input variables. However, this would assume that the relationships between the input and output variables, through the hidden layers of the CANYON model, will remain unchanged in the future, and the sensitivity of retrieved carbonate system parameters to departures from this assumption remains to be explored.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.