E ﬀ ects of Model Coupling on Typhoon Kalmaegi (2014) Simulation in the South China Sea

: Typhoon Kalmaegi (2014) in the South China Sea (SCS) is simulated using a fully coupled atmosphere–ocean–wave model (COAWST). A set of sensitivity experiments are conducted to investigate the e ﬀ ects of di ﬀ erent model coupling combinations on the typhoon simulation. Model results are validated by employing in-situ data at four locations in the SCS, and best-track and satellite data. Correlation and root-mean-square di ﬀ erence are used to assess the simulation quality. A skill score system is deﬁned from these two statistical criteria to evaluate the performance of model experiments relative to a baseline. Atmosphere–ocean feedback is crucial for accurate simulations. Our baseline experiment successfully reconstructs the atmospheric and oceanic conditions during Typhoon Kalmaegi. Typhoon-induced sea surface cooling that weakens the system due to less heat and moisture availability is captured best in a Regional Ocean Modeling System (ROMS)-coupled run. The Simulated Wave Nearshore (SWAN)-coupled run has demonstrated the ability to estimate sea surface roughness better. Intense winds lead to a larger surface roughness where more heat and momentum are exchanged, while the rougher surface causes more friction, slowing down surface winds. From our experiments, we show that these intricate interactions require a fully coupled Weather Research and Forecasting (WRF)–ROMS–SWAN model to best reproduce the environment during a typhoon. and H.L.; visualization, K.T.C.L.K.S.; data curation, H.Z.; supervision, C.D. and H.L.; writing—review and editing, C.D., H.L., H.Z. and R.W.


Introduction
Typhoons (also called tropical cyclones (TC) or hurricanes, depending on location; hereafter TC for simplicity) are one of Earth's most destructive natural disasters, which induce rapid and severe modulations in both the atmospheric and marine environments through heat and momentum energy transfers across the air-sea interface. TCs cause huge loss of lives and billions of dollars worth of damage every year. Accurate TC forecast is the crux of the matter, as the last few decades have witnessed the growth of coastal populations [1], thus increasing the number of people at risk and their

Materials and Methods
This section introduces the observation, reanalysis and model data that were employed in this study, and describes the model settings, experimental design, sensitivity tests and statistical criteria that were used to assess the performance of the different model experiments.

Best-Track Data
Best-track data for Typhoon Kalmaegi are taken from the Japan Meteorological Agency (JMA [47]) and the Joint Typhoon Warning Center (JTWC [48]). The data provided at 6-h intervals includes TC center coordinates (latitude and longitude), P min and V max .

In-Situ Data
Observational data (10-m wind, sea level pressure (SLP), SST, significant wave height (SWH) and ocean current) from 4 moored buoys that formed part of an array (yellow triangles in Figure 1) deployed in the middle of the SCS is employed to validate our simulation results. The location coordinates and the water depth at the observation points are listed in Table 1. More information about the mooring specifications is found in Zhang et al. [49]. Remotely sensed SST data from the optimally interpolated SST product using both microwave (MW) and infrared (IR) SST measurements (MW-IR product version 5, distributed by Remote Sensing Systems [50]) is used to validate the modeled SST. The daily data has a spatial resolution of 0.09 • × 0.09 • .

Reanalysis and Model Data
The National Centers for Environmental Protection (NCEP) Final (FNL) Operational Global Analysis data [51] is used as initial and boundary conditions for the atmospheric model. The data has a temporal resolution of 6 h and a spatial resolution of 1 • × 1 • . The initial (currents, current depth, water level, temperature and salinity) and open-boundary conditions (currents, temperature and salinity) for the ocean model are derived from the HYbrid Coordinate Ocean Model (HYCOM) data [52]. WaveWatch III (WW3) Global Wave Model (NWW3_Global_Best [53]), which is produced by the University of Hawaii, is used as the boundary conditions for the wave model. The hourly data is provided at a spatial resolution of 0.5 • × 0.5 • .

Model and Experimental Design
The model settings introduced in this subsection are the best results from a series of uncoupled model tuning experiments for factors such as initial conditions, model settings and parameterization schemes. It is assumed that the configuration that yields the best results for the uncoupled simulations will also generate the best coupled simulations.

COAWST
This study uses the COAWST Modeling System [41], which has been developed for atmosphere, ocean and wave coupling studies. The coupling system is comprised of three models: the Weather Research and Forecasting (WRF) atmospheric model [54], Regional Ocean Modeling System (ROMS) [55] and Simulated Wave Nearshore (SWAN) spectral model [56]. COAWST allows the inclusion or exclusion of the coupling components. Communication among these components is accomplished through the Model Coupling Toolkit (MCT) [57]. The Spherical Coordinate Remapping and Interpolation Package (SCRIP) [58] is used to generate interpolation weights between the individual model grids for data transmission in the MCT. Previous research employing the COAWST modeling system (e.g., [7,17,46,[59][60][61]) provides useful references for the successful configuration of the current case study. No nesting is implemented. In the coupled model simulations, the time interval for data exchange among the models is set to 600 s. The models are run from 00:00 on 12 September 2014, to 00:00 20 September (8 days) but only results from 12 September, when the depression started to deepen, to 17 September, when the TC made landfall, are analyzed in the present study.

Atmospheric Model
The WRF model is a non-hydrostatic, quasi-compressible, terrain-following, three-dimensional atmospheric model which offers several advanced physical parameterizations for simulating atmospheric mesoscale and microscale motions. In the COAWST modeling system, the Advanced Research WRF (ARW) solver [54] is applied. The bottom roughness in the WRF version of the COAWST modeling system has been modified to consider the effects of ocean surface waves [41]. Without wave coupling, WRF uses Charnock's [62] roughness length to wind stress to parameterize sea surface roughness: where z 0 is the roughness length, α is the Charnock parameter (value set to 0.011), u * is the friction velocity (m s −1 ) and g is the acceleration due to gravity (9.81 m s −2 ). With the addition of wave coupling, surface roughness is computed from Taylor and Yelland's [63] wave-steepness-dependent ocean roughness relation: where H s is the significant wave height (m) and L p is the wave period (s). 1200 and 4.5 are constants calculated in Taylor and Yelland [63]. In this study, the WRF model version 3.9.1 runs on a 6 km × 6 km grid. A compromise is made between high model resolution and computation time. The WRF domain extends from 2-30 • N and 99-140 • E (Figure 1). The SCS spans from the equator to 23 • N and 99-121 • E. Since Typhoon Kalmaegi simulation is initialized on 12 September when it was located at about 135 • E, the model domain is therefore designed to cover this area. The domain is then extended further east and north to give the TC enough space to develop and intensify. Although the current study focuses on the surface level, 58 sigma levels are applied in the vertical. The model is forced at its lateral and surface boundaries by FNL data. In the uncoupled WRF and coupled WRF-SWAN run, SST data from the WRF-ROMS experiment is used to force the atmospheric model every 6 h, while in the ROMS-coupled simulations, SST is directly obtained from the ocean model. The WRF model is initialized at 00:00 on 12 September 2014, and is run until 00:00 20 September (8 days) with an integration time step of 30 s and hourly output.
Cloud detrainment and subgrid-scale convection are computed using the Kain-Fritsch cumulus convection scheme [64]. Grid-scale precipitation processes are simulated using the WRF single-moment six-class moisture microphysics [65]. The Rapid Radiative Transfer Model [66] is used to simulate long-wave radiation, while the Dudhia scheme [67] is used for short-wave radiation simulation. At the planetary boundary layer, the MYNN2.5 turbulence scheme [68] is used. In the uncoupled WRF run, the atmospheric model uses the Charnock parameterization (Equation (1)) for the surface roughness, while the Taylor and Yelland formulation (Equation (2)) is used to compute the wave-induced enhanced bottom roughness [63] in the wave-coupled runs. Spectral nudging [69], with a nudging coefficient of 0.0003 s −1 , is applied to temperature, horizontal zonal and meridional wind, and geopotential height at layers k > 14 every 6 h during the whole simulation. A wavenumber of 3 is used in the zonal and meridional direction, which implies that waves longer than 2000 km are nudged. Simulations without activating spectral nudging were reported to have a poor TC track simulation (e.g., [17,46]).

Oceanic Model
The ROMS [55] model is a free-surface, three-dimensional ocean model that solves the hydrostatic primitive equations for momentum by using a split-explicit time-stepping scheme. A stretched terrain-following coordinate is used in ROMS to discretize the primitive equations over variable topography in the vertical direction. These coordinates allow higher resolutions over the study area. ROMS contains a range of options for pressure gradient, horizontal and vertical advection, and subgrid-scale parameterizations. Additional details can be found in Shchepetkin and McWilliams [55].
This study employs the ROMS version 3.7. The ROMS domain spans 5-24 • N and 99-135 • E (white box in Figure 1) with a grid resolution of 1/12 • × 1/12 • . The domain is slightly smaller than the WRF domain for the WRF domain to fully incorporate ROMS, thus increasing stability at the open-boundaries during model coupling. However, it is computationally more expensive, as transmission time of re-gridded data among model components is increased. Thirty-two (32) levels are used in the vertical for the stretched terrain-following coordinate with the following vertical stretching parameters: Vtransform = 2, Vstretching = 4, θ s = 10, θ b = 0.5, Tcline = 150 and h = 5000 (where θ s is the surface stretching parameter, θ b is the bottom stretching parameter, Tcline is thermocline depth in m and h is the maximum depth in m). These settings ensure a higher vertical resolution of the upper ocean layers, with 4 m intervals between levels from 0-20 m, 5 m intervals from 20-40 m, then 6-25 m intervals between 40 m and 150 m. The northern, eastern and southern lateral boundaries are set as open boundaries. Eight tidal constituents (M 2 , K 1 , O 1 , S 2 , N 2 , P 1 , K 2 and Q 1 ) from the TOPEX/Poseidon global tidal model 8 (TPXO8) [70] are also applied along the open boundaries. The K-Profile Parameterization (KPP) vertical turbulent mixing scheme is employed [71]. The surface-wave-breaking-induced turbulent kinetic energy is incorporated by the surface boundary layer KPP mixing. Sea surface forcing includes wind stress and heat fluxes from the atmosphere. In ROMS, wind stress is calculated from Equations (3)-(5): where C D is the drag coefficient, ρ a is the air density, τ is the surface wind stress and U 10 is 10-m wind.
Friction velocity is calculated using Charnock's relationship [62]: where U 10 is the 10-m wind speed and κ is the Von Kármán constant (value set to 0.4). While wind stress depends on surface roughness, heat flux is regulated by the difference between SST and air temperature. The heat exchange coefficient C K is calculated from Equations (6) and (7): where H K = ρ a u * k is the moist enthalpy flux and k is the specific moist enthalpy, defined as: where q is the specific humidity, c p is the specific heat of air, T is air temperature, L v is the latent heat of vaporization of water, and c pv is the specific heat of water vapor.
In the atmosphere-ocean coupled run, wind stress and surface fluxes are directly provided by the WRF model, and the SST exchanged with the atmospheric model is updated in a dynamically consistent way, unlike in the uncoupled WRF run where SST is fed to WRF every 6 h. Furthermore, a flux-conservative remapping scheme is applied in the WRF-ROMS coupled runs to ensure that ROMS and WRF use the same fluxes at the atmosphere-ocean interface. ROMS is first run from 00:00 1 September to 00:00 12 September (12 days). Typhoon Kalmaegi simulation is performed from 00:00 12 September to 00:00 20 September. The time step for the simulation is 60 s, and the output is written every 1 h.

Wave Model
The SWAN model [56], a third-generation wave model, is used to simulate the ocean surface waves in COAWST. The model version is 41.20. The SWAN model domain co-locates with the ROMS domain (white box in Figure 1). The SWAN model requires U 10 wind but, dynamically, the wave model is driven by the wind stress calculated according to Zijlema [72]: where U ref is the reference wind speed (31.5 m s −1 ). The boundary conditions are taken from the WW3 Global Wave Model. The SWAN model runs in nonstationary mode. It considers the effects of changes in currents [73] and free surface elevations, dissipation processes (such as whitecapping, bottom friction and diffraction) [74], depth-induced wave breaking, and quadruplet nonlinear wave-wave interactions [56]. ROMS provides currents and free surface elevations. Wave growth is computed based on the Komen formulation [75]. SWAN is run from 00:00 1 September to 00:00 12 September. The model is then run from 00:00 12 September to simulate surface waves during the passage of Typhoon Kalmaegi. An integration time step of 60 s and hourly output is applied during the simulation.

Sensitivity Tests
In this subsection, we introduce the sensitivity tests to investigate the effects of model coupling on Typhoon Kalmaegi simulation (Table 1).
From Table 2, W, R and S are the uncoupled model runs for WRF, ROMS and SWAN, respectively; WR is the WRF-ROMS coupled run, WS is the WRF-SWAN coupled model run and WRS is WRF-ROMS-SWAN coupled model run.

Evaluation of Model Performance
To evaluate the performance of the various experiments, we apply two statistical criteria: correlation (R 2 ; Equation (9)) and root-mean-square difference (RMSD; Equation (10)). Correlation estimates the relationship between simulation results and a reference dataset (i.e., observed data or baseline experiment). RMSD shows the goodness of fit by determining the average distance of a data point from a fitted line and is directly interpretable as it has the same units as the variable under consideration. Using this method, we can immediately assess the performance of each experiment compared to the baseline experiment (WRS). R 2 and RMSD are calculated (1) between simulated and observed data to validate the simulation results, and (2) between sensitivity and baseline experiment results to assess the performance skill of each experiment.
Correlation is given by: where var is the variable under consideration; N var is the total number of data points for variable var used in the calculation; i is the data point; e is the experiment index; and m var and r var are the model/experiment data and reference data (observations/baseline), respectively, for variable var. RMSD is calculated using: We further introduce a statistically based modeling skill assessment, a skill score (SS) [76,77], to objectively quantify the performance of the different numerical experiments relative to the baseline Atmosphere 2020, 11, 432 9 of 27 experiment. Using RMSD and R 2 , skill scores SS r (Equation (11)) and SS c (Equation (12)) are respectively calculated for each experiment with reference to the baseline experiment WRS as follows: where b is the baseline experiment and e is the sensitivity experiment. From Equation (11), SS r > 0 when RMSD(e) < RMSD(b) which means an improved simulation. SS r < 0 when RMSD(e) > RMSD(b), implying a poorer model performance. SS r = 0 if RMSD(e) = RMSD(b) (no change in model performance) and SS r = 1 if RMSD(e) = 0 (perfect simulation). From Equation (12), SS c > 0 when r(e) > r(b) which means a better simulation of the variability. SS c < 0 when r(e) < r(b), that is, a poorer model performance. SS c = 0 if r(e) = r(b) (the variabilities in the two experiments are identical).

Results
In this section, the TC simulation sensitivity of different model couplings is presented. The results are taken at the nearest grid point when comparisons are made at buoy locations. JMA and JTWC observation data are used concurrently to demonstrate that observation data from different typhoon centers contain deviations possibly attributable to different estimation methods.

Typhoon Track and Intensity
We first evaluate the simulated TC track, minimum central SLP (P min ) and maximum central sustained wind speed (V max ).
Comparing modeled and observed TC track ( Figure 2), it can be seen that all the experiments accurately simulated the TC track, with only minor departures from observations. A common feature of the simulated tracks is that they meander slightly during the first two days of simulation before stabilizing and following a path more consistent with observations. The WRF uncoupled experiment exhibits a slight northward deviation upon approaching the Philippines. Nevertheless, all the tracks display almost similar patterns, with nearly the same landfall time and location. The track geometric distance difference from the JMA track and TC translation speed given in Figure 3 quantify the results presented above. The simulated tracks' distances from the JMA track are relatively small (Figure 3a), with the largest departure exhibited by JTWC on 18:00 15 September. The simulated translation speed trends are consistent with observations. From 12:00 12 September, the translation speeds steadily increase until 12:00 on 14 September, when Typhoon Kalmaegi passed over the Philippines. Afterward, the TC traveled at about 8.5 m s −1 until it made landfall at Zhanjiang. Surprisingly, WR presents the poorest translation speed simulation (RMSD = 1.68 and R 2 = 0.62) while WS has the best performance (RMSD = 1.14, R 2 = 0.82). This result suggests that translation speed is influenced more significantly by surface roughness than surface temperature. Analyzing WRS, it is seen that this experiment performs a little poorer than WS, probably due to the wave-current interactions occurring in the ROMS-SWAN coupling.
Atmosphere 2020, 11, x FOR PEER REVIEW 10 of 28  . The same colors are used for RMSD and R 2 , which are calculated relative to JMA data. The dots represent data between 00:00 12 September and 00:00 17 September, with 6-h intervals.
Analysis of Vmax shows that, although the model can reproduce its variability (Figure 4a), all the experiments underestimate Vmax after 06:00 14 September when Typhoon Kalmaegi passed over the Philippines (12:00 14 September) and entered the SCS (18:00 14 September). The sudden drop in Vmax between 06:00 and 12:00 14 September is caused by the TC interacting with the land surface. However, our simulations were not able to match the wind speed of the TC in the SCS. The difference between the model and best-track data might be related to the conversion from instantaneous wind to sustained wind speed [78,79]. The WRF uncoupled experiment (W) presents the best model  Analysis of Vmax shows that, although the model can reproduce its variability (Figure 4a), all the experiments underestimate Vmax after 06:00 14 September when Typhoon Kalmaegi passed over the Philippines (12:00 14 September) and entered the SCS (18:00 14 September). The sudden drop in Vmax between 06:00 and 12:00 14 September is caused by the TC interacting with the land surface. However, our simulations were not able to match the wind speed of the TC in the SCS. The difference between the model and best-track data might be related to the conversion from instantaneous wind to sustained wind speed [78,79]. The WRF uncoupled experiment (W) presents the best model Analysis of V max shows that, although the model can reproduce its variability (Figure 4a), all the experiments underestimate V max after 06:00 14 September when Typhoon Kalmaegi passed over the Philippines (12:00 14 September) and entered the SCS (18:00 14 September). The sudden drop in V max between 06:00 and 12:00 14 September is caused by the TC interacting with the land surface. However, our simulations were not able to match the wind speed of the TC in the SCS. The difference between the model and best-track data might be related to the conversion from instantaneous wind to sustained wind speed [78,79]. The WRF uncoupled experiment (W) presents the best model performance (RMSD = 10.35; R 2 = 0.81). It is observed that model coupling decreases V max simulation skills. The more accurate surface roughness estimation in WS might contribute to the slightly smaller MSD. The ROMS-coupled run presents a poorer performance as compared to the SWAN-coupled run, from which it can be inferred that a more accurate representation of the ocean surface roughness in WS improves surface wind simulations.
A comparison of P min from simulations and observations (Figure 4b) shows that the models better simulate P min than V max with R 2 above 0.98 and RMSD less than 5.77. However, a TC with a slightly lower P min as compared to the observations is obtained from the experiments. From our simulations, Typhoon Kalmaegi suddenly intensified upon reaching the eastern Philippines coast on 06:00 14 September and weakened while passing over the Philippines island, before re-intensifying in the SCS. This sharp pressure decrease is explained by the warmer waters east of the Philippines island (~1.5 • C; figure not shown), which provided more energy to the TC. The TC weakened upon passing over the land before re-intensifying in the warm waters of the SCS. Moreover, coupling to ROMS shows an improved simulation (RMSD = 4.10), while coupling to SWAN decreases the simulation performance (RMSD = 4.66). While coupled simulation does not significantly improve the modeling skill for V max , model coupling, especially coupling to ROMS, improves the simulation P min performance.
Having compared the simulation results with best-track data from JMA and JTWC, experimental results are now compared with in-situ (buoy) data (yellow triangles in Figure 1).
Atmosphere 2020, 11, x FOR PEER REVIEW 11 of 28 performance (RMSD = 10.35; R 2 = 0.81). It is observed that model coupling decreases Vmax simulation skills. The more accurate surface roughness estimation in WS might contribute to the slightly smaller MSD. The ROMS-coupled run presents a poorer performance as compared to the SWAN-coupled run, from which it can be inferred that a more accurate representation of the ocean surface roughness in WS improves surface wind simulations. A comparison of Pmin from simulations and observations (Figure 4b) shows that the models better simulate Pmin than Vmax with R 2 above 0.98 and RMSD less than 5.77. However, a TC with a slightly lower Pmin as compared to the observations is obtained from the experiments. From our simulations, Typhoon Kalmaegi suddenly intensified upon reaching the eastern Philippines coast on 06:00 14 September and weakened while passing over the Philippines island, before re-intensifying in the SCS. This sharp pressure decrease is explained by the warmer waters east of the Philippines island (~1.5 °C; figure not shown), which provided more energy to the TC. The TC weakened upon passing over the land before re-intensifying in the warm waters of the SCS. Moreover, coupling to ROMS shows an improved simulation (RMSD = 4.10), while coupling to SWAN decreases the simulation performance (RMSD = 4.66). While coupled simulation does not significantly improve the modeling skill for Vmax, model coupling, especially coupling to ROMS, improves the simulation Pmin performance.
Having compared the simulation results with best-track data from JMA and JTWC, experimental results are now compared with in-situ (buoy) data (yellow triangles in Figure 1). . The same colors are used for RMSD and R 2 , which are calculated relative to JMA data. The dots represent data between 00:00 12 September and 00:00 17 September, with intervals of 6 h.

Sea Level Pressure
The observed and simulated SLP at the four buoy locations ( Figure 5) shows that the simulation generally agrees well with observations, albeit suffering from a slightly lower SLP at all four

Sea Level Pressure
The observed and simulated SLP at the four buoy locations ( Figure 5) shows that the simulation generally agrees well with observations, albeit suffering from a slightly lower SLP at all four locations. A lower SLP can be seen at Buoys 2 and 5, since the TC passes closer to these two buoys. Coupling to ROMS improves the simulation accuracy on the right side of the TC but not on the left side. Mixing from currents, and the additional mixing when the wave component is activated in WRS, result in a slightly less intense TC due to enhanced surface cooling. On the temporal scale, the R 2 of the different experiments are relatively similar.
Atmosphere 2020, 11, x FOR PEER REVIEW 12 of 28 locations. A lower SLP can be seen at Buoys 2 and 5, since the TC passes closer to these two buoys.
Coupling to ROMS improves the simulation accuracy on the right side of the TC but not on the left side. Mixing from currents, and the additional mixing when the wave component is activated in WRS, result in a slightly less intense TC due to enhanced surface cooling. On the temporal scale, the R 2 of the different experiments are relatively similar.

10-m Wind
The spatial distribution of 10-m wind from WRS ( Figure 6) shows faster winds along the TC track. Moreover, a stronger wind magnitude is displayed on the right side of the track.
From the analysis at the buoy locations, the 10-m wind from simulations ( Figure 7) shows good agreement with the observations. The intensification of the TC is evident at Buoys 1 and 2, which displays a faster wind. Moreover, since the TC's center was closer to Buoys 2 and 5, the wind speed at these two locations is slightly weaker than at Buoys 1 and 4. The simulated wind is slightly stronger than observations as the TC passes over the mooring array, except at Buoy 1, where the wind speed peak is almost similar to that of observation. Contrary to Vmax, where coupling to SWAN increases the wind simulation accuracy, this is not the case here. The results also show that the inclusion of the waves in the calculation of surface roughness did not contribute to a more accurate result on the left side of the track but displayed better results on the right side. Interestingly, coupling to ROMS indicates that SST has a degree of influence on the wind speed by modulating the TC intensity. The sudden wind speed increase or decrease from the buoy data, that is not captured in our simulations, might be due to the limitations of our model.
The comparison of wind direction from observation and simulation (Figure 8) shows that all experiments successfully capture the shifts in wind direction. The wind direction is generally clockwise on the right side (Buoys 1 and 4) of the TC track and anticlockwise on the left side (Buoys 2 and 5). The same colors are used for RMSD and R 2 , which are calculated relative to observations.

10-m Wind
The spatial distribution of 10-m wind from WRS ( Figure 6) shows faster winds along the TC track. Moreover, a stronger wind magnitude is displayed on the right side of the track.
From the analysis at the buoy locations, the 10-m wind from simulations ( Figure 7) shows good agreement with the observations. The intensification of the TC is evident at Buoys 1 and 2, which displays a faster wind. Moreover, since the TC's center was closer to Buoys 2 and 5, the wind speed at these two locations is slightly weaker than at Buoys 1 and 4. The simulated wind is slightly stronger than observations as the TC passes over the mooring array, except at Buoy 1, where the wind speed peak is almost similar to that of observation. Contrary to V max , where coupling to SWAN increases the wind simulation accuracy, this is not the case here. The results also show that the inclusion of the waves in the calculation of surface roughness did not contribute to a more accurate result on the left side of the track but displayed better results on the right side. Interestingly, coupling to ROMS indicates that SST has a degree of influence on the wind speed by modulating the TC intensity. The sudden wind speed increase or decrease from the buoy data, that is not captured in our simulations, might be due to the limitations of our model.
The comparison of wind direction from observation and simulation (Figure 8) shows that all experiments successfully capture the shifts in wind direction. The wind direction is generally clockwise on the right side (Buoys 1 and 4) of the TC track and anticlockwise on the left side (Buoys 2 and 5).
Atmosphere 2020, 11, x FOR PEER REVIEW 13 of 28   Atmosphere 2020, 11, x FOR PEER REVIEW 13 of 28   Overall, SLP is better simulated than wind speed, with smaller RMSD and higher R 2 . Having achieved good agreements between the observed and simulated wind and SLP fields, an analysis of ocean parameters can now be performed confidently.
The drag coefficient CD is displayed and compared at the four buoy locations (Figure 9). The differences in CD point to the different feedback of sea surface roughness on the wind. As the TC intensifies, W and WR generally present larger CD, which corresponds to the larger wind speed, as observed in Figure 7. However, the peak CDs for all experiments are relatively similar. A larger CD corresponds to a larger surface roughness, which acts to reduce surface wind while at the same time providing a larger surface area for more heat exchange to the atmosphere.  Overall, SLP is better simulated than wind speed, with smaller RMSD and higher R 2 . Having achieved good agreements between the observed and simulated wind and SLP fields, an analysis of ocean parameters can now be performed confidently.
The drag coefficient C D is displayed and compared at the four buoy locations (Figure 9). The differences in C D point to the different feedback of sea surface roughness on the wind. As the TC intensifies, W and WR generally present larger C D , which corresponds to the larger wind speed, as observed in Figure 7. However, the peak C D s for all experiments are relatively similar. A larger C D corresponds to a larger surface roughness, which acts to reduce surface wind while at the same time providing a larger surface area for more heat exchange to the atmosphere. Overall, SLP is better simulated than wind speed, with smaller RMSD and higher R 2 . Having achieved good agreements between the observed and simulated wind and SLP fields, an analysis of ocean parameters can now be performed confidently.
The drag coefficient CD is displayed and compared at the four buoy locations (Figure 9). The differences in CD point to the different feedback of sea surface roughness on the wind. As the TC intensifies, W and WR generally present larger CD, which corresponds to the larger wind speed, as observed in Figure 7. However, the peak CDs for all experiments are relatively similar. A larger CD corresponds to a larger surface roughness, which acts to reduce surface wind while at the same time providing a larger surface area for more heat exchange to the atmosphere.

Oceanic Parameters
Here, the performance of the models in simulating ocean parameters, namely, ocean surface, SST current and SWH, is evaluated. We assume that changes in these parameters are mainly induced by atmospheric forcing.

Ocean Current
As Typhoon Kalmaegi moves over the ocean, wind stress caused by strong surface winds results in the adjustment of ocean currents to flow typically in the direction of the wind forcing. Following the TC's passing, ambient ocean surface currents slowly recover to their initial states. The spatial distribution of surface currents ( Figure 10) shows a corresponding increased current speed as Typhoon Kalmaegi passes over a region. As the TC intensifies, stronger currents are observed on the right side of the TC track (15 and 16 September), which agrees with past studies (e.g., [3]). The near-surface currents are mostly dependent on the wind stress, with maximum current speed to the right of the TC [4] caused by the resonance effect of the wind stress on the right side of the track. On 12 September, the region of strongest current appears at the rear of the TC. Due to the proximity of the TC to the Kuroshio at that time, it is possible that the surface current is influenced by the strong Kuroshio current. It is also noted that the strong currents east of the Philippines are still present on 15 September but have nearly recovered on 16 September. After the passage of the TC, the current on the right side of the track rotates clockwise, because the wind stress changes to a clockwise rotation.
Atmosphere 2020, 11, x FOR PEER REVIEW 16 of 28  From Figures 11 and 12, the simulated current direction matches that from observation, albeit experiencing a higher current speed, especially on the right side of the TC (Buoys 1 and  4). R generally presents the fastest currents during the TC event while WR has the slowest current. The difference between observed and simulated currents is possibly caused by the limitations of the model in high-wind conditions. The surface current direction follows the Ekman flow, which is a shift of 45 • clockwise (in the northern hemisphere) from the direction of wind flow (Figure 8). Moreover, the current direction is observed to rotate clockwise at all three buoy locations.

Sea Surface Temperature
A comparison of the SST difference from observations (MW-IR) and WRS simulation, calculated by subtracting the daily mean of 12 September from the daily mean of 16 September, shows that the WRS experiment reasonably reproduces the observed ocean surface cooling along the TC track (Figures 13a and b). Stronger surface cooling is observed on the right side of the TC track in both cases, due to the faster wind speed on the right-hand side of the TC track. The maximum cooling in

Sea Surface Temperature
A comparison of the SST difference from observations (MW-IR) and WRS simulation, calculated by subtracting the daily mean of 12 September from the daily mean of 16 September, shows that the WRS experiment reasonably reproduces the observed ocean surface cooling along the TC track (Figure 13a,b). Stronger surface cooling is observed on the right side of the TC track in both cases, due to the faster wind speed on the right-hand side of the TC track. The maximum cooling in MW-IR is about 4 • C and extends toward the continental slope. In the model, by contrast, it is about 3.5 • C and is mainly confined to the continental shelf. Figure 13c,d are the differences between WRS and WR, and R, respectively. The SST cooling in WR is slightly weaker than that in WRS. However, for R, the cooling along the TC track is much pronounced, which can be attributed to the stronger surface wind, thus resulting in an enhanced mixing.
WRS is expected to have an improved SST simulation, as mixing due to currents and waves are both included in a fully coupled model. Comparing the simulated SST with mooring data (Figure 14), R 2 is greater than 0.92, while RMSD is less than 0.63, which points to good model performance. As Typhoon Kalmaegi passes over the moorings, a drop of~2-3 • C is observed on the right side of the TC (Buoys 1 and 4) and a drop of~1.5 • C is noted at Buoys 2 and 5, which are on the left side of the TC track. However, the models do not capture the cooling on the TC right-hand side between 00:00 15 September and 00:00 16 September. Due to the thick mixed layer, the deep cooler water is not brought to the surface. The SST cooling on the left-hand side is more in line with observations. Additionally, Buoys 1 and 4 show a slight SST recovery after the passage of Typhoon Kalmaegi, while at Buoys 2 and 5, the SST continues to decrease steadily. According to the literature, SST recovers within a few days to weeks [26,80,81]. It is worth noting that there is a difference in SST at the beginning of the run, which results from the initial conditions applied to ROMS. That is, the comparison of water temperature from the initial condition and buoy data shows a small difference between the two data.
Atmosphere 2020, 11, x FOR PEER REVIEW 18 of 28    The comparison of MLD from the three experiments ( Figure 15) shows a larger departure of R from WRS as compared to WR from WRS. Moreover, WRS generally presents a deeper MLD than WR, indicating a possible added turbulent mixing from the waves. The comparison of MLD from the three experiments ( Figure 15) shows a larger departure of R from WRS as compared to WR from WRS. Moreover, WRS generally presents a deeper MLD than WR, indicating a possible added turbulent mixing from the waves.

Significant Wave Height
SWH includes a complex interaction between waves and currents, modulated by surface wind forcing. Higher waves are expected in areas of higher wind speed. In the simulation of wave conditions during Typhoon Kalmaegi, the spatial distribution of SWH ( Figure 16) shows the largest SWH distributions along the TC track, but more specifically, is higher on its right-hand side. The wind direction on the right-hand side of the TC is nearly in the same direction as the TC propagation direction, while on the left side, the wind direction and TC propagation direction are in opposite directions. Therefore, due to a resonance effect, the waves on the right-hand side are higher as a result

Significant Wave Height
SWH includes a complex interaction between waves and currents, modulated by surface wind forcing. Higher waves are expected in areas of higher wind speed. In the simulation of wave conditions during Typhoon Kalmaegi, the spatial distribution of SWH ( Figure 16) shows the largest SWH distributions along the TC track, but more specifically, is higher on its right-hand side. The wind direction on the right-hand side of the TC is nearly in the same direction as the TC propagation direction, while on the left side, the wind direction and TC propagation direction are in opposite directions. Therefore, due to a resonance effect, the waves on the right-hand side are higher as a result of longer exposure to the TC's strong winds as the system moves forward. By contrast, the waves on the left side are influenced by the TC winds for a shorter period. Figure 17 compares SWH from observation and simulations at Buoys 1 and 4, which are both on the right side of the TC track. Wind speed at Buoy 1 is faster, thus matching in a slightly higher SWH at that location. By comparing all the experiments, it is observed that the 3-way coupling improves the SWH simulation. The lower RMSD of WS as compared to S indicates that atmosphere-wave coupling indeed influences the wave simulation, while the best performance of WRS in simulating SWH, demonstrates the importance of wave-current interactions in the correct simulation of SWH. The comparison of observed and simulated wave directions ( Figure 18) shows only minor differences between the two as the TC approaches the moorings. During and after the passage of Typhoon Kalmaegi, the simulated wave direction is more consistent with observations.

Skill Score
This study also aims to identify the best model coupling combination for Typhoon Kalmaegi simulation. In this section, the model experiments are evaluated using the skill score (SS) method introduced in Section 2.4. The R 2 , RMSD, SS c and SS r values are averages calculated from the individual criteria at the different observation points.
The SS for wind and SLP is presented in Table 3. In most cases, the baseline experiment shows the most accurate results, while WS displays the worst results. As we observed in the analyses above, to have a better representation of the SST and thus an improved TC intensity simulation, coupling to the ocean is crucial [12][13][14]. Table 3. Correlation coefficient (R 2 ), root-mean-square difference (RMSD) and skill score (SS c and SS r ) for 10-m wind and SLP from the different experiments.

Experiments
Wind SLP The SS for SST and surface current quantify the performance of the models in simulating these variables (Table 4). For SST, the uncoupled ROMS performs best as compared to the baseline experiment, and WR performs worst. The reason why R performs best in terms of RMSD might be due to the wind forcing for the ROMS uncoupled model. The stronger unabated wind might have caused heightened mixing, thus bringing the results close to the observations as compared to WRS. However, wave-induced mixing is not included in WR, hence the slightly weaker cooling. Analyzing ocean surface current, R performs best for R 2 . Although coupling to surface waves in the WRS experiment allows for a better estimation of the surface roughness, the results in terms of R 2 and RMSD are slightly poorer. Accurately simulating the oceanic conditions still proves to be a challenge. Among the three sensitivity tests involving SWAN (Table 5), the baseline experiment presents the best results. The WS and S experiments exhibit the lowest performance for R 2 and RMSD, respectively. The results suggest that wave simulation requires coupling not only to the atmosphere but also to the ocean. In both S and WS, SWAN is forced by a wind field that might not correspondingly weaken due to surface cooling, since it is not coupled to ROMS. Therefore, the stronger winds generate higher waves and contribute to the observed results.
The results from this section show the performance of the different experiments in simulating Typhoon Kalmaegi, relative to the WRS baseline experiment. WRS generally presents the best simulation performance for wind, SLP and SWH but not for SST and ocean current. However, considering that air-sea exchanges and wave-current interactions are crucial for a consistent physical and dynamic TC system, WRS is the best coupling combination.

Discussion and Conclusions
Due to the difficulty and high expense of obtaining in-situ data in the middle of the ocean during TC events, numerical simulations are widely applied in TC air-sea interactions studies. Previous studies significantly contributed to the knowledge of the processes occurring during a TC event. The current study complements past research by investigating the effects of model coupling on the simulation of Typhoon Kalmaegi, determining the model coupling combination that produces the most accurate results and contributing to discussions about the physical and dynamic processes during a TC. An extensive evaluation of the model experiment results against observations is first conducted to ensure consistent results for the atmospheric and oceanic variables. A sophisticated statistical method, the skill score, is then applied to compare the simulation performance of the different model combinations against a baseline experiment.
We first start by running the uncoupled atmosphere, ocean and wave models separately, followed by WRF-ROMS and WRF-SWAN coupling to assess the possible contribution of these two components to the simulation accuracy. Finally, a WRF-ROMS-SWAN coupling is run to investigate the role of both ocean and wave coupling in TC simulation. The effects of parameterization schemes, model configurations and initial conditions are not investigated in this study.
As presented in the previous sections, the numerical experiments all show favorable comparisons of the atmosphere and ocean with in-situ data during Typhoon Kalmaegi. Since the differences between the tracks' geometric distances are not very large (Figures 2 and 3a), it is suggested that model coupling does not significantly influence track simulation. Instead, the TC track is mainly dependent on large-scale circulation [6]. Translation speed, on the other hand, is moderated by surface roughness; a larger surface roughness leads to more friction, which ultimately slows down the TC. As can be seen in our results (Figure 3b), WS has a better representation of surface waves, thus a better estimation of surface roughness, which produces an improved translation speed simulation. The poorer translation speed simulation of WRS as compared to WS is attributed to the uncertainties brought by wave-current interactions into the surface roughness estimation. Figure 4 evidences that V max is more significantly influenced by surface roughness while SST more strongly modulates P min . Moreover, it is observed that the deepening of the TC produces a lower P min , but V max does not intensify accordingly. The difference between pressure and wind is explained as follows: a larger surface roughness acts to increase the surface area of the ocean where more heat can be exchanged, thus strengthening the TC (lower P min ). On the other hand, a rougher sea surface also means more friction, which works to slow down the surface wind and might even influence the TC track. Therefore, increased heat fluxes intensify the TC, while increased friction weakens the TC, thus creating a balance so that the TC does not over-intensify or over-weaken.
Sea surface cooling, mainly from entrainment [82], has been demonstrated to modulate TC intensity [83] significantly. A clear feedback mechanism exists between atmosphere and ocean: as the TC approaches and strengthens, upper ocean cooling occurs. The TC thus weakens due to less heat and moisture availability, and upper ocean cooling decreases. WR and WRS experiments simulate the SST cooling reasonably well, while R displayed a weaker cooling ( Figure 13). Since SST corresponds to the amount of energy transferred to the TC, an accurate SST field is essential to ensure a TC with the correct core pressure [46]. The W and WS experiments simulate a slightly more intense TC due to the SST field