As the labeled data, three measurements are used in this study: the mooring system, high-resolution profiler, and shipboard survey cruises. A time series of temperature profiles over 9 months (from 22 July 2019 to 15 May 2020), recorded by a mooring system (named M1) which deployed in the SYS, near the western boundary of SYSCWM (35.18°N,122.26°E, Figures 1A, B ). The M1 data has 244 temperature profiles after quality control, including 17 depth levels (from 1 m to 55 m), covering the maturation to disappearance of the SYSCWM. The moored high-resolution profiler (named H1), which was deployed at the same location as M1 from 3 June 2022 to 4 July 2022, provides a fine temperature profiles time series. This profiler recorded vertical temperature profiles from 1 m to 50 m during the growth to maturity of the SYSCWM. The sample interval of H1 is 30 min and the vertical resolution is 0.1 m. In this study, the spatiotemporal resolution of the H1 data is averaged to daily and 1 m. In addition, the 55 m depth level of H1 data is extrapolated from several adjacent temperatures for their similarity. Cruise observations were carried out with 1 m vertical resolution in the western SYS in April, July and October 2019. The cruise covered the sea west of 124°E, from 33°N to 37°N, and a total of 5 latitude sections were used in this study. The five temperature latitude sections obtained by CTD castings during the cruise survey along different latitudes (33°N, 34°N, 35°N, 36°N, 37°N), named S33-S37 ( Figure 1B ).

Multisource satellite observations are used as input data, including absolute dynamical topography (ADT), SST, SSS, and SSW. The SSW contains u and v components (USSW, VSSW). The ADT data are provided by SSALTO/Data Unification and Altimeter Combination System (DUACS) and were available through the Copernicus Marine Environment Monitoring Service (CMEMS, https://marine.copernicus.eu/). The product merged multiple L3 along-track measurements and conducted the tidal corrections (Taburet et al., 2019). The SST data are obtained from Daily Optimum Interpolation Sea Surface Temperature (DOISST, https://psl.noaa.gov/), developed by National Oceanic and Atmospheric Administration Physical Sciences Laboratory (NOAA PSL). It is a blend of in situ SST with satellite SST derived from the Advanced Very High Resolution Radiometer (Banzon et al., 2016; Huang et al., 2021). The SSS data are obtained from SMOS L3OS 2Q Debiased daily valid ocean salinity values product (https://sextant.ifremer.fr/), which are distributed by Centre Aval de Traitement des Données SMOS (CATDS) and corrected the offshore SSS through various in situ observations (Boutin et al., 2018). The SSW data are provided by the Cross-Calibrated Multi-Platform (CCMP; https://rda.ucar.edu/datasets/ds745.1/). The CCMP uses a variational analysis method to smoothly fuse multisource surface wind data into the gridded data at 6 hours intervals (Atlas et al., 2011). The temporal resolution of the CCMP data is 6 hourly while the rest is daily, and the spatial resolution of all these data is 0.25°×0.25°.

We coupled the tide model data into the inputs of machine learning model. The tide model data, including surface tidal elevation and tidal currents, are estimated by the TPXO7 global tidal model provided by Oregon State University, which was built hourly on a 0.25°×0.25° grid. The tide model is based on the hydrodynamic equation and uses the generalized inversion method to assimilate the measured data, including satellite altimetry data and tide observations. Furthermore, it was recently used for the hydrographic study in the YS (Bi et al., 2021; Lin et al., 2021; Sun et al., 2022). The M₂ tide is the most dominant tidal component in the SYS, having stronger tidal current ( Figure 1A ). The tides have complex structures in the SYS, which is detrimental to temperature inversion. In this study, the tidal time series of eight basic tidal components (M₂, S₂, N₂, K₂, K₁, O₁, P₁, and M₄) are extracted by the Matlab Tide Model Driver toolbox (https://www.esr.org/research/polar-tide-models/tmd-software/). The tide model data and satellite observations, which have the same spatial resolution, were co-located with the temperature profiles by the nearest neighbour method, and the temporal resolution is unified to daily.

We also validate the temperature estimation with the ARMOR3D dataset (Guinehut et al., 2012), which was obtained through CMEMS. The ARMOR3D used multiple linear regression and optimal interpolation, providing the weekly temperature and salt fields at 0.25° × 0.25° resolution over 15 regularly spaced vertical levels between surface and 80 m depth. The weekly averaged three-dimensional temperature field in April, July and October 2019 from ARMOR3D is used to compare. The YSCWM below the thermocline is clearly visible in the observation-based ARMOR3D data ( Figure 2A ). In addition, the M1 temperature data are used to evaluate ARMOR3D. In order to match the temporal resolution, the M1 data are first calculated as weekly average and then compared to the nearest neighboring grid in ARMOR3D. As shown in Figure 2B , most of the data points are distributed along the equal line with low bias, absolute error and high Pearson’s correlation coefficient. The evident seasonal temperature variations in ARMOR3D are well simulated compared to the M1 observations ( Figures 2C, D ). Even though ARMOR3D presents a shallower mixed layer and a more durable YSCWM which lasts until October, it well reproduces the vertical thermal structure at the M1 station and is worth to refer for the thermal structure of SYS.

Considering the large seasonal variation of the thermal structure in the SYS, we use unsupervised GMM clustering techniques to shrink the sample space and improve the accuracy (Landschutzer et al., 2013; Parard et al., 2015). As a probabilistic model, GMM is often used for data clustering (Attal et al., 2015). First, the GMM randomly initializes the Gaussian distribution parameters of each cluster. Then the posterior probability of each sample is calculated and used to compute the new Gaussian distribution parameters. The process is repeated until the expectation function is maximized. Compared with the K-means method, GMM is more suitable for non-spherical clusters with different sizes and densities (Wang et al., 2019; Askari, 2021). Therefore, it is appropriate for the classification of ocean temperature profiles (Maze et al., 2017; Sambe and Suga, 2022). GMM requires the number of classes (K) as an input parameter. Therefore, the Davies-Bouldin index (DBI) is used to determine the appropriate number of classes in this study. The number of classes having the minimized DBI is considered the optimal result. Since the initial values of the Expectation-Maximization algorithm are randomized, the GMM clustering was applied 20 times, and 80% of the data were randomly selected from the M1 and H1 data each time to stabilize the clustering results. Figure 3 shows the DBI from clustering results with different K. As a result, we judge that stable and good clustering results could be obtained if K = 3. The clustering results are shown in Figure 4 . Although the YSCWM temperature structure from H1 data is still growing, it is approaching maturity. Therefore, they are named after a specific stage of YSCWM: the maturity stage, the declining stage, and the disappearance stage. During the maturity stage of YSCWM with weaker wind, the sea surface is subjected to strong thermal radiation, forming a stable upper mixed layer and a strong thermocline, which prevents heat transfer, so the bottom water stays cold (Lee et al., 2016). It leads to a multi-layer temperature structure in the SYS, with a large temperature difference between the sea surface and the bottom ( Figure 4A ). In the YSCWM declining stage, the cooling at the sea surface and stronger mixing lead to a thicker and colder upper mixed layer and the subsequent weakening and deepening of the thermocline ( Figure 4B ). Meanwhile, critical tidal currents raise the temperature at the bottom layer then decline the YSCWM (Li et al., 2016). Thermal forcing at the air-ocean interface and agitation by strong winds together cause strong vertical mixing, forming a well-mixed low temperature structure ( Figure 4C ) from the sea surface to the bottom in the YSCWM disappearance stage (Chu et al., 1997).

To tackle the limitations of small data and complex computations, we adopt the LGBM algorithm to predict the temperature by taking advantage of its lightweight. LGBM is a gradient boosting framework based on decision trees, which has been well used in the marine field and shown a faster training speed and higher accuracy for small data (Su et al., 2021; Dong et al., 2022). Same as the other boosting algorithms, it sums the results of multiple decision trees as the final prediction output. Gradient-based One-Side Sampling (GOSS) and Exclusive Feature Bundling (EFB) are two important features of LGBM. The GOSS excludes most of the samples with small gradients and calculates the precise information gain by the remaining samples. The EFB approach integrates many mutually exclusive features and reduces the data dimension. To build a better model, the Bayesian optimization strategy is used to optimize several important parameters of LGBM. The optimization method is a Gaussian process with a faster speed. According to previous studies, three essential hyperparameters need to be adjusted: the number of leaf nodes (num_leaves), the learning rate, and the number of iterations (n_estimators). The bounds of n_estimators were set 100 and 1000, and the best n_estimators is 400 without overfitting. It improves the accuracy by 16.6% compared to n_estimators=100. However, the accuracy at n_estimators=1000 is only increased by 0.1% compared to the best n_estimators. When the learning_rate is increased to 0.01 from 0.001, the performance is improved by 21% compared to the starting learning_rate=0.001, but the effect does not enhance when it is increased further until 0.1. The test range of num_leveas is from 5 to 30. The performance of the model at the best num_leveas=5 improved by 3.4% over num_leveas=30. The optimal parameters are shown in Table 1 . In addition, the max depth is set to 5, to prevent overfitting due to excessive complexity of the model. The other parameters are set to default values.

First, we input the eight harmonic components (M₂, S₂, N₂, K₂, K₁, O₁, P₁, and Q₁), geographic location and time parameters into the TPXO7.2 global tidal model, to extract tidal elevation and tidal currents data. TPXO7.2 fits best the Laplace tidal equation in the least squares sense. Second, the datasets consisting of tide model data, satellite observations and in situ temperature profiles are divided into three different stages by GMM clustering. The surface parameters (ADT, SST, SSS, SSW, tidal elevation, and tidal currents) are used as independent input variables and the temperature time series are used as labels to prepare the training and test data. To ensure that the training and test sets have a similar seasonal distribution, all samples at the location of M1 are normalized and randomly sampled into the training set (60%) and the test set (40%) by month ( Figure 5 ). Third, the model is tuned and trained using the Bayesian optimization method to obtain suitable temperature estimators at 17 depth levels. Figure 6 shows the technique flowchart of one stage at a certain depth. We use a total of 162 samples to train and 114 samples to test when using mooring observations for validation. Finally, temperature predictions are applied to a larger horizontal space and verified with cruise observations in the SYS where the number of training data and evaluating data are 276 and 78, respectively.

To evaluate the tide model coupled temperature inversion method, we designed comparative trials named GLGBM-tides and GLGBM. They both use the LGBM method with pre-clustering process but the former couples the tide model while the latter does not. Additionally, we compared other reconstruction methods. Case SVR and Case ANN use Support Vector Regression (SVR) model and Artificial Neural Network (ANN) model, respectively. Table 1 summarizes the different trials. These are optimized by the Bayesian optimization strategy, and the parameters of different models are shown in Table 1 . The ARMOR3D dataset is also used for comparison.

Previous studies have shown that strong tidal mixing has an important effect on the temperature structure and enhances vertical heat exchange in the water column during summer in the YS (Lü et al., 2010; Yao et al., 2012; Li et al., 2016; Yu et al., 2016). Here, we first compared GLGBM-tides and GLGBM to investigate how tides affect temperature estimation in this study. Figure 7 shows the comparison between the temperature profiles obtained by the two models and in situ observations. The profiles are randomly selected according to spring tide and neap tide in the maturity stage of YSCWM. In this stage, bottom vertical disturbances are stronger (Li et al., 2016), which affects the heat transfer and thermal structure significantly. Besides, the air-sea heat flux and the cooling process of the previous winter strongly influences the intensity of YSCWM (Zhu et al., 2018). This leads to machine learning models having more difficulty accessing these temperature variations and more considerable differences between in situ and estimated temperature ( Figure 7 ). However, it can be seen that the temperature profiles obtained from GLGBM-tides are more consistent with the measured profiles, especially deeper than 30 m. This confirms that the method coupled with tide model can effectively improve the structure of the predicted temperature profiles during the maturity stage. To further validate the above results, several evaluation indicators metrics are used to assess the two models. Except for root mean square error (RMSE), coefficient of determination (R²) and absolute difference, the error (defined as the proportion of RMSE in the actual mean temperature observations) is also used to evaluate the accuracy and reliability of the model. The evaluation indicators are computed as follows:

Here, T _i denotes the observed temperature while T i ' is the estimated temperature by models. The T ¯ is the mean values of T _i over the whole observation. N is the number of test samples.

From Figures 8A-C , the evaluation indicators of the two methods are similar within the 1-28 m depth layer. However, in the 40 m depth level, the RMSEs of the two are 0.806 and 0.863, respectively. Meanwhile, the accuracy of other layers has been improved by different degrees from 30 m to 55 m. Figure 8D shows that the smaller absolute errors occupy a larger proportion in the GLGBM-tides model. In addition, the enhancement is mainly manifested during the maturity stage ( Figure 9 ). It may be attributed to the tidal mixing primarily influencing the range up to 30 m from the bottom during summer (Qiao et al., 2004b). In this trial, GLGBM-tides coupled the tide model while GLGBM not. Meanwhile, strong tides affect the heat transfer and thermal structure of the profile, especially the bottom layer. As a result, GLGBM-tides better learn the temperature variation affected by tidal mixing, and it presents a more consistent vertical thermal structure with in situ observations ( Figure 7 ) and better performance than GLGBM ( Figure 8 ).

Furthermore, we analyzed the accuracy of the models at three specific stages from 30 m to 50 m ( Figure 9 ). The averaged R² and RMSE are significantly different in the maturity stage of YSCWM and similar in the decline and disappearance stages. It performs less well in the maturity stage than the other two stages in the YSCWM deep. Strong stratification leads to a large difference in temperature between YSCWM and the upper layer. Besides, YSCWM is influenced not only by the air-sea heat flux but also by the cooling process of the previous winter (Zhu et al., 2018). It means that the thermal structure of YSCWM is more difficult to be described by sea surface parameters in the machine learning models hence lower R² and higher RMSE. The averaged R² of GLGBM-tides and GLGBM are 0.614/0.547, with approximately 12% improvement. It results from the stronger influence of tidal mixing on the temperature structure in summer. Therefore, tides are worth considering in the offshore temperature field reconstruction.

Overall, the GLGBM-tides has good accuracy with errors of less than 8% at all depth layers and most absolute difference of less than 2°C ( Figures 8C, D ). It is worth noting that a bump appears above 30 m in Figure 8B . This phenomenon may be related to the depth of the mixed layer. According to previous research, the depth of the mixed layer in SYSCWM is about 5-25 m (Qiao et al., 2004b). The temperature does not vary significantly within the mixed layer, which causes the lower RMSE and higher R². The tidal mixing primarily influences the range up to 30 m from the bottom, enhancing the vertical temperature variability (Qiao et al., 2004a) and the particular structure of the YSCWM makes it difficult for the model to accurately describe the temperature variations. Therefore, the accuracy of reconstruction at these depths will be worse ( Figures 8A-C ).

It helps to understand the different effects of each sea surface parameter on the ST, by analyzing the importance of sea surface parameters at different depths. The LGBM reflects the importance of different features by calculating the number of times the sea surface parameters are used to segment the data across all trees. The relative importance of each parameter is calculated by summing and normalizing the feature importance from the LGBM. Figure 10A shows the relative importance of each sea surface parameter from GLGBM-tides. According to previous studies, the vertical thermal structure in the Yellow Sea (YS) is influenced by air-sea heat flux, the wind, tidal vertical mixing, and freshwater input (Chu et al., 1997). The temperature in the mixed layer is vertically quasi-uniform due to the mixing of multiple dynamic processes, such as wave motion and wind. Meanwhile, the mixed layer gradually thickens from the maturity stage to the disappearance stage of YSCWM, which means that the sea surface temperature (SST) can explain more subsurface temperature variations. Consequently, SST is the main driver of the model, with a more than 30% contribution at 17 depth levels ( Figure 10A ). However, below the mixed layer, the heat transfer is blocked, and it is difficult to explain the temperature change by relying on SST alone. Therefore, the trend of SST contribution decreases with deepening ( Figure 10A ).

Warming or cooling mainly drives density changes, causing sea level changes since salinity variation is not significant in the SYS. There is a close correlation between ADT and subsurface thermal structure. The sea level variations are influenced more significantly by those depths where temperature sharply changes, such as the thermocline. Therefore, the ADT contribution is higher at those depths where the temperature fluctuates drastically ( Figure 10A ), such as the thermocline in the maturity and declining stages and the bottom layer affected by tides. It leads to an average relative importance of 10% and 16% for ADT above and below 15 m depth, respectively.

SSS and SSW are also important parameters (Wu et al., 2012; Klemas and Yan, 2014; Su et al., 2015). The SSS is related to freshwater input (Nieves et al., 2014), which causes density anomalies and then affects the dynamics. The contribution of SSS is less variable from surface to 40 m depth but increases at the bottom ( Figure 10A ). This may be related to the Yellow Sea Warm Current (YSWC) in the winter, which brings a more salty and warmer water mass, especially at the bottom and manifests in the SSS. Wind forcing changes sea level and also affects ocean mixing, intensifying heat exchange between layers. Southerly winds prevail in summer and northerly winds during winter in the SYS, which causes VSSW to contribute more than USSW ( Figure 10A ). The vertical distribution of the wind (USSW and VSSW) contribution is roughly same but increases slightly at the bottom ( Figure 10A ), which is due to the mixed layer deepening during the declining stage of YSCWM.

Tide-induced mixing causes changes in the ocean heat vertical distribution. Even though the overall tidal contribution is weak and less variable, it may be important for a particular stage. During the maturity stage of YSCWM, the tides contribution (u, v, and z in Figure 10A ) is about 15% within the mixed layer but can exceed 30% below the mixed layer ( Figure 10B ) causing the tidal-induced mixing mainly affects the bottom and above 30 m range (Qiao et al., 2004b). It is comparable to the SST contribution ( Figure 10B ).

We compared other temperature prediction methods. The SVR and ANN methods have no pre-clustering process and tides. The overall R² of SVR and ANN are 0.862/0.888 with the RMSE of 1.506/1.22°C, respectively on the time series. It shows that the GLGBM coupled tides have better accuracy from Figures 11A-C . However, the ANN has similar accuracy above 20 m compared to GLGBM-tides, which may be related to the dominance of SST in this depth range. Additionally, GLGBM-tides allows errors to be smaller and more concentrated, effectively improving model performance, as revealed by the error density distribution ( Figure 11D ).

We choose H1 data to demonstrate the performance of different methods for fine and continuous data. Since deep learning is more applicable to large data, ANN performs unstable. We implement ANN 20 times to obtain the average temperature estimation. Figure 12 shows the observation from H1 and the reconstructed temperature structure from different methods. The seasonal warming in the upper mixed layer has been reproduced by all methods. Here we adopt the upper boundary of the thermocline as the mixed layer depth (MLD) to further evaluate the performance of models. The reconstructed temperature fields are interpolated to 1 m vertical resolution before calculating MLD. The results show that the MLD is maintained around 10-15 m in the H1 observations ( Figure 12A ). For the reconstructed temperature field by GLGBM-tides ( Figure 12B ), the MLD changed generally consistent with the H1 observation. Influenced by atmospheric processes, the MLD becomes shallower from 17 June to 1 July. This process is well reproduced by GLGBM-tides. Reconstructions from other methods failed to capture this variation. The MLD from reconstructed temperature by ANN is stabilized at about 15 m ( Figure 12C ) while the MLD reconstructed by SVR ( Figure 12D ) is too deep. The reconstructed temperature from ANN can indicate the trend of YSCWM but has large noises ( Figure 12C ). The temperature field estimated from SVR fails to reproduce the strong thermocline and YSCWM ( Figure 12D ). GLGBM-tides can reproduce the vertical temperature structure well compared to the observations. However, the overall estimate of the YSCWM by GLGBM-tides is slightly warmer than the observations from surface to bottom. Hence, the intensity of YSCWM from estimation is weaker. It is noticeable that the reconstruction of the thermocline is well, which assists in predicting the depth of the YSCWM.

We attempt to apply the temperature estimation at the locations of the cruise observations by training the samples from H1 and M1 and use S33-S37 data for verification. The ARMOR3D reanalysis data is used to compare as well. The temperature estimation beyond the topography is deleted. Figure 13 shows the temperature structure of 35°N and 36°N sections (S35 and S36) in three stages of YSCWM. The overall RMSE by all samples of GLGBM-tides and ARMOR3D is 1.781/2.133°C, respectively. It is higher than above due to the spatial heterogeneity of the thermal structure in SYS but the reconstructed vertical temperature structure is still in general agreement with the observations. In the mixed layer, the reconstructed temperature was colder than observation while the ARMOR3D is warmer and the reconstruction has a small zonal variation. It is the result of the training data containing inadequate spatial features. In contrast, the reanalysis data shows a clear spatial difference for fully considering spatial features during production but shows a shallower mixed layer, such as Figures 13B, D . In the declining and disappearance stages, the temperature reconstruction is better for the strong mixing but the ARMOR3D still shows a significant temperature gradient from surface to 35 m depth ( Figures 13C, D ). The estimates provide a better reconstruction of the thermocline than ARMOR3D ( Figures 13A, B ). The intensity of the thermocline in the ARMOR3D data is strong ( Figures 13A, B ) in maturity stage while it is weak in declining stage ( Figures 13C, D ). The estimated temperature of YSCWM by the GLGBM-tides is slightly warmer especially in the declining stage ( Figures 13C, D ), but consistent in terms of depth and spatial distribution. The ARMOR3D have the shallower upper boundary of the YSCWM so the temperature of YSCWM is cold as observations (see Figure 13B ). Both have good reconstruction of well-mixed temperature structure in the disappearance stage ( Figures 13E, F ). However, the cold cores in S36 could be observed ( Figures 13B, D, F ) by cruise data but this special structure is difficult to reproduce. Figure 14 shows the spatial distribution of RMSE in three stages. The accuracy of the proposed method is good from 121°E to 123.5°E. From Figure 14A , the RMSE increases from the center (location of M1) along longitude towards the sides, but with larger differences in farther regions, which may stem from the sparseness of the offshore observations. On the contrary, the RMSE of ARMOR3D decreases gradually from the center to the outside but is similar on the west side of the study area. However, the GLGBM-tides and ARMOR3D have close overall R², which are 0.927 and 0.884, respectively. Generally, our reconstruction results are reliable through comparison with ARMOR3D data.