Estimating the mixed layer depth of the global ocean by combining multisource remote sensing and spatiotemporal deep learning

ABSTRACT

Estimating the ocean mixed layer depth (MLD) is crucial for studying the atmosphere-ocean interaction and global climate change. Satellite observations can accurately estimate the MLD over large scales, effectively overcoming the limitation of sparse in situ observations and reducing uncertainty caused by estimation based on in situ and reanalysis data. However, combining multisource satellite observations to accurately estimate the global MLD is still extremely challenging. This study proposed a novel Residual Convolutional Gate Recurrent Unit (ResConvGRU) neural networks, to accurately estimate global MLD along with multisource remote sensing data and Argo gridded data. With the inherent spatiotemporal nonlinearity and dependence of the ocean dynamic process, the proposed method is effective in spatiotemporal feature learning by considering temporal dependence and capturing more spatial features of the ocean observation data. The performance metrics show that the proposed ResConvGRU outperforms other well-used machine learning models, with a global determination coefficient (R²) and a global root mean squared error (RMSE) of 0.886 and 17.83 m, respectively. Overall, the new deep learning approach proposed is more robust and advantageous in data-driven spatiotemporal modeling for retrieving ocean MLD at the global scale, and significantly improves the estimation accuracy of MLD from remote sensing observations.

KEYWORDS:

• Mixed layer depth
• remote sensing observations
• residual convolutional gate recurrent unit
• global ocean

1. Introduction

One of the three fundamental layers that make up the ocean’s vertical structure is the upper ocean mixed layer. Almost vertically uniform density, temperature, and salinity characterize this quasi-homogeneous ocean layer (Kara, Rochford, and Hurlburt 2003; de Boyer Montégut et al. 2004; Sallée et al. 2021). The mixed layer depth (MLD) is a significant physical parameter to describe the dynamic upper ocean mixed layer. It is crucial for controlling the uptake and subduction of anthropogenic carbon (Karleskind, Lévy, and Memery 2011; Graven et al. 2012; Bopp et al. 2015), modulating heat transport and storage between the atmosphere and the ocean (Yu et al. 2019; Siegelman et al. 2020). Its variability thus highly determines the heat content, and affects the thermal inertia of the surface ocean, as well as the heat waves and ocean memory (Shi et al. 2022; Taves et al. 2022; Yao, Wang, and Fu 2022). In addition, the MLD shapes marine ecosystems by hosting the most primary production and providing the oxygen of the deep ocean (Diaz et al. 2021; Sugimoto 2022). It also plays a vital role in the ocean’s submesoscale phenomena, mode water formation, and the Antarctic intermediate water formation (Dong et al. 2008; Wang et al. 2022). Therefore, accurate and effective estimation of the MLD is an indispensable study that significantly improves the understanding of the atmosphere-ocean interaction and the ocean’s physical processes, which help assess and predict the changes in the climate system under global warming. However, the technical framework for accurate estimation of MLD is not fully established and remains challenging globally.

Two main methods, classified as indirect and direct, have been employed to estimate the MLD. The direct estimation is based on the 3D subsurface temperature or density using different criteria including threshold, gradient, and their combination (Holte and Talley 2009). The direct MLD estimation from 3D thermohaline structure generally relies on in situ observations (de Boyer Montégut et al. 2004; Holte et al. 2017), data interpolation (Li et al. 2017; Zhang et al. 2022), numerical simulation and reanalysis (Saha et al. 2014; Storto, Masina, and Dobricic 2014), and remote sensing retrieval (Guinehut et al. 2012; Jeong et al. 2019). However, in situ observations cannot adequately provide spatiotemporal continuous data, and data interpolation mainly relies on profile data but has large deviations without sufficient observation data before the Argo era. In addition, numerical simulation and reanalysis are restricted by the prior assumption and complicated parameters, and the subsurface thermohaline reconstruction has uncertainties from remote sensing data. Overall, the above of these methods are direct estimates with limitations caused by the sparsity of 3D ocean observations.

Previous studies have demonstrated the advantages of the deep ocean remote sensing (DORS) technique in indirectly estimating the ocean interior variables from satellite observations combined with Argo observations (Ali, Swain, and Weller 2004; Su et al. 2015; Lu et al. 2019; Su et al. 2021; Su et al. 2021), These advantages include the reduction of the complexity and subjectivity, and filling the data gaps before the Argo era. Moreover, indirect methods for estimating MLD based on DORS techniques have been developed, including dynamic theoretical models and machine learning approaches. The complexity of dynamic models makes modeling and estimation accuracy for large-scale subsurface oceans uncertain (Yan, Schubel, and Pritchard 1990; Zervakis, Kokkini, and Potiris 2017). However, machine learning-based methods have gained popularity and are widely used to estimate MLD. These methods include the bootstrap multilayer perceptron to estimate the MLD in the Gulf Stream (Pauthenet et al. 2022), the clustering neural network method to estimate the MLD in the tropical Indian Ocean (Gu et al. 2022), and the convolutional neural network and variational autoencoder to estimate the MLD in the equatorial Pacific and Southern Indian Ocean (Foster, Gagne, and Whitt 2021). However, these machine learning methods to estimate the MLD have not yet been extended to the global ocean and do not explicitly consider the inherent spatiotemporal dependence of the ocean dynamic process.

Spatiotemporal deep learning methods have been well employed to capture and analyze the complex and nonlinear spatiotemporal features of the ocean processes (Song et al. 2022; Su et al. 2022). Among them, the convolutional gate recurrent unit (ConvGRU) algorithm has exhibited a powerful ability to efficiently learn spatiotemporal features while considering training efficiency to make accurate estimates (Ballas et al. 2016). However, the ConvGRU algorithm mainly focuses on local information and has limited capability to fully exploit the spatial features of ocean data. This limitation may lead to inaccurate estimations when spatial information is essential to the phenomena being studied. Therefore, developing new approaches that can effectively exploit both local and global spatial features and time dependence is important to increase the estimation precision of the ocean interior state and change.

This study proposes a novel approach that combines ConvGRU with residual network and spatial attention blocks to construct a spatiotemporal deep learning model to accurately estimate the global MLD from satellite observations, along with a comparison of various methods. The Argo MLD gridded data are used for training and validation. The performance and robustness of the ResConvGRU approach and the contribution of each block to MLD estimation are also evaluated. By showcasing the potential and feasibility of combining satellite remote sensing and spatiotemporal deep learning in global MLD reconstruction, this study serves as a proof of concept. This study aims to provide a new technical framework based on satellite observations and artificial intelligence for the estimation of the MLD and to achieve accurate and robust global-scale MLD mapping from a new remote sensing perspective. In the future, this model is expected to contribute to the reconstruction of a long-time series of global MLD before the Argo era and enable MLD super-resolution reconstruction from satellite remote sensing.

The remaining sections of this paper are organized as follows: Section 2 offers an overview of the study area and the data used. Section 3 outlines the research methods and experimental setup employed in this study. Section 4 presents the estimated results and discussions. Finally, Section 5 concludes this study.

2. Study area and data

This study focuses on the Pacific, Atlantic, Indian, and Southern Oceans, primarily situated in the region of the global ocean that spans 180°W-180°E and 60°S-60°N. Among them, the Southern Ocean is defined according to the Coupled Model Intercomparison Project standard as the region located south of 30°S. In this study, multisource remote sensing data are employed as input data, including sea surface temperature (SST), sea surface salinity (SSS), significant wave height (SWH), sea surface density (SSD), northward and eastward components of sea surface wind (USSW and VSSW), absolute dynamics topography (ADT), zonal and meridian absolute geostrophic velocities (UGOS and VGOS), longitude (LON), and latitude (LAT). The SST sets use Daily Optimum Interpolation produced based on radiometer satellite and ship observations from NOAA, with a 0.25° horizontal resolution (Huang et al. 2021). The SSS and SSD sets were from the Copernicus Marine Environment Monitoring Service (CMEMS), with a 0.25° horizontal resolution covering 1993 to the present (CMEMS, 2022). The monthly SWH and ADT were provided from the AVISO multimission product, with a 1° horizontal resolution (AVSIO 2019; CMEMS, 2019). The USSW and VSSW monthly data were from the Cross-Calibrated Multi-Platform (CCMP) with 0.25° horizontal resolution (Wentz et al., 2016). All variables were transformed into monthly averages and interpolated into a 1° grid using the nearest neighbor algorithm to maintain the uniform spatiotemporal resolution. Besides, in this study, the Argo MLD global monthly 1° gridded data were employed as labels during both training and validation, which were provided by the International Pacific Research Center (2020). The Argo MLD was defined as the depth at which density increases from 10 m to the value equivalent to the temperature drop of 0.2°.

3. Methods

3.1. Residual convolutional gate recurrent unit

This study proposed a novel spatiotemporal deep learning approach, i.e. ResConvGRU, to estimate the global MLD. ResConvGRU consists of three ResNet blocks (residual blocks), two ConvGRU cells, and one spatial attention block, as shown in Figure 1. ConvGRU has been proposed as an expansion of convolutional long short-term memory, which is more suitable for the era of ocean big data due to its higher training efficiency (Ballas et al. 2016). However, capturing broader context information directly in ConvGRU presents a challenge because of the locality of the convolution operation. In addition, ConvGRU encounters the issue of gradient disappearance when modeling long sequences, resulting in the loss of essential features. To tackle the above problems, the integration of ResNet with dilated convolutions has been employed. ResNet has been proposed to overcome gradient disappearance (He et al. 2016), while dilated convolutions aim to expand the receptive field. As shown in Figure 1 (a), three parallel ResNet blocks, each employing distinct dilated convolutions. This design significantly amplifies the receptive field, addresses information loss due to gradient disappearance, and extracts multiscale spatial features of ocean observation data. In addition, to bolster the model’s focus on distinct spatial positions within the satellite observation data, thereby better capturing valuable information embedded in the data, the spatial attention mechanism has been adopted. As shown in Figure 1(c), the spatial attention block is attached to the back of the ConvGRU, which improves the model’s attention to spatial information and integrates local features by dynamically adjusting feature weights. The approach in this study combines these algorithms to build a spatiotemporal series model that can effectively capture the spatiotemporal information of ocean data and extract broader spatial features to estimate the global MLD.

Figure 1. Schematic of ResConvGRU formed from three residual blocks (a) with two ConvGRU cells (b) and a spatial attention block (c).

As depicted in Figure 1, each residual block includes dilated convolution layers with a kernel size of 3 × 3, a batch normalization layer, and a rectified linear activation function (ReLU). After testing the combination of different dilation rates, the optimal dilation rates of the three residual blocks finally were set to 1, 2, and 2, respectively. The first dilated convolution layer output channel is set to 64, and the remaining dilated convolution layer output channel is set to 128 in each residual block. The spatial features are extracted from the three residual blocks and then concatenated for input into the ConvGRU cell. After testing one, two, and three ConvGRU cells and combined with the same number of dropout layers, the network finally consists of two ConvGRU cells with hidden dense layers of size 128 and 64, respectively, and two dropout layers with rates of 0.3 and 0.2, respectively. In the spatial attention block, the two convolution layers’ output channels are set to 1 and 1, respectively.

A grid search strategy was used for selecting optimal hyperparameters to achieve accurate estimation results with the ResConvGRU model. The mean squared error loss function was used. The Adam optimizer with an initial learning rate of 0.1 and the batch size was set to 4. Figure 2 (a) illustrates the loss as a function of the number of epochs in the ResConvGRU model. It can be observed that the best optimal epoch is 120. In addition, the normalized root mean squared error (NRMSE) and determination coefficient (R²), averaged over five repetitions of the training and validation processes, were applied to determine the optimal timestep. It is worth noting that the NRMSE is obtained by calculating the RMSE and then divided by the value range of the MLD. The timestep, an important parameter in spatiotemporal neural networks, controls the model’s memory for processing sequence data. As depicted in Figure 2 (b), the results reveal that the model achieves optimal performance when the timestep is set to 6. Accordingly, the timestep was configured as 6.

Figure 2. Hyperparameter Settings for the model. The (a) represents the value of the loss function varies with the number of epochs in the ResConvGRU for MLD estimation. The (b) represents the relationship between the residual convolutional gate recurrent unit’s (ResConvGRU) performance and timestep. The bar and error bar indicate the uncertainty estimated through five repetitions of the training and validation processes while computing the mean and standard deviation (SD) values for the average of R² and NRMSE, respectively. Similarly, the below error bars are computed following the same methodology.

3.2. Comparative methods

Comparisons of the ConvGRU, GRU, and Random Forest (RF) models are conducted to demonstrate that the proposed ResConvGRU is a preferred strategy. The ConvGRU model consists of two stacked ConvGRU cells with hidden dense layers of size 128 and 64, two dropout layers with rates of 0.3 and 0.2, and a 2D convolution with an output channel of 1. RF has been extensively used in satellite observation information extraction (Su, Li, and Yan 2018) since it was proposed in 2001 (Breiman 2001). Utilizing the grid search method, two critical parameters within the RF, the number of decision trees and the count of the subsets of features were set to 150 and 2, respectively. GRU is an improved long short-term memory neural network (Cho et al. 2014), which simplifies its internal structure to make it execute faster with higher accuracy in some fields. After testing various combinations of GRU layers and fully connected layers, the model consisting of three GRU layers with neuron counts of 157, 100, and 100 and one fully connected layer was finally chosen as a comparison experiment. Meanwhile, there are dropouts of 0.2, 0.1, and 0.2 between each layer.

3.3. Experimental setup

Figure 3 shows the technical flowchart of this study, which is divided into four steps: data preprocessing, model training, accuracy testing, and MLD estimation. The multisource satellite surface data and Argo MLD data were preprocessed to ensure uniform spatiotemporal resolution and coverage. After that, for the training of the ResConvGRU and ConvGRU models, the input data was processed to a 5D tensor due to the characteristics of the spatiotemporal model. For the GRU and RF models, all input data were normalized and scaled using maximum and minimum normalization before training to narrow the gap between units and accelerate the training process.

Figure 3. Schematic of the experiment modeling procedure. The MLD estimation is divided into four steps: data preprocessing, model training, accuracy testing, and MLD estimation.

Following the data preprocessing step, the model was trained using data from January 2010 to December 2016. The remaining data from January 2017 to December 2019 were set aside for validation. All experiments were implemented on an NVIDIA Geforce RTX A4000 16-GB GPU using the Pytorch library as a backend. Ten experimental cases were designed (as presented in Table 1) to investigate the impact of varied variable combinations and building blocks. Cases 1∼6A used different combinations of variables as inputs to the ResConvGRU and comparative methods. Based on previous studies (Nardelli et al. 2017; Foster, Gagne, and Whitt 2021; Gu et al. 2022; Pauthenet et al. 2022), the SSS, SST, ADT, VSSW, and USSW were chosen as the baseline variable combination, and new variables (SWH, LAT, LON, UGOS, VGOS) are added based on them. Cases 1∼3B represent the ResConvGRU conditions with no residual block, one residual block, and two residual blocks, respectively. Case 4B represents the ResConvGRU model without the spatial attention block. Moreover, as presented in Table 2, the cross-validation strategy was applied to the ResConvGRU model using optimal input parameters to prove the robustness of the proposed approach. Subsequently, performance metrics including NRMSE, RMSE, R², mean absolute error (MAE), and Pearson correlation coefficient (ρ) were employed to evaluate the estimation accuracy. The performance metrics used in this study represent a global result derived from the validation set, rather than being an average of monthly values.

Table 1. Experimental design.

Case	Training methods
Case 1A	MLD = ResConvGRU/ConvGRU/RF/GRU (SSS, SST, ADT, VSSW, USSW)
Case 2A	MLD = ResConvGRU/ConvGRU/RF/GRU (*5, LAT, LON)
Case 3A	MLD = ResConvGRU/ConvGRU/RF/GRU (*5, SSD)
Case 4A	MLD = ResConvGRU/ConvGRU/RF/GRU (*5, UGOS, VGOS)
Case 5A	MLD = ResConvGRU/ConvGRU/RF/GRU (*5, SWH)
Case6A	MLD = ResConvGRU/ConvGRU/RF/GRU (*5, LAT, LON, SWH)
Case 1B	MLD = ResConvGRU without residual block (*5, SWH)
Case 2B	MLD = ResConvGRU with one residual block (*5, SWH)
Case 3B	MLD = ResConvGRU with two residual blocks (*5, SWH)
Case 4B	MLD = ResConvGRU without spatial attention (*5, SWH)

*5 variables: SSS, SST, ADT, VSSW, and USSW

Table 2. Cross-validation. Training and validation datasets for the ResConvGRU model.

Time spans	Training dataset	Validation dataset
Time span 1	2010.01–2016.12	2017.01–2019.12
Time span 2	2010.06–2017.06	2017.07–2019.12, 2010.01–2010.06
Time span 3	2011.01–2017.12	2018.01–2019.12, 2010.01–2010.12
Time span 4	2011.06–2018.06	2018.07–2019.12, 2010.01–2011.06
Time span 5	2012.01–2018.12	2019.01–2019.12, 2010.01–2011.12
Time span 6	2012.06–2019.06	2019.07–2019.12, 2010.01–2012.06
Time span 7	2013.01–2019.12	2010.01–2012.12

4. Results and discussions

4.1. Variable impact evaluation

To examine the effect of various combinations of variables on the global MLD estimation, ResConvGRU, ConvGRU, GRU, and RF were used to establish training methods with Cases 1A-6A to estimate MLD. Figure 4 shows the performance of four methods in all cases on the validation dataset and Figure 5 shows the NRMSE spatial distribution of ResConvGRU in different Cases.

Figure 4. Effect of different variable combinations on the MLD estimation by ResConvGRU, ConvGRU, GRU, and RF models.

Figure 5. Spatial distribution of the NRMSE for ResConvGRU from Cases 1A to 6A (a–f). The computation is based on the average of five-time estimations and the Argo MLD from 2017 to 2019 (a total of 36 months).

As depicted in Figure 4, the performance metrics of ResConvGRU consistently outperformed other methods across all cases, with an average R² of 0.876 and RMSE of 18.54 m. In comparison, ConvGRU, GRU, and RF achieved average R² values of 0.823, 0.765, and 0.729, and average RMSE values of 22.18, 25.56, and 27.40 m, respectively. For ResConvGRU, Case 5 including SWH showed the best performance in terms of NRMSE spatial distribution and evaluation accuracy. Notably, the SD of R² in Case 5A accounted for a mere 0.14% of the total variation, which was significantly lower than in other cases. Furthermore, Figure 5 illustrates that Case 5A had a smaller NRMSE value on a global scale, particularly in the eastern Equatorial Pacific Ocean. These findings suggest that SWH has a positive impact on MLD estimation within the ResConvGRU. Conversely, when comparing Case 1A with Cases 3A and 4A, the RMSE was relatively lower, indicating that the contribution of UGOS, VGOS, and SSD to MLD estimation was insignificant.

For the other methods (ConvGRU/GRU/RF), Case 6A including SWH, Lat, and Lon demonstrated the best performance, with R² values of 0.833, 0.805, 0.801, and RMSE values of 21.55, 23.34, and 23.55 m, respectively. Compared with Case 1A and Case 5A, the addition of SWH also improved the R² by 1.4%, 5.6%, and 11.5%, respectively. Similarly, when comparing Cases 1A with 3A and 4A, it can be concluded that the contribution of UGOS, VGOS, and SSD to MLD estimation is negligible. Overall, these results strongly suggest that SWH plays an important role in the estimation of MLD.

4.2. Model comparisons

To fully assess the performance of the proposed approach, Case 5A of the ResConvGRU model, and Case 6A of the ConvGRU/GRU/RF model were used as comparative experiments in this section.

Figure 6 displays the temporal variation of the metric results, in which the ResConvGRU model exhibited relatively high R² and ρ values and lower MAE and RMSE values. The R² and RMSE values of ResConvGRU, ConvGRU, GRU, and RF were 0.886, 0.833, 0.805, and 0.801 and 17.83 m, 21.55, 23.34, and 23.55 m, respectively. These metric results indicate that the ResConvGRU approach outperformed the other methods in global MLD estimation. Upon examining the metric results of the temporal variation across all methods, it was observed that the R², ρ, RMSE, and MAE exhibited a similar pattern of change over time. It is suggested that the performance of these methods is primarily influenced by the properties of the MLD or input data rather than the approaches being used.

Figure 6. Performance metrics of R², ρ, RMSE (unit: m), and MAE (unit: m) values for ResConvGRU-based (red dot), ConvGRU-based (green dot), GRU-based (blue dot), and RF-based (gray dot) MLD estimation. The computation was based on optimal estimates and the Argo MLD from 2017 to 2019 (a total of 36 months).

Figure 7 shows the boxplot distributions of R² and RMSE values for the major ocean basins. The ResConvGRU model showed a better performance compared to other models with average R² and RMSE values of 0.854, 0.870, 0.765, 0.495 and 9.64, 12.45, 6.23, 26.45 m in the Atlantic, Pacific, Indian, and Southern Oceans, respectively. Nevertheless, the Southern Ocean revealed a higher bias in each model. This might be a result of the higher intrinsic variability of the MLD in the Southern Ocean (Sallée, Speer, and Rintoul 2010), The degradation of model performance in the Southern Ocean indicates that additional sea surface variables can be incorporated to simulate relationships comparable in strength to those in other ocean basins. For example, the freshwater flux, which influences the MLD in the Antarctic below 60°S, can be added as a sea surface variable (Kara, Rochford, and Hurlburt 2003). Moreover, the sparse Argo profiles in the Southern Ocean introduce uncertainty in the relationship between Argo MLD grid data and true data, making it difficult to establish a true relationship between satellite observations data and Argo MLD data.

Figure 7. Distributions of R² and RMSE (unit: m) values of various regions of all models. The middle black line indicates the median R² or RMSE values, while the white star symbol indicates the average R² or RMSE values. The boxes capture 25–75% of the monthly R² or RMSE values. Outliers are defined as dots outside the box with values 1.5 × lower/upper quantile.

Figure 8 shows the spatial distribution of R² of the ResConvGRU, ConvGRU, GRU, and RF models. It should be noted that the spatial distribution of R² is roughly similar across models, that is, low R² values can be found in the equatorial Pacific and Indian Ocean, the west coast of Australia, and Southern Ocean regions, while high R² values can be found in midlatitude regions. Figure 8 also implies that the accuracy of estimation is affected not only by the methods but also by the variability of target MLD in different regions. For example, the intraseasonal oscillation of the MLD in the equatorial regions is substantial (Whitt, Nicholson, and Carranza 2019). However, the ResConvGRU model substantially improved the representation of the MLD in these regions.

Figure 8. Spatial distribution of the R² of the ResConvGRU, ConvGRU, GRU, and RF models (a–d). The computation is based on the average of five-time estimation and Argo MLD from 2017 to 2019.

4.3. Spatial distribution and pattern

In this section, the optimal ResConvGRU model was employed to estimate and evaluate the global and local MLD. Figures 9 and 10 display the spatial distribution results of the global ocean and Indian Ocean MLD estimated by ResConvGRU modeling with Argo MLD in Case 5A. The residual error is calculated from Argo gridded data minus the ResConvGRU-estimated data.

Figure 9. Spatial distribution of the ResConvGRU-estimated MLD in Case 5A, Argo MLD data, and the residual error of the ResConvGRU-estimated MLD (unit: m) at different seasons in 2018.

Figure 10. Spatial distribution of the ResConvGRU-estimated MLD in Case 5A, Argo MLD data, and the residual error of the ResConvGRU-estimated MLD (unit: m) in the Indian Ocean at different seasons in 2018.

As shown in Figure 9, the spatial distribution of the MLD shows a strong seasonal variation, that is, the mixed layer reaches minimum values in summer versus the deepest one in winter in each hemisphere. In Figure 9 (a and b), a small zone is observed in the highest MLD that surpasses 150 m over the south and north of the Kuroshio, and its bifurcation coincides with the strong surface cooling and wind stress (Oka, Talley, and Suga 2007). Notably, the MLD of the North Atlantic Ocean surpasses that of the Pacific Ocean, due to the absence of a halocline-formed barrier layer (Kara, Rochford, and Hurlburt 2000). In May (Figure 9 (d and e)), the MLD sharply decreased in both the North Atlantic and the North Pacific Oceans. On the contrary, the mixed layer gradually deepens in the Southern Ocean. In August (Figure 9 (g and h)), the MLD of the Southern Ocean is characterized by a deep and large zonal region, where it occurs in the northern Antarctic Circumpolar Current, which relates to the strong wind in the subpolar westerlies, making the surface water shear, leading to the mixed layer’s deepening. After the southern hemisphere winter, the enhancement of buoyancy and warming sharply shallows the mixed layer (Figure 9 (j and k)), reaching its shallowest in summer. Analyzing the ResConvGRU-estimated MLD and the residual error, it becomes evident that the North Pacific Ocean was affected by Kuroshio, leading to obvious estimation errors (Figure 9 c). Similarly, the Gulf Stream impacts the North Atlantic resulting in an evident estimation error (Figure 9 (c and f)). Moreover, obvious estimation errors are found in the Southern Ocean, which is influenced by the strong Antarctic Circumpolar Current (Figure 9 (i and l)). Overall, the spatial pattern of the MLD is well captured by the ResConvGRU model in the global ocean despite some biases over the western boundary current and the Antarctic Circumpolar Current regions.

Compared with other oceans, the Indian Ocean’s MLD variability is mainly influenced by the monsoons, as shown in Figure 10. During the summer monsoon (Figure 10 (d), (e), (g), and (h)), the spatial pattern of MLD stands out with high values, extending from 40°S to 30°S associated with the northward shift of strong winds in the subpolar westerlies. During winter monsoon (Figure 10 (a), (b), (j) and (k)), the Arabian Sea’s mixed layer deepens due to the cold and dry winds from the Tibetan Plateau (Smith and Madhupratap 2005; de Boyer Montegut et al. 2007). The MLD of the southeastern equatorial Indian Ocean exhibits shallower depths with convection suppression and light winds caused by MJO forcing (Figure 10 (j and k)). According to the ResConvGRU-estimated MLD and the residual error, both an overestimation and an underestimation are observed in the south of the Indian Ocean, which are affected by subpolar westerlies with obvious estimation errors (Figure 10 (i) and (l)). However, the spatial pattern of the estimated MLD is comparable to that of the Argo MLD, displaying a remarkable semiannual signal in the Indian Ocean.

The seasonal cycle of Argo MLD and ResConvGRU-estimated in the global ocean and various basins, as shown in Figure11. Notable, there is a little bit of underestimation in the Southern Ocean from February to March. These estimation errors can be attributed to the pronounced seasonal cycle of MLD in the Southern Ocean. Overall, the estimated MLD demonstrates a similar time-series change with the Argo MLD gridded data on both global and basin scales.

Figure 11. The monthly mean variation of Argo MLD (dotted line) and ResConvGRU-estimated (line) at the global ocean (GO), Atlantic Ocean (AO), Pacific Ocean (PO), Indian Ocean (IO), and Southern Ocean (SO).

4.4. Ablation study and cross-validation

In this section, the effect of residual block and spatial attention mechanisms on the performance of ResConvGRU models was investigated in estimating MLD. To do this, a set of experimental control groups, named Cases 1B–4B, was designed with varying degrees of ablation of these building blocks. The baseline model is Case 5A. Meanwhile, the robustness of the ResConvGRU model also is explored based on cross-validation.

Figure 12 shows the results of the experiments, illustrating the effect of different building blocks on the MLD estimation by the ResConvGRU models. The experiment was repeated five times, and the average R² values across Cases 1B∼4B and Case 5A were 0.793, 0.872, 0.880, 0.875, and 0.886, while the average NRMSE values were 1.7%, 1.9%, 1.30%, 1.32%, and 1.27%, respectively. The comparison of Cases 1B and 2B indicates that the residual block plays a critical role in improving the MLD estimation accuracy, with R² increasing by 7.9% from 0.793–0.872. This highlights the significance of using residuals as key building blocks for the ResConvGRU models to perform well. Furthermore, the comparison of Cases 2B, 3B, and 5A shows that the three residual block structures of ResConvGRU are optimal. Finally, the comparison of Cases 5A and 4B also indicated that the spatial attention mechanism can contribute to improving the MLD estimation, albeit by a relatively small amount.

Figure 12. Effect of different building blocks on the MLD estimation by ResConvGRU models.

The R² distribution of Cases 1∼4B and 5A were analyzed to better understand the spatial distribution of MLD estimation accuracy, as shown in Figure 13. As illustrated in Figure 13, regions within 30° N have a relatively low R², while the south and north of the Pacific and Atlantic Oceans have a relatively high R². The comparison of Cases 1B and 2B (shown in Figure 13(a) and (b)) highlights the significant improvement in MLD estimation accuracy with the help of residuals, especially in the equatorial, Indian, and Southern Oceans. Moreover, the design of three residual blocks in the ResConvGRU models further improves the accuracy of spatial MLD estimation, as seen in the comparison of Cases 2B, 3B, and 5A (shown in Figure 13(b), (c), and (e)). The same holds for the spatial attention mechanisms, as shown in the comparison of Cases 4B and 5A.

Figure 13. Spatial distribution of the R² of Cases 1B∼4B (a–d) and 5A (e). The computation is based on the average of five times estimation and Argo MLD from 2017 to 2019 (a total of 36 months).

A cross-validation strategy was used to estimate the MLD for various periods beyond the original 2017–2019 dataset to evaluate the robustness of the ResConvGRU model. Figure 14 shows the robustness of the ResConvGRU model based on the cross-validation strategy. The MLD estimation was repeated five times. The R² values were 0.886, 0.887, 0.891, 0.892, 0.893, 0.893, and 0.891, the NRMSE values were 1.27%, 1.25%, 1.26%, 1.81%, 1.82%, 1.71%, and 1.72%. These findings confirm that the ResConvGRU model performs robustly in estimating the MLD for different periods.

Figure 14. Robustness verification of the ResConvGRU model based on the cross-validation strategy.

4.5. Discussions

Previous studies for the MLD reconstruction have achieved low bias with Argo MLD, with RMSE of 3.79 m in the Indian Ocean (Gu et al. 2022). The proposed approach has a similar bias, with an average RMSE of 3.94 m in the same region, which suggests that ResConvGRU performs well on MLD reconstruction on the basin scale. However, it is crucial to emphasize that the ResConvGRU maintains a high level of reconstruction accuracy on a global scale.

As with previous studies, basic sea surface variables (SST, SSS, ADT, USSW, VSSW) can be used to estimate MLD. However, this study found that including SWH as an input improved the performance of the models, especially RF and GRU. The SWH provides information on ocean waves but has never been considered by previous studies. The surface geostrophic currents (UGOS, VGOS) were found to provide smaller contributions to the estimation, similar to previous findings (Pauthenet et al. 2022). It may be possible to include additional training variables, such as surface heat flux, wind stress, and surface freshwater flux, to provide information about the heat exchange, turbulent kinetic energy, and precipitation that dominate the mixing in tropical regions (Kara, Rochford, and Hurlburt 2000; 2003; Keerthi et al. 2016).

In this study, several machine learning methods were compared for MLD estimation. RF uses decision trees by randomly selecting samples and features, but they cannot capture spatiotemporal relationships in the data. The GRU method uses convolutional operations to obtain feature sequences and establishes a time sequence model using a memory unit, but it fails to fully capture the spatial learning of geographical features. The ConvGRU method is based on spatiotemporal characteristics of data, but it does not fully exploit the geospatial features of ocean data. In contrast, the proposed method considers the spatiotemporal dependence and full mining of the spatial features of ocean data, resulting in a more comprehensive and accurate estimation model.

In other related work, various machine-learning approaches have been employed to estimate or predict ocean subsurface variables, such as subsurface thermohaline structures, dissolved oxygen, and chlorophyll profiles (Su et al. 2021; Su et al. 2021; Tiyasha et al. 2021). However, the lack of physical interpretation and constraints is a common issue in most methods. To overcome this limitation, physics-informed deep learning has been proposed as an interesting approach (Raissi, Perdikaris, and Karniadakis 2019). Thus, applying physics-informed deep learning to estimate MLD holds promise for future study.

5. Conclusions

This study proposed a new spatiotemporal deep learning approach, i.e. ResConvGRU, to challenge the global MLD estimation from multisource remote sensing data (SST, SSS, SSW, SWH, ADT) and spatial information (LON and LAT), combined with Argo MLD gridded data. It incorporates multiple parallel residual blocks and spatial attention to enhance spatial feature extraction and to obtain spatiotemporal memory by ConvGRU cell. Through experiments involving various input variables and methods, the proposed approach demonstrates remarkable improvements in MLD estimation. The performance metrics demonstrate that the ResConvGRU model has effectively achieved a high global R² of 0.886 and a low global RMSE of 17.83 m, indicating strong estimated capability across the entire global ocean. Furthermore, it outperforms the ConvGRU, GRU, and RF models, achieving accuracy improvements of 5.3%, 8.4%, and 8.5%, respectively. Meanwhile, the ablation study suggests that the residual blocks with dilated convolution effectively extracted relevant information from ocean data to further improve MLD estimation. Cross-validation results suggest the proposed approach is robust and reasonable for estimating the MLD.

Overall, ResConvGRU has more superiority and robustness than ConvGRU, RF, and GRU in modeling for estimating MLD due to considering the spatiotemporal dependence of global ocean processes. The SWH is an important variable in estimating the MLD. Furthermore, seasonal variation in the MLD at a global scale can be well captured and depicted from remote sensing estimation. Future work will aim to combine deep learning with physical information, which will provide realistic physical constraints to further improve the estimation accuracy, and apply in higher resolution subsurface variables reconstruction and prediction.

The primary highlight of this study is the estimation of global ocean MLD using the novel ResConvGRU approach to provide a new technical framework based on a remote sensing perspective. It will help the reconstruction of time series subsurface dynamic parameters hidden in satellite remote sensing data. Meanwhile, this study helps to support the atmosphere-ocean interaction studies and advances our comprehension of the mixed layer’s contribution to the carbon cycle and the redistribution of heat within the climate system.

Acknowledgments

We thank the International Pacific Research Center (IPRC) for the Argo gridded data (http://apdrc.soest.hawaii.edu/datadoc/argo_iprc_gridded.php), the Copernicus Marine Environment Monitoring Service (CMEMS) for the SSS and SSD data (https://resources.marine.copernicus.eu/product-detail/MULTIOBS_GLO_PHY_S_SURFACE_MYNRT_015_013/INFORMATION), the AVISO altimetry for the SWH and ADT data (https://www.aviso.altimetry.fr/en/data/products/wind/wave-products/mswh/mwind.html; https://www.aviso.altimetry.fr/en/data/products/sea-surface-height-products/global/gridded-sea-level-heights-and-derived-variables.html), the National Oceanic and Atmospheric Administration (NOAA) for the OISST SST data (https://www.ncei.noaa.gov/products/optimum-interpolation-sst), and the Research Data Archive at the NCAR for the CCMP SSW data (https://rda.ucar.edu/datasets/ds745.1/), which are freely accessible for the public.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data supporting the findings of this study are available upon request from the corresponding author.

Additional information

Funding

This study was supported by the National Natural Science Foundation of China [grant no 41971384], and the Natural Science Foundation for Distinguished Young Scholars of Fujian Province of China [grant no 2021J06014].