The SCS is a semi-enclosed maritime area in the western Pacific marginal sea that runs from northeast to southwest. This region has abundant natural resources, making it one of the three largest marginal seas in Asia. Resembling an irregular rhombus in shape, a stepped deepening trend emanates from the peripheral coasts toward the center, with the maximum depth reaching 5559 m and an average depth of approximately 1212 m (Zhu et al., 2019; Chen et al., 2021). The northern boundaries of the SCS begin at the line linking Nanao Island in Guangdong Province and E-luanbi on Taiwan Island. They extend southward to Sumatra Island and Kalimantan Island. Its western boundary is defined by the Indochina Peninsula and the Malay Peninsula, whereas its eastern boundary extends to the Philippine archipelago. The northeast SCS is connected to the East China Sea through the Taiwan Strait, its eastern part connects with the Pacific Ocean through the Bashi Strait and the Bahrain Strait, and its southern portion is adjacent to the Indian Ocean through the Malacca Strait (Jenner and Thuy, 2016; Sun, 2016). All inland areas adjacent to the SCS are economically developed regions in China. Coastal cities are significantly affected by any abnormal changes in sea level. Therefore, the SCS was selected as the study area. The location and topographic characteristics of the SCS are illustrated in Figure 1 .

From the Archiving, Validation, and Interpretation of Satellite Oceanographic (AVISO) data center, we obtained monthly SLA data gridded at 0.25°for 1993–2022, which were compiled from various satellite altimetric observations (e.g., T/P, Jason-1/2, ERS-1/2, and Envisat). To comprehensively analyze each component of SLA, monthly and 1°-gridded data products of the steric sea level (SSL) at 0–300 m depths from 1993 to 2022 were also employed. The SSL data products were calculated by integrating temperature and salinity data by Cheng et al. (2017). The datasets are publicly accessible at http://www.ocean.iap.ac.cn. Another component of SLA is the non-steric sea level (NSSL), which is calculated by subtracting the inverted barometer correction from gravity recovery and climate experiment (GRACE)-derived ocean bottom pressure data (Willis et al., 2008). For the calculation of ocean bottom pressure, the GRACE Level-2 Release 06 (RL06) Stokes’ coefficients from the Center for Space Research (CSR), the GeoForschungsZentrum Potsdam (GFZ), and the Jet Propulsion Laboratory (JPL) from 2003 to 2015 were employed. After the low-degree coefficients were corrected via satellite laser ranging data [cf. Swenson et al. (2008)] for detailed procedures), the ocean bottom pressure was destriped and Gaussian spatially-smoothed at 500 km to minimize uncertainties (Swenson and Wahr, 2006). The three datasets were then averaged to minimize noise further and resampled into 1°×1° gridded solutions for calculating the NSSL.

The sea surface temperature (SST) dataset provided by the Met Office Hadley Centre, known as HadISST, which has a 1°×1° spatial resolution from 1993 to 2022, was also utilized. These datasets underwent reconstruction via a two-stage reduced-space optimal interpolation procedure, followed by the superposition of quality-improved gridded observations onto the reconstructions to restore local detail (see Rayner et al. (2003) for details). The data can be accessed at https://www.metoffice.gov.uk/hadobs/hadisst/.

In this study, the satellite-altimetric SLA was resampled to 1°×1° to unify the spatial resolution with other datasets, with grid points within the SCS shown in Figure 1 .

EOF, or principal component analysis (PCA), is a common method for extracting dominant spatiotemporal modes in gridded meteorological and climate datasets (Zou et al., 2021). To prevent the reordering of EOF modes caused by trends and seasonal variability, the trend and seasonal signals for each equatorial Pacific variable (SST, SLA, SSL, and NSSL) were removed beforehand via least squares (LS) fitting. As a result, the PC1 time series for SST, SLA, SSL, and NSSL in the equatorial Pacific were obtained, reflecting their ENSO representative characteristics (Ashok et al., 2007). The process of EOF decomposition involves two steps (Hannachi et al., 2007):

First, the spatiotemporal gridded data are represented in matrix form as follows:

Where w m n represents the observed values of the SST, SLA, SSL, and NSSL at the corresponding position m and time n. Then, the anomaly field, W x , can be defined as:

where P is the identity matrix composed of N ones in the diagonals.

The covariance matrix D is then calculated as:

Finally, the spatial modes and time coefficients are obtained from the eigenvalues and eigenvectors of covariance matrix D.

WCO is the covariance strength of two time series in the time‒frequency domain calculated by wavelet coefficients. It provides a clear and intuitive representation of the linear or nonlinear correlation between two time series over time and frequency (Grinsted et al., 2004). In this study, we applied WCO to explore the relationship between SLA/SSL/NSSL PC1 (x) and SST PC1 (y). The wavelet coherence between two time series x and y is represented as:

where s is the scale; S is the smoothing operator; W n x ( s ) and W n y ( s ) are the continuous wavelet transforms of time series x and y, respectively. W n x y ( s ) is a cross-spectrum of the x and y time series. A more detailed description can be found in Grinsted et al. (2004).

To reveal the lagged effect and strength between the SCS SLA and equatorial Pacific oceanic variables (i.e., SST, SLA, SSL, and NSSL), the cross-correlation coefficient (Yang et al., 2024), R ( τ ) , between the SCS SLA (X) and equatorial Pacific oceanic variables (Y) is formulated as:

where δ 12 is the cross-covariance of X and Y; δ 11 and δ 22 are the variances of X and Y; and τ is the time lag between X and Y, ranging from –12 to 12 months in our study.

WALS is a relatively novel model averaging method that is designed to address the model uncertainty resulting from model selection introduced in 2010 (Magnus and De Luca, 2016). Compared with the Akaike information criterion (AIC), the Bayesian information criterion (BIC), the combined criterion (ABIC) (Fok and Liu, 2019; Liu et al., 2019) or frequency theory alternatives, employing semiorthogonal transformations of regression equations (Magnus et al., 2010; Rahman et al., 2020) substantially reduces the computational burden. The estimation process of WALS includes the following steps (detailed in Magnus et al. (2010) and Magnus and De Luca (2016)):

1. Set up a linear regression model with the target data, y ( n × 1 ), expressed as:

where X 1 is the n × p focus regressor matrix, X 2 is the n × q auxiliary regressor, ϵ ∼ N ( 0 ,   σ 2 I n ) is the n × 1 random error vector, σ 2 represents the error variance, I n is an n × n identity matrix, and β 1 ( p × 1 ) and β 2 ( q × 1 ) are unknown parameter vectors to be estimated in the following steps.

2. Calculate a projection matrix M as follows:

where M is a symmetric matrix. Thus, we can define M X 2 whose columns are orthogonal to the columns of X 1 . The columns of M X 2 can be orthogonalized through ( M X 2 ) ′ ( M X 2 ) = X ′ 2 M X 2 , followed by diagonalizing X ′ 2 M X 2 through

where P is an orthogonal q × q matrix, and Λ is a diagonal matrix whose diagonal elements are the corresponding eigenvalues. Through P and Λ , we can define a new orthogonalized auxiliary regressor X 2 * and a new regression parameter vector β 2 * as in Equation 9.

where X 2 * β 2 * = X 2 β 2 and X 2 * ′ M X 2 * = I q . Thus, Equation 6 can also be expressed as y = X 1 β 1 + X 2 * β 2 * + ϵ , in which β 1 and β 2 * are to be estimated.

3. Select the k-th model selection for the auxiliary variables such that m k is expressed as

where X 2 k * = X 2 * V k and k = 1 , ⋯ 2 q . V k = d i a g ( v k 1 , ⋯ , v k j ) is an q × q diagonal matrix with diagonal elements v k j ∈ { 0 , 1 } ,   j = 1 , ⋯ , q . Then, the least squares estimators of β 1 and β 2 * under model m k are given by Equation 11.

where Q * = ( X 1 ′ X 1 ) − 1 X 1 ′ X 2 * . Because ϵ ∼ N ( 0 ,   σ 2 I n ) , the joint distribution of β ^ 1 k and β ^ 2 k * is expressed as follows:

4. Define the WALS estimator of β 1 k as

where the weight function λ k ≥ 0 is assumed to be satisfied with the following conditions:

where s 2 is the estimator of σ 2 in the unrestricted model (all auxiliary variables are used). β ^ 2 * = X 2  ' * M y . Under this assumption, the WALS estimator of β 1 k can then be written as

5. Let η ^ = β ^ 2 * / σ , which implicitly implies that η ^ j ∼ N ( η j , 1 ) ,     j = 1 , ⋯ , q . The Laplace estimator η ¯ j is subsequently calculated via Equation 16 to determine how η j must be estimated.

where h ( η ^ j ) = e − c η ^ j ϕ ( η ^ j − c ) − e c η ^ j ϕ ( − η ^ j − c ) e − c η ^ j ϕ ( η ^ j − c ) + e c η ^ j ϕ ( − η ^ j − c ) and ϕ denote the cumulative distribution functions of the standard normal distribution. c = log 2 (Seya et al., 2014). Define N ¯ = ( η ¯ 1 , ⋯ , η ¯ q ) ′ .

6. According to the above processing, the WALS estimators of β 2 and β 1 are presented as follows:

In this study, the focus and auxiliary variables are chosen to be the lagged equatorial Pacific ENSO-related quantities and the SCS interannual SLA in three adjacent months, forming a Lag-WALS model for the parameter estimation (c.f., Section 5.1).

Following previous studies, we performed an EOF analysis of SST in the equatorial Pacific, revealing that the first EOF mode represents the conventional ENSO phenomenon (Ashok et al., 2007; Li et al., 2019). The same procedure is applied to the SLA, SSL, and NSSL ( Figures 3A–E ). The first spatial model (EOF1) of SLA and SSL exhibits strong similarities with SST, and a typical ENSO horseshoe spatial pattern in SST can be revealed in EOF1 of SLA and SSL ( Figures 3A–C ). Furthermore, the changes in the first EOF temporal mode (i.e., PC1) also confirm the consistent variability of SLA/SSL with SST, with a high correlation coefficient of 0.93 ( Figure 3E ). In contrast to SLA and SSL, NSSL PC1 displays large disparities, yielding a low correlation coefficient of 0.39 with SST in addition to the spatial pattern manifested from NSSL EOF1 ( Figures 3D, E ). This suggests weak or even no ENSO-related signals in NSSL.

To confirm the above speculation, WCO spectra between the PC1 time series of SST and SLA (SSL or NSSL) were computed ( Figures 4A–C ). Strong (weak) coherences between SST and SLA/SSL (NSSL) are exhibited in the WCO spectra with statistical significance higher than 95% at timescales between 1 and 4 years. This can be explained by the consistent phase (i.e., direction) between the SST and SLA/SSL vectors, which has common interannual variability ( Figures 4A, B ). The NSSL and SST, however, differ by approximately –90° in phase, implying that both are orthogonal to (i.e., independent of) each other ( Figure 4C ). This confirms that weak or even no ENSO-related signals are present in NSSL.

The mechanism of the influence of ENSO on the SCS SLA is believed to involve atmospheric circulation and the intrusion of the Kuroshio Current (Rong et al., 2007, 2009; Xi et al., 2019b; Xiong et al., 2023; Yang et al., 2024). Given the above coherence between the PC1 time series of the equatorial Pacific SLA/SSL/NSSL and SST, we further tested the consistency via cross-correlation analysis between the equatorial Pacific SLA/SSL/NSSL/SST and the interannual SLA in the SCS. We found that the spatial pattern of the maximum cross-correlation between the SCS SLA and SLA/SSL PC1 (NSSL PC1) was similar (different) to (from) that between SST PC1 ( Figure 5 ). This finding indicates that the time series of the equatorial Pacific SLA and SSL PC1 can affect the interannual SCS SLA, as can the SST PC1. The resulting spatial distribution of the maximum cross-correlation between the SST/SLA/SSL PC1 and the SCS SLA indicates a significant impact of ENSO on the southeastern coastal boundaries of the SCS, which is consistent with the findings of Soumya et al. (2015). This characteristic may be attributed to the geographical setting in which the eastern boundary of the SCS is adjacent to the western Pacific, where ENSO-related coastal Kelvin waves can propagate to the southeastern SCS through the Mindoro Strait, inducing negative (positive) SLA variations during El Niño (La Niña) events (Liu et al., 2011). A relatively weak correlation was found between the SST/SLA/SSL PC1 and the SCS SLA in the interior SCS, which may be related to the wind-stress curl anomalies associated with ENSO and local wind patterns that can dampen the SLA variability in the central SCS (Cheng et al., 2016). During the mature phase of El Niño in winter, an anomalous anticyclone dominates the SCS (Wu et al., 2003), leading to downwelling in the central basin where the SLA is expected to increase. The opposite situation occurs during La Niña events (Wang et al., 2018). Additionally, the topography of the SCS, characterized by an abyssal plain with a mean depth of 4300m, may also influence the spatial distribution of this correlation (Wang et al., 2018).

Figure 6 also shows the spatial patterns of the time lag corresponding to the maximum cross-correlation. Again, we find that the time lag pattern between SLA/SSL PC1 (NSSL PC1) and the SCS SLA is consistent (inconsistent) with that of SST PC1. This further demonstrates that the PC1 time series of the SLA and SSL almost synchronized with the SST response before affecting the SCS SLA. This is consistent with the wavelet coherence analysis results illustrated in Figure 4 . This suggests that the SLA and SSL could be indices equivalent to SST in revealing the time lag between ENSO-related signals and the interannual SCS SLA.

In summary, a significant correlation between the SCS SLA and SST/SLA/SSL PC1 is found, with most grid points showing the maximum cross-correlation when the SCS SLA lags the SST, SLA, and SSL by several months. It is anticipated that interannual sea level variations in the SCS can be predicted in advance on the basis of ENSO-related SST/SLA/SSL variations. Therefore, the SST, SLA, and SSL PC1 time series were used to establish the SCS interannual SLA prediction model in Section 5.1.

The interannual SCS SLA also affects itself by the previous months (Fu et al., 2019; Zhao et al., 2021). In other words, an autoregressive time lag pattern exists for the interannual SLA in the SCS. To illustrate this, the cosine (latitude) weighting scheme was employed to calculate the spatially averaged interannual SCS SLA time series followed by calculating its time autocorrelation. Figure 7 depicts the correlation matrix of the average interannual SCS SLA in the adjacent months. We found that the interannual SCS SLA yields high time autocorrelation coefficients for the adjacent three months (i.e., > 0.8), indicating a moderately strong relationship between the current interannual SLA and the previous three months. The shorter the time interval, the greater the autocorrelation. This implies the potential use of the interannual SLA from the previous month for forecasting the future interannual SLA. Therefore, the interannual SLA for the past three months was used as an auxiliary variable to constrain the forecasted interannual SCS SLA in the Lag-WALS model.

In this study, we incorporated the autoregressive property of the SCS SLA and the time lag between each equatorial Pacific ENSO-related SST/SLA/SSL PC1 and SCS SLA into a Lag-WALS model based on the WALS method to forecast the interannual SCS SLA. The specific steps of establishing the Lag-WALS model are as follows:

where β 1 , 1 ,   β 1 , 2 ,   β 1 , 3 and β 2 , 1 ,   β 2 , 2 ,   β 2 , 3 are the coefficients of the Lag-WALS model; SST_PC1 , SLA_PC1 , and SSL_PC1 are the PC1 time series of the SST, SLA, and SSL in the equatorial Pacific; S L A i n t e r is the interannual SLA in the SCS; t is the t-th epoch; Δ t is the weighted average time lag; and ϵ is a random vector of unobservable disturbances.

Then, the forecast S L A i n t e r in the epoch, t, can be obtained via the Lag-WALS model coefficient, β , solved via Equation 18. Note that the S L A i n t e r prediction values must be denormalized to restore their data magnitude.

To evaluate the performance of the Lag-WALS model, both internal and external accuracies are assessed in terms of the standard deviation (STD), root-mean-square error (RMSE), mean absolute error (MAE), and Bias (Fok et al., 2011; He et al., 2018; Ma et al., 2022). The STD, MAE, and Bias for the internal accuracy of the Lag-WALS model from 1993 to 2019 in the SCS are relatively small ( Figure 8 ). The resulting performance of the SCS boundaries is more accurate (i.e., a lower STD and MAE) than that of the central region. Their spatial patterns are largely similar to the maximum cross-correlation spatial pattern, as displayed in Figure 5 . Figure 9 shows the RMSE, MAE, and Bias for the model’s external accuracy from 2020 to 2022. Similar conclusions are drawn, except for the outliers at particular locations. This finding indicates that the impact of ENSO might be a potential reason for the higher accuracy in coastal areas, as the maximum cross-correlation between the interannual SLA and ENSO is greater in coastal areas than in the central SCS ( Figure 5 ). This finding is consistent with the previous literature discussed in Section 4.2.

Another reason for the higher accuracy in coastal areas is that SSL significantly alters the SLA over deep water (Dangendorf et al., 2021), whereas NSSL substantially impacts the coastal SLA within the SCS (Cheng and Qi, 2010). Note that SSL and NSSL are two major SLA components. In this study, SSL data products, which were calculated on the basis of the integration of temperature and salinity data by Cheng et al. (2017), were downloaded at http://www.ocean.iap.ac.cn. The temperature and salinity data used to calculate the SSL are mostly reanalysis data. Different results are present among different SSL data products owing to different datasets and processing strategies (Camargo et al., 2020), thus resulting in unreliable SSL estimates over deep water. Coastal regions are shallower than SSL-dominated regions, such that fewer vertical depth structures are present. This makes the impact of NSSL dominate the coastal SLA.

Two notable regions with apparent error features are found. One is the abnormally high positive values in Box A of Figure 9 , which is attributable to significant bottom topography changes within the region (c.f., Figure 1 ). This observation is consistent with that of Xu et al. (2012), who highlighted that significant changes in the mean dynamic topography have significant impacts on the SLA in the central SCS. Additionally, tidally induced energy dissipation amplified over the rough slope topography with complex seamounts and canyons in the northern SCS (Chang et al., 2006) significantly contributes to extreme sea level variability in the SCS (Menéndez and Woodworth, 2010; Zhang and Sheng, 2015). The spatial distribution of sea level maxima in the SCS is also influenced by tropical cyclones (Feng and Tsimplis, 2014; Zhang and Sheng, 2015).

Box B should be attributed to our interannual SLA prediction model, which merely considers ENSO-related oceanic signals. In essence, ENSO contributes 30% to the SLA, whereas the Indian Ocean Dipole (IOD) contributes 40% to the SLA in the southwestern SCS (Soumya et al., 2015; Li J. et al., 2023). Additionally, decadal processes, such as the Pacific Decadal Oscillation (PDO) and North Pacific Gyre Oscillation (NPGO), can indirectly modulate the interannual SLA in the SCS via subsurface temperature and salinity (Zhou et al., 2012; Deng et al., 2013; Wang et al., 2021; Li J. et al., 2023). In the southeastern SCS and the northern Gulf of Thailand, the PDO accounts for ~30% of the interannual SLA (Soumya et al., 2015). Specifically, the SCS, located between the western Pacific and the Indian Oceans, is significantly affected by interannual and decadal sea level fluctuations linked to ENSO and the PDO in the Pacific Ocean, as well as the IOD in the Indian Ocean (Mohan and Vethamony, 2018). The PDO is the predominant source of interdecadal climate variability in the Northwest Pacific, defined via monthly SST anomalies in the Pacific poleward of 20°N (Mantua et al., 1997). The decadal variability in the SLA is governed primarily by the PDO, particularly in the northwestern region of the SCS (Cheng et al., 2016; Mohan and Vethamony, 2018). In contrast, the IOD and ENSO are more influential in explaining the interannual sea level variability in the eastern Indian Ocean and SCS. While the IOD predominantly drives interannual SLA variations in the eastern Indian Ocean, the ENSO is the primary influence on the SCS (Mohan and Vethamony, 2018). Nonetheless, the interannual SLA variability in certain areas of the SCS is also affected by the IOD and PDO. Therefore, the interannual SLA prediction accuracy in the SCS should be substantially improved when the combined effects of ENSO, the IOD, and the PDO are incorporated into Equation 18, albeit with potential overfitting. This represents the major limitation of this study. Another possibility is that the correction of glacial isostatic adjustment to satellite altimetry has been ignored potentially affecting the long-term trend, although its impact is relatively small (Huang et al., 2013), at approximately –0.3 mm/yr over the SCS (Feng et al., 2012).

Overall, the prediction accuracy (RMSE) of the Lag-WALS model in the SCS ranges from 0.49 to 5.93 cm, with average STD/MAE/Bias (RMSE/MAE/Bias) values of the model’s internal (external) accuracy in the SCS of 1.01/0.80/–0.002 cm (1.39/1.13/–0.08 cm), respectively ( Table 1 ). Table 2 displays the prediction accuracy of all existing SLA prediction methods. The results obtained by the Lag-WALS method are comparable to those of the other methods. In particular, our resulting accuracy closely agrees with the best result of the hybrid prediction model proposed by Fu et al. (2019) over the SCS (i.e., RMSE: 1.32 cm; MAE: 1.03 cm). This finding indicates that our proposed model yields good reliability for interannual SLA estimates, whereas the calculation steps are simpler and computationally more efficient than those of Fu et al. (2019). Note that Fu et al. (2019) is a millimeter more accurate than our prediction, which is attributed to the greater number of involved steps and time-consuming data processing procedures.

To assess the performance of our proposed model seasonally, the forecasted interannual SLAs were compared with the altimetric observations from 2020 to 2022 across four seasons: spring (March to May), summer (June to August), autumn (September to November), and winter (December to February). The spatial patterns of the altimetric-observed mean interannual SCS SLA and the model-forecasted SLA for different seasons are shown in Figures 10A, B . Both patterns basically agree for each individual season. Notably, the model’s estimates in the central region of the SCS during summer exhibit a more diverse spatial pattern than those in the other three seasons. The residuals, which are the difference between the observed and model values ( Figure 10C ), are relatively small regardless of season. The regions with significant outliers are largely consistent with the findings shown in Figure 9 , as the potential reasons discussed above. Overall, the statistics of the residuals for the four seasons are –0.14, 0.17, –0.06, and –0.28 cm, respectively ( Table 3 ).

In addition to the comparison with the observed values in different seasons, changes in regionally averaged time series were also analyzed. Figure 11 shows the changes in the regionally averaged time series of the observed and model values, along with their differences from 1993 to 2022. The comparison reveals a strong consistency between the observed and modeled time series changes, with differences mostly within ±1 cm. The RMSE and MAE of the training (testing) set for the difference between the observed and model values are 0.36 cm and 0.28 cm (0.49 cm and 0.39 cm), indicating that the relevant parameters obtained from the training set have good performance in the testing set (i.e., 2020–2022) and that the model did not exhibit overfitting issues in this study. Furthermore, we set the significance level at 0.05 during model training stage. This study focuses only on the interannual variation in the SLA, eliminating the influence of other signals, followed by normalizing the data during the model construction stage to avoid the impact of differences in data magnitude. These strategies effectively prevented the occurrence of overfitting issues in this study. The Lag-WALS model generally yields good reliability and higher accuracy in predicting the SCS SLA at the interannual scale.