The WRF model is a quasi-compressible, non-hydrostatic model. A single-domain setup was used for the study area, spanning 110°E to 135°E and 5°N to 35°N, with a horizontal resolution of 27 km and a temporal resolution of 3 hours. The model consisted of 38 vertical levels, with the top boundary pressure set to 50 hPa.

The ROMS is a high-resolution stretched terrain-following hydrostatic model (Haidvogel et al., 2000; Shchepetkin and McWilliams, 2005). The model was run using WRF simulation results. Bathymetric data for the model were obtained from GEBCO. The vertical grid configuration was setting to 35 layers with increased resolution in the surface layer, thermocline, and deep water areas.

The SWAN model is a model established for nearshore wave simulation (Booij et al., 1999), and subsequently improved in later versions. SWAN is a non-hydrostatic wave model designed to simulate wave characteristics across complex terrains and varying water depths. In this research, the SWAN model grid was configured identically to that of ROMS. SWAN was run both independently and in coupled mode, using wind fields from the WRF model and ocean currents from the ROMS model.

The COAWST model (Warner et al., 2010) integrates the WRF model, the ROMS, and the SWAN model, coupled through the MCT, enabling simultaneous simulation of atmosphere-ocean-wave interactions. In this coupling process, the WRF model supplies ROMS with wind speed and atmospheric pressure, while receiving sea surface temperature from ROMS. Additionally, WRF supplies the wind speed to drive the SWAN model, influencing wave dynamics, and receives Hs and wavelength from SWAN. The SWAN model supplies ROMS with data on Hs, wavelength, direction, and period, and receives current velocity and depth. Data exchange between the component models occurs at 10-minute intervals.

The simulation duration was 120 hours (12 UTC 23 July 2023 to 12 UTC 28 July 2023) during which Super Typhoon Doksuri passed through the coastal waters of China.

To investigate the spatiotemporal evolution characteristics of Stokes drift induced by typhoon waves, a total of four experiments were designed. The first three simulations used the WRF, ROMS, and SWAN models individually. Subsequently, we compared these three models with COAWST. Table 1 summarizes the configurations of these component models.

Two primary methods currently exist for calculating Stokes drift. One approach uses wave spectra to calculate Stokes drift (Breivik et al., 2014), while the other uses Hs, period, and wavenumber. In this study, the second method was adopted for calculating Stokes drift. The second method was chosen for its simplicity and practicality, allowing for easier implementation using Hs, wave period, and wavenumber. This approach is also consistent with previous studies, where it has been widely used to model Stokes drift and provides a reasonable approximation. Given the large-scale nature of this study, a simplified approach was preferred to ensure computational efficiency without sacrificing significant accuracy. The vertical influence of Stokes drift reaches a maximum depth, referred to as the depth of influence. When the water depth exceeds this influence depth, Stokes drift and its effects become negligible. The equation is as follows:

where d s represents the depth of influence of Stokes drift; k is the wave number; λ is the wavelength; c is the phase velocity; T is the wave period.

For a monochromatic deep-water gravity wave under conditions of homogeneous and irrotational flow, Stokes drift is expressed by the following equation:

where u → s t   represents wave-induced Stokes drift velocity; U → s  represents the surface Stokes drift velocity; z represents depth; d_s  represents the depth of Stokes influence. During the influence of Super Typhoon Doksuri, the large typhoon waves triggered significant Stokes drift. Assuming the fluid is irrotational and inviscid, Stokes drift velocity can be determined using Equation 5.

The movement of water particles driven by Stokes drift is known as Stokes transport, which is essential for mass and energy exchange in the upper ocean. Classical coastal hydrodynamics theory suggests that greater Hs and shorter wave periods lead to enhanced mass transport via Stokes drift. Stokes transport is determined by vertically integrating Stokes drift, as represented by the following equation:

where T s represents Stokes transport; u → s t represents Stokes drift.

Stokes transport facilitates the mass and energy transfer in the upper ocean, playing a significant role in global ocean transport processes. To compare the magnitude of Stokes transport with total transport (McWilliams and Restrepo, 1999), the Ekman-Stokes number (Es) is defined, which represents the proportion of Stokes transport to net Ekman transport. The specific equation is as follows:

where   T → E represents net Ekman transport; T → s = τ ρ w f , τ represents wind stress magnitude; ρ w represents seawater density; f represents the Coriolis parameter; U 10 represents the wind speed at 10 m above sea level; κ is the von Karman constant, set at 0.4; C D is the drag coefficient (Wu, 1980):

Verification and analysis of wind and wave field simulation accuracy confirm the reliability of COAWST in simulating Stokes drift. Stokes drift velocity characteristics under typhoon conditions are analyzed at four typical times ( Figures 7A–D ). In the early stage of typhoon formation ( Figure 7A ), low wind and wave intensities resulted in a Stokes drift of about 0.35 m/s. As waves intensified ( Figure 7B ), Stokes drift increased, showing a rightward bias with a maximum near the typhoon center, reaching 0.67 m/s in crescent-shaped contours along the right of the path. The drift field’s vortex structure also intensified, indicating strong wave energy transfer. Due to Taiwan Island’s obstruction ( Figure 7C ), the vortex structure fragmented, with Stokes drift peaking at 0.73 m/s southeast of the island and extending over a larger area. Before landfall ( Figure 7D ), Stokes drift decreased to about 0.4 m/s, while areas outside the typhoon’s influence had minimal drift, around 0.1 m/s or less. Under typhoon conditions, Stokes drift is markedly enhanced relative to non-typhoon conditions.

These results show that Stokes drift velocity changes significantly with the movement and intensity of the typhoon, especially in the intense core region where it peaks, enhancing surface material drift in the ocean. Stokes drift distribution is generally positively correlated with wind field intensity and Hs: stronger wind increases energy input, Hs, and wavelength, intensifying wave asymmetry and Stokes drift. This causes peak Stokes drift near the typhoon center. The spatiotemporal variation of Stokes drift velocity, influenced by typhoon-induced wave asymmetry, exhibits a rightward bias with crescent-shaped contours, reflecting the nonlinear wave characteristics and net drift effect.

EOF decomposition of Stokes drift velocity further analyzes the spatiotemporal distribution of Stokes drift ( Figures 8A–D ). The first spatial pattern (EOF1) represents an average field characteristic, showing the primary spatial distribution pattern, with the corresponding temporal mode time coefficients generally matching the periods of Stokes drift velocity enhancement or reduction. EOF1 indicates that Stokes drift in this region shows significant spatial non-uniformity under the influence of the typhoon ( Figure 8A ), explaining 62.40% of the total variance. EOF1 of the entire Stokes drift velocity field is positive, suggesting that wave energy during the typhoon tends to concentrate in specific areas. The values are relatively highest in the southeastern Taiwan sea and coastal areas of mainland China, ranging from approximately 0.15 to 0.2, indicating a significant impact from the typhoon.

The first principal component (PC1) starts with a rapid increase from a large negative value (around -2), reaching a peak (approximately 1), and then gradually decreases, stabilizing around zero ( Figure 8C ). The trend of the temporal mode reflects the temporal variation of the typhoon’s impact on Stokes drift velocity in the study area. The rapid increase phase of the curve corresponds to the strengthening of the typhoon, while the plateau phase corresponds to the peak stage of the typhoon. The subsequent decrease indicates a weakening of the typhoon’s influence, with the drift velocity gradually returning to normal levels.

Combining PC1 with EOF1 allows for an analysis of Stokes drift velocity under typhoon conditions. At the onset of the typhoon, the PC1 is negative, indicating that the corresponding Stokes drift velocity at these moments is in an opposite phase to the spatial mode associated with the temporal coefficient at these times, with a generally lower Stokes drift velocity across the field. Over time, as the typhoon’s impact reaches the sea area near Taiwan, PC1 rapidly rises to a peak, corresponding to the positive center near Taiwan. When PC1 declines again, it indicates that the typhoon has moved away, and the Stokes drift velocity across the study area gradually returns to near-normal levels.

Figure 8B displays the second spatial pattern (EOF2) of the EOF analysis of Stokes drift velocity under typhoon conditions, which accounts for 25.64% of the total variance. EOF2 represents the primary pattern of temporal changes in Stokes drift velocity, rapidly progressing from southeast to northwest roughly along the typhoon path, it first increases and then decreases. The increasing (red) regions are mainly distributed in the northwestern Pacific east of the Philippines, while the decreasing (blue) regions are concentrated in the South China Sea and East China Sea. The increasing and decreasing regions exhibit a northwest-southeast directional distribution. As the typhoon progresses, the wind field distribution typically generates different wind stresses along both sides of the track, leading to an asymmetric distribution of Stokes drift velocity in the ocean surface layer. Topographic differences lead to significant variations in Stokes drift response across these regions. The terrain and current characteristics between Taiwan Island and Luzon Island may also modulate the distribution of Stokes drift velocity, such as the weakening of Stokes drift velocity in coastal areas due to the blocking effect of terrain.

The second principal component (PC2) is shown in Figure 8D . Combined with EOF2, the second mode is positively correlated with the overall Stokes drift velocity at the beginning and end phases, while it is negatively correlated at the intermediate times.

The first and second modes account for 88.04% of all the modes. The first mode primarily reflects the overall impact of the typhoon on Stokes drift velocity across the field, while the second mode captures the asymmetric changes across both sides of the typhoon track as well as the dynamic temporal evolution of drift velocity. Further analysis of EOF modes for different years and seasons could help distinguish the varied impacts of typhoons on Stokes drift in this region under different conditions.

Stokes drift drives transport phenomena, playing a crucial role in mass and energy exchange in the upper ocean. Stokes transport during the typhoon, calculated using Equation 6, is shown in Figure 9 . Consistent with the Stokes drift distribution in Figure 7 , high-value areas of Stokes transport align with those of drift, mainly along the typhoon track with the maximum on the right side. Stronger wind and wave fields lead to greater Stokes drift, enhancing mass and energy transport. In the early typhoon stage ( Figure 9A ), Stokes drift and transport are relatively weak, around 0.2 m²/s. As typhoon intensity increases ( Figure 9B ), Stokes transport intensifies, reaching 5-6 m²/s near the typhoon eye with a distinct crescent-shaped structure, indicating significant wave energy within the strong wind region. During northward movement ( Figure 9C ), transport further increases to 10-13 m²/s, concentrated in the central and northern sea areas, driven by higher wind speeds and Hs. Upon landfall, Stokes transport weakens, with peak values reduced to 5-6 m²/s in the Taiwan Strait, reflecting the typhoon’s decline and decreasing wind and wave forces.

These phenomena are primarily driven by the strong impact of the typhoon wind field on ocean waves. During intensification, high waves from strong winds greatly increased Stokes transport intensity. Wind fields asymmetry and water particle movement created a distinct spatial pattern of Stokes transport along the typhoon track. As the typhoon weakened, decreasing wind and wave strength led to a reduction in Stokes transport intensity.

In ocean regions with strong currents below the surface, the shear effect of the current may have a significant impact on Stokes drift and the associated mass transport. In the case of shear flow propagating obliquely to the wave direction, the maximum induced mass transport no longer occurs at the still water surface but at a certain depth below (Li et al., 2024). The shear of the current can alter momentum transfer and vortex structures between water layers, thus affecting the spatial distribution and intensity of Stokes drift. Particularly under extreme weather conditions, such as typhoons, the interaction between currents and wind-waves may lead to significant ocean dynamic effects in localized areas, thereby altering the transport processes of matter and energy. Therefore, further consideration of the vertical shear effects of currents is crucial for improving the accuracy and predictive capability of Stokes drift models, especially in ocean regions with complex flow fields.

The influence depth of Stokes drift during Typhoon Doksuri, calculated from Equation 1 and displayed in Figure 10 , reflects the spatial distribution of Stokes drift, emphasizing dynamic changes as the typhoon traversed sea areas. In the early stage ( Figure 10A ), Stokes drift depth was relatively shallow and uniformly distributed, reflecting weaker wind field impacts on the surface layer. In the mid-stage ( Figure 10B ), depth increased significantly, especially in the right-front (northeastern) region, due to enhanced wind speed and pressure gradients. At the peak stage ( Figure 10C ), Stokes drift depth reached its maximum, exceeding 20 m, around the typhoon center, indicating intense wave energy transfer and vertical mixing over a broader area. In the late stage ( Figure 10D ), drift depth decreased, as the influence of wind waves and currents weakened. These results show that Stokes drift depth varies significantly with typhoon intensity and position, particularly in the core region, promoting mass and energy exchange, enhancing vertical turbulence, and cooling the upper ocean, thus improving the realism of typhoon simulations.

To investigate the variation of the Ekman-Stokes number during the typhoon process, the analysis was conducted based on the track of Typhoon Doksuri, focusing on the changes of the Ekman-Stokes number along the latitude of the typhoon track. The Ekman-Stokes number at typical times was calculated using Equation 7, as shown in Figure 11 . From the analysis of the Ekman-Stokes number ( Figures 11A–D ), it can be observed that it exhibited significant fluctuations across different latitudinal ranges over time. As the dates progressed, the Ekman-Stokes number generally showed an increasing and then decreasing trend in the north-south direction (latitude range 15°N to 25°N).

Overall, the Ekman-Stokes number decreases sharply with increasing latitude, primarily due to the enhancement of the Coriolis force at higher latitudes, which significantly strengthens the Ekman flow, while the Stokes drift becomes relatively weaker. The relative interplay between these two factors results in the observed sharp decline of the Ekman-Stokes number with increasing latitude.

From July 24 to July 26, the latitude of the first minimum in each figure corresponds to the location of the typhoon’s high-value area at that time, with significant peaks on both sides indicating strong Stokes transport. Subsequently, the Ekman-Stokes number decreases along the latitude. In the early stages of the typhoon at lower latitudes, Stokes transport comprises approximately 27% of the net transport. As the typhoon grows and moves, the proportion of Stokes transport gradually increases. On July 27, in the region from 22°N to 24°N, the Ekman-Stokes number reaches its peak (close to 0.35). Throughout the typhoon’s passage, Stokes drift remains relatively weak at lower latitudes, while wave contributions in higher latitude areas gradually intensify.