2.1 Temperature and salinity in situ profiles

We use the in situ temperature and salinity (T−S) vertical profiles sampled by Argo floats (Argo, 2022) and ships from the CMEMS quality-controlled Coriolis Ocean Dataset for Reanalysis (CORA) database (Cabanes et al., 2013; Szekely et al., 2019). We keep only the profiles extending at least from 25 to 1000 m for the period 1993 to 2019, totaling 67 767 T−S profiles for the region 80 to 30^∘ W and 23 to 50^∘ N. All the profiles that do not reach 1000 m or start deeper than 25 m are discarded. Profiles are interpolated on an uneven vertical grid with 51 levels, with spacing increasing with depth (27 levels are within the first 100 m leading to a vertical resolution of 1 m in the upper levels and 450 m at 1000 m depth). Profiles without data at the surface are extrapolated by repeating the shallowest observation point. There is little seasonal bias in the distribution of data with 5647 profiles by month on average, a minimum of 5027 in February and a maximum of 6257 in October. The spatial distribution of profiles kept in the analysis is shown in Fig. 1b. It reveals a general lack of data in the center of the subtropical gyre compared to the Gulf Stream region west of 60^∘ W. The temporal distribution (Fig. 1a) reveals a significant increase in sampling after 2000 thanks to the Argo program (Wong et al., 2020). After 2012, the number of T−S profiles stabilizes at ∼ 400 profiles per month.

Figure 1Count of temperature and salinity profiles extending from 25 to 1000 m for the region. Profiles are counted for $\mathrm{1}/\mathrm{4}$^∘ bins (b) and represented where they exceed 0 with a truncated color bar at 20 profiles; the maximum number of profiles for a bin is 185 profiles. The count in (a) is by months.

2.2 Input data

Input data (Table 1) are surface satellite data and bathymetry. Mean dynamic topography (MDT) is from the Centre National d'Etude Spatiale (CNES-CLS18) (Mulet et al., 2021), and the bathymetry is the ETOPO1 bedrock, distributed by NOAA (2009). The SLA is the level-4 daily product from CNES-CLS (6.2_DUACS_DT2018). It has a $\mathrm{1}/\mathrm{4}$^∘ horizontal grid resolution. We also use geostrophic surface velocities derived from the SLA and distributed by CNES-CLS (Table 1). MDT is calculated by merging information from altimeter data, GRACE, and GOCE gravity field and oceanographic in situ measurements (drifting buoy velocities, hydrological profiles) (Mulet et al., 2021), while SLA is from altimeter data only. Keeping MDT and SLA separated in the inputs allows us to determine their respective importance in the prediction (see Sect. 4.4). The SST is from the European Space Agency (ESA) Climate Change Initiative (CCI) and Copernicus Climate Change Service (C3S), v2.3 and level-4 product. It provides gap-free maps of daily average SST at 20 cm depth and 0.05^∘ × 0.05 ^∘ horizontal grid resolution, using satellite data from the Along Track Scanning Radiometer (ATSR), Sea and Land Surface Temperature Radiometer (SLSTR) and Advanced Very High Resolution Radiometer (AVHRR) series of sensors (Merchant et al., 2019).

Table 1List of input variables for the neural network and additional datasets to compare with. ESA CCI: European Space Agency Climate Change Initiative. C3S: Copernicus Climate Change Service. NOAA: National Oceanic and Atmospheric Administration. DUACS: Data Unification and Altimeter Combination System. CMEMS: Copernicus Marine Environment Monitoring Service.

n/a: not applicable.

Download Print Version | Download XLSX

2.3 Additional datasets

We validate our results against other observational and synthetic datasets. The sea surface salinity (SSS) CCI dataset is distributed at $\mathrm{1}/\mathrm{4}$^∘ horizontal grid resolution from 2010 to 2019 (Boutin et al., 2021). We do not use SSS as an input variable for several reasons. SSS satellite observations only cover the period 2010–2019, and its quality is questionable at high latitudes and in cold water (Boutin et al., 2021). As our prediction only depends on the input variables, it is risky to rely on data with systematic errors. Moreover, we tested an architecture with SSS as input, and results were not improved significantly: only the surface salinity was slightly better. The relative importance algorithm also showed that SSS was not used significantly in predictions. It has the same order of importance than geostrophic currents (see Sect. 3.4 on the explainability of the NN). However, SSS is a useful product to compare with, and we further discuss this in Sect. 4.5 and Appendix B.

We compare OSnet gridded fields to Armor3D and Glorys12 because they are the only ocean products, to our knowledge, that extend from 1993 to today with a spatial resolution of at least $\mathrm{1}/\mathrm{4}$^∘ and a frequency under the month (weekly for Armor3D and daily for Glorys12). The global eddy-resolving reanalysis Glorys12 (Lellouche et al., 2021) is based on the physical model NEMO (Madec, 2015) and ocean observations assimilated by means of a reduced-order Kalman filter. It is provided at $\mathrm{1}/\mathrm{12}{}^{\circ }$ horizontal grid resolution and a daily mean. We also compare our results to the observation-based Armor3D weekly product (Guinehut et al., 2012). This later product is built in two steps: (i) prediction of synthetic T−S fields by multiple linear regression from SST and SLA and (ii) optimal interpolation combining synthetic and in situ T−S profiles. A section of OSnet is compared to hydrographic section AT20 sampled along 52.3^∘ W by the research vessel Atlantis from 1 to 11 May 2012 (McCartney, 2012). Finally, we use T−S data sampled at moorings of the Line W array, deployed in April 2004 between Cape Cod and Bermuda. We use profiles from the third mooring located at 69.11^∘ W, 38.51^∘ N.

3.1 The multilayer perceptron

The neural network used is a multilayer perceptron (MLP), which is a class of feed-forward artificial neural networks (Rosenblatt, 1961; Gardner and Dorling, 1998). The MLP guesses the non-linear relation between inputs and outputs, through one or more hidden layers with many neurons stacked together. The learning mechanism that allows the MLP to iteratively minimize the loss function is called backpropagation. We keep the architecture simple with only two layers of 256 neurones each. Dropout is used as a regularization method to reduce overfitting and improve generalization (Srivastava et al., 2014). Activation functions are a rectified linear activation function (ReLU) for the hidden layers, linear for T and S output, and a sigmoid for the MLD output. We tested more complex architectures (additional layers, convolutional layers, bottleneck architecture) but could not improve the accuracy of the results with them. A simple architecture is advantageous for its lower computation time. The input consists of 12 values, listed in Table 1: latitude, longitude, day of the year (cosine and sine), bathymetry, MDT, SST, SLA, four geostrophic velocities U, V and both anomalies. Inputs are linearly interpolated at the in situ location of each profile. The outputs are predictions of three vectors of 51 depth levels (temperature, salinity and the MLD mask). The depth levels are the ones presented in the data section, on which the CORA profiles are interpolated.

The dataset is split (randomly with no replacement) as follows: 20 % of the profiles are set aside for testing. In the remaining 80 %, we use 80 % as training data and 20 % as validation data. Validation data are used to avoid overfitting by assessing the performance of the trained model after each epoch (one epoch sees all the training data). Be aware that the train and test data are not truly independent: the selection is random without accounting for spatial and temporal autocorrelation. Once the training of the NN is finished, we select the model with the best performance on the validation dataset. We then run this model on the test dataset (data not seen in training or validation) to confirm the good generalization of the model (i.e., training and test errors are similar, Fig. 3). Given a NN architecture with good generalization properties, we retrain a NN using all 67 767 profiles (Fig. 3, orange) as training data. The training of one model takes ∼20 min on an eight-core CPU with 32 Go of RAM.

Figure 3Normalized root mean square error (nRMSE) between temperature (a) and salinity (b) observed (CORA) and predicted profiles (Glorys12, Armor3D, NN and OSnet). The normalization is done with the standard deviation of the observed temperature and salinity by depth. The upper panels are a zoom of the first 100 m of the full-depth lower panels. The Glorys12 (blue) and Armor3D (purple) profiles are collocated with the CORA profiles, and the error is calculated between these subsamples. The NN profiles are only predicted with the NN, without adjustment of MLD, for 15 trained datasets (dark grey) and 15 test datasets (light grey). NN full (orange dotted) corresponds to the predictions using the full dataset (test + train) and is averaged for 15 models (bootstrap). Finally, the OSnet profiles (orange) are predicted with a NN bootstrapped 15 times and the MLD adjustment is performed, which slightly increases the error at the surface.

Download

To further improve the prediction performance and assess the associated confidence intervals, we exploit a bootstrapping scheme (Breiman, 1996). More precisely, we bootstrap the training procedure 15 times using a different initialization and training dataset each time. Indeed, because of the instability of the prediction method, bootstrapping can give substantial gains in accuracy. Overall, given 15 trained models, we compute the mean T, S and K profiles for all input data and their standard deviation (Fig. 3, grey). The latter delivers an estimate for the confidence interval. This bootstrap method reduces the estimation bias. Finally, the T−S prediction is generalized on a daily $\mathrm{1}/\mathrm{4}$^∘ horizontal grid. The spatial resolution of the input data (Table 1) is unified to $\mathrm{1}/\mathrm{4}$^∘ in longitude and latitude by a nearest neighbor interpolation method. This produces T and S fields with 51 depth levels from the surface to 1000 m for each day between 1 January 1993 and 31 December 2019. It is freely available (Pauthenet et al., 2022a).

3.2 Prediction of the mixed-layer depth

We define the mixed-layer depth H with a density deviation from the surface method, as proposed by de Boyer Montegut et al. (2004). It is the depth at which the potential density referenced to the surface, σ₀, exceeds by a threshold of 0.03 kg m⁻³ the density of the water at 10 m, σ₀ ($z=-H$) =σ₀ ($z=-\mathrm{10}$ m) +0.03 kg m⁻³. This definition is chosen for its simplicity of application but tends to overestimate deep winter MLDs compared to the more sophisticated hybrid algorithm of Holte et al. (2017). These regions of deep winter MLD are rarely observed in our dataset (1145 profiles with MLD deeper than 300 m, i.e., 1.7 % of the dataset), comforting our choice of using a simple density threshold. For the NN, the mixed layer is represented in the form of a unitless profile K of size Z=51 that is filled with zeros between the surface and the mixed-layer depth H and with ones from H to −1000 m.

This formulation allows the NN to give an estimation of the gradients around the MLD instead of a single depth value (Fig. 4d). The resulting mask K is also convenient for the MLD adjustment performed on predicted profiles.

Figure 4Example of a profile (orange) sampled on 3 March 2012 at 35.55^∘ W and 23.53^∘ N, truncated at 300 m deep. We display temperature (a), salinity (b), potential density (c) and profile K that is a mask of MLD (see Eq. 1). The T−S profile predicted by the NN is in red and the adjusted using the ${\widehat{K}}^{*}$ profile (OSnet) in green (see Sect. 3.3). The green and red bands are the confidence interval for each profiles, i.e., the standard deviation of the 15 bootstrapped models. MLDs observed, predicted and adjusted are shown with dotted horizontal lines. SST is added with a purple dot and a horizontal bar for its mapping uncertainty.

Download

3.3 Adjustment of the mixed layer

We identified two types of error in the direct prediction of T−S profiles. Firstly, the MLD predicted by the $\widehat{K}$ profile has a better accuracy (MLD RMSE of 40 m) than the MLD computed from the T−S profiles directly (MLD RMSE of 50 m). The latter is systematically too shallow due to unrealistic T−S excursions on the vertical in the mixed layer, causing the density threshold to be reached too close to the surface. These sharp variations of T and S in the mixed layer also create density inversions. Secondly, the gradients of the layer under the MLD are systematically underestimated compared to the observed profiles. The mean and standard deviation of gradients of σ₀ over a 200 m-thick layer under the MLD is 1.33±0.9 kg m⁻³ for the observed profiles and 1.24±0.8 kg m⁻³ for the predicted ones. The presence of strong gradients under the MLD has been documented (Johnston and Rudnick, 2009), and the profile $\widehat{K}$ seems to contain this information. Indeed, $\widehat{K}$ is a sigmoid-like profile (Fig. 4d), and its vertical gradients are proportional to the T−S gradients around the MLD. The summer $\widehat{K}$ profiles have sharp vertical gradients compared to the winter ones (not shown), which is coherent with the seasonal cycle of the transition layer thickness (Johnston and Rudnick, 2009).

We thus choose to apply an MLD adjustment to the predicted profiles, in the same spirit as convective adjustment schemes are used in numerical hydrostatic models (Madec, 2015). We want to weight the vertical gradients of $\widehat{T}$ and $\widehat{S}$ by $\widehat{K}$ in order to reduce the gradients of $\widehat{T}$ and $\widehat{S}$ in the mixed layer and increase the gradients of $\widehat{T}$ and $\widehat{S}$ just below the MLD ($\widehat{K}=\mathit{\lambda }$) while keeping the deeper gradients unchanged. We first modify the $\widehat{K}$ into a new mask ${\widehat{K}}^{*}$ as follows:

with λ=0.57 the value of K corresponding to the MLD. The calibration of λ is done by a cross-validation procedure according to the estimation bias between $\widehat{T}$ at the sea surface and the SST value (Fig. 5). In other words, λ=0.57 allows us to adjust the MLD while keeping null the mean difference between $\widehat{T}$ and SST (green in Fig. 5). We expect this value (λ) to be specific to our region and NN parameterization. It would likely require a new calibration for another study.

Figure 5Density distribution of the difference between the in situ surface temperature and salinity and the remote sensing SST and SSS. The predicted profiles (blue) correspond to the profiles produced by the NN alone. The adjusted distribution with $\widehat{K}$ is in red (also red in Fig. 4d). The adjusted distribution with ${\widehat{K}}^{*}$ is in green and corresponds to the OSnet product (i.e., NN + adjustment), also shown in green in Fig. 4d.

Download

After computing the ${\widehat{K}}^{*}$ profiles, we reconstruct iteratively T−S profiles with the gradients multiplied by ${\widehat{K}}^{*}$, starting from the bottom value (at z=1000 m) where ${\widehat{K}}^{*}=\mathrm{1}$ (because the deepest MLD never reaches 1000 m, the maximum observed in the CORA dataset is H=628 m for our region). On a predicted temperature profile $\widehat{T}$ (the same is applied to salinity), the adjusted temperature profile ${\widehat{T}}^{*}$ is computed as follows:

and we retrieve the temperature profiles iteratively along depth by starting from the bottom $z+\mathrm{1}=\mathrm{1000}$ m, where ${\widehat{K}}_{z+\mathrm{1}}^{*}=\mathrm{1}$ and ${\widehat{T}}_{z+\mathrm{1}}^{*}={\widehat{T}}_{z+\mathrm{1}}$:

Figure 5 presents the overall bias of surface $\widehat{T}$ and $\widehat{S}$ relative to SST and SSS. If we adjust gradients with a $\widehat{K}$ profile, the surface temperature is systematically too warm compared to SST (Fig. 5a, red). Now, if the adjustment is also increasing the gradients under the MLD (${\widehat{K}}^{*}$), the surface temperature bias is null (Fig. 5a, green). This supports our choice to amplify gradients just under the MLD (${\widehat{K}}^{*}>\mathrm{1}$) and to reduce them in the mixed layer (${\widehat{K}}^{*}<\mathrm{1}$). Note that the direct prediction of temperature at the surface (Fig. 5a, blue) is more accurate compared to SST than in situ observations, because OSnet learns from SST. The salinity difference relative to SSS is too large for the adjustment to cause a significant issue (Fig. 5b). Still, the adjusted salinity profiles with K predicted create a fresh bias, and the use of ${\widehat{K}}^{*}$ corrects that.

A good example of profile prediction and adjustment is presented in Fig. 4. In this case the adjustment corrects perfectly the MLD estimate, from 30 m predicted (red) to 133 m adjusted (green). It reduces the T−S gradients above the MLD estimated from the K profile (K=0.5) and increases the T−S gradients just under this MLD. The decrease in T−S gradients in the mixed layer also removed the density inversion (Fig. 4c).

The adjustment proposed here reduces the variance of $\widehat{T}$ and $\widehat{S}$ above the MLD, reduces the number of density inversions, improves the predictions of the MLD (see Table 2) and increases the gradients in the transition layer under the MLD. For the adjusted profiles, the mean and standard deviation of gradients of σ₀ over a 200 m-thick layer under the MLD are 1.38±0.9 kg m⁻³, which is closer to the observed values of $-\mathrm{1.33}\pm \mathrm{0.9}$ than the predicted values of 1.24±0.8 kg m⁻³. The large summer gradients are especially well retrieved (not shown).

Table 2Metrics of accuracy for predictions of Armor3D, Glorys12, a neural network (NN, i.e., OSnet without MLD adjustment) and OSnet compared to the in situ CORA profiles. The σ₀ inversions larger than 0.01 kg m⁻³ are counted. The Armor3D and Glorys12 statistics are computed on the subsampled products at the locations of the profiles of CORA.

^* All the MLDs are computed with the density criterion of 0.03, except for Armor3D, for which a different criterion is used to bypass their density inversion issues.

Download Print Version | Download XLSX

An alternative idea to solve the density inversion issue is to constrain the NN to predict profiles that are hydrostatically stable. The physical relationship can be implemented in the NN to enforce consistency on the predictions (Karpatne et al., 2017). This can be done by modifying the loss of the NN and penalizing the predictions with density inversions (Appendix A). This solution is elegant and allows us to predict directly profiles without density inversions. We provide this alternative approach here for the record of a negative result with regard to the design of such a NN. Indeed, profiles predicted with this custom loss still have poor MLD estimates compared to $\widehat{K}$ profiles and hence still need a posteriori adjustment. The modified loss (Appendix A) is not needed in this case because the MLD adjustment presented above happens to remove density inversions efficiently.

3.4 Explainability of the neural network

Explaining the predictive skills of the neural network is key to interpreting the prediction and strengthening trust in the model. It is also a useful tool in the development phase. Here, it gives insights into the relationship between surface data and in situ profiles. In this section, we use a game-theoretic approach to retrospectively estimate the relative importance of each input for each output. The algorithm, called SHapley Additive exPlanations (SHAP) (Lundberg and Lee, 2017), is a unified framework combining state-of-the-art methods to explain deep neural networks. It is based on a method called Deep Learning Important FeaTures (DeepLIFT) (Shrikumar et al., 2017) and Shapley values. DeepLIFT is a method for computing the effect of changing the original input to a reference value (uninformative background value for the input). The change in the output is representative of the importance of the input for predicting the output. Shapley regression values (Shapley et al., 1953) represent the impact of an input on the output by removing it from the input set and retraining this model with the subset of inputs. This being computationally expensive, it is possible to obtain an approximation of the effect of removing a variable from the model by integrating over samples from the training dataset using the Shapley sampling values (Štrumbelj and Kononenko, 2014). It produces an “importance” value for each particular prediction. The importance value is positive or negative, indicating the direction in which the input influences the output relative to the averaged output. The SHAP algorithm being computationally expensive, it was not possible to run it over the full dataset. After some tests, we found that 300 random samples were representative enough to obtain stable results for the average feature importance across the entire dataset.

4.1 Accuracy of the predictions compared to observations

Table 2 presents several metrics to evaluate OSnet, Armor3D and Glorys12 relative to CORA. Each dataset is predicted or subsampled at the locations of CORA profiles (in longitude, latitude, depth and time). The first two lines indicate the number and size of vertical density inversions. They indicate that the MLD adjustment of OSnet suppresses almost all density inversions, from 17.3 % to 0.32 % of profiles. Meanwhile, Armor3D has about 50 % of profiles with density inversions and Glorys12 almost none (0.01 %). Regarding the amplitude of density inversions, the MLD adjustment suppresses sufficiently large inversions and decreases by 1 order of magnitude the mean amplitude of inversions, from 10⁻³ to 10⁻⁴ kg m⁻³.

The root mean square error (RMSE) of the MLD (Table 2) indicates that the MLD adjustment improves the MLD RMSE of OSnet from 50 to 40 m, which is the same order of magnitude as Glorys12 (38.6 m) and Armor3D (39.4 m). Note that the MLD of Armor3D is computed with a different criterion to bypass density inversions. Armor3D uses the minimum of temperature and density threshold equivalent to a 0.2 ^∘C decrease from the surface. The MLD of Armor3D computed with a density criterion of 0.03 kg m⁻³ yields a RMSE of 62.6 m. Finally, the global errors of temperature, salinity and density indicate that Armor3D profiles are the closest to the observed profiles. The metric is the normalized RMSE (nRMSE), which is the RMSE between predicted and observed profiles divided by the standard deviation of the observed profiles. It is a ratio of error compared to the variability observed. OSnet has a smaller nRMSE than Glorys12, and the MLD adjustment slightly increases the temperature nRMSE at the surface. Figure 3 reveals the vertical distribution of nRMSE. OSnet gives a more accurate prediction at the surface compared to both Armor3D and Glorys12, but Armor3D is closer to observations for the rest of the water column. Overall, the nRMSE of OSnet predictions is the same order of magnitude compared to other products, and it does not contain significant density inversions.

4.2 Temperature and salinity maps

Let us examine a daily map of temperature and salinity at $\mathrm{1}/\mathrm{4}$^∘ resolution (Fig. 6). We chose a date in the pre-Argo era to illustrate the generalization skill of the OSnet product. All the maps reveal coherent horizontal structures. At the surface, the warm Gulf Stream detaches from Cape Hatteras and meanders further east, transforming into the North Atlantic Current (Fig. 6a, b). The surface confidence intervals are maximum for the cold and fresh waters near the edge of the continental slope and inside the cold and fresh core eddies and meanders (Fig. 6c, d). On average, confidence intervals highlight cold waters north of the Gulf Stream (Fig. 7a, b), which is consistent with the error of prediction maps presented in Buongiorno Nardelli (2020). This could be due to the lack of profiles containing these cold waters in our dataset. At depth (1000 m in Fig. 6e, f), the signature of large eddies is visible, associated with a maximum of the confidence interval again (Fig. 6g, h). The salty and warm Mediterranean Overflow Waters are seen in the southeast of the region. The average confidence interval at 1000 m is maximum along the Gulf Stream and its meanders, rather than in waters north of the Gulf Stream like at the surface (Fig. 7c, d). It corresponds to areas with the largest variability (Gaillard et al., 2016; Forget and Wunsch, 2007). Note the different color scale: the maximum confidence interval values at 1000 m depth are twice as small for temperature and 5 times smaller for salinity compared to the surface.

Figure 6OSnet temperature and salinity maps for 7 January 1993, at the surface (a, b) and at 1000 m (e, f). Their respective confidence intervals are displayed too, i.e., the standard deviation of the 15 bootstrapped models (c, d, g, h).

Figure 7Time average maps of the confidence intervals for surface and 1000 m temperature and salinity, i.e., the standard deviation of the 15 bootstrapped models.

4.3 MLD maps

To illustrate the quality of the predicted MLD of OSnet, we show MLD maps for a given day (5 January 2018) in Fig. 8. We picked a winter day to emphasize deep MLD areas. The direct prediction of the NN (Fig. 8a) has shallow patches in a few places that are due to density inversions. The density threshold is met too shallow due to these artifacts in the water column (see the profile example in Fig. 4). The MLD adjustment corrects these shallow patches, and the MLD field of OSnet looks more consistent and more similar to Glorys12. The OSnet MLD does not exhibit any very deep patch (MLD >300 m). These deep MLD events are rarely observed in CORA (1.7 % of the profiles have a MLD deeper than 300 m) but are often present in the MLD fields from Glorys12 (Fig. 8). They occur in warm meanders and eddies of the Gulf Stream. The MLD field of Armor3D (Fig. 8d) is for the week that contains 5 January 2018. It has several patches of either shallow or deep MLD (i.e., 28^∘ N, 52^∘ W or 48^∘ N, 32^∘ W) which look very sharp compared to OSnet and Glorys12. These patches might be regions around observed profiles for the given week, and the optimal interpolation of Armor3D is overfitting the profile at the expense of the general property field coherence. The Armor3D MLD is computed with a different criterion to bypass their density inversion issues. Still, some patches have no values where the criterion could not be matched (Fig. 8d).

Figure 8MLD maps for 5 January 2018. Panel (a) shows the result of the NN alone (step 1 from the schematic Fig. 2) and panel (b) is the final OSnet product. Panels (c) and (d) are the MLDs of the TS profiles of Glorys12 and Armor3D. Armor3D is a weekly product, so the profiles here are for the first week of January 2018. All the MLDs are computed with the density criterion of 0.03, except for Armor3D, for which a different criterion is provided to bypass their density inversion issues. The shelf break is traced in black with the bathymetry contour of 1000 m. The maximum of the color bar is set by the maximum of OSnet.

Monthly MLD averages are presented for March and August in Fig. 9. The averages in 1^∘ × 1^∘ boxes for the in situ profiles (Fig. 9a, d) are compared to OSnet (b, e) and Glorys12 (c, f). The three estimates are in good agreement with each others and with climatologies not shown here (Holte et al., 2017; Sallée et al., 2021). The main structures are well respected with a large winter patch of deeper MLD extending between the Gulf Stream and the subtropical gyre. In winter, vigorous air–sea fluxes extract heat from the ocean and erode the superficial stratification. This process activates convective mixing and deepens the mixed layer, ventilating and creating the Eighteen Degree Mode Water (Speer and Forget, 2013; Maze and Marshall, 2011). In summer the near-surface water warms and caps the mode water layer. The summer MLD is shallower everywhere with a slightly deeper signature in the core of the Gulf Stream, as it separates from the coast at Cape Hatteras. A deeper summer MLD is also found south of ∼30^∘ N along the equatorward edge of the subtropical gyre (Fig. 9d, e, f), a feature also observed in the climatology of de Boyer Montegut et al. (2004). This tropical summer MLD deeper than 30 m is marked by the trade winds (Stramma and Schott, 1999). The large permanent anticyclonic “Mann eddy” (Mann, 1967; Rossby, 1996) is clearly visible as a deep mixed-layer patch at 43^∘ N, 42^∘ W (Fig. 9). A region of the deeper MLD is also visible along the North Atlantic Current, deeper in OSnet than in Glorys12.

Figure 9Maps of the monthly mean of the MLD defined with a density threshold of 0.03 kg m⁻³ for (a, d) CORA T−S profiles averaged by bins of 1^∘ (b, e) OSnet and (c, f) Glorys12. The months of March and August are representative of the seasonal range of the MLD for the region. Note the different color-bar range for each month. The maximum of the color bar is set by the maximum of OSnet. The shelf break is traced in black with the bathymetry contour of 1000 m.

4.4 Importance of each input for the reconstruction

In Figs. 10 and 11 we present the absolute values of the relative importance of each input, on each output, averaged over depth, over the 300 test profiles and over the 15 bootstrapped models. Error bars correspond to the standard deviation of the 15 models. To be comparable, the importance values of the inputs are normalized so that the sum per output is equal to 1. Figures 10 and 11 give a general overview of what the NN uses for the predictions. These importance values can also be displayed for a specific profile or by depth levels, seasons, or geographical regions, providing insights to elucidate the behavior of the NN. The main result here is that SST is the main driver for estimating T, S and MLD profiles (Fig. 10). As expected, it is especially important for predicting surface temperature. MDT is the second-most important variable, and it is the most important deeper than ∼200 m (Fig. 11). The latitude, longitude, SLA and day of the year equally concur to explain the predictions. The rest of the input variables, i.e., bathymetry and the surface geostrophic currents derived from SLA, are smaller contributions to the predictions. Even if they have small contributions on average, they can be important for a specific profile. The cosine of the day of the year is very important for the prediction of the MLD (Fig. 10), probably because it is in phase with the MLD seasonal cycle, while temperature and salinity cycles are in phase with the sine of the days (not shown). Still, it means that the day of the year alone drives ∼29 % of the MLD predictions, which is equivalent to the importance of SST for MLD predictions (∼29 % too).

Figure 10Relative importance of each input for each output, averaged by depth. The inputs (x axis) are sorted by importance for the temperature to have the largest importance on the left of the plot. The cosine of the day of the year is more important than the sine for the MLD prediction because the cosine is in phase with the seasonal cycle of the MLD.

Download

Figure 11Relative importance of each input for each output, averaged by depth. The areas are sorted by variance to have the input with the largest difference of impact by depth to the right of each panel.

Download

4.5 Time series of surface properties

To assess the accuracy of OSnet through time, we analyze the time series of the spatially averaged surface temperature and compare it to the observed SST time series as well as the Glorys12 and Armor3D products. We average data over the region after removing values on the shelf (bathymetry >1000 m). The long-term trends are obtained by applying a seasonal-trend decomposition based on local regression (STL) (Cleveland et al., 1990). STL is a filtering procedure that extracts three components: (i) the variations in the data at the seasonal frequency (Fig. 13), (ii) the low-frequency variation together with nonstationary and long-term changes (Fig. 12) and (iii) a remaining high-frequency component.

Figure 12Nonseasonal low-frequency time series of surface temperature (a) and salinity (b) averaged over the region excluding the shelf shallower than 1000 m. It is extracted with a seasonal-trend decomposition. The linear trends of temperature are displayed with dashed lines, and their slope and R² are in the legend of panel (a). The grey shaded areas are the SST mapping error in (a) and the SSS random error (b). The green shaded areas are the OSnet confidence intervals, i.e., the standard deviation of the 15 bootstrapped models.

Download

Figure 13Seasonal variation of the mean surface temperature (a), surface salinity (b) and MLD (d) for OSnet, Armor3D and Glorys12 compared to remote sensing SST and SSS (black). The difference between SST and the surface temperature is displayed in panel (c), with the direct prediction of the NN hence without MLD adjustment in red. It is averaged over the period 1993–2019 excluding the shelf shallower than 1000 m, except for the surface salinity because SSS only ranges from 2010 to 2019. The bands or errors are the standard deviation over the time period.

Download

OSnet follows best the long-term trend of SST with a linear trend of 0.197 ^∘C/decade, close to the SST trend of 0.190 ^∘C/decade (Fig. 12a). For comparison, the surface-averaged temperature of Armor3D is warmer in the pre-Argo years, probably due to their background global average that is mostly composed of Argo floats. This is a validation that OSnet generalizes well and is not biased towards the recent years of in situ observations. Regarding surface salinity long-term trends, no significant trend is captured over the 1993–2019 period (Fig. 12b). There is no clear agreement between the different datasets, except in the last period, 2010 to 2019, where all averages increase like the SSS signal. Armor3D mean surface salinity drops significantly during the last 2 years, 2018 and 2019, out of the SSS error margin. Note that OSnet does not include the areas over the continental shelf and does not predict deeper than 1000 m. A significant part of the climatic signal takes place in the coastal regions (Ezer et al., 2013; Davis et al., 2017), and an improvement of OSnet would be to deal with profiles of different lengths in order to include these regions.

The mean seasonal variation of surface temperature, salinity and MLD is presented in Fig. 13. We compute it by averaging the signal by day of the year (week of the year for Armor3D). The surface temperature seasonal signal is well reproduced for each dataset, as expected considering that SST is included in the input of OSnet and used to produce both Armor3D and Glorys12 (Fig. 13a). A close observation of the curves in Fig. 13a shows that the seasonal cycle of OSnet surface temperature is too cold between May and September (Fig. 13c). This is due to the MLD adjustment, since the direct prediction of surface temperature gives a very precise seasonal cycle of SST (red line in Fig. 13c). Armor3D seasonal temperature is warmer from August to April, which could be a bias caused by the undersampling of the pre-Argo, colder years (Fig. 13c). This bias is also present in the temperature time series (Fig. 12a).

The surface salinity seasonal signal is noisier, in part because it presents inter-annual variations that are the same order of magnitude as the seasonal variations (Fig. 13b). We compare 2010–2019 because it is the only period available for the SSS product. OSnet surface salinity is the closest to SSS compared to the predictions of Armor3D and Glorys12, which are both fresher than the observed SSS. We observe a delay in the SSS seasonal variations: it is fresher by almost 0.1 psu from January to March. We discuss this in Appendix B and Fig. A1.

Finally, the seasonal cycle of MLD (Fig. 13d) in the region is asymmetrical, with slow deepening from summer to winter and fast shoaling during the early spring. This asymmetry is expected since buoyancy loss at the surface leads to convective mixing (hence deepening the mixed layer requires buoyancy loss over an ever-increasing water column depth), while buoyancy gain directly leads to a stable stratification and shallow mixed layer (e.g., Taylor and Ferrari, 2010). OSnet MLD compares well to the MLD computed on Glorys12. We do not present the MLD of Armor3D here because it is computed with a different criterion. The winter MLD variance is larger in Glorys than OSnet, which is also observed on daily maps of MLD. Events of deep MLDs are not represented in OSnet (Fig. 8).

5.1 Comparison to observed data

One major feature of OSnet is the possibility of estimating T and S profiles at any location, given that surface data are available. Here we predict T−S at the location of a mooring of the line W (Fig. 14) and along the hydrographic section AT20 (Fig. 15). Temperature and salinity at 1000 m are plotted in comparison to mooring data for a period of 3 years (10 May 2004 to 11 March 2007). OSnet corresponds well to observations but with a slight warm and fresh shift in the first year (Fig. 14). A warm-core eddy crosses the location of the mooring in 2006 (see the SLA map in Fig. 14a), and its warm and salty deep signature is well captured by OSnet. Smaller warm and salty spikes appear in October 2004 and July and November 2005 but are not visible in the mooring data. They correspond to warm meanders of the Gulf Stream revealed by the SST and SLA at these three periods (not shown), causing deep T−S changes in OSnet.

Figure 14OSnet prediction of temperature and salinity at 1000 m (blue) at the location of a mooring of line W3 (black, 69.11^∘ W, 38.51^∘ N). The map of SLA in the temperature panel is for 25 April 2006, when a warm-core eddy went through the mooring location. These mooring data are not included in the learning dataset of OSnet.

Figure 15Hydrographic section AT20 (a, b, c) sampled along the 52.3^∘W meridian by the research vessel Atlantis from 1 to 5 May 2012. Temperature and salinity profiles are estimated with OSnet (e, f) at the exact locations of the sampled CTDs by interpolating the input data at those locations. Glorys12 profiles (h, i) are collocated in time and space with the CTD profiles. The SSH in panel (a) is averaged over the sampling period of the section. The MLDs in panel (d) are computed with the density threshold of 0.03. Salinity segments at 200 and 800 m are plotted along latitude in panel (g) with confidence intervals in blue bands around the mean value for OSnet, i.e., the standard deviation of the 15 bootstrapped models. The contour intervals for the plot section are 0.5 psu from 33 to 37 psu for salinity and 2 ^∘C from 4 to 24 ^∘C for temperature.

We compare the OSnet T−S structure along the hydrographic section AT20 sampled in May 2012 by the research vessel Atlantis (Fig. 15). The OSnet prediction is done at the exact locations of the conductivity–temperature–depth (CTD) profiles, interpolated linearly on the maps of input data (Table 1). The comparison is also done for Glorys12 but by collocating the profiles with a nearest neighbor method. Both Glorys12 and OSnet predictions are coherent, but we note two specific biases in Glorys12. First, Glorys12 displays a deep patch of MLD around 41^∘ N, north of the Gulf Stream, that is not observed by the CTDs nor predicted by OSnet (Fig. 15d). This deep MLD could be due to the nearest neighbor selection of the profile that is not exact in the case of Glorys12, or it could be an artifact of their model. Indeed, deep patches of MLD are also visible on daily MLD maps of Glorys12 but are absent on the OSnet daily MLD maps (Fig. 8). Second, the salinity reconstruction of Glorys12 differs more from the observed data north of the Gulf Stream (Fig. 15g). OSnet performs very well along this section in comparison to Glorys12.

5.2 OSnet to explore theoretical inputs

Since OSnet is very easy and fast to manipulate to make predictions, it can be used to make predictions using theoretical inputs. To illustrate this, we seek the interior signature of eddies detected with altimetry. We make two predictions of temperature and salinity for a section across two particular eddies observed on 6 October 2006 (Fig. 16a): one prediction is based on all observational inputs (Fig. 16b, c), and the second prediction is made by removing the eddy signature in SLA; we simply set it to 0 (Fig. 16e, f, g). The interior temperature and salinity anomalies associated with the eddies are obtained by difference (Fig. 16h, i). Anomalies are the largest at depth around $\mathrm{400}/\mathrm{500}$ m with amplitudes of a few degrees per meter of SLA. These are reasonable amplitudes and structures for the region (Castelao, 2014) and illustrate how OSnet could easily be used to extract more knowledge than the standard realistic T−S predictions on a grid.

Figure 16Prediction of T−S profiles for a section across two meanders of the Gulf Stream (b, c) and for a simulated SLA flattened to zero (e, f). The meanders are visualized with maps of SLA (a) and SSH d, which are computed by adding MDT and SLA for 6 October 2006. The anomaly between the T−S sections for the observed SLA and the simulated SLA is displayed in panels (h) and (i). The SSH contour in panel (d) is 0.1 m to represent the Gulf Stream.

5.3 Potential improvement of the method

Several improvements of the method could be made in the future. One important limitation is the current NN architecture that prevents predictions of profiles of different lengths. It constrains the analysis to use a fixed maximum depth and thus prevents profiles from being kept shorter than that maximum depth. A solution would be to develop a custom loss that can deal with empty variables. Another way would be to add depth as an input variable like in Buongiorno Nardelli (2020). In that case, we could predict properties over shallow bathymetry and also deeper than our rather arbitrary 1000 m limit.

OSnet produces coherent horizontal and temporal patterns even though each profile is predicted independently. Yet we wonder how using horizontal and temporal surface gradients as inputs could improve predictions, especially in frontal regions. To test this, we would need to work with three-dimensional (long–lat–time) patches of input data for each profile location and to build a different NN architecture (e.g., Ouala et al., 2018; Jouini et al., 2013; Tandeo et al., 2013; Lguensat et al., 2018) that takes patches of data as input and profiles as output. Convolutional neural architectures accounting for irregularly sampled space–time observations might also be appealing (Fablet et al., 2021). The expected result would be sharper fronts, as we observe that fronts from OSnet are smoother than in observations, e.g., Fig. 15. We wonder whether the NN could learn from the temporal surface gradient to anticipate vertical changes in stratification.

Prediction intervals (PIs), also called “coverage probability”, could be computed as a complement of confidence intervals (Khosravi et al., 2011). While the confidence intervals (Fig. 7) give the range of variation of a set of NNs with different initializations and training datasets, the PI gives the probability of the observation associated with the prediction being in a range of values. The PI can be obtained by adding output variables to the proposed architecture. Here, the PI could be represented as temperature and salinity profiles around the prediction, representing the 95 % interval. It means there would be a 95 % probability of the true value associated with the prediction lying within the interval.

Given the very promising results of our study, obvious future work would be to apply OSnet in other more challenging regions, with fewer data or more complicated vertical structures and different dynamics. Also, the 3D geostrophic velocities of the OSnet gridded product could be estimated using the thermal wind equation combined with surface altimeter geostrophic currents (Mulet et al., 2012). Finally, since we found that OSnet correctly captures the SST warming trend (Fig. 12a) and mesoscale structures, it would be interesting to apply OSnet to other boundary currents and to compare the resulting OHC estimates to previous reconstructions (e.g., Cheng et al., 2017). Other global ocean indicators such as ocean freshwater content or steric sea level could be investigated as well.