The Argo profiles are obtained from Global Sea Ocean Argo Scatter Dataset (V3.0) (Liu et al., 2021) provided by the China Argo Real-time Data Center (ftp://ftp.argo.org.cn/pub/ARGO/global/). This dataset collects more than 2.3 million temperature and salinity depth profiles observed by over 15,000 automated profiling buoys deployed in the global ocean by international Argo member countries during July 2000 through June 2020. Data inside the region of 95° W-135° W and 5° S-45° N for the period 2011-2021 are considered here. Profiles with depths exceeding 700 m are selected, and the profiles are interpolated through a spline into the same 71 vertical levels extending from the surface (5 m) to 705 m (vertical step is 10 m).

The SST data of this study are created by the OSTIA (Operational SST and Ice Analysis) system, using re-processed ESA SST CCI, C3S EUMETSAT and REMSS satellite data and in situ data from the HadIOD dataset, distributed through the Copernicus Marine Environment Monitoring Service (CMEMS, http://marine.copernicus.eu/services-portfolio/access-to-products/, product_id = SST_GLO_SST_L4_REP_OBSERVATIONS_010_011). This product provides daily maps of the SST and SST uncertainty on a global regular grid at 0.05° resolution, which are stored using the netCDF format using the Group for High Resolution SST specification.

The SSS data is a Level 4 products on a 0.25 degree spatial and 4-day temporal grid produced by the International Pacific Research Center (IPRC) at the University of Hawaii at Manoa in collaboration with the Santa Rosa Remote Sensing System (RSS) in California in conjunction with observations from NASA’s Aquarius/SAC-D and Soil Moisture Active (SMAP) satellite missions. The product is a continuous, consistent multi-satellite SSS data obtained by optimal interpolation with a 7-day decorrelation time scale (Melnichenko et al., 2016). Their mean root mean squared difference from globally synchronized in situ data is about 0.19 psu and the product bias is about zero.

The altimeter sea level anomalies with daily and 0.25°×0.25° resolutions are provided by Sea Level TAC (Thematic Assembly Centre, https://resources.marine.copernicus.eu/product-detail/SEALEVEL_GLO_PHY_L4_MY_008_047/DATA-ACCESS). The data produced in the frame of this TAC are generated by the processing system including data from all altimeter Copernicus missions (Sentinel-6A, Sentinel-3A/B) and other collaborative or opportunity missions (e.g.: Jason-3, Saral[-DP]/AltiKa, Cryosat-2, OSTM/Jason-2, Jason-1, Topex/Poseidon, Envisat, GFO, ERS-1/2, Haiyang-2A/B/C).

The climate data used in this study are from World Ocean Atlas 2018 (Boyer et al., 2018) (https://www.ncei.noaa.gov/products/world-ocean-atlas), which is provided by the National Oceanographic Data Center (now the National Centers for Environmental Information - NCEI)(Boyer et al., 2018) (https://www.ncei.noaa.gov/products/world-ocean-atlas). The atlas is an objectively analyzed, quality-controlled collection of temperature, salinity, oxygen, phosphate, silicate, and nitrate averages based on profile data from the World Ocean Database and distributed online by NCEI. Monthly climatology fields of temperature and salinity at standard depth levels at a spatial resolution of 0.25°×0.25° were used in this study and interpolated by cubic spline onto regularly spaced vertical grids (10 m apart).

Note that all satellite data and climatology data are interpolated to the same spatial distribution as the Argo observations by bilinear interpolation in the present study. And anomalies are defined as the observation data (SST,SSS,Argo) subtracted by the monthly WOA18 data.

In this section, GPR is used to estimate the vertical profiles of TA and SA. The SAs and TAs in each level are considered as a collection of some random variables and obey a joint Gaussian distribution, defined as y = f ( x ) + ϵ , where x is input vector of predicted parameters, which are sea surface temperature anomaly (SSTA), sea surface salinity anomaly (SSSA), sea surface height anomaly(SSHA), longitude (LON), latitude(LAT) and the day of the year projected on a circle (JULD). The prior distribution of f ( x ) is assumed to be a GP given by f ( x ) ∼ G P ( 0 , κ ( x , x ) ) , and κ is the semi-positive definite kernel. ϵ ∼ G P ( 0 , χ 2 ) denotes the Gaussian noise term, here χ is the standard deviation. Based on the historical Argo profiles and remote sensing data, we can collect n TA observations or SA observations at depth z to form the observation set y = { y 1 , y 2 , … , y n } . Corresponding to each observation, we have also collected a set of input parameters X g = [ x 1 ∣ x 2 | … | x n ] ∈ ℝ d × n , here x i = ( S S T A i , S S S A i , S S H A i , L O N i , L A T i , Then the prior distribution of y can be given as:

GPR estimates the posterior distribution of an unknown quantity under the assumptions of a Gaussian process and the likelihood of a normal distribution. Specifically, from the above assumptions, it is known that the joint probability distribution of the estimated value f * at the new input parameter x * and the existing observation y is a joint normal distribution of the following form:

where K * = κ ( X g , x * ) , K * * = κ ( x * , x * ) . Combining the prior assumption and the likelihood, the posterior distribution of f * can be derived from the maximum likelihood method (Williams and Rasmussen, 1995),

The reconstruction of TA and SA profiles using GPR is carried out in the python library (GPy, 2012), and the output is a vector consisting of the TA profile and SA profile in series at the same location. The kernel function chosen in this paper is the radial basis function. The minimum and maximum of the input/output can be utilized to scale the input/output to [0,1] to eliminate the magnitude effect before the model is trained.

In summary, the derivation of the vertical structure of TA and SA using GPR involves two processes: online and offline stages. In the offline stage, historical observations are collected to build the corresponding input/output training set, from which the GPR is trained to learn the mapping of the input parameters to the TA and SA profiles. In the online stage, the corresponding TA and SA profiles are estimated from the already trained GPR for the new input parameters. The flowchart of TA and SA estimation using the GPR is shown in Figure 1 .

LSTM is an extension of the traditional recurrent neural network (RNN), which is mainly used to deal with the case of traditional RNN failure. It solves the problem of gradient disappearance and gradient explosion of traditional RNNs to a certain extent to learn long-term dependent information. Same as general RNN, it also consists of a series of repeating cells with sequential connections, but the structure of the cells is more complex. Crucially, the LSTM adds a state vector updated over time to the cells, which can selectively record the information of the system and preserve the long-term state of the system. Each cell has three inputs, the input x i of the network at the current time, the output value h i − 1 at the previous time, and the cell state C i − 1 at the last time. The input variables are passed through three gates with different functions, namely forget gate, input gate and output gate, to forget and add information to the cell state C_i at the current moment and update the output state h i at the current moment. The structure of each cell is shown in Figure 2 (Hochreiter and Schmidhuber, 1997), and the tasks of the different gates are implemented using the following equations (Hochreiter and Schmidhuber, 1996; Gers et al., 2000):

where σ denotes the sigmoid activation function, W f , W I , W C , W O are weight matrices, b f , b I , b C , b O are model biases, and * represents dot product operation.

Referring to (Buongiorno Nardelli, 2020), we also use LSTM to estimate the TA and SA profiles. TAs and SAs from different depths are considered the output states of different cells. The input to each cell is the same, a multivariate vector consisting of SSTA, SSSA, SSHA, LON, LAT, and JULD corresponding to the current position. The structure of the LSTM we used is the same as that in (Buongiorno Nardelli, 2020), i.e., a 2-layer stacked network. Each layer contains 35 hidden units, and the optimization algorithm is Adam (Kingma and Ba, 2014). Similarly, we also use max−min normalization to preprocess the data before the LSTM training.

POD, which has a wide range of applications in various fields (Liang et al., 2002; Pinnau, 2008; Chapelle et al., 2012; Singler, 2014), provides a means to obtain a low-dimensional description of the system. This method extracts a small number of modes from historical TA and SA profiles that can represent the main features of the field, which simplifies and reduces the dimensionality of the data. The temperature-salinity anomaly profiles at the same location are placed in a multivariate matrix X

where m is the number of vertical levels, n is the number of vertical profiles. In order to calculate the k modes, we only need to do the singular value decomposition of X as follows:

where U ∈ ℝ 2 m × r and Z ∈ ℝ n × r are orthogonal matrices. The matrix Σ = d i a g ( σ 1 , σ 2 , … , σ r ) contains the singular value σ 1 ≥ σ 2 ≥ … ≥ σ r , where r is the rank of X . The first k columns of U are the k modes to be computed, which we denote as

The dimensions of L i = [ L 1 i , L 2 i , … , L m i ] ⊤ and M i = [ M 1 i , M 2 i , … , M m i ] ⊤ is m × 1 . Then the corresponding TA profile can be represented by a linear combination of L i   ( i = 1 , 2 , … , k ) , and the SA profile is represented by a linear combination of M i   ( i = 1 , 2 , … , k ) , where the unknown coefficients are shared. A common method for solving the coefficients is to solve a linear system that is obtained by equating the corresponding sea surface elements inthe reconstructed anomaly profiles to sea surface observations. In this case, k = 2 is required. To increase k , other physical quantities of the ocean need to be added to X to impose constraints on the combined coefficients.

Different from the above methods, in this paper, two methods, GPR and LSTM, are utilized to learn the relationship from the sea surface parameters to the coefficients. In particular, U k is the solution of the following minimization problem

where Y k = { W ∈ ℝ 2 m × k : W T W = I k } , ‖ · ‖ F is the Frobenius norm, and the error is estimated as

From this, we can naturally compute the coefficient vector α i = U k T u i = Δ [ α i 1 , α i 2 , … , α i k ] ⊤ corresponding to the historical temperature and salinity anomaly profile u i . And, we can determine k by the following equation according to the required accuracy

where δ is user defined cut-off threshold. To retrieve the vertical structure of TA and SA from the sea surface parameters, we collect the input parameters x i corresponding to the obtained state vector

u i   ( i = 1 , 2 , … , n ) . Here x i is still a vector consisting of SSTA, SSSA, SSHA, LON, LAT and JULD. Then based on the collected dataset D t r = { ( x 1 , α 1 ) , ( x 2 , α 2 ) , … , ( x n , α n ) } , we can train GPR or LSTM to build the mapping of input parameters x to coefficients α , then the TA and SA are reconstructed as (Hesthaven et al., 2016; Guo and Hesthaven, 2019)

The corresponding flow chart is given in Figure 3 .

Further, we point out that L i and M i are actually functions of depth z . To estimate the TA and SA at any depth z * , we use cubic spline interpolation to interpolate L i and M i with respect to depth z to obtain the corresponding interpolation functions L ˜ i ( z ) and M ˜ i ( z ) satisfying L ˜ i ( z j ) = L j i and M ˜ i ( z j ) = M j i   ( i = 1 , 2 , … , k , j = 1 , 2 , … , , m ) , where z i   ( i = 1 , 2 , … , m ) denotes the vertical levels, so that the TA and SA at depth z * can be reconstructed as

To evaluate the performance of different models, we randomly selected 20% of the 11,374 Argo profiles as the test set and the remaining 80% as the training set. Different models are then employed to retrieve TA and SA profiles with the same inputs, i.e., SSTA, SSSA, SSHA, LON, LAT, and JULD. The RMSE between TA and SA profiles obtained from these models and the Argo profiles is calculated to evaluate the performance of the different methods. The RMSE of the WOA is obtained by comparing the temperature and salinity profiles obtained after interpolating the monthly mean of the climatology of WOA18 with the Argo profiles. Let δ=0.04% and we choose k=14 modes, in fact, before performing the POD, we similarly scaled the data between 0 and 1. Figure 4 shows the vertical distribution of the RMSEs estimated by the different models. The RMSEs of all proposed methods are smaller than the RMSE of the WOA18, and a similar vertical structure can be seen in the different models. The RMSEs are small at the sea surface but increase rapidly with depth, decrease rapidly after reaching the maximum, and stabilized. This may be related to the complex dynamical processes in the ocean’s upper layers and the perturbations in the mixed and thermocline layers, while the seawater is relatively stable in the deeper layers. The RMSEs of temperature reach their maximum at 105 - 115 m, and only near this depth (where the temperature variation is relatively large), LSTM shows superior performance in predicting TA profiles compared to other methods. This suggests that LSTM may have greater potential for approximating strongly nonlinear functions since, at depths where temperature changes rapidly, the relationship between temperature and depth is more complex and may have stronger nonlinearity.

However, the performance of the four methods is comparable from an overall perspective. The RMSEs of salinity reach their maximum between 55 - 65 m. Although the four reconstruction methods have comparable accuracy, GPR performs a little better than LSTM in the reconstruction of upper salinity, which may arise because the rate of change of the salinity profile is not as large as that of the temperature profile, GPR can produce sufficiently accurate approximations, while LSTM has a more complex structure and numerous parameters, making them reliant on large amounts of training data to obtain more accurate predictions. In addition, the POD combined with the regression methods does not cause a large loss in the accuracy of the estimated profiles, which is better reflected in the RMSEs of GPR and GPR-POD, as the RMSEs of both are almost the same. The selected modes are sufficient to characterize the vast majority of the temperature and salinity fields. To further measure the uncertainty of the estimates of GPR-based methods, Figure 5 shows the standard deviation of the posterior distribution of the GPR and GPR-POD predictions. We can see that both methods display similar distribution patterns, with comparable uncertainty estimates at most points, but the uncertainty of the model prediction is higher at the western boundary of the region, which can be attributed to the limited training data available in this region. This also shows a drawback of the GPR-based approach, which is better suited for interpolated predictions rather than extrapolation. Notably, the predictions of GPR-POD exhibit higher levels of uncertainty, despite comparable RMSEs between GPR-POD and GPR. This discrepancy may be attributed to the fact that more information is available in the training data for GPR, making it more confident in its predictions.

Considering the inhomogeneity and sparsity of spatio-temporal distribution of thermohaline data, all the thermohaline relationships in the test set are drawn in the same graph ( Figure 6 ), the purpose of which is to compare the predictions of different models, in which the different colors of points are to better distinguish the differences between points. The main distribution structure of T-S graph reconstructed by all methods is similar to that of Argo field, but the distribution range of points is more concentrated than that of Argo field. From the T-S graphs, it can be observed that the reconstructed results are generally weaker for the points with larger absolute values of salinity anomalies. To better visualize the predicted distribution of temperature and salinity anomalies, a histogram of the number of salinity anomalies in different intervals is displayed at the top of the T-S graph, and a histogram of temperature anomalies is displayed on the right side of the T-S graph. The histogram of salinity anomalies reveals that the distribution of LSTM-POD is more consistent with that of Argo, while LSTM underestimates the number of points of salinity anomalies in the interval [-0.5,0], which is overestimated by GPR and GPR-POD. The histograms of temperature anomalies from different models indicate that all four methods overestimate the number of temperature anomalies in the interval [0,2.5], but their overall distribution is very similar to that of Argo. In conclusion, the reconstructed methods can predict the main distribution structure of the T-S graph. Although there are slight differences in the distribution range and there may be problems of underestimating the thermohaline state, the performance of the proposed reconstructed methods is generally satisfactory. Furthermore, Figure 7 displays the scatter plots of the LSTM estimated and Argo’s temperature and salinity anomalies at a depth of 105 m. This depth is significantly impacted by the variation of the upper mixed layer, and the RMSE estimated for temperature and salinity at this depth is also large. The results indicate that the overall discrepancy between them is minimal, as the LSTM prediction preserves the fundamental characteristics of the temperature and salinity anomalies, and its pattern (positive or negative) is accurately reconstructed.

To further compare the improvement of POD on the computational efficiency of different models, we compare the running times of different models. The CPU times of LSTM and LSTM-POD are 11216 s and 1488 s, respectively, computed using a laptop with 4 Intel(R) Core(TM) i7-6770 CPU @ 3.40GHz and 8G Memory. Since GPR requires higher memory for matrix operations, GPR and GPR-POD are run using a single node with a single core in Tianhe-2, which adopts an Intel Xeon E5-2692 v2 CPU @ 2.20GHz and 64GB Memory. Furthermore, the running times of GPR and GPR-POD are 891 s and 788 s, respectively. We can see that the training time of LSTM can be significantly reduced by combining POD, as LSTM takes 7.5 times longer than LSTM-POD. However, this reduction is not apparent in GPR, this is attributed to the fact that the computational complexity of GPR-POD is Ο ( k n 3 ) , while the computational complexity of GPR is Ο ( 2 m n 3 ) . In this case, n=9099,k=14,m=71, the improvement in training time is insignificant when the number of training samples is much larger than the number of output dimensions.

To test the sensitivity of the estimated TA and SA to different input parameters, we compare GPR-POD estimates with different training inputs. Figure 8 . shows the RMSEs of temperature with various inputs, and as well as the correlation between test input parameters, the correlation between test input parameters, and test TA profiles , and the correlation between test input parameters and test SA profiles. By observing the RMSEs of temperature, we can find that the RMSE of the model without SSHA in the input is large at all depths, which indicates that SSHA plays a vital role in ocean motion and processes. In fact, changes in the subsurface layer may give rise to changes in SSH resulting from the interaction of several factors, such as heat exchange, internal thermal expansion, and ocean circulation. A rise in sea temperature leads to an increase in SSH, while a decrease in sea temperature leads to a reduction in SSH. The correlation coefficients between SSHA and TA profiles also affirm the excellent association between SSHA and TA. Notably, the correlation coefficient between SSHA and TA is the highest among all input parameters except SSTA. In addition, SSH changes under the influence of wave and wind shear. Incorporating SSHA into the model input can provide more information for predicting TA. Latitude, longitude, and time all improve the TA profile reconstruction to different degrees, showing that geographic and temporal information is helpful in predicting TA. In particular, incorporating latitude data yields a more substantial reduction in RMSE than incorporating longitude data, suggesting that latitude has a greater impact on TA reconstruction. This is because the discrepancy in solar radiation across different latitudes results in significant variations in sea temperature in the north-south direction. Therefore, latitude information is essential to improve the performance of subsurface reconstruction. Furthermore, the correlation between latitude and TA profiles is higher compared to that between longitude and TA profiles, highlighting the importance of considering latitude when predicting TA. The model without SSTA in the input has the largest RMSE in the upper layer (<100 m ), but when the depth is greater than 100 m, there is almost no difference between the model without SSTA and GPR-POD ( all parameters as input). This suggests that SSTA data mainly improve the reconstruction of upper layer TA, while deeper temperature variations are difficult to be interpreted from satellite measurements. In addition to the fact that SSSA has no significant relationship with ocean interior temperature, which can also be verified by observing the correlation between SSSA and TA profiles. Similar results can be observed in the RMSEs of salinity, where adding both latitude and SSHA data reduces the RMSEs of salinity at all depths. The relationship between ocean salinity and temperature on density is a fundamental aspect that affects changes in SSH, causing SSHA and SA profiles to have some correlation. Latitude also influences changes in ocean salinity as a result of a combination of factors, including precipitation, evaporation, and mixing of water masses. Consequently, incorporating latitude and SSHA into the prediction of subsurface salinity can significantly improve model accuracy. SSSA data play an important role in the reconstruction of upper ocean SA profiles, while they do not help much in the reconstruction of deep SA profiles, and the correlation between SSSA and SA gradually decreased with increasing depth. The RMSEs of several other models show that longitude, time, and SSTA have somewhat limited improvements on the models, where longitude and time can reduce the maximum of salinity RMSE, while SSTA has no significant improvement on the reconstruction of SA profiles.

To test the influence of climatology on the temperature and salinity reconstruction models, we considered three different inputs and outputs for LSTM and GPR. The inputs are 1) X: SSTA, SSSA, SSHA, LAT, LON, JULD; 2) X-Without-WOA: SST, SSS, SSHA, LAT, LON, JULD, and 3) X-WOA: SST, SSS, SSHA, LAT, LON, JULD, SS (climatology), ST (climatology) and the outputs are 1) X: STA, SSA; 2) X-Without-WOA: SS, ST, and 3) X-WOA: SS, ST. The RMSEs estimated by different models are shown in Figure 9 . Observing the RMSEs of the two regression models with various inputs and outputs, we can find a similar behavior, especially when the temperature at depths of less than 100 m and the salinity at depths of less than 200 m, training the models with the temperature and salinity anomaly fields can significantly improve the prediction accuracy than directly using temperature and salinity as the input and output of the model. This is because the seasonal variation signal disappears after calculating the anomaly field using the monthly mean of climatology, thus making the SS and SA easier to predict. On the other hand, incorporating monthly temperature and salinity climatology fields into the model input gives the most accurate estimates. Climatological information added to the input can effectively reduce the impact of the sea surface parameter errors. In summary, seasonal variations in seawater temperature and salinity play a crucial role in oceanic processes. The information provided by seasonal variations in temperature and salinity can help us better estimate temperature and salinity structure within the ocean.

Comparing LSTM and GPR with the same input and output, it can be found that the accuracy of GPR with temperature and salinity climatology as input is higher than that of LSTM. This is due to the inclusion of monthly climate fields of temperature and salinity in the input, which makes the model less complex. Similar to the findings in Section 4.1, GPR can provide more accurate predictions when the amount of training data is not abundant. In contrast, the performance of GPR without considering climatology is the worst, expressly, by up to 20% (∼245 m) over the GPR-WOA estimate. When only sea surface parameters are included in the input without climatology as an aid, it leads to a more complex relationship between the input and output in the model. The LSTM, on the other hand, can better approximate the strongly nonlinear function, so the GPR-Without-WOA does not perform as well as the LSTM-Without-WOA.

Since the SSS data is the 4-day temporal grid, in order to make full use of the other daily satellite and in situ observations, we use the SSTA, SSHA, LAT, LON, and JULD as model inputs in this section to train the LSTM model (named LSTM-Daily) to estimate the TA and SA profiles. Thus, there are 45,031 Argo observation profiles available, and again, 80% are randomly selected as the training set, and the remaining 20% are used to test the model’s performance. The RMSEs of the temperature and salinity estimated by LSTM-Daily are shown in Figure 10 , respectively, besides the RMSE reduction rate of LSTM relative to climatological estimates also shown in Figure 10 . It is clear that the RMSE of the estimated TA profiles is reduced by 20%−70% compared to climatology. By increasing the training data size, a more accurate reconstruction of the temperature profiles can be obtained. On the other hand, the prediction accuracy decreases for salinity at depths less than 100 due to the lack of SSS information, which again validates the importance of SSSA for SA reconstruction at a depth of less than 100 m. However, when the depth is greater than 100 m, as the effect of SSSA data disappears, higher accuracy SA estimation can be achieved again due to increased training data size.