The WW3 model (Tolman, 2009) was developed based on the third-generation wave model WAM, with its governing equations modeled by solving the action balance equation over the wave number-direction spectrum. The model uses the global digital elevation model (DEM) dataset released by the General Bathymetric Chart of the Oceans (GEBCO), with a resolution of 15″×15″, and wind field data at a height of 10 meters from the ERA5 reanalysis data by the European Centre for Medium-Range Weather Forecasts (ECMWF), with a resolution of 0.25°×0.25°, from January 1, 2017, 00:00 to December 31, 2021, 23:00. The extent of the wind field should be greater than or equal to the extent of the WW3 and SWAN numerical simulations. No additional data is input into the boundary conditions of the WW3 model. The wave spectrum grid of the WW3 model is set to 32×24, with a frequency range from 0.0373 Hz to 0.7159 Hz, divided into 32 bands, and wave direction divided into 24 directions. The calculation area of the model covers the longitude range of 110°E to 130°E and the latitude range of 10°N to 30°N, with a spatial resolution of 0.25°×0.25°. The layout of the model region is shown in Figure 1 .

The SWAN model was modified and improved by Booij et al. (1996) from Delft University of Technology based on the third-generation wave model WAM. The model discretizes the governing equations using an implicit method, taking into account wave-wave interactions and the breaking effects caused by depth changes during wave propagation, making it effective in simulating the evolution of nearshore waves. The computational range of the SWAN model is from 115.59°E to 117.71°E in longitude and from 21.78°N to 23.66°N in latitude, using an unstructured grid, as shown in Figure 1B . Bathymetric data comes from the GEBCO dataset, wind field data uses ERA5 reanalysis data, and the wave spectrum data at open boundary points is obtained from the wave spectrum output of the WW3 model. The simulation time range is from 00:00:00 on January 1, 2017, to 23:00:00 on December 31, 2021, with an output time interval of one hour.

To verify the accuracy of the WW3-SWAN numerical model, a MARK III Wave Rider instrument was deployed in the waters off the Stone Tablet Mountain Cape, at the coordinate position (22°55.7046′N, 116°31.4034′E), as shown in Figure 2 . The Wave Rider instrument has a wave height measurement range of ±20m. The measured data were processed by the instrument’s built-in software, which then statistically generated hourly wave height observation data. The observation period was from 00:00 on April 1, 2021, to 23:00 on November 18, 2021.

Figure 3 shows a comparison between the significant wave heights from the numerical model and the measured values. The significant wave height values from the WW3-SWAN numerical simulation are consistent with the observed values in terms of the overall trend. However, due to the ERA5 reanalysis data underestimating the intensity of typhoons in the Northwest Pacific (Li and Hu, 2021), the numerical simulation slightly underestimates the peak values of the significant wave heights.

Figure 4 shows the situation of some typhoons in the Western Pacific in 2021, with longitude on the horizontal axis, latitude on the vertical axis, and wind speed represented by the color scale. As shown in Figure 4A , during the spring season, Typhoon Surigae formed on April 10, 2021, with wind speeds rapidly increasing from 28 m/s to 60 m/s, and was upgraded to a super typhoon on April 17-18, 2021. The typhoon’s center was located approximately 1,280 km southeast of Manila, Philippines, in the Northwest Pacific Ocean (10.3°N, 131.9°E), with maximum winds near the center reaching 15 on the Beaufort scale (50 m/s). It transitioned into an extratropical cyclone on April 25. At 12:00 on April 18, 2021, the South China Sea was affected by the typhoon, with wind speeds around 10 m/s in the area of the wave monitoring site, leading to higher waves. Therefore, during the typhoon period, the average significant wave height measured in Figure 3A was 1.4m, slightly higher than the numerical simulation value. In May, with no typhoon influence, the average significant wave height at the wave monitoring site was 0.66m, with relatively calm sea conditions, and the numerical simulation values were closer to the measured values at this time.

As shown in Figure 4B , during the summer season, Typhoon Choi-Wan entered the South China Sea on June 3, 2021, with maximum sustained winds near the center reaching 65 km/h. At 9:00 on June 4, 2021, wind speeds at the wave monitoring site reached 6-10 m/s, with a peak significant wave height of around 1.5 m. As shown in Figure 4C , Typhoon Lupit formed in Zhanjiang, Guangdong, on August 2, 2021, and gradually approached the coasts of Fujian and Guangdong. By 15:00 on August 6, 2021, wind speeds from Typhoon Lupit along the Fujian-Guangdong coast reached around 10 m/s, causing the wave height at the monitoring site to reach a maximum of approximately 2.5 m. Therefore, during the typhoon periods shown in Figure 4B , the observed values of significant wave height were consistently higher than the values simulated by the numerical model.

As shown in Figure 4D , during the autumn season, Typhoon Kompasu formed in the Philippine Sea on October 8, 2021, and steadily moved westward after entering the South China Sea, with its center approaching the coastal areas of the South China Sea. Therefore, as seen in Figure 3C , the measured wave heights increased significantly during mid-October 2021, while the simulated wave heights were slightly lower.

To further validate the accuracy of the numerical simulation results in this study against the measured data, Mean Absolute Error (MAE), Root Mean Square Error (RMSE), and the coefficient of determination (R²) were used to quantitatively evaluate the accuracy of the numerical results. The calculation formulas are shown in Equations 1–3.

In the formulas, x_i represents the numerical simulation values, y_i represents the measured values, n is the total number of samples, and x ¯ and y ¯ are the mean values of the numerical simulation and measured values, respectively.

To evaluate the WW3-SWAN numerical simulation model, Table 1 uses MAE, RMSE, and R² for a quantitative assessment of model performance. Statistical analysis shows that the WW3-SWAN model performs well across different seasons. The MAE ranges from 0.1413 m to 0.2130 m, indicating that the average deviation between the simulated and observed values is quite small. RMSE, which is more sensitive to larger errors, is slightly higher, ranging from 0.1828 m to 0.2844 m. This is mainly due to the impact of extreme weather conditions like typhoons, which cause deviations in significant wave height at peak values. The R² values are notably high, between 0.7801 and 0.8493, indicating a strong linear relationship between the simulated and observed significant wave heights. The model performs best in the spring, with the highest R² value of 0.8493. In summer and autumn, frequent typhoons lead to reduced accuracy in the numerical simulation. Overall, the WW3-SWAN model reliably reflects the significant wave height in the study area, capturing the magnitude and temporal variation, and can serve as input data for the VMD-CNN-BiLSTM model.

VMD (Variational Mode Decomposition) is an adaptive, fully non-recursive signal processing technique that combines Wiener filtering, Hilbert transform, and the Alternating Direction Method of Multipliers (ADMM). As a non-stationary time series, significant wave height is well-suited for decomposition using VMD. The VMD decomposition process effectively transforms into an optimization process. The two main components of VMD are constructing the variational problem and solving it. Variational modes refer to the modes obtained by solving the variational problem. VMD iteratively searches for the optimal solution of the variational modes, adaptively updating the optimal center frequency and bandwidth for each Intrinsic Mode Function (IMF). VMD redefines the intrinsic mode function, as shown in Equation 4. Compared to other decomposition methods like Empirical Mode Decomposition (EMD) or Wavelet Transform, VMD was chosen for its superior ability to reduce mode mixing and provide more stable component separation under complex wave conditions.

Where k represents the mode number, A_k (t) is the amplitude of the k-th mode, φk(t) is the phase of the k-th mode, and uk(t) is the k-th mode function.

At this point, the variational problem constructed by VMD is shown in Equation 5:

Where u_k represents the corresponding mode function, and w ¯ k represents the center frequency of the corresponding mode.

By introducing Lagrange multipliers, the constrained optimization problem above is transformed into an unconstrained problem, as shown in Equation 6:

Where α represents the variance regularization parameter, and λ represents the Lagrange multiplier.

To solve this problem, the Alternating Direction Method of Multipliers (ADMM) is used. The specific solving steps are as follows:

By constructing and solving the variational problem, VMD can effectively decompose non-stationary data. However, the number of modes after VMD decomposition needs to be manually selected. Multiple tests are required to find the most appropriate number of modes.

CNN are an effective deep learning model widely used for feature extraction in image processing and spatio-temporal data. Through mechanisms like local receptive fields and weight sharing, CNNs can effectively capture local spatial features in the data. In this forecasting model, CNN is used to extract the spatial features of wind fields and wave heights, which will serve as inputs for subsequent time series modeling. CNN architecture is constructed by stacking three main types of layers: convolutional layers, pooling layers, and fully connected (FC) layers. Each convolutional layer contains a set of learnable filters, which aim to automatically extract local features from the input matrix. These filters perform convolution operations based on two important concepts: weight sharing and local connections, which help reduce computational complexity and enhance model performance. The pooling layer follows the convolutional layer, performing down-sampling. A notable feature of the pooling layer is its ability to reduce the dimensionality of feature maps, thus preventing overfitting. Typically, FC layers are used in the final layers of CNN architecture, and their role is to learn nonlinear combinations of features extracted by convolutional layers, generating the final output. Since wave height and wind field data usually exhibit significant spatio-temporal dependencies, CNN can effectively extract local features and patterns from this data through its receptive fields. Therefore, CNN is selected to extract features from wind fields and wave heights in this study.

Figure 5 illustrates the specific process of wind field and wave height data processed through a one-dimensional Convolutional Neural Network (1D-CNN). The input data, representing a sample at a certain time from the dataset, is preprocessed and fed into the convolutional layer of the CNN in sequence form. In the convolutional layer, multiple filters (also known as convolutional kernels) slide over the input sequence, extracting local temporal features through local connectivity and weight sharing. After the convolution operation, the data moves to the pooling layer for downsampling. By selecting the maximum value (max pooling) or the average value (average pooling) within a window, the dimensionality of the feature map is reduced. This not only decreases the computational complexity of the model but also effectively prevents overfitting. After processing by the pooling layer, the dimensionality of the feature map is significantly reduced, preserving key features while lowering computation costs. Finally, these processed feature maps are flattened into a one-dimensional vector, which serves as input for subsequent fully connected layers or other models (such as LSTM or BiLSTM) for the final prediction task.

Additionally, since the current wave height is not only related to the current wind field conditions but also influenced by historical wind field and wave height changes, traditional neural networks struggle to capture this long-term dependency. LSTM, with its special architecture, can effectively retain and utilize information from long-term time series, allowing it to capture complex temporal patterns in the data. Moreover, LSTM can solve the vanishing gradient problem found in conventional Recurrent Neural Networks (RNN), making it more stable and accurate in predicting long sequences. Therefore, in wave height forecasting tasks, LSTM becomes a natural choice to better model the temporal dependency and dynamic changes in the data.

A typical LSTM unit contains three types of gates: the input gate i_t , forget gate f_t , and output gate o_t , as shown in Figure 6 . In each gate, the state of the memory cell is controlled through element-wise multiplication and the Sigmoid function. The inputs to the LSTM model are the input data at the current state x_t and the output of the hidden state from the previous layer h_{t -1}.

The input data first passes through the forget gate, which determines which information should be discarded or retained. The equation for the forget gate is as follows:

Here, σ represents the Sigmoid activation function, and W_f and b_f represent the weights and biases of the forget gate, respectively. The current input x_t and the previous hidden state h_{t -1} are fed into the Sigmoid function. By transforming values between 0 and 1, the forget gate determines which information needs to be updated, where 0 represents unimportant information and 1 represents important information.

Next, the data passes through the input gate, with the calculation formula as follows:

Then, the current input x_t and the hidden state h_{t -1} are fed into the hyperbolic tangent function (tanh). At this point, the cell state is calculated and updated to the new cell state. The formula is as follows:

Here, tanh is the hyperbolic tangent activation function, and ⊙ denotes the element-wise multiplication operation, with C_t being the new cell state.

Finally, the output gate selects the next hidden state. The new cell state C_t and the new hidden state h_t are passed to the next time step. The formula for the output gate is as follows:

A unidirectional LSTM can only process information flow in one direction, whereas a bidirectional LSTM (BiLSTM) enhances the model’s ability to understand wave height and wind field temporal evolution by analyzing both forward and backward information in parallel. BiLSTM consists of two LSTM layers operating in opposite directions, as illustrated in Figure 7 . The horizontal dashed line represents the time axis flow of the time series data, while the vertical slanted lines depict the information transmission paths between network layers.

The VMD-CNN-BiLSTM model is shown in Figure 8 . VMD decomposes the wave height data into several Intrinsic Mode Functions (IMFs), breaking down the non-stationary wave height time series into relatively stationary subcomponents. The CNN network extracts local features from wind speed and IMFs, while the BiLSTM network models the wave time series data to accurately predict future wave heights. The detailed process is as follows:

Before being input into the prediction model, the wave height dataset was normalized to a range between 0 and 1, which accelerates the model’s convergence and improves prediction accuracy.

Due to the influence of the monsoon climate and typhoons in this sea area, the effective wave height sequence fluctuates greatly, requiring data processing. This paper uses Variational Mode Decomposition (VMD) to decompose the original sequences of wind fields and effective wave heights into several relatively smooth components. Taking the spring period from 2017 to 2021 as an example, with a data time interval of 1 hour, the VMD decomposition results are shown in Figure 10 .

From the decomposition, we can observe that the wind field ( Figure 10A ) and the effective wave height ( Figure 10B ) sequences are decomposed into five components (IMF1 to IMF5), transitioning from high-frequency to low-frequency components.

To compare the impact of VMD decomposition of wind fields and effective wave heights on wave height prediction, two cases were designed: Case 1 includes seven vectors, namely the wind field, IMF1 to IMF5 of the effective wave height, and the original effective wave height; Case 2 includes 12 vectors, specifically IMF1 to IMF5 of the wind field, the wind field, IMF1 to IMF5 of the effective wave height, and the original effective wave height.

According to the data in Table 2 , the prediction results of Case 1 and Case 2 show significant differences across different seasons. In all seasons, the errors in Case 1 are generally smaller than those in Case 2, indicating that decomposing only the effective wave height better captures its intrinsic features, while introducing the IMF components of the wind field increases the model’s complexity, leading to greater errors. Notably, the computation time for Case 1 is significantly shorter than for Case 2, especially in winter, where the CPU time for Case 2 is more than 10 times that of Case 1. This further suggests that introducing the IMF components of the wind field not only increases the model’s computational complexity but also significantly prolongs the computation time. Therefore, in the subsequent predictions, to simplify the computation and improve model efficiency, only the effective wave height data, which has a more significant impact on the predictions, will be decomposed.

This experiment used data from spring, summer, autumn, and winter between 2017 and 2021 as model driving data, with data from 2017 to 2020 used as the training set and data from 2021 as the test set. In the univariate model, only significant wave height is used as the input parameter for the BiLSTM, CNN-BiLSTM, and VMD-CNN-BiLSTM models.

Figure 11 compares the univariate predictions of BiLSTM, CNN-BiLSTM, and VMD-CNN-BiLSTM models with the WW3-SWAN simulation values. As shown in Figure 11 , the bidirectional LSTM (BiLSTM) is capable of considering both past and future information and performs well in predicting the overall trend, especially in periods with smaller fluctuations. However, in regions of sharp changes in wave peaks and troughs (as indicated by the black boxes in Figure 11 ), BiLSTM shows significant errors compared to the WW3-SWAN values. This may be due to BiLSTM’s tendency to over-smooth the predictions during periods of sharp fluctuations. In contrast, the CNN-BiLSTM model is more effective at capturing the short-term fluctuations of wave peaks, particularly in areas of peak changes, outperforming BiLSTM. However, CNN-BiLSTM is less effective at capturing troughs, possibly due to limitations in its ability to extract local features. By decomposing the significant wave height data using VMD, the model can effectively extract important frequency components, and combined with CNN’s ability to extract local features, it significantly improves prediction accuracy in areas of sharp changes in wave peaks and troughs. Overall, the VMD-CNN-BiLSTM model performs best in capturing changes in wave peaks and troughs.

Figure 12 and Table 3 compare the error metrics of the three models (BiLSTM, CNN-BiLSTM, and VMD-CNN-BiLSTM) in univariate significant wave height prediction, including mean absolute error (MAE), root mean square error (RMSE), and coefficient of determination (R²). The left y-axis of Figure 11 represents the specific values of MAE, RMSE, and R² for each model, while the right y-axis shows the relative values of each model compared to BiLSTM. Negative values of MAE and RMSE indicate that the model performs better than BiLSTM, while positive values indicate poorer performance; for R², larger positive values indicate better prediction accuracy. The results show that the VMD-CNN-BiLSTM model’s error is significantly lower than the other two models, especially in regions of sharp changes in wave peaks and troughs. Across all seasons, the VMD-CNN-BiLSTM model demonstrates the best prediction performance, particularly in the autumn and winter seasons, where complex wave height changes caused by typhoons and strong monsoons are present. For example, in the spring season, the MAE of the VMD-CNN-BiLSTM is 0.0159 meters, a 51.23% reduction compared to BiLSTM; across the entire year, the RMSE of the VMD-CNN-BiLSTM is 0.0256 meters, a 62.30% reduction compared to BiLSTM. Furthermore, the R² of the VMD-CNN-BiLSTM is the highest across all seasons and in annual statistics, reaching 0.9979 in the spring, a 1.04% improvement compared to BiLSTM. This indicates that the VMD-CNN-BiLSTM model has a stronger correlation between the predicted results and the actual observations, reflecting the actual wave height changes more accurately. Therefore, the MAE and RMSE of the VMD-CNN-BiLSTM model are significantly lower than those of other models across different seasons, indicating superior performance in capturing wave peaks and troughs. The higher R² value further demonstrates the model’s advantage in trend prediction, particularly in handling complex fluctuations.

Due to the influence of the monsoon climate in the region, this study conducted a multivariate prediction research to further improve prediction accuracy. In the experiment, both wind speed and significant wave height were used as parameters for the prediction model, considering the impact of the wind field. Similarly, data from spring, summer, autumn, and winter between 2017 and 2021 were used as model driving data, with data from 2017 to 2020 as the training set and data from 2021 as the test set.

Figure 13 presents the comparison curves of multivariate predictions from BiLSTM, CNN-BiLSTM, and VMD-CNN-BiLSTM models with WW3-SWAN simulation values. The figure shows that the monsoon climate significantly affects the wind field and wave height variations in the Southeast China Sea, especially in summer and autumn, where the frequent occurrence of typhoons exacerbates the complexity of wave height changes. Therefore, considering multivariate factors such as the wind field is crucial for improving the accuracy of wave height predictions. Spring is a transitional period from the winter to summer wind directions, with complex wind field changes, especially during the impact of typhoon “Shuriki,” where wind speed and wave height fluctuations significantly increase. Figure 13A shows that, compared to univariate predictions, multivariate predictions more accurately capture the overall trend of spring wave heights. Notably, the VMD-CNN-BiLSTM model, by effectively integrating instantaneous changes in the wind field, can accurately predict wave peak and trough changes, with its prediction curve highly aligning with WW3-SWAN simulation values, demonstrating high prediction accuracy. In summer, the prevailing southeast monsoon leads to a relatively stable wind field. Figure 13B indicates that, compared to univariate predictions, the VMD-CNN-BiLSTM model performs particularly well when considering wind field factors, with its prediction curve closely matching the WW3-SWAN simulation values, especially during August 17 to August 20, when VMD-CNN-BiLSTM accurately captures the characteristics of wave troughs. Autumn is a transitional period from summer to winter winds, with frequent typhoons and significant wave height changes. Figure 13C shows that the VMD-CNN-BiLSTM model better utilizes the intense changes in wind field data to accurately capture extreme wave peak values. The model performs excellently under extreme weather conditions such as typhoons, with its prediction curve closest to the WW3-SWAN simulation values. In winter, the Northeast monsoon prevails in the Southeast China Sea, with strong winds and long durations, resulting in higher overall wave height levels and frequent fluctuations. Figure 13D shows that in winter, the VMD-CNN-BiLSTM model effectively captures the overall trend and local fluctuations of wave heights, with its prediction curve highly consistent with the WW3-SWAN simulation values, demonstrating the best prediction performance.

Figure 14 and Table 4 show the error performance of multivariate significant wave height prediction models (BiLSTM, CNN-BiLSTM, and VMD-CNN-BiLSTM) in different seasons and annual statistics, including mean absolute error (MAE), root mean square error (RMSE), and coefficient of determination (R²), as well as the ratios of each model relative to BiLSTM. The influence of the Southeast China Sea monsoon climate and typhoons was considered to evaluate each model’s performance under complex meteorological conditions. The data in Figure 14 and Table 3 indicate that, within the same season, the VMD-CNN-BiLSTM model has significantly lower errors than the other two models. Particularly, after considering multivariate factors such as the wind field, the prediction performance of VMD-CNN-BiLSTM has significantly improved. Seasonal differences in prediction performance indicate that VMD-CNN-BiLSTM performs exceptionally well in autumn and winter, accurately capturing the drastic wave height changes brought by typhoons and strong monsoons. VMD-CNN-BiLSTM shows optimal performance in MAE and RMSE across all seasons, indicating that this model significantly outperforms BiLSTM and CNN-BiLSTM in multivariate prediction accuracy. For example, in spring, the MAE of VMD-CNN-BiLSTM is 0.0147 meters, a 51.48% reduction compared to BiLSTM; in annual statistics, the RMSE of VMD-CNN-BiLSTM is 0.0244 meters, a 61.81% reduction compared to BiLSTM. Additionally, the R² value of VMD-CNN-BiLSTM is the highest across all seasons and annual statistics, reaching 0.9981 in spring, a 0.79% improvement compared to BiLSTM, indicating stronger correlation and consistency in multivariate predictions. The inclusion of wind speed significantly improved the predictive performance of the VMD-CNN-BiLSTM model, particularly under complex meteorological conditions in autumn and winter, resulting in lower MAE and RMSE, as well as higher R². This indicates that the model is more effective at capturing the complex relationship between wind fields and wave heights, thereby enhancing the accuracy and stability of wave height predictions.