The Coupled Ocean-Atmosphere-Wave-Sediment Transport (COAWST) modeling system, originally detailed by Warner et al. (2010), integrates multiple numerical models to consider interactions among various environmental processes ( Figure 1 ). In this study, the latest version of the COAWST system is employed, comprising the following model components: 1) the Weather Research and Forecasting (WRF) model, a nonhydrostatic, quasi-compressible atmospheric model utilizing the Advanced Research WRF (ARW) core; 2) the Regional Ocean Modeling System (ROMS), a free-surface, terrain-following ocean model based on hydrostatic and Boussinesq approximations, and 3) the Simulating Waves Nearshore (SWAN) model, a spectral wave model that solves the wave spectral density evolution equation. The Model Coupling Toolkit (MCT), a fully parallelized coupler leveraging MPI, was employed to facilitate the exchange of state variables between each model component. In addition, while COAWST includes hydrological and sediment transport models, they were not activated in current study of TC Choi-wan (2021). Further detailed model descriptions can be referred to: http://woodshole.er.usgs.gov/operations/modeling/COAWST/index.html.

This study carried out coupled numerical simulations of the Tropical Storm Choi-wan (2021) in the western North Pacific and the South China Sea ( Figure 2A ). It originated as a tropical depression on May 30, 2021, about 1,700 km east-southeast of Manila. The storm moved west northwestward, intensifying into a tropical storm by May 31, with maximum sustained winds of 75 km/h at its peak. Choi-wan crossed the Philippines, causing significant rainfall and flooding, before entering the South China Sea on June 2. Over the following days, it turned northward and gradually weakened, eventually transitioning into an extratropical cyclone near the Ryukyu Islands on June 5.

To investigate the effects of wave coupling during storm propagation, four experimental scenarios were designed. The fully coupled configuration (Case 1, Figure 2B ) represents complete integration with two-way exchanges of model variables among WRF, ROMS, and SWAN. WRF provided hydrodynamic models with atmospheric forcing fields, including surface wind vectors (U₁₀, V₁₀), relative humidity (R_h), air temperature (T_air), cloud fraction, precipitation, atmospheric pressure (P_atm), and radiative fluxes, enabling ROMS to calculate heat and momentum fluxes via the COARE algorithm. ROMS, in turn, supplied WRF with SST. SWAN simulated wave parameters, such as significant wave height (H_wave) and wave period (T_surf), which were incorporated into WRF through a wave-dependent roughness scheme (e.g., Oost et al., 2002). The WRF-ROMS coupling (Case 2, Figure 2B ) excluded SWAN to isolate atmospheric-oceanic interactions. In the WRFSWAN one-way coupling (Case 3, Figure 2B ), wave fields simulated by SWAN were not fed back into WRF, with SWAN driven solely by WRF-generated winds to avoid inconsistencies in external wind data. The WRF-ROMS-SWAN coupling without wave-current feedback (Case 4, Figure 2B ) included all three models but excluded exchanges between SWAN and ROMS. These configurations systematically evaluate wave coupling effects, providing critical insights into wave-atmosphere-ocean interactions in TC modeling.

Specifically, in all model scenarios, the WRF model was configured with a single grid at a spatial resolution of 15 km, featuring 50 vertical layers extending from the sea surface to 50 hPa. Initial and boundary conditions were sourced from the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis version 5 dataset (Hersbach et al., 2020). To resolve mesoscale physical processes, the Single-Moment 6-class scheme (Hong et al., 2006) was employed for microphysics, and the Kain-Fritsch scheme (Kain, 2004) was used for cumulus convection. The RRTM model (Mlawer et al., 1997) and the Dudhia model (Dudhia, 1989) were applied for longwave and shortwave radiation, respectively. For ocean responses, the ROMS and SWAN models utilized the same curvilinear grid, with a horizontal resolution of 2 km. Bathymetry was interpolated from a combination of the General Bathymetric Chart of the Oceans gridded bathymetry dataset (https://www.gebco.net/) and nautical charts provided by maritime authorities. To better resolve vertical mixing during tropical cyclones, ROMS employed 16 terrain-following vertical layers, with vertical mixing parameterized using the generic-length-scale turbulence closure model (Warner et al., 2005). In SWAN, the computational grid was discretized into 36 directional bins and 24 frequency bins spanning 0.04 Hz to 1.0 Hz. Wave growth due to wind was modeled using the exponential formulation of Komen et al. (1984). Bottom friction was parameterized with the Madsen et al. (1988) model, applying a friction coefficient of 0.035. White-capping dissipation was represented using the approach of Komen et al. (1984), with a dissipation rate of 2.36×10^-5. Additionally, in scenarios incorporating wave coupling, surface wind stress was computed using two distinct parameterization schemes. The first, referred to as CHNK, is a wind-speed-dependent approach based on the Charnock relation (Charnock, 1955). The second, referred to as OOST, employs a wave-dependent surface roughness scheme as proposed by Oost et al. (2002).

The fully coupled model scenario incorporating the wave-dependent surface roughness scheme (Case 1) was evaluated against observational data. The simulated track of tropical storm Choi-wan was validated against the best-track data provided by the China Meteorological Administration (CMA, https://tcdata.typhoon.org.cn/en/zjljsjj.html), as illustrated in Figure 3A . Additionally, modeled wind speeds were compared with the Global Ocean Hourly Reprocessed Sea Surface Wind and Stress dataset from the European Centre for Medium-Range Weather Forecasts (ECMWF, https://data.marine.copernicus.eu/product/WIND_GLO_PHY_L4_MY_012_006/description), as shown in Figure 3B . A Pearson correlation coefficient (R) of 0.84 and a Root Mean Square Error (RMSE) of 2.06 m/s were obtained for the wind field. Similarly, the wave field simulations were validated using ECMWF’s Global Ocean Waves Reanalysis data (https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_WAV_001_032/description), as depicted in Figure 3C , yielding R and RMSE values of 0.79 and 0.35 m, respectively. Additionally, the modeled water level time series were compared with in-situ observations from the Currimao station in the Philippines ( Figure 3D ). During the initial 12 hours (00:00–12:00 on June 3), discrepancies in water levels were observed, potentially attributed to errors in the initial conditions during the model’s warming-up period. These inaccuracies can prevent the model from accurately capturing water level dynamics. Additionally, the atmospheric model’s difficulty in resolving the track of a mild vortex may lead to inaccuracies in the wind field, contributing to an underestimation of water levels. However, from 12:00 on June 3 to 12:00 on June 4, the model accurately reproduced the observed water levels. Overall, the modeled tropical cyclone characteristics and associated oceanic responses aligned well with observational data. Nonetheless, discrepancies were noted, particularly in the meteorological fields, which could stem from insufficient boundary value accuracy for a relatively weak tropical cyclone. Biases in the wave field may also be linked to errors in tracking the tropical cyclone’s trajectory.

Modeled water temperature was compared with in-situ observations at point S3 ( Figures 4B ). Akima interpolation was employed to align the simulation data with observation times and depths. As shown in Figure 4A , the modeled temperatures during tropical storm Choi-wan were generally lower than the observations, with the average error across the seven experimental cases remaining below 0.3°C throughout the 36-hour simulation period. Despite the limitations of Akima interpolation, this accuracy underscores the reliability and stability of the simulations. The performance of the experiments was assessed by analyzing normalized temperature results during two distinct periods: before the tropical storm approached point S3 ( Figure 4C ) and at its closest proximity to S3 ( Figure 4D ). On average, the modeled sea temperatures were lower than the observations, indicating that higher modeled temperatures correspond to better performance. During the first period, Case 2 produced the highest sea temperatures, while during the second period, maximum temperatures were observed in Case 1. The results in Figure 4D also reveal that the OOST parameterization scheme outperformed the CHNK scheme across both coupling configurations (as described in Figure 2 ). Overall, the fully coupled OOST parameterization (Case 1-OOST) yielded the most accurate temperature results, effectively capturing observed SST variations. Conversely, simulations involving only wave-current coupling (Case 3) produced comparatively lower temperatures, highlighting the importance of fully coupled parameterizations in accurately representing tropical cyclone induced thermal responses.

Furthermore, Case 2, which excluded wave effects, served as the control group. Figure 5 illustrates the differences in normalized vertical sea temperature over time at point S3 during tropical storm Choi-wan (from 00:00 on June 3 to 12:00 on June 4) across all experimental groups compared to Case 2. When Choi-wan was closest to point S3 (03:00–04:00 on June 4, marked by the red dashed rectangular box), temperature differences were primarily concentrated in the upper layers, reflecting the strong wind-driven mixing that entrains cooler water from below the thermocline into the surface layer. Subsequently, bottom-layer temperatures consistently increased relative to Case 2, suggesting that wave effects enhanced vertical mixing and facilitated the downward transport of heat. Taking the TC lag effects into account, the OOST parameterization showed superior performance in promoting vertical mixing compared to the CHNK scheme, evidenced by larger temperature differences across all coupling configurations ( Figures 5D-F ). The wave-dependent surface roughness scheme in OOST effectively modified surface stress, enhancing momentum transfer into the water column and leading to more pronounced thermal adjustments. Moreover, in Case 3 with the OOST parameterization scheme ( Figure 5E ), stronger vertical mixing caused a noticeable temperature drop after the storm passage, with larger areas showing differences exceeding 0.01 relative to Case 1 in the middle and bottom layers. This implied the impacts of wave-induced turbulence in facilitating upwelling, which redistributed heat by cooling the surface while transporting it into deeper layers. Such processes are essential in shaping the post-storm thermal structure. However, neglecting wave effects (Case 2) underestimated vertical mixing and heat redistribution, leading to less accurate simulations of storm-induced cooling and recovery. Fully coupled models incorporating advanced wave parameterizations, such as OOST, demonstrated improved capabilities in capturing these dynamics. By accurately representing interactions between waves, currents, and atmospheric forcing, these models enhance predictions of SST changes during extreme events. These changes are critical for understanding tropical storm intensity, feedback mechanisms, and marine biogeochemical processes. The results underscore the necessity of coupling schemes that account for wave dynamics to accurately simulate the complex thermal responses of the ocean during tropical cyclones.

To further explore the role of wave coupling in tropical cyclone simulations, the vertical temperature distribution of Case 2 along the tropical storm’s track (at 3-hour intervals, Figure 6A ) is shown in Figure 6B . The results clearly demonstrate the effect of vertical mixing on sea temperature and highlight a lag in the oceanic response to tropical storm intensity. The deepest downward extension of the temperature anomaly occurred approximately 8 hours after the maximum wind speed at the tropical storm’s center began to decline. However, this depth persisted for only about 2 hours before gradually diminishing. This delay in the response underscores the complex interaction between atmospheric forcing and ocean dynamics, where momentum transfer and mixing processes continue to affect the thermal structure even after the storm’s peak intensity has passed.

Figure 7 presents the differences in vertical sea temperature distribution over time between the various experiments and Case 2. The results show that the impact of wave coupling on sea temperature was minimal when simulating a tropical storm of weak intensity. Temperature variations across all experiments and vertical levels were generally below 0.2°C, with most of the changes concentrated in the upper layers. However, the OOST parameterization (lower panel of Figure 7 ) consistently produced stronger vertical mixing compared to the CHNK parameterization scheme (upper panel of Figure 7 ). The temperature variations in the OOST configuration were more widely distributed both temporally and spatially, with higher values observed.

In Case 1 and Case 4 ( Figures 7A, D , 7C, F ), the differences between the two parameterization schemes were primarily observed at the deepest temperature intrusions, with smaller temperature changes in the mid-layers due to OOST mixing. The most significant differences occurred in Case 3 ( Figures 7B, E ), where negative variations exceeding 0.04°C were observed in the CHNK parameterization at upper levels, while the OOST parameterization showed predominantly positive changes. As the tropical storm intensified, vertical mixing became more pronounced, enhancing the influence of wave coupling on the simulation results. Under these conditions, the advantages of the OOST parameterization became increasingly evident, implying its superior performance in capturing wave-induced vertical mixing during relatively intense cyclonic events.

Comparisons of the experiments at point S3 and along the tropical storm’s track revealed positive outcomes of wave coupling, with the OOST parameterization exhibiting notable advantages. To evaluate model performance more comprehensively, the overall response of different schemes at a representative time step (03:00 on 4 June) was selected for comparison and analysis.