Analysis of Wave-Induced Stokes Transport Effects on Sea Surface Temperature Simulations in the Western Paciﬁc Ocean

: This study investigated the performance of two ocean wave models, that is, Simulation Wave Nearshore (SWAN) and WAVEWATCH-III (WW3), and the interannual and seasonal variability of transport induced by Stokes drift during the period from 1989 to 2019. Three types of sea surface wind products were used for wave simulation: the European Centre for Medium-Range Weather Forecasts (ECMWF) Interim, the Cross Calibrated Multi-Platform Version 2.0 (CCMP V2.0) from Remote Sensing Systems (RSS), and the National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS). The modeling was validated against wave measurements from the Jason-2 altimeter in 2015. The analysis found that the root mean square error (RMSE) of signiﬁcant wave height (SWH) from the WW3 model using CCMP wind data was 0.17 m, which is less than the ~0.6-m RMSE of SWH from the SWAN model using the other types of wind data. The simulations from the WW3 model using CCMP wind data indicated that the Stokes transport is up to 2 m 2 /s higher in the South China Sea and Japan Sea than that at other ocean regions in January. The interannual variation showed that the Stokes transport generally increased from 0.25 m 2 /s in 1989 to 0.35 m 2 /s in 2018. We also found that the accuracy of the sea surface temperature (SST) simulation using the Stony Brook Parallel Ocean Model (sbPOM) is improved by as much as 0.5 ◦ C when Stokes transport is considered to validate the sbPOM-simulated SST against the measurements from Argo in 2012–2015. In particular, the Stokes transport has a negative effect on Summer (March to June) and has a positive effect in Autumn (July to September), which is probably caused by the tropical cyclones.


Introduction
It is well known that ocean waves play an important role in the ocean dynamic environment, especially in energy exchange in the air-sea interaction layer. Therefore, waves are an important aspect of physical oceanography. Ocean waves are also of great significance to offshore human activities and for monitoring of the global marine environment. The generation and evolution mechanism of ocean waves has been well studied over the past decades [1][2][3][4][5]. Ocean waves are small-scale wind-generated gravity waves classified as wind-sea and swell. The long-term variability of ocean waves has been analyzed in global [6] and regional seas, that is, the Northeast Pacific Ocean [7], the Northeast Atlantic Ocean [8], the tropical Indian Ocean [9], and the Irish Sea coasts [10].
Techniques for monitoring ocean waves include moored buoys, numeric model predictions, and remote sensing. Wave buoy observation is currently recognized to be an accurate method. The U.S. National Oceanic and Atmospheric Administration (NOAA) National Data Buoy Center (NDBC) ocean buoys are a valuable open-access resource for ocean wave research. Remotely sensed ocean wave parameters derived from a satellite altimeter [11], a synthetic aperture radar [12], and the Chinese-French Oceanic SATellite (CFOSAT) are also commonly used. However, these methods have obvious limitations, such as small spatial coverage and lack of data for coastal waters. Numerical simulation is a key part of the study of ocean waves [13], especially in the development of numeric wave models. Currently, the so-called third-generation model is the most advanced kind of numeric model for wave simulation [14]. It is based on the principle of energy density balance in terms of wave propagation, that is, the ocean wave model (WAM) [15], with its successor WAVEWATCH-III (WW3) (Version 5.16) [16], developed by the National Centers for Environmental Prediction (NCEP) [17,18], and Simulating Waves Nearshore (SWAN), developed by the Delft University of Technology [19,20]. The WAM model was the first operational model. The WW3 model has good performance at modelling waves [21] considering the wave-current interactions in typhoons, the depth-induced wave breaking [22], the sea water elevation [23], and sea ice [24]. These models use a rectangular grid for global wave analysis, whereas the SWAN model uses an unstructured triangular grid that performs well in near-shore regions [25]. In the previous study [26], it was found that switch ST2 in the WW3 model has a better capability for simulating waves, which is an input/dissipation source term established by Tolman and Chalikov [27] and is labelled as TC96. The validation of simulated results against measurements from moored buoys indicated a 0.79-m root mean square error (RMSE) of SWH for typhoon Fung-wong (2014) and a 1.12-m RMSE of SWH for typhoon Chan-hom (2015).
Stokes drift is the mean Lagrangian velocity of a water particle traveling in the direction of wave propagation [28]. The volume of Stokes transport is obtained by integrating the Stokes drift over all wave directions. Stokes transport is an important component of the ocean circulation system [29]. By quantitatively comparing the global Stokes transport with Ekman transport, in high-wind dominant seas, Stokes transport reaches a magnitude equal to that of wind-driven circulation transport [30]. A recent study also concluded that the ratio of Stokes transport to Ekman transport can be as high as 50% [31]. When wave-induced Stokes drift is calculated using simulated results from wave models, a suitable model and optimal wind forcing field should be selected to achieve reliable wave simulation.
Ocean circulation models are important research tools for the ocean dynamics system, especially for global climate change [32]. The Princeton Ocean Model (POM) is an ocean circulation model that has been widely adopted for operational applications [33,34]. In recent studies, the POM model is successfully used for the water level field simulation in the Bohai sea and the analysis of summer circulation structure in the Beibu Gulf [35,36]. The POM also could be applied for the simulation in typhoon conditions [37]. The new parallel version of POM, denoted as Stony Brook Parallel Ocean Model (sbPOM) [38], is now available to public researchers. It has good scalability as compared with its predecessors. Global sea surface temperature (SST) can be simulated using the sbPOM model [39]. This model was used to study SST cooling in the China Seas during the typhoon events [40]. Stokes drift enhances the shear instability of the upper ocean and strengthens the mixing effect at the air-sea layer, indicating that the Stokes drift affects the SST. Some studies have focused on the effects of Stokes drift on the mixed layer in the upper ocean [39]. The influence of Stokes drift is supposedly included in the sea surface dynamics simulation in order to improve the accuracy of SST.
In this study, we analyzed the long-term variability of Stokes transport induced by ocean waves in the Western Pacific Ocean during the period from 1989 to 2019. The data are briefly described in the "Materials and Methods" section, as are the settings for the SWAN, WW3, and sbPOM models. The calculation of Stokes transport is also presented in this section. The accuracy of the ocean wave modeling and the effect of Stokes transport on SST simulation using the sbPOM model is discussed in the "Results" section, and the conclusions are summarized in the "Conclusions" section.

Materials and Methods
Three types of atmospheric-marine data are open-access for worldwide research, including ECMWF reanalysis [41], Cross Calibrated Multi-Platform Version 2.0 (CCMP V2.0) from Remote Sensing Systems (RSS) [42], and the NCEP Global Forecast System (GFS). Although winds in the above datasets are commonly selected as forcing fields for ocean wave modeling that aims to analyze Stokes drift and Stokes transport over global seas, the performances of these resources necessitate the further study of regional seas using two well-known numeric wave models.

Wind and Wave Data
The ECMWF reanalysis is a global atmosphere-ocean forecast system that has provided operational products since 1958. We analyzed the global wave distributions simulated from SWAN and WW3 [43] models using ECMWF wind data. We used 6-h ECMWF Interim winds with a typical spatial resolution of a 0.125 • grid in both longitude and latitude directions. The CCMP wind data are a Level-3 ocean vector wind analysis product, produced using satellite observations, the measurements from moored buoys, and the simulations from the numeric model. Using 6-h CCMP wind data with a spatial resolution of a 0.5 • grid as a forcing field, the waves in the China Sea were simulated using WW3 and SWAN models. The NCEP Global Forecast System (GFS) system provides wind fields at 6-h intervals from 1948 to the date on a grid with a spatial resolution of 1.0 • . Collectively, the CCMP, the ECMWF, and the NCEP GFS winds ensure the same spatial resolution, which is suitable for studying the performances of wave simulations. Figure 1 shows the wind maps from CCMP, ECMWF, and NCEP GFS at 18:00 UTC on 1 July 2015. Figure 2 shows the temporal evolution of the averaged wind stress magnitude (a), wind stress curl (b), and wind speed (c) for the three kind of products in 2015. The calculation of wind stress is described in [44]. There are apparent differences at wind speeds greater than 10 m/s, especially for CCMP.   The advantage of satellite altimeters is to measure sea surface height and other characteristics of the ocean surface, such as wind and tide. These characteristics are useful for understanding the global ocean and deriving the ocean bathymetry. With global coverage between 66 • N and 66 • S latitude and a 10-day repeat of the ground track, the Jason-2 altimeter maps more than 90% of the global ocean every 10 days. The geophysical data record (GDR) wave data, from the Jason-2 altimeter proceeded to Level-2 (L-2) swath product, are available for 2009-2016. The wave product derived from Jason-2 altimeter data has been used for validating the model simulations, as conducted in our previous studies [23,45]. The accuracy of the sea surface height measured by the Jason-2 altimeter is currently within 3.4 cm. Under these circumstances, the global wave product of the Jason-2 altimeter has sufficient accuracy for studying the large-scale ocean waves, such as swell pools [46]. Measurements from the Jason-2 altimeter in 2015 were used to validate the model-simulated results. Figure 3 shows the significant wave height (SWH) in the satellite's path from July 1 to July 30 2015, overlaid with water depth.

Settings of SWAN and WW3 Models
Generally, the numeric wave models, SWAN and WW3, follow the unique wave propagation balance equation: where N is the wave action density spectrum, S represents the net sources and sinks from the two-dimensional spectrum at wave numbers k and directions θ, σ is the intrinsic frequency, and x and t are the space and time coordinates, respectively. S comprises the input and dissipation source terms, which can be expressed as follows: where S in includes an atmosphere-wave interaction term of the function for the windinduced wave growth, S bot is the friction induced by the wave-bottom interaction, S db is the wave decay due to white-capping and depth-induced wave breaking, S nl is the nonlinear wave-wave interaction term, and S tq is the three-(triad) and four-wave components (quadruplets) of the wave-wave interactions. The complete description of these terms is provided in the technical manual of the SWAN and WW3 models [47]. The default settings of the SWAN model include [23]: WESTHuysen (non-linear saturation-based white-capping combined with wind) for the propagation scheme and switches QUADrupl, TRIad, and BREakinge for bottom friction. The unstructured grid of the simulated area for the SWAN model is shown in Figure 4. The four locations, A-D, were selected for further analysis. The WW3 model provides three kinds of packages for solving the non-linear term of quadruplet wave-wave interactions, that is, Discrete Interaction Approximation (DIA), Full Boltzmann Integral (WRT), and Generalized Multiple DIA (GMD). GMD has two kinds of coefficients, herein called GMD1 and GMD2. Shao et al. (2018) [45] concluded that the GMD2 package is the best option because minimum error is achieved when comparing the simulated SWH from the WW3 model with that from the ECMWF and the Jason-2 altimeter. Therefore, the switch ST2 and the GMD2 package for the WW3 model are selected in the setting of the WW3 model in this study. The outputs from the two models (SWH, dominant wave propagation direction, and mean wave period) have a spatial resolution of a 0.25 • grid at a 1-h interval.

Settings of the sbPOM Model
The sbPOM model is based on an improvement of the well-known Parallel Ocean Model (POM model). It is a three-dimensional, primitive-equation. In the vertical, the sbPOM uses a bottom-following sigma coordinate system: where H( in which U ss is the drift rate of Stokes on the surface, k is the unit wavenumber vector of the fluctuation, H s is the significant wave height, T is the mean wave period, g = 9.8 m/s 2 , and z is the water depth where the sea surface is assumed to be 0 with upwards being positive. The Simple Ocean Data Assimilation (SODA) dataset was produced by an analysis system developed by the University of Maryland in the 1990s, which provides an improved estimation of ocean data from those based solely on observations and numeric simulations. The monthly average SST and sea surface salinity from SODA data with a spatial resolution of a 0.5 • grid are set as the initial field of the model. The temperature difference between day and night is significant in July due to long sunshine time and strong solar radiation. Figure 5 shows the monthly average sea surface temperature and sea surface salinity in July 2015 derived from SODA data. The NCEP/NCAR Reanalysis total heat flux (latent heat flux, sensible heat flux, long-wave radiation, and solar radiation) is treated as the upper boundary forcing field of the model, which has a 6-h time resolution and a spatial resolution of 1.875 • × 1.905 • in the longitude × latitude direction. The total net heat flux is directly calculated by NCEP heat data according to the formula: in which Q s is the solar radiation, Q b is the long-wave radiation, Q s is the sensible heat radiation and Q h is the latent heat flux. It is observed that the relative larger gradient of SST appears at near-shore waters, for example, in the Bohai Sea, East China Sea and Japan Sea. In order to ensure the stability of the simulation results from the sbPOM model, the time of the outer mode (positive compression mold) was set to 20 s and the time of the inner mode (oblique compression mold) was set to 600 s. The simulation period is from January 2012 to December 2015. The bathymetric topography was based on the General Bathymetry Chart of the Oceans (GEBCO), which has a minimum water depth of 10 m. The maximum water depth was set to 5000 m, and water depths exceeding 5000 m were replaced by 5000 m to match the maximum water depth in the SODA data. The spatial resolution of the simulated sea surface temperature from the sbPOM model is a 0.25 • grid at a 1-h interval. The shore boundary adopts the solid wall boundary, and the lateral open boundary adopts the definite boundary conditions, which are directly given by the SODA data. The wind field and heat flux are used as the upper boundary conditions of the model, and the monthly average temperature, salinity and flow field data of SODA are used as the initial field of the model. The basic settings for sbPOM are listed in Table 1.  A global array of temperature/salinity profiling floats, known as Argo, is a major component of the ocean observation system that began in 2000. In this study, high-quality sea surface temperature measurements from the Argo data were used to validate the simulations from the sbPOM model. Figure 6 shows a map of NCEP total heat flux at 18:00 UTC on 1 July 2015. The triangles represent the geographic locations of data collected by Argo [48].

Calculation of Stokes transport
The Stokes drift profile u st (z) based on a one-dimensional wave spectrum is as follows [49]: where g is the gravitational acceleration, ω is the angular frequency, k is the wave number, z is the water depth at which the sea surface is assumed to be 0, and S(ω) is the onedimensional frequency spectrum. The sea surface Stokes drift velocity u 0 is as follows: Therefore, the Stokes transport U st is stated as follows: For deep water, z → −∞; the profile of Stokes drift becomes: As mentioned above, the Stokes drift can be calculated using the one-dimensional wave spectrum S(ω), which is obtained from numeric wave modeling or a parametric wave function. Alternately, Stokes drift superimposed with the environment flow field vector and the calculation formula can be expressed as follows: u st (z) = u 0 e 2kz (13) u 0 = a 2 ωk· D (14) where u st (z) represents the Stokes drift caused by waves, u 0 represents the surface Stokes drift, ω is the frequency, k is the wave number, z is the ocean depth, D represents the direction of wave propagation, and a is the amplitude of waves. Using the diffusion relationship σ 2 = kgtanh(kh), the Stokes drift can be derived from the significant wave height Hs and the mean wave period T m . The specific expression form is as follows: The existence of Stokes drift causes a non-zero net transport, which is called the Stokes transport. The vertical integration of the Stokes drift with depth provides the Stokes transport T s . The specific expression is as follows: We used the equation above to calculate the Stokes transport T s from the simulated wave height H s , dominant wave propagation, and the mean wave period T m derived from the numeric wave model. Figure 7 shows six SWH maps of the results from the SWAN and WW3 models for 1 July 2015 at 12:00 UTC using three forcing winds: ECMWF (Figure 7a,d), CCMP ( Figure 7b,e), and NCEP (Figure 7c,f). Due to the low resolution of the NCEP wind field, there are few data for the Sea of Japan. In general, simulations using ECMWF wind data produce SWHs that are relatively higher than those produced using CCMP and NCEP wind data. This is particularly true for the East Japan Sea, although the patterns of simulated wave fields are almost consistent. There are more than 2000 match-ups, which are available for the statistical analysis of model-simulated SWH during the period from January to December 2015. The results indicate a 0.17-m RMSE (Figure 8e), which is less than thẽ 0.6-m RMSE of SWH from the SWAN model using three types of winds (Figure 8a-c) and the~0.2-m RMSE of SWH from the WW3 model using ECMWF wind data (Figure 8d) and NCEP wind data (Figure 8f). The large error between the two models is probably caused by the low spatial resolution of the NCEP wind data. The ECMWF Interim data contain numerous of uncertainties and errors. This is probable explanation that the accuracy of wind speed for ECMWF is slightly worse than that for CCMP. The performance of wave simulations could be improved by taking re-analysis hourly ERA-5 winds. In this study, CCMP wind data are recommended for wave simulation using the WW3 model.     Figure 9 shows the monthly average Stokes transport in the China Sea and the Japan Sea in January, April, July and December of 2015, based on simulations from the WW3 model using CCMP wind data. It is clearly observed that the Stokes transport caused by the waves is higher in the South China Sea and the Japan Sea (up to 2 m 2 /s in January) than that in other regions. The Stokes transport is relatively weak (<0.8 m 2 /s) in July, which dominates the east end of Taiwan in December (>1 m 2 /s). The CCMP wind data were applied for the WW3 model in 1989-2018. Interannual variation was studied at four spots in the China Sea (Figure 9). Spot A is located in the Bohai Sea, spot B in the Yellow Sea, spot C in the East China Sea, and spot D in the South China Sea. As shown in Figure 10, the Stokes transport generally increases at all spots from 0.25 m 2 /s in 1989 to 0.35 m 2 /s in 2018. However, this kind of behavior is relatively weak at spot A, located in the semi-enclosed Bohai Sea. Winds and waves are internally correlated. The local winds determine the development of the wind-sea. In other words, the generation of waves affects the Stokes transport. With respect to Spot A, located in the inner oceanic region (Bohai Sea), the Stokes transport should not change significantly. It is not surprising that the Stokes transport generated by swell in the inland sea is restricted by the fetch. This is probably the reason that Stokes transport in this area has had a steady pattern over the past two decades. As mentioned in the Introduction, Stokes transport is an important component of global wind-driven circulation. It also affects the distribution of SST. Figure 11 shows a time series of SST from Argo (ID: 2901578) for the period from 30 June to 31 July 2015. The blue and red lines represent sbPOM simulations with Stokes transport and without Stokes transport, respectively. Again, the large difference in temperature in July is suitable for analyzing the SST change induced by Stokes transport. Moreover, the 'V'-type tendency of SST (black line) measured by Argo (ID: 2901578) is observed. A small change in the SST magnitude from the sbPOM model without Stokes transport is present. Beginning on 10 June, Stokes transport is consistent with that of the Argo (ID: 2901578). Based on the seasonal analysis, the RMSE of the SST are listed in the Table 2. It is found that the Stokes transport plays an important role in Summer (March to June) and Autumn (July to September). For example, there is about a 0.5 • C difference between the SST simulations with Stokes transport and SST simulations without Stokes transport. In particular, the difference between RMSE without Stokes transport and that with Stokes transport is about 0.7 • C in Summer, while the difference in RMSE is about −1.0 • C in Autumn. This kind of behavior is probably caused by the tropical cyclones in Autumn [40]. Figure 12 shows a comparison of SST from the Argo and the sbPOM in 2012-2015. The statistical analysis indicates a 3.25 • C RMSE with Stokes transport, which is less than the 3.64 • C RMSE without Stokes transport. From this, we concluded that Stokes transport should be included in the SST simulation from the sbPOM model. In doing so, the accuracy of SST is improved by as much as 0.5 • C.

Conclusions
Because of the non-closed trajectory of a sea surface water particle, the net displacement of particles in the direction of wave propagation has a nonlinear effect, called the Stokes drift. The Stokes drift causes non-zero mass transport called the Stokes transport, although the Stokes transport is relatively small [49]. Besides, the Stokes drift produces the Coriolis-Stokes force in the background of planetary voracity, which affects the distribution of the sea surface current field. In other words, wave-induced Stokes drift should be included in SST simulations when using the sbPOM ocean circulation model. In this study, the interannual and seasonal variability of transport induced by Stokes drift during the period from 1989 to 2019 is investigated. Specifically, the wave-induced Stokes transport is included in the SST simulation using sbPOM.
The WW3 and SWAN are well-known numeric models for wave simulation. The WW3 model is commonly used for global wave analysis. The SWAN model has an unstructured triangular grid that makes it particularly useful in near-shore conditions. The applicability of the two models in regional seas using various types of wind products has not been well studied. In this study, the performances of the two numeric wave models were investigated in the Western Pacific Ocean using ECMWF, CCMP and NCEP GFS wind data. The simulated SWHs were validated against measurements from the Jason-2 altimeter in 2015. The validation shows a 0.17-m RMSE of SWH from the WW3 model using CCMP wind data, which is less than the~0.6-m RMSE of SWH from the SWAN model using three types of winds. Relatively high errors occur when using NCEP wind data, that is, 0.7-m RMSE from the SWAN model and about 0.3-m RMSE from the WW3 model. These errors were probably caused by the coarse resolution of the NCEP wind data. Therefore, the CCMP wind data are employed for wave simulation using the WW3 model in our work.
Sea surface waves were simulated with the WW3 model using 2015 CCMP wind data. The Stokes transport was high in the South China Sea and the Japan Sea (up to 2 m 2 /s in January). Four spots in the China Seas were selected for analyzing interannual variation during the period from 1989 to 2019. The analysis showed that the Stokes transport generally increased from 0.25 to 0.35 m 2 /s. SST was simulated continuously using the sbPOM model. The qualitative analysis confirmed that the simulated SST was closer to the Argo measurements when Stokes transport was considered. Furthermore, the validation of sbPOM-simulated SST in 2012-2015 showed a 3.25 • C RMSE with Stokes transport, which is less than the 3.64 • C RMSE without Stokes transport. Therefore, we concluded that the inclusion of the Stokes transport could improve the accuracy of the SST simulation with the sbPOM. The seasonal analysis yields that the Stokes transport plays an important role in SST simulations, for example, the difference in RMSE is about −0.5 • C in Summer and −1.0 • C in Autumn, respectively. In particular, the Stokes transport has a positive effect in Autumn due to the tropical cyclones.
The Stokes drift is supposedly strong at a high state, which may enhance the energy exchange in the air-sea layer [50]. The characteristics of wave-induced Stokes drift in specific regions, for example, the Arctic Ocean, in tropical cyclones, and the effect of the monsoon will be further studied in the future.

Informed Consent Statement: Not applicable.
Data Availability Statement: Due to the nature of this research, the participants of this study did not agree to their data being shared publicly; so, supporting data are not available.