Waves in SKRIPS: WAVEWATCH III coupling implementation and a case study of Tropical Cyclone Mekunu
In this work, we integrated the WAVEWATCH III model into the regional coupled model SKRIPS (Scripps–KAUST Regional Integrated Prediction System). The WAVEWATCH III model is implemented with flexibility, meaning the coupled system can run with or without the wave component. In our implementations, we considered the effect of Stokes drift, Langmuir turbulence, sea surface roughness, and wave-induced momentum fluxes. To demonstrate the impact of coupling we performed a case study using a series of coupled and uncoupled simulations of Tropical Cyclone Mekunu, which occurred in the Arabian Sea in May 2018. We examined the model skill in these simulations and further investigated the impact of Langmuir turbulence in the coupled system. Because of the chaotic nature of the atmosphere, we ran an ensemble of 20 members for each coupled and uncoupled experiment. We found that the characteristics of the tropical cyclone are not significantly different due to the effect of surface waves when using different parameterizations, but the coupled models better capture the minimum pressure and maximum wind speed compared with the benchmark stand-alone Weather Research and Forecasting (WRF) model. Moreover, in the region of the cold wake, when Langmuir turbulence is considered in the coupled system, the sea surface temperature is about 0.5 ∘C colder, and the mixed layer is about 20 m deeper. This indicates the ocean model is sensitive to the parameterization of Langmuir turbulence in the coupled simulations.
Ocean surface waves play a key role in mediating exchanges of momentum, heat, and gases across the air–sea boundary (Fan et al., 2009; D'Asaro et al., 2014). The importance of surface waves in mediating air–sea interactions has been studied for decades (Fairall et al., 2003; Chen et al., 2007b). Surface waves can enhance upper-ocean mixing through Langmuir turbulence, and neglecting the Langmuir mixing process may contribute to a shallow bias in mixed-layer depth (MLD) (Li et al., 2016). In addition, waves determine the sea surface roughness, which affects wind stress that is important for short-term forecasting of tropical and sub-tropical cyclones (Olabarrieta et al., 2012).
Several regional coupled model studies have considered the effect of waves in air–sea interactions (e.g., Chen et al., 2007a; Warner et al., 2010; Liu et al., 2011; Wu et al., 2019a; Lewis et al., 2019; Sauvage et al., 2022). Because of the importance of air–sea heat fluxes on the energy budget of a tropical cyclone (TC) (Emanuel, 1991), many of these studies have focused on TCs and demonstrated increased accuracy in simulated intensity of TCs when coupled (e.g., Bender and Ginis, 2000; Chen et al., 2007b; Warner et al., 2010; Wu et al., 2019a; Lewis et al., 2019; Li et al., 2022). Studies have also shown a strong coupled feedback in conditions where the heat content of the upper-ocean layer is low and a weak feedback when the ocean has a thick mixed layer (Mogensen et al., 2017). Saxby et al. (2021) highlight outstanding challenges, with high-resolution convection-permitting atmosphere-only and coupled configurations both accurately simulating TCs in the Bay of Bengal, suggesting that many of the deficiencies originate in the atmospheric model, but improvements could also be gained by coupling this to a wave model.
The sea state is highly complex and variable in TC conditions, with Langmuir turbulence playing an important role in the upper-ocean mixing (Rabe et al., 2015; Reichl et al., 2016a, b). This turbulence is associated with coherent Langmuir circulation structures that exist and evolve over a range of spatial and temporal scales in the surface ocean (Langmuir, 1938; McWilliams et al., 1997; Thorpe, 2004). These structures arise through an interaction between ocean surface waves and the background Eulerian current. Langmuir turbulence enhances turbulent entrainment, deepening the mixed layer and leading to sea surface cooling, which in turn affects the air–sea heat fluxes that modulate the development of TCs. Studies of idealized TCs suggest including Langmuir turbulence in model simulations may cool the sea surface temperature (SST) by 0.5–0.7 ∘C and increase the mixed-layer depth by up to 20 m (Reichl et al., 2016b; Blair et al., 2017).
Because of the importance of ocean surface waves in air–sea interaction, we implemented a regional coupled ocean–wave–atmosphere model, with the capability of investigating the impact of surface waves on air–sea interaction. The goal of this work is twofold. First, we demonstrate the integration of the wave model WAVEWATCH III to the Scripps–KAUST Regional Integrated Prediction System (SKRIPS, Sun et al., 2019), which is a regional coupled ocean–atmosphere model that has been used to investigate extreme heat wave events on the shore of the Red Sea (Sun et al., 2019), North Pacific atmospheric rivers (Sun et al., 2021), and sea ice evolution in the Southern Ocean (Cerovečki et al., 2022). The second goal is to evaluate the implementations of ocean surface waves in the coupled system, especially for Langmuir turbulence that alleviates the model bias (Li et al., 2016; Li and Fox-Kemper, 2017). The coupled model is also sensitive to the parameterization of Langmuir turbulence because it increases ocean mixing and cools down the SST during the simulation. Here, we perform a series of coupled and uncoupled numerical simulations of Tropical Cyclone Mekunu in the Arabian Sea. Because of the chaotic nature of the atmosphere, we ran an ensemble of 20 members for each coupled and uncoupled experiment. The Arabian Sea is investigated in this work because of its rich and diverse ecosystem, its economic impact on the surrounding countries, and its important role in international trade. Continued climate warming is expected to further amplify the risk of cyclones in the Arabian Sea (Dube et al., 1997; Evan et al., 2011; Evan and Camargo, 2011) and increase socio-economic implications for coastal communities in that region (Henderson-Sellers et al., 1998; Murakami et al., 2017; Bhatia et al., 2018). We investigated Tropical Cyclone Mekunu because it was the strongest tropical cyclone in the northern Indian Ocean in 2018 (Government of India, 2018). It had a clear signature of SST cooling and MLD deepening that can be used for testing the model. We investigated the sensitivities of the coupled model to the parameterizations of surface-wave-driven mixing to examine the effect of surface waves on air–sea interactions.
The rest of this paper is organized as follows. The implementation of the coupled model is described in Sect. 2. An overview of Tropical Cyclone Mekunu, the design of the experiments, and the validation data are presented in Sect. 3. Section 4 details the numerical simulation results, and Sect. 5 discusses the sensitivity of the simulation results to parameterizing Langmuir turbulence based on the evolving wave state. Section 6 concludes the paper with a summary of the main findings.
2.1 Model description
In this work, version 6.07 (Tolman, 1991; WW3DG, 2019) of the Wave Height, Water Depth And Current Hindcasting Third-Generation Wave Model (WAVEWATCH III, hereinafter, WW3) is integrated into the SKRIPS model. The SKRIPS model (Sun et al., 2019) is a regional coupled ocean–atmosphere model: the oceanic model component is the MIT general circulation model (MITgcm) (Marshall et al., 1997; Campin et al., 2019), and the atmospheric model component is the Weather Research and Forecasting (WRF) model (Skamarock et al., 2019). The Earth System Modeling Framework (ESMF) (Hill et al., 2004) is used as the coupler to drive the coupled simulation. The National United Operational Prediction Capability (NUOPC) layer in the ESMF is used to simplify the implementations of component synchronization, execution, and other common tasks in the coupling (Hill et al., 2004).
The schematic description of the coupled model is shown in Fig. 1. In the coupling process, all model components send data to ESMF: MITgcm sends SST and ocean surface velocity; WRF sends surface atmosphere fields, including (1) net surface longwave and shortwave radiative fluxes, (2) surface latent and sensible heat fluxes, (3) 10 m wind speed, (4) precipitation, and (5) evaporation; and WW3 sends the wave variables to ESMF, including the (1) bulk wave parameters (i.e., significant wave height, peak wavelength, and mean wavenumber), (2) surface Stokes drift, (3) Langmuir turbulence parameters (i.e., Langmuir number and enhancement factor), and (4) momentum flux terms due to surface waves. Then all model components read the data they need from ESMF: MITgcm reads surface atmospheric variables and wave variables; WRF reads SST, ocean surface velocity, and wave variables; and WW3 reads wind speed and surface current velocity. The surface current velocity sent to WRF and WW3 is consistent in our model, using the current velocity in the first layer of MITgcm. We used the surface current based on previous literature (Warner et al., 2008, 2010; Couvelard et al., 2020), but this may overestimate the strength of surface currents impacting the wave model, as suggested by Fan et al. (2009), who used the current velocity at (L is the mean wavelength). The wind speed sent to WW3 and MITgcm is the relative 10 m wind speed from WRF based on the Monin–Obukhov similarity theory (Monin and Obukhov, 1954; Renault et al., 2020), and WW3 and MITgcm use the relative 10 m wind speed without correcting the current velocity in the simulations.
Similar to our previous work Sun et al. (2019), the MITgcm model uses the surface atmospheric variables received from ESMF to prescribe surface forcing, including (1) total net surface heat flux, (2) surface wind stress, and (3) freshwater flux. The total net surface heat flux is computed by adding surface latent heat flux, sensible heat flux, net shortwave radiation flux, and net longwave radiation flux. The surface latent and sensible heat fluxes are computed using the COARE 3.0 bulk algorithm in WRF (Fairall et al., 2003). The implementations of the wave effects are discussed in the following sections: the Stokes forces, the Langmuir turbulence parameters, and the momentum fluxes are detailed in Sect. 2.2, 2.3, and 2.4, respectively; the sea surface roughness parameterizations are summarized in Sect. 2.5. The Stokes forces, the Langmuir turbulence parameters, and the momentum fluxes are detailed in Sect. 2.2, 2.3, and 2.4, respectively. The sea surface roughness parameterizations are summarized in Sect. 2.5.
To implement the WW3–ESMF interface, we followed our previous implementations (Sun et al., 2019) and the WW3–ESMF interface (WW3DG, 2019) for the prototype case. We separated the WW3 main program into three subroutines that handle initialization, execution, and finalization. These subroutines are used by the ESMF–NUOPC coupler that controls the wave component in the coupled run. During the simulation, WW3 receives and sends boundary fields via subroutine calls by the WW3–ESMF interface, shown in Fig. 1. In addition, WW3 grid information is provided to the coupler in the initialization subroutine. To carry out the coupled simulation on HPC (high-performance computing) clusters, the WW3–ESMF interface runs in parallel via MPI (message passing interface) communications. We have also updated the MITgcm–ESMF and WRF–ESMF interfaces by including the inputs and outputs associated with the wave model. The wave component is implemented with flexibility, meaning the coupled system can run with or without the wave component. It is noted that ESMF online regridding options are also implemented for the wave component when exchanging boundary fields, but these are not used in this work because we aim to present the implementation of wave components. The online regridding option will be used when using a higher-resolution ocean model for the Arabian Sea operational model.
2.2 Stokes forces in MITgcm
where t is time; ; ρw is the density of water; p is the pressure; is the buoyancy term; is the vertical unit vector; Du is the diffusion; is the Coriolis parameter; is the Stokes shear force; is the wave-filtered Eulerian velocity (Eulerian velocity of the flow solved in MITgcm); uS is the Stokes drift; and is the wave-filtered Lagrangian velocity. Here, the Einstein summation convention is used (e.g., ), although vector notation is used when it is unambiguous.
where c is a scalar quantity, such as potential temperature and salinity, and Dc is the diffusion.
The Stokes advection and Stokes Coriolis terms are implemented in MITgcm by modifying the source term of the governing equations. The profiles of Stokes velocity are determined based on Breivik et al. (2014). Considering the effect of Langmuir turbulence, the Stokes shear term in Eq. (1) is parameterized according to the literature (Li et al., 2016; Li and Fox-Kemper, 2017; Li et al., 2017), detailed in Sect. 2.3. It is noted that our implementations and tests aim to demonstrate the impact of Langmuir turbulence on the ocean, and thus the divergence of the Stokes drift is not considered in our governing equations as discussed in Wu et al. (2019a, b). There are also other options to better approximate the Stokes velocity profiles (e.g., Breivik et al., 2016; Romero et al., 2021) that remain to be tested in future work.
2.3 Parameterization of Langmuir turbulence
Considering the impact of the surface waves, the Stokes drift provides a source of the turbulent kinetic energy (TKE) through the vortex force and modified pressure (Craik and Leibovich, 1976) or, more cleanly, the Stokes shear force (Suzuki and Fox-Kemper, 2016; Li and Fox-Kemper, 2017) as mentioned in Eq. (1). Evidence of this enhanced vertical mixing has been documented from observations and large-eddy simulations (McWilliams et al., 1997; D'Asaro, 2001; Van Roekel et al., 2012; D'Asaro et al., 2014). In this work, we aim to implement the Stokes shear force in the coupled model and investigate its effect on the coupled system. Although it is not explicitly accounted for in KPP (K‐profile parameterization) (Large et al., 1994), KPP might have implicitly incorporated some effects by tuning the parameters to ocean observations (Reichl et al., 2016b). In addition, our implemented model is about 8 km resolution (0.075∘), and the horizontal gradients of the Stokes drift are several orders of magnitude smaller than vertical gradients. Following Suzuki and Fox-Kemper (2016), we only consider the effects of Stokes shear force due to Langmuir turbulence because of this scale separation.
Although there are many unknowns about the role of Langmuir mixing in ocean modeling, there are many parameterizations that aim to represent these processes and alleviate model bias. Within the KPP scheme, we implemented three Langmuir turbulence parameterizations (Van Roekel et al., 2012; Li et al., 2016; Li and Fox-Kemper, 2017; Li et al., 2017): (1) VR12-MA, (2) LF17, and (3) LF17-ST. Both VR12-MA and LF17 parameterize the Langmuir turbulence based on the parameters computed from WW3: in VR12-MA the KPP turbulent velocity scale is multiplied by an enhancement factor, while in LF17 the KPP turbulent velocity scale is treated in the same way as VR12-MA, and the entrainment buoyancy flux is also considered. On the other hand, LF17-ST parameterizes the Langmuir turbulence similarly to LF17, but parameters are computed using the 10 m winds instead of using the output from WW3. VR12-MA and LF17 are implemented because they are used in a variety of case studies and substantially improve the shallow biases of mixed-layer depth (Li et al., 2016; Li and Fox-Kemper, 2017; Li et al., 2019). We aim to compare the performance of VR12-MA and LF17 to demonstrate the impact of entrainment on the simulations. We also used the well-validated LF17-ST implementation by Schultz et al. (2020) to validate LF17 in the coupled simulations due to the similarity of these two options. Because LF17-ST does not need bulk wave parameters as input, it can be also used in uncoupled MITgcm simulations (Schultz et al., 2020) or coupled simulations without waves to parameterize the Langmuir turbulence.
For all the implemented schemes, the turbulent velocity scale ws is modified by multiplying with an enhancement factor ε as follows:
where α is the angle between wind and Langmuir cells and LaSLP is the Langmuir number (Van Roekel et al., 2012; Li et al., 2016, 2019). Here we used the projected Langmuir number defined in Eq. (6) of Li et al. (2019) based on the surface averaged Stokes velocity. More details of the projected Langmuir number can be found in Appendix A.
The enhanced turbulent velocity scale affects the vertical viscosity, tracer diffusivity, and nonlocal flux in KPP. In particular, the KPP eddy diffusivity profile κv is
where h is the boundary layer depth, h is the KPP boundary layer depth, G(σ) is the shape function, and is the normalized depth. The boundary layer depth is determined based on the bulk Richardson number:
where Br and ur are the buoyancy and velocity averaged in the surface layer, B(z) and u(z) are the local buoyancy and local velocity, and Vt is the unresolved vertical shear proportional to ws (Large et al., 1994; Li et al., 2016). The KPP boundary layer depth is defined at the smallest KPP boundary layer depth h that reaches the critical bulk Richardson number Ricr=0.3.
Different from VR12-MA, LF17 parameterized the entrainment flux due to Langmuir turbulence by revising the bulk Richardson number:
the turbulence shear term is defined as follows:
where dimensionless coefficient , N(z) is the local buoyancy frequency, Ric is the critical Richardson number, is the convective velocity scale, and B0 is the surface buoyancy flux. Here LaSL is the Langmuir number defined in Eq. (5) of Li et al. (2019) and is detailed in Appendix A.
In LF17-ST, the enhancement coefficient and entrainment flux are calculated similarly to LF17, but the Langmuir turbulence coefficient La is determined by the Stokes velocity parameterized from 10 m wind and mixed-layer depth. The details can be found in Eq. (25) in Li et al. (2017). In this work, the implementations of LF17-ST in MITgcm followed the code provided by Schultz et al. (2020).
2.4 Momentum flux in MITgcm
The surface boundary condition is also modified by waves. The surface wind stress τa is modified by subtracting the part that goes into wave growth τaw and adding the wave-to-ocean momentum flux due to wave breaking τow. Hence, the momentum flux in MITgcm is (Jenkins, 1989; Weber et al., 2006; Janssen, 2012)
Here, Sin and Sds are the wind input and wave dissipation source terms, respectively; k is the wave number; ω is the angular wave frequency; and θ is the wave direction. In the coupled model, τa is calculated in MITgcm (Large and Yeager, 2004) because WRF does not directly output the momentum flux terms. The parts that go into wave growth τaw and wave breaking τow are calculated in WW3 (Tolman, 1995; Ardhuin et al., 2003).
2.5 Ocean roughness closures
The parameterization of ocean roughness in WRF is also important. When WRF is not coupled to WW3, the bottom roughness length z0 is computed with the formulation proposed by Smith (1988), which is a combination of the formulae described by Liu et al. (1979) and Charnock (1955). When coupled with WW3, we parameterize the surface roughness based on the Charnock coefficient calculated from WW3 to make the surface roughness consistent. We have also implemented a few other ocean roughness closure models that have been used in COAWST Olabarrieta et al. (2012): DGHQ (Drennan et al., 2003, which is based on wave age), TY2001 (Taylor and Yelland, 2001, which is based on wave steepness), and OOST (Oost et al., 2002, which considers both the effects of wave age and steepness). These models parameterize z0 using the bulk wave parameters from WW3. We implemented these options in the WRF Mellor–Yamada–Nakanishi–Niino (MYNN) surface layer scheme. More detailed descriptions of these closure models and sensitivity analysis are presented in Appendices B and C.
3.1 Overview of the event
Tropical Cyclone Mekunu formed in the southeastern Arabian Sea on 20 May 2018 and then propagated northwest before making landfall in southwestern Oman on 26 May. Categorized as ESCS (extremely severe cyclonic storm), Tropical Cyclone Mekunu was the second cyclonic storm over the Arabian Sea in 2018 and the strongest tropical cyclone in the northern Indian Ocean that year. The peak maximum sustained surface wind speed was 170–180 km h−1, gusting to 200 km h−1 (95 knots), and the lowest estimated central pressure was 960 hPa on 25 May (Government of India, 2018). Salalah, the capital city of southern Oman's Dhofar province, received 278.2 mm (10.95 inches) of rain in just 24 h ending around 10:30 LT (local time) on 26 May, over double the city's average annual rainfall of about 127 mm (5 in.), with a total of 617 mm of rainfall during 23–27 May (Government of India, 2018). In addition to the extremely heavy rainfall in Oman and Yemen, Tropical Cyclone Mekunu caused heavy rainfall that created desert lakes over the “Empty Quarter” in Saudi Arabia. The warm, sandy, and wet soil was the perfect environment for the outbreak of desert locusts, posing a serious risk to food security and livelihoods (Salih et al., 2020).
3.2 Model setup
To illustrate the coupled model capabilities, we perform the following types of model runs.
CPL.AOW. Coupled ocean–wave–atmosphere (MITgcm–WW3–WRF) simulations.
CPL.AO. Coupled ocean–atmosphere (MITgcm–WRF) simulations, where the ocean–atmosphere model is not coupled to the wave model, aiming to demonstrate the impact of the wave model on the simulation results.
ATM.DYN. Stand-alone atmosphere (WRF) simulations, where the atmosphere model is not coupled to the wave or ocean models. The SST forcing is from the HYCOM/NCODA daily global analysis data (the Global Ocean Forecast System, Version 3.1 (Chassignet et al., 2007), hereinafter, HYCOM). Compared with CPL.AO and CPL.AOW, this run serves as a benchmark that aims to demonstrate the impact of waves and coupled air–sea interactions on the simulation results.
The grid spacing and computational domain are outlined in Table 1. To generate the grids, we choose the latitude–longitude (cylindrical equidistant) map projection for MITgcm, WW3, and WRF. The model domain extends from 0 to 30.6∘ N and from 30 to 78∘ E. The horizontal grid has 408 × 640 (lat × long) cells, and the spacing is 0.08∘ in both directions. We use identical horizontal grids for all model components to eliminate the complication of regridding winds near steep orography and complex coastlines, although the regridding capability is implemented in SKRIPS. There are 40σ layers in the atmosphere model (top pressure is 50 hPa) and 50 z layers in the ocean model (dz=4 m in the top). The wave model has a spectral grid of 48 directions (7.5∘ resolution) and 32 frequencies exponentially spaced from 0.0343 to 1.1 Hz. Because of the chaotic nature of the atmosphere, we generated 20-member ensembles for each run by adding small random perturbations to the initial SST (<0.01 ∘C) at every grid point in the coupled model. The random perturbations are added without any spatial or temporal correlation, aiming to demonstrate the internal variability of the model.
The initial conditions, boundary conditions, and forcing terms of the simulations are also outlined in Table 1. In the coupled runs, the ocean model uses the HYCOM data as initial and boundary conditions for ocean temperature, salinity, and horizontal velocities (Chassignet et al., 2007). At each time step, the boundary conditions for the ocean are updated by linearly interpolating between the daily HYCOM data. A restoring layer with a width of 13 grid cells is applied at the lateral boundaries to enforce the boundary conditions. The inner and outer boundary relaxation timescales are 10 and 0.5 d, respectively. The atmosphere is initialized using the NCEP (National Centers for Environmental Prediction) FNL (Final) Operational Global Analysis data. The same data also provide the boundary conditions for air temperature, wind speed, and humidity. The atmospheric boundary conditions are updated based on linearly interpolating between 6-hourly NCEP FNL data. The “specified” zone in WRF prescribes the lateral boundary values, and the “relaxation” zone is used to nudge the solution from the domain interior toward the boundary condition value. Here we use the default width of one point for the specified zone and four points for the relaxation zone. In the wave model, the wave spectra at the offshore boundary come from the global wave modeling system described by Rascle and Ardhuin (2013). To initialize the wave model, we allowed the wave field to spin up for 19 d from 1 May 2018, and we then analyze the period from 20 May 2018. In contrast, we did not spin up MITgcm or WRF, trying to initialize the coupled model using the analysis data directly. This may cause an initial shock in the coupled simulation, but large initial shocks are not observed in the simulations. We performed downscaled hindcasts in this work, which allows us to focus on the impacts of air–sea interactions during the tropical cyclone event by minimizing the boundary errors.
The time step of the ocean model is 120 s. The horizontal sub-grid mixing is parameterized using nonlinear Smagorinsky viscosities, and the K-profile parameterization is used for vertical mixing processes (Large et al., 1994) with modifications accounting for Langmuir mixing (Van Roekel et al., 2012; Li et al., 2016; Li and Fox-Kemper, 2017; Li et al., 2017). The time step of the atmosphere model is 30 s. The Morrison two-moment scheme (Morrison et al., 2009) is used to resolve the microphysics, the updated version of the Kain–Fritsch convection scheme (Kain, 2004) is used for cumulus parameterization, the Mellor–Yamada–Nakanishi–Niino 2.5-order closure scheme (Nakanishi and Niino, 2004, 2009) is used for the planetary boundary layer (PBL) and the surface layer (SL), the Rapid Radiation Transfer Model for General Circulation Models (RRTMG; Iacono et al., 2008) is used for longwave and shortwave radiation transfer through the atmosphere, and the Noah land surface model is used for the land surface processes (Tewari et al., 2004). The selection of WRF physics schemes is the same as our previous work (Sun et al., 2021). The wave model uses a global integration time step of 600 s, spatial advection time step of 60 s, spectral advection time step of 60 s, and minimum source term time step of 10 s. In CPL.AOW and CPL.AO, the coupling interval is 120 s to allow for capturing the diurnal cycle of air–sea fluxes. We output the simulation results every 3 h to demonstrate the tropical cyclone evolution.
When the effects of surface waves are considered in CPL.AOW, the model setup is as follows. The Stokes Coriolis and the Stokes advection in Eq. (1) are considered, the impact of Langmuir turbulence is parameterized in the same way as Li and Fox-Kemper (2017), the ocean surface roughness is parameterized using the Charnock coefficient (CHNK) from WW3, and the wind stress in the ocean model is treated as mentioned in Eq. (9). We have compared the coupled model with and without wave effects in Sect. 4, and we further illustrate the sensitivity of the coupled model to the wave effects in Sect. 5 and Appendix C. It is noted that when the Stokes Coriolis and the Stokes advection are not explicitly considered in the experiments, the model setups are consistent with Li et al. (2016), assuming the simulated velocity is Lagrangian.
3.3 Validation data
To evaluate the performance of the simulations, the model outputs are compared with available data. The track of the tropical cyclone, the tropical central pressure, and the maximum wind speed are validated against IBTrACs data (Knapp et al., 2010, 2018). Here we use the IBTrACS–World Meteorological Organization version. IBTrACS provides a compilation of historical TC data as recorded by meteorological centers and/or forecast agencies, and in this case the data from Indian Meteorological Department (IMD) are used.
The simulated SST fields are validated against the HYCOM data. Because the HYCOM data are the initial and boundary conditions in the coupled model, this aims to show the error increase from the initial condition. We used bilinear interpolation to map the validation data onto the model grid to compare the results in a consistent way. When interpolating SST, only the values on ocean points are used. The SST is also validated by using in situ observations from the satellite-tracked drifters of the Global Drifter Program (GDP, from https://www.aoml.noaa.gov/envids/gld/index.php, last access: 9 June 2023) (Lumpkin and Pazos, 2007).
In this section, the ensemble coupled simulation results (i.e., CPL.AOW and CPL.AO) are compared with the results from the uncoupled runs (i.e., ATM.DYN) to assess the performance of the models, the impact of coupled feedbacks, and the effect of the waves. We compared the ensemble-averaged characteristics of the tropical cyclone (e.g., track, intensity, and wind speed), the changes in the ocean (e.g., sea surface temperature and mixed-layer depth), and the waves generated by the tropical cyclone. By comparing the coupled run with uncoupled runs, we aim to (1) demonstrate the capability of the coupled model and (2) illustrate the impact of including ocean–wave–atmosphere interactions on simulating this tropical cyclone event.
4.1 Cyclone track, intensity, and wind speed
First, we examine the characteristics of Tropical Cyclone Mekunu obtained from CPL.AOW, CPL.AO, and ATM.DYN to demonstrate the capability of the coupled model. The tracks of the tropical cyclone, defined by the positions of the low-pressure center, are presented in Fig. 2, where it can be seen that all models can qualitatively match the observed evolution and track. Although the translation speed of the tropical cyclone from CPL.AOW is somewhat slower (CPL.AOW: 245 km d−1; IBTrACS: 254 km d−1), the distances between the cyclone centers for all model runs and IBTrACS data are less than 250 km until 26 May, as shown in Fig. 3a.
The characteristics of Tropical Cyclone Mekunu (i.e., cyclone central pressure and maximum wind speed) obtained in the ensemble simulations are compared quantitatively with IBTrACS data in Fig. 3b and c. From 22 to 27 May, the root-mean-square errors (RMSEs) of the cyclone low-pressure center are 9.53, 9.25, and 10.55 hPa for CPL.AOW, CPL.AO, and ATM.DYN, respectively (ensemble standard deviations are 2.85, 2.67, and 2.66 hPa). The RMSEs of the maximum wind speed are 9.06, 8.81, and 9.81 m s−1 for CPL.AOW, CPL.AO, and ATM.DYN (ensemble standard deviations are 3.13, 3.00, and 2.80 m s−1). respectively. In addition, the ensemble mean lowest pressures in CPL.AOW, CPL.AO, and ATM.DYN are higher than the IBTrACS data by 4.9, 5.3, and 17.3 hPa, respectively. The overestimation in ATM.DYN is more significant than 1 standard deviation between 25 to 26 May. For the maximum wind speed, CPL.AOW and CPL.AO underestimate the maximum wind speed by about 6.2 and 5.9 m s−1, while in ATM.DYN the underestimation is as large as 14.1 m s−1. ATM.DYN does not capture the intensification of the TC between 24 and 26 May that is present in the IBTrACS observations and in the coupled simulations.
To highlight the surface fluxes from the simulation, we show the snapshots of the wind speed and latent heat fluxes (LHFs) in Fig. 4. Instead of plotting the entire computational domain, we highlight the region around the center of the tropical cyclone (from 7 to 22∘ N and from 46 to 62∘ E). The 10 m wind speeds obtained in CPL.AOW and CPL.AO are generally consistent, except for the region near the center of the cyclone. Figure 4c shows weaker wind speed and a smaller area of high wind speeds (indicated by the contour of 15 m s−1) in ATM.DYN because the uncoupled run underestimates the intensity of the cyclone. In addition, we present the LHFs because they are the major component of the net surface heat fluxes and they are associated with the water vapor uptake. It can be seen in Fig. 4d–f that the LHFs are weak along the cyclone tracks due to the cold wake (shown in Fig. 5) but are generally consistent in CPL.AOW and CPL.AO. Near the center of the tropical cyclone in Fig. 4, the LHFs in ATM.DYN are weaker than the coupled runs (by a few hundred W m−2). This is because the tropical cyclone is weaker in ATM.DYN, which will be further discussed in Sect. 4.2. The largest differences in cyclone track, intensity, and wind speed are between the uncoupled (ATM.DYN) and coupled (CPL.AO, CPL.AOW) simulations, with the latter being closer to IBTrACS observations.
In summary, both CPL.AOW and CPL.AO runs better simulate the tropical cyclone characteristics than ATM.DYN in comparison with the IBTrACS data. For the cyclone central pressure and wind speed, CPL.AOW is better for extreme conditions but is outperformed by CPL.AO for their RMSEs throughout the event. CPL.AO also better simulates the track of the tropical cyclone. The ATM.DYN also underestimates the latent heat loss from the ocean compared with CPL.AOW and CPL.AO. The differences between the tropical cyclones simulated in the coupled and uncoupled simulations are associated with the SST cooling in the simulations, which are further discussed in Sect. 4.2. The differences between CPL.AOW and CPL.AO are further investigated in Sect. 5.
4.2 SST and mixed-layer depth
To highlight the impact of the tropical cyclone on the ocean, we plot the evolution of SST and ocean MLD from CPL.AOW and ATM.DYN in Fig. 5. The results obtained from CPL.AO are not presented here, but in Sect. 5 we investigate the effect of wave coupling on SST and MLD in the coupled system. Figure 5b and c show the differences of ensemble-averaged SST between 20 and 27 May (27 May SST minus 20 May SST). This aims to highlight the development of an SST cold wake during this event. It can be seen that the SST cools down by a maximum of 4 ∘C along the TC track, which can impact the ocean and air–sea interactions (Price, 1981; Stramma et al., 1986; Pasquero et al., 2021). It can be seen in Fig. 5 that the SST cooling in CPL.AOW is weaker than that in HYCOM, indicating the SST is warmer throughout the simulation in CPL.AOW. Contributed by the warmer SST, the intensity of the tropical cyclone is also stronger in CPL.AOW than ATM.DYN. Due to the chaotic nature of the atmosphere, it is still unknown why the warmer SST in CPL.AOW improves the simulation of tropical cyclone characteristics. In addition, although we compared the simulation results using available data, the lack of in situ observations in this region makes it challenging to validate the SST used the simulations. It should be noted that CPL.AOW also captures the SST warming in the Arabian Gulf and the Gulf of Aden compared with HYCOM data.
The evolution of the mixed-layer depth during the event is shown in Fig. 6 to demonstrate the impact of the tropical cyclone on ocean mixing. Here we are using Δρ=0.03 kg m−3 to define the MLD. The initial MLD is about 30–40 m along the track of the tropical cyclone. To highlight the evolution of the mixed layer, the differences of ensemble-averaged MLD between 20 and 27 May (27 May MLD minus 20 May MLD) are plotted in Fig. 6b and c. It can be seen that the MLD deepens by approximately 30–40 m (the standard deviation is about 10 m) along the track of the tropical cyclone, which is almost a 100 % increase compared to its initial value in Fig. 6a. It is noted that CPL.AOW has stronger MLD deepening than HYCOM but weaker SST cooling. We hypothesize that this is because (1) the parameterization of the ocean mixing layer is different when the effects of Langmuir turbulence are considered in CPL.AOW and (2) the atmosphere forcing used in the coupled model has a higher spatial and temporal resolution that makes the SST and MLD different.
Ocean surface waves are expected to affect air–sea interactions. The ensemble mean significant wave height (Hs) and the ensemble standard deviation obtained from the coupled simulation are shown in Fig. 7. The snapshots of the ensemble-averaged Hs are presented in Fig. 7a–c. On 26 May, the ensemble-averaged Hs is as high as 8 m near the eye wall of the tropical cyclone. Figure 7b and c shows that alternating regions of high and low waves can be observed between 12 and 24∘ N and between 60 and 75∘ E. The spatial pattern of high and low beams of Hs is due to surface wave refraction by ocean currents. We have performed uncoupled WW3 simulations to investigate these beams, and more details can be found in Sun et al. (2022).
Near the eye wall of the tropical cyclone, the standard deviation of Hs from the ensembles is approximately 3 m, showing greater variance near the eye wall (Fig. 7d–f), while the spatial variability of Hs, with alternating high and low beams, is consistent throughout the ensemble. Although CPL.AOW captures the overall spatial variability of Hs, the exact location of the beams deviates from the altimetry observations, since the central location of the tropical cyclone is not well captured by CPL.AOW. The comparison of the modeled Hs with altimeter data is shown in Appendix D.
To explore the effects of the surface waves, coupled simulations were run using three recent parameterizations of Langmuir turbulence. We compared the characteristics of the tropical cyclone (e.g., track, intensity, and wind speed) and the changes in the ocean (e.g., SST and MLD). In the sensitivity analysis, we compare the simulation results without Langmuir turbulence (NoLT) and those with Langmuir turbulence (LF17, VR12-MA, and LF17-ST). This aims to illustrate the sensitivity of the coupled model to Langmuir turbulence, which may deepen the MLD and cool the SST during a tropical cyclone event. Note that the coupled run using LF17 is identical to CPL.AOW in Sect. 4. To evaluate the effect of Langmuir turbulence on the ocean, we also performed the simulations using spectral nudging in WRF in addition to the “free runs” (simulations without spectral nudging). The spectral nudging is performed because of the uncertainties of the atmosphere model, especially for the tracks of the cyclones. By restraining the uncertainty of the atmosphere using spectral nudging, we are able to highlight the impact of Langmuir turbulence on the ocean. In the present study, WRF nudges the model fields to NCEP FNL data and the wavenumber for spectral nudging is about 600 km. Although the simulations with spectral nudging have smaller internal variability, they underestimate the intensity of the tropical cyclone in the simulations. Similar to the simulation results shown in Sect. 4, the characteristics of the cyclones (i.e., tracks, cyclone central pressure, wind speed) are not significantly different from the simulations using different parameterizations. The cyclone characteristics in the sensitivity analysis are detailed in Appendix C. In this section, we only highlight the sensitivity of SST, ocean mixed layer, and other surface fluxes to the parameterization of Langmuir turbulence. We also performed a similar sensitivity analysis for different surface roughness parameters that may impact the atmosphere surface variables. The sensitivity analyses of Stokes advection, Stokes Coriolis, and wave-induced momentum fluxes are also performed. However, these results are summarized in Appendix C because they are not significantly different.
5.1 SST and ocean mixed layer
To highlight the SST differences between the simulations, we plotted the SST differences between the runs with and without Langmuir turbulence in Fig. 8, with regions of significant SST changes (P<0.05) highlighted. It can be seen that in CPL.LF17 and CPL.LF17-ST the SST cooling near the tracks of the cyclone is stronger by about 0.5 ∘C in comparison with CPL.NoLT due to the effect of Langmuir turbulence. When the spectral nudging is added to reduce the randomness of the atmosphere model, the cyclones within the ensemble simulations are more similar and thus the SST cooling is more significant. The results obtained using LF17 and LF17-ST are generally consistent because they use a similar way to calculate the enhancement coefficient and entrainment flux. Their differences are because of different options to parameterize the Langmuir number La. On the other hand, when using VR12-MA, we observed weaker SST cooling compared with the simulation without Langmuir turbulence. Though it is demonstrated in Reichl et al. (2016b) that VR12-MA is not adequate to parameterize the Langmuir turbulence, this non-intuitive SST change needs to be documented and discussed. After we examine the vertical profiles of the ocean, we hypothesize that too much diffusivity is added to the ocean current velocity when using VR12-MA in this case. The reduction of vertical gradient in ocean current velocity reduces the ocean mixing in KPP and contributes to a weaker SST cooling. More details on the SST cooling are presented in Sect. 5.2.
The comparison between the MLDs obtained from the simulations is shown in Fig. 9. Again, we highlight the regions with significant MLD changes (P<0.05) for both “free run” and those with spectral nudging. Due to Langmuir turbulence, the MLDs in CPL.LF17 and CPL.LF17-ST are deeper by a maximum of about 20 m than that of CPL.NoLT. When using VR12-MA, the MLD is shallower again due to the reduction of turbulent shear from the parameterization of Langmuir turbulence in this case. It is noted that the largest SST and MLD changes are not centered on the location of the tropical cyclone. We hypothesize that this is because (1) the SST and MLD changes need some time to develop and (2) the winds on the right-front quadrant of the cyclone are strongest (Moon et al., 2004; Fan et al., 2009).
5.2 The vertical profiles
To investigate the SST warming and cooling in the wake zone due to Langmuir turbulence, we examined the vertical profiles of the ocean aiming to illustrate the impact of different Langmuir turbulence options. Here we averaged the quantities of interest in the region between 11 to 14∘ N and from 55 to 57∘ E because large SST and MLD changes are observed in this region. Due to the internal variability of the atmosphere, we only compare the simulation results when spectral nudging is used.
To analyze the impact of different Langmuir turbulence options, in Fig. 10a–d we plotted the domain-averaged bulk Richardson number, buoyancy difference, vertical density gradient, and current velocity, which are the dominant terms in Eq. (6). The bulk Richardson number is plotted because it is used in the MITgcm KPP scheme to determine the boundary layer depth, which is crucial to parameterize vertical mixing. The critical Richardson number Ricr is 0.3, meaning the ocean is assumed dynamically unstable and turbulent when Ri<0.3. It can be seen that when VR12-MA is applied, the Richardson number increases compared with CPL.NoLT (no Langmuir turbulence), indicating the parameterized turbulence is getting weaker. Examining each component of the Richardson number in Eq. (6), it can be seen that buoyancy and vertical density gradient terms do not change significantly, shown in Fig. 10b and c, while the changes in horizontal current speed can be seen in Fig. 10d near the surface.
When VR12-MA is applied, the Langmuir enhancement factor (see Eq. 3) is used to amplify the KPP diffusivity term (see Eq. 5). This reduces the vertical gradient of horizontal current, shown in Fig. 10d. When the velocity gradient is reduced in VR12-MA, the term in Eq. (6) decreases, and thus the Richardson number increases. This Richardson number increase results in a decrease in estimated MLD, and thus the SST cooling is weaker than the simulation without Langmuir turbulence (CPL.NoLT). To verify this, we run an ensemble of simulations (CPL.VR12-MA-NoU) that do not enhance the KPP diffusivity for horizontal currents, then we observed cooler SST due to enhanced mixing by Langmuir turbulence. The simulation results of this verification test are detailed in Appendix D. It is noted that the drawbacks of VR12-MA are also discussed in Reichl et al. (2016b), showing that using Lagrangian currents uL on parameterizing Langmuir turbulence can also alleviate the bias when using VR12-MA.
On the other hand, when using LF17 in the simulations, the same enhancement factor ε as VR12-MA is added, but the VtL(z) term is used in Eq. (7) for parameterizing the Richardson number. Although the velocity gradient is also smaller, shown in Fig. 10d, the entrainment flux Vtl(z) decreases the Richardson number. This implies stronger vertical mixing due to the Langmuir entrainment by the tropical cyclone. Hence, the SST cooling in the near wake region of the tropical cyclone is stronger when LF17 is used than VR12-MA. This shows parameterizing the Langmuir turbulence using LF17 gives more realistic results than VR12-MA.
This work described the integration of WAVEWATCH III into the SKRIPS regional coupled model. The implementation allows for the use or disuse of the wave model in the SKRIPS model. The parameterizations of the surface waves are implemented in MITgcm to account for the impact of waves on the ocean and the Stokes advection and Stokes Coriolis forces in the coupled model.
To test the coupled ocean–wave–atmosphere model, we performed a series of simulations of Tropical Cyclone Mekunu in the Arabian Sea, which is a representative tropical cyclone case. In order to model the uncertainty due to the atmospheric internal variabilities, we added small perturbations to the initial SST to generate 20 ensembles for each test case. We then compared the fully coupled simulations (CPL.AOW) with coupled runs without the wave model (CPL.AO) and stand-alone atmosphere simulations (ATM.DYN). The characteristics of the tropical cyclone (e.g., track, intensity, and wind speed) obtained in the two coupled simulations are similar within uncertainty. However, the stand-alone atmosphere model sees SST cooling from the HYCOM analysis that is stronger than in the coupled runs. Compared with the coupled simulations, in ATM.DYN the tropical cyclone has higher pressure and lower wind speed, making it less consistent with the observations.
We further tested the sensitivity of simulated characteristics of cyclone to wave coupling. The simulation results show that the characteristics of the tropical cyclone (e.g., intensity, pressure, maximum wind speed) are not sensitive to wave coupling compared with the internal variability of the model as resolved by the ensemble simulations. When the effect of Langmuir turbulence is parameterized using the LF17 and LF17-ST options that account for the Langmuir entrainment, the maximum SST cooling is about 0.5 ∘C cooler and the maximum mixed-layer deepening is about 20 m along the track of the tropical cyclone, indicating the surface waves play an important role in modulating the response of the upper ocean to tropical cyclone surface forcing. On the other hand, when the effect of Langmuir turbulence is parameterized using the VR12-MA, the SST cooling and MLD deepening are weaker due to the changes of current shear in the coupled simulation.
The results presented here motivate further studies to evaluate and improve this and other regional or high-resolution coupled models for investigating dynamical processes and forecasting applications, especially the interaction between the ocean and waves and their feedback with the atmosphere.
In this work, we used the Langmuir number summarized in Li et al. (2019). The Langmuir number LaLP is defined from the surface layer (a fraction of the mixed layer, which is the upper 20 % in their definition) averaged Stokes drift 〈uS〉SL and a reference Stokes drift near the base of the mixed layer:
The Langmuir number LaLP is used in the scaling of Langmuir-enhanced entrainment in LF17 in Eq. (8).
The projected Langmuir number LaSLP considers the misalignment between wind and waves:
where θww is the misalignment between wind and waves and θwl is the misalignment between wind and Langmuir cells. The Langmuir number LaSLP is used in VR12-MA, and the scaling of Langmuir-enhanced diffusivity is used in LF17 in Eq. (4).
In this work, we implemented three different ocean roughness closure models in (Olabarrieta et al., 2012): (1) a DGHQ model based on wave age (Drennan et al., 2003), (2) a TY2001 model based on wave steepness (Taylor and Yelland, 2001), and (3) a OOST model that considers both the effects of wave age and steepness (Oost et al., 2002). These options are implemented in the MYNN surface layer scheme in WRF.
where z0 is the ocean roughness, Lp is the wavelength at the peak of the wave spectrum, and Hs is the significant wave height.
where u* is the wind friction velocity and Cp is the wave phase speed at the peak frequency.
We also implemented another option that uses the Charnock coefficient (CHNK) calculated from WAVEWATCH III. In the present study, we used the ST4 option in WW3, and thus the Charnock coefficients are calculated based on Ardhuin et al. (2010). In this paper we used WAVEWATCH III version 6.0.7 compiled with the following switches:
F90 NOGRB NOPA LRB4 SCRIP SCRIPNC NC4 TRKNC DIST MPI PR3 UQ FLX0 LN1 ST4 STAB0 NL1 BT4 DB1 MLIM TR0 BS0 IC2 IS2 REF1 IG0 XX0 WNT2 WNX1 CRT1 CRX1 TIDE O0 O1 O2 O2a O2b O2c O3 O4 O5 O6 O7
The comparison between the simulations of Tropical Cyclone Mekunu as resulting from the coupled models using different setups. Here we tested different setups in parameterizing Langmuir turbulence and sea surface roughness. For the Langmuir turbulence we test VR12-MA, LF17, LF17-ST, and no Langmuir turbulence; for the sea surface roughness we test CHNK, TY2001, DGHQ, and OOST. The simulation using CHNK is the same as CPL.AOW in Sect. 4. Finally, we examine the impact of Stokes forces and wind stress.
Figure C1 shows that the tracks of tropical cyclones from the coupled simulations are generally consistent within the ensemble spread, although the track from CPL.AOW (CHNK) is slightly closer to IBTrACS than the other simulations. The distance error, simulated cyclone central pressure, and maximum wind speed shown in Figs. C2 and C3 are also close. Note that the differences between the ensemble-mean pressure and wind speed in the coupled models are smaller than the standard deviations shown in Figs. C2 and C3. The snapshots of the 10 m wind speed and latent heat fluxes in Fig. C5 aim to illustrate the sensitivity of the surface atmosphere to parameterizing surface roughness. It can be seen that the 10 m wind speed and latent heat loss are different when using TY2001, DGHQ, and OOST in comparison with using CHNK in CPL.AOW. It can be seen that the simulations using TY2001, DGHQ, and OOST has stronger wind obtained from the simulations, similar to the findings in Olabarrieta et al. (2012).
The impact of Stokes forces and wind stress on the characteristics of the tropical cyclone is shown in Fig. C4. The coupled simulation results without using Stokes advection and Stokes Coriolis are shown using “NoStokes”; the results without correcting the wind stress due to the waves are shown using “NoStress”. It can be seen in the simulation results that the effects of Stokes forces and wind stress are not significant. Note that the NoStokes experiment is consistent with the implicit scheme in Li et al. (2016).
In addition to the experiment performed without Stokes advection and Stokes Coriolis, we performed a test to further illustrate the difference when Langmuir turbulence options are applied. Simulations are performed to test VR12-MA, LF17, and LF17-ST options in comparison with no Langmuir turbulence (NoLT). The SST differences for these simulations are shown in Fig. C6. It can be seen that SST differences in Fig. C6a–c are generally consistent with those shown in Fig. 8, except for the regions near the track of the tropical cyclone where the uncertainty is large. The SST cooling and warming patterns obtained from the spectral nudging experiments in Fig. C6d–f are generally consistent with those shown in Fig. 8. This demonstrates that the simulation results are not sensitive to the implicit or explicit options for Stokes advection and Stokes Coriolis.
To illustrate the impact of velocity shear in KPP, we performed a 20-member ensemble simulation using VR12-MA scheme but do not enhance the velocity scale when calculating the KPP diffusivity. The simulation results of this experiment (CPL.VR12-MA-NoU) are shown in Fig. D1. In this experiment, It can be seen here that when the diffusivity is not enhanced, the MLD deepens by about 20 m and SST cools down by about 0.5 ∘C. This indicates that the diffusivity of the velocity scale in VR12-MA makes an impact on the mixing layer and SST when applied to parameterize the Langmuir turbulence.
To evaluate the simulation performance of surface waves, we compared the modeled Hs with along-track Hs measurements from the Jason‐3 and SARAL/AltiKa altimeters. We use quality‐controlled, unfiltered, and not resampled along‐track Hs measurements provided by the Institut Français de Recherche pour l'Exploitation de la MER (IFREMER; ftp://ftp.ifremer.fr/ifremer/cersat/products/swath/altimeters/waves/, last access: 9 June 2023, Queffeulou and Croizé-Fillon, 2013).
The comparison of the significant wave height Hs with the altimeter data is shown in Fig. E1. We also plotted the simulation results obtained in Sun et al. (2022). WAV.WND indicates the simulation using a stand-alone WAVEWATCH III model driven by ERA5 wind only; WAV.CUR indicates the simulation using a stand-alone WAVEWATCH III model driven by ERA5 wind and HYCOM currents. It can be seen from Fig. E1 that the coupled model captures the focusing and defocusing of the waves. However, because of the error in the location of the tropical cyclone, the patterns of the Hs are not completely consistent with the observational data.
To illustrate the impact of Langmuir turbulence on the upper ocean, a sensitivity analysis of the SST and the ocean mixed layer is performed. First we compare the SST cooling obtained in the simulations with the observational data from Global Drifter Program (Lumpkin and Centurioni, 2019). In Fig. F1 we show the SST changes throughout the event along the drifter tracks. It can be seen that a drifter closest to the track of the TC (highlighted in Fig. F1) recorded a cooling of 3.81 ∘C, while the SST cooling in both CPL.NoLT and CPL.LF17 is not significantly different (CPL.NoLT: 2.18 ∘C; CPL.LF17: 2.17 ∘C; standard deviation: 0.37 ∘C). It is noted that the SST cooling in HYCOM is 3.23 ∘C, which is closer to the drifter data than the coupled simulations. However, because the in situ observations are sparse in this region, future work still needs to be done to evaluate the performance of the coupled model.
The source code of the coupled model is maintained on GitHub (https://github.com/iurnus/scripps_kaust_model, last access: 9 June 2023) and Zenodo (https://doi.org/10.5281/zenodo.7972577, Sun, 2020). The code documentation is available at https://github.com/iurnus/scripps_kaust_model_doc, last access: 9 June 2023.
All validation data used in this work are publicly available. The HYCOM/NCODA daily global analysis data are available at https://www.hycom.org/dataserver/gofs-3pt0 (Lumpkin and Centurioni, 2019). The drifter data are provided by Global Drifter Program at https://www.aoml.noaa.gov/envids/gld/index.php (Naval Research Laboratory, 2014–2021). The wave height data are provided by the Institut Français de Recherche pour l'Exploitation de la MER at ftp://ftp.ifremer.fr/ifremer/cersat/products/swath/altimeters/waves/ (Queffeulou and Croizé-Fillon, 2013). The IBTrACs data are obtained from https://doi.org/10.25921/82ty-9e16 (Knapp et al., 2018). The source code of the coupled model is maintained on GitHub https://github.com/iurnus/scripps_kaust_model and Zenodo https://doi.org/10.5281/zenodo.7972577 (Sun, 2020).
RS worked on the coding tasks for integrating WW3 to the coupled system, wrote the code documentation, and performed the ensemble simulation. RS and AC drafted the initial manuscript. RS, ABVB, and SL implemented the WW3 wave model. All authors designed the computational framework and the numerical experiments. All authors discussed the results and contributed to the writing of the final manuscript.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We gratefully acknowledge the research funding (grant no. OSR-2016-RPP-3268.02) from KAUST (King Abdullah University of Science and Technology). We also appreciate the computational resources of supercomputer Shaheen II and the assistance provided by KAUST Supercomputer Laboratory. Rui Sun and Aneesh C. Subramanian were supported by ONR ASTRAL research initiative (N00014-23-1-2092). Ana B. Villas Bôas was supported by NASA award no. 80NSSC19K1004 through the S-MODE program. Aneesh C. Subramanian was supported by NOAA grant no. NA18OAR4310405 and ONR MISOBOB research initiative (N00014-17-S-B001). Bruce D. Cornuelle and Matthew R. Mazloff were supported by NOAA grant nos. NA21OAR4310257, NA18OAR4310403, and NA22OAR4310597. Arthur J. Miller was partly supported by the National Science Foundation (OCE-2022868). We thank Baylor Fox-Kemper, Qing Li, and Zhihua Zheng for discussing the implementation of the Langmuir turbulence parameterizations. Rui Sun, Ana B. Villas Bôas, and Matthew R. Mazloff were supported by NASA award no. 80NSSC23K0979.
This research has been supported by the King Abdullah University of Science and Technology (grant no. OSR-2016-RPP-3268.02), the National Oceanic and Atmospheric Administration (grant nos. NA21OAR4310257, NA21OAR4310403, NA21OAR4310405, and NA22OAR4310597), the National Science Foundation (grant no. OCE-2022868), the National Aeronautics and Space Administration (grant nos. 80NSSC19K1004, 80NSSC23K0979), and the Office of Naval Research (grant no. N00014-17-S-B001).
This paper was edited by Deepak Subramani and reviewed by three anonymous referees.
Ardhuin, F., O'reilly, W., Herbers, T., and Jessen, P.: Swell transformation across the continental shelf. Part I: Attenuation and directional broadening, J. Phys. Oceanogr., 33, 1921–1939, 2003. a
Ardhuin, F., Rogers, E., Babanin, A. V., Filipot, J.-F., Magne, R., Roland, A., Van Der Westhuysen, A., Queffeulou, P., Lefevre, J.-M., Aouf, L., and Collard, F.: Semiempirical dissipation source functions for ocean waves. Part I: Definition, calibration, and validation, J. Phys. Oceanogr., 40, 1917–1941, 2010. a
Bender, M. A. and Ginis, I.: Real-case simulations of hurricane–ocean interaction using a high-resolution coupled model: Effects on hurricane intensity, Mon. Weather Rev., 128, 917–946, 2000. a
Bhatia, K., Vecchi, G., Murakami, H., Underwood, S., and Kossin, J.: Projected response of tropical cyclone intensity and intensification in a global climate model, J. Climate, 31, 8281–8303, 2018. a
Blair, A., Ginis, I., Hara, T., and Ulhorn, E.: Impact of Langmuir turbulence on upper ocean response to Hurricane Edouard: Model and observations, J. Geophys. Res.-Oceans, 122, 9712–9724, 2017. a
Breivik, Ø., Janssen, P. A., and Bidlot, J.-R.: Approximate Stokes drift profiles in deep water, J. Phys. Oceanogr., 44, 2433–2445, 2014. a
Breivik, Ø., Bidlot, J.-R., and Janssen, P. A.: A Stokes drift approximation based on the Phillips spectrum, Ocean Model., 100, 49–56, 2016. a
Campin, J.-M., Heimbach, P., Losch, M., Forget, G., edhill3, Adcroft, A., amolod, Menemenlis, D., dfer22, Hill, C., Jahn, O., Scott, J., stephdut, Mazloff, M., baylorfk, antnguyen13, Doddridge, E., Fenty, I., Bates, M., Martin, T., Abernathey, R., samarkhatiwala, Smith, T., Lauderdale, J., hongandyan, Deremble, B., raphael dussin, Bourgault, P., dngoldberg, and T., A. T.: MITgcm/MITgcm: checkpoint67m, Zenodo [code], https://doi.org/10.5281/zenodo.3492298, 2019. a
Cerovečki, I., Sun, R., Bromwich, D. H., Zou, X., Mazloff, M. R., and Wang, S.-H.: Impact of downward longwave radiative deficits on Antarctic sea-ice extent predictability during the sea ice growth period, Enviro. Res. Lett., 17, 084008, https://doi.org/10.1088/1748-9326/ac7d66, 2022. a
Charnock, H.: Wind stress on a water surface, Q. J. Roy. Meteor. Soc., 81, 639–640, 1955. a
Chassignet, E. P., Hurlburt, H. E., Smedstad, O. M., Halliwell, G. R., Hogan, P. J., Wallcraft, A. J., Baraille, R., and Bleck, R.: The HYCOM (hybrid coordinate ocean model) data assimilative system, J. Marine Syst., 65, 60–83, 2007. a, b
Chen, S. S., Price, J. F., Zhao, W., Donelan, M. A., and Walsh, E. J.: The CBLAST-Hurricane program and the next-generation fully coupled atmosphere–wave–ocean models for hurricane research and prediction, B. Am. Meteorol. Soc., 88, 311–318, 2007a. a
Chen, S. S., Price, J. F., Zhao, W., Donelan, M. A., and Walsh, E. J.: The CBLAST-Hurricane program and the next-generation fully coupled atmosphere–wave–ocean models for hurricane research and prediction, B. Am. Meteorol. Soc., 88, 311–318, 2007b. a, b
Couvelard, X., Lemarié, F., Samson, G., Redelsperger, J.-L., Ardhuin, F., Benshila, R., and Madec, G.: Development of a two-way-coupled ocean–wave model: assessment on a global NEMO(v3.6)–WW3(v6.02) coupled configuration, Geosci. Model Dev., 13, 3067–3090, https://doi.org/10.5194/gmd-13-3067-2020, 2020. a
Craik, A. D. and Leibovich, S.: A rational model for Langmuir circulations, J. Fluid Mech., 73, 401–426, 1976. a
D'Asaro, E. A.: Turbulent vertical kinetic energy in the ocean mixed layer, J. Phys. Oceanogr., 31, 3530–3537, 2001. a
D'Asaro, E. A., Thomson, J., Shcherbina, A., Harcourt, R., Cronin, M., Hemer, M., and Fox-Kemper, B.: Quantifying upper ocean turbulence driven by surface waves, Geophys. Res. Lett., 41, 102–107, 2014. a, b
Drennan, W. M., Graber, H. C., Hauser, D., and Quentin, C.: On the wave age dependence of wind stress over pure wind seas, J. Geophys. Res.-Oceans, 108, https://doi.org/10.1029/2000JC000715, 2003. a, b, c
Dube, S. K., Rao, A. D., Sinha, P. C., Murty, T. S., and Bahulayan, N.: Storm surge in the Bay of Bengal and Arabian Sea the problem and its prediction, Mausam, 48, 283–304, 1997. a
Emanuel, K. A.: The theory of hurricanes, Annu. Rev. Fluid. Mech., 23, 179–196, 1991. a
Evan, A. T. and Camargo, S. J.: A climatology of Arabian Sea cyclonic storms, J. Climate, 24, 140–158, 2011. a
Evan, A. T., Kossin, J. P., and Ramanathan, V.: Arabian Sea tropical cyclones intensified by emissions of black carbon and other aerosols, Nature, 479, 94–97, 2011. a
Government of India, Ministry of Earth Sciences, I. M. D.: Extremely severe cyclonic storm, “MEKUNU” over the Arabian Sea (21–27 May 2018): A report., Tech. rep., New Delhi: India Meteorological Department, https://rsmcnewdelhi.imd.gov.in/uploads/report/26/26_e580e2_mekunu.pdf (last access: 9 June 2023), 2018. a, b, c
Henderson-Sellers, A., Zhang, H., Berz, G., Emanuel, K., Gray, W., Landsea, C., Holland, G., Lighthill, J., Shieh, S.-L., Webster, P., and McGuffie, K.: Tropical cyclones and global climate change: A post-IPCC assessment, B. Am. Meteorol. Soc., 79, 19–38, 1998. a
Iacono, M. J., Delamere, J. S., Mlawer, E. J., Shephard, M. W., Clough, S. A., and Collins, W. D.: Radiative forcing by long-lived greenhouse gases: Calculations with the AER radiative transfer models, J. Geophys. Res.-Atmos., 113, https://doi.org/10.1029/2008JD009944, 2008. a
Jenkins, A. D.: The use of a wave prediction model for driving a near-surface current model, Deutsche Hydrografische Zeitschrift, 42, 133–149, 1989. a
Kain, J. S.: The Kain–Fritsch convective parameterization: an update, J. Appl. Meteorol., 43, 170–181, 2004. a
Knapp, K. R., Kruk, M. C., Levinson, D. H., Diamond, H. J., and Neumann, C. J.: The international best track archive for climate stewardship (IBTrACS) unifying tropical cyclone data, B. Am. Meteoro. Soc., 91, 363–376, 2010. a
Knapp, K. R., Diamond, H. J., Kossin, J. P., Kruk, M. C., and Schreck, C. J.: The international best track archive for climate stewardship (IBTrACS) unifying tropical cyclone data, version 4, Tech. rep., [NI – North Indian], NOAA National Centers for Environmental Information, https://doi.org/10.25921/82ty-9e16, 2018. a, b
Langmuir, I.: Surface motion of water induced by wind, Science, 87, 119–123, 1938. a
Large, W. G. and Yeager, S. G.: Diurnal to decadal global forcing for ocean and sea-ice models: the data sets and flux climatologies, Tech. rep., NCAR Technical Note: NCAR/TN-460+STR, CGD Division of the National Center for Atmospheric Research, 2004. a
Lewis, H. W., Castillo Sanchez, J. M., Arnold, A., Fallmann, J., Saulter, A., Graham, J., Bush, M., Siddorn, J., Palmer, T., Lock, A., Edwards, J., Bricheno, L., Martínez-de la Torre, A., and Clark, J.: The UKC3 regional coupled environmental prediction system, Geosci. Model Dev., 12, 2357–2400, https://doi.org/10.5194/gmd-12-2357-2019, 2019. a, b
Li, Q. and Fox-Kemper, B.: Assessing the effects of Langmuir turbulence on the entrainment buoyancy flux in the ocean surface boundary layer, J. Phys. Oceanogr., 47, 2863–2886, 2017. a, b, c, d, e, f, g
Li, Q., Webb, A., Fox-Kemper, B., Craig, A., Danabasoglu, G., Large, W. G., and Vertenstein, M.: Langmuir mixing effects on global climate: WAVEWATCH III in CESM, Ocean Model., 103, 145–160, 2016. a, b, c, d, e, f, g, h, i, j
Li, Q., Reichl, B. G., Fox-Kemper, B., Adcroft, A. J., Belcher, S. E., Danabasoglu, G., Grant, A. L., Griffies, S. M., Hallberg, R., Hara, T., Harcourt, R. R., Kukulka, T., Large, W. G., McWilliams, J. C., Pearson, B., Sullivan, P. P., Roekel, L. V., Wang, P., and Zheng, Z.: Comparing ocean surface boundary vertical mixing schemes including Langmuir turbulence, J. Adv. Model. Earth Sy., 11, 3545–3592, 2019. a, b, c, d, e
Li, Z., Tam, C.-Y., Li, Y., Lau, N.-C., Chen, J., Chan, S., Lau, D.-S. D., and Huang, Y.: How does air-sea wave interaction affect tropical cyclone intensity? An atmosphere–wave–ocean coupled model study based on super typhoon Mangkhut (2018), Earth Space Sci., 9, e2021EA002136, https://doi.org/10.1029/2021EA002136, 2022. a
Liu, B., Liu, H., Xie, L., Guan, C., and Zhao, D.: A coupled atmosphere–wave–ocean modeling system: Simulation of the intensity of an idealized tropical cyclone, Mon. Weather Rev., 139, 132–152, 2011. a
Liu, W., Katsaros, K., and Businger, J.: Bulk parameterization of air-sea exchanges of heat and water vapor including the molecular constraints at the interface, J. Atmos. Sci., 36, 1722–1735, 1979. a
Lumpkin, R. and Centurioni, L.: Global Drifter Program quality-controlled 6-hour interpolated data from ocean surface drifting buoys, Tech. rep., NOAA National Centers for Environmental Information [data set], https://doi.org/10.25921/7ntx-z961, 2019. a, b
Lumpkin, R. and Pazos, M.: Measuring surface currents with Surface Velocity Program drifters: the instrument, its data, and some recent results, Lagrangian analysis and prediction of coastal and ocean dynamics, 39, 2007. a
Marshall, J., Adcroft, A., Hill, C., Perelman, L., and Heisey, C.: A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers, J. Geophys. Res.-Oceans, 102, 5753–5766, 1997. a
Mogensen, K. S., Magnusson, L., and Bidlot, J.-R.: Tropical cyclone sensitivity to ocean coupling in the ECMWF coupled model, J. Geophys. Res.-Oceans, 122, 4392–4412, 2017. a
Monin, A. S. and Obukhov, A. M.: Basic laws of turbulent mixing in the surface layer of the atmosphere, Contrib. Geophys. Inst. Acad. Sci. USSR, 151, e187, 163–187, 1954. a
Moon, I.-J., Ginis, I., and Hara, T.: Effect of surface waves on air–sea momentum exchange. Part II: Behavior of drag coefficient under tropical cyclones, J. Atmos. Sci., 61, 2334–2348, 2004. a
Morrison, H., Thompson, G., and Tatarskii, V.: Impact of cloud microphysics on the development of trailing stratiform precipitation in a simulated squall line: Comparison of one-and two-moment schemes, Mon. Weather Rev., 137, 991–1007, 2009. a
Murakami, H., Vecchi, G. A., and Underwood, S.: Increasing frequency of extremely severe cyclonic storms over the Arabian Sea, Nat. Clim. Change, 7, 885–889, 2017. a
Nakanishi, M. and Niino, H.: An improved Mellor–Yamada level-3 model with condensation physics: Its design and verification, Bound.-Lay. Meteorol., 112, 1–31, 2004. a
Nakanishi, M. and Niino, H.: Development of an improved turbulence closure model for the atmospheric boundary layer, J. Meteorol. Soc. Jpn. Ser. II, 87, 895–912, 2009. a
Naval Research Laboratory: Global Ocean Forecast System (GOFS) 3.1, [data set], https://www.hycom.org/dataserver/gofs-3pt1/analysis, 2014–2021. a
Olabarrieta, M., Warner, J. C., Armstrong, B., Zambon, J. B., and He, R.: Ocean–atmosphere dynamics during Hurricane Ida and Nor'Ida: an application of the coupled ocean–atmosphere–wave–sediment transport (COAWST) modeling system, Ocean Model., 43, 112–137, 2012. a, b, c, d
Price, J. F.: Upper ocean response to a hurricane, J. Phys. Oceanogr., 11, 153–175, 1981. a
Queffeulou, P. and Croizé-Fillon, D.: Global altimeter SWH data set, version 10., Tech. rep., Laboratoire d'Ocánographie Spatiale, IFREMER, Plouzané, France, ftp://ftp.ifremer.fr/ifremer/cersat/products/swath/altimeters/waves/ (last access: 9 June 2023), 2013. a, b
Rabe, T. J., Kukulka, T., Ginis, I., Hara, T., Reichl, B. G., D’Asaro, E. A., Harcourt, R. R., and Sullivan, P. P.: Langmuir turbulence under hurricane Gustav (2008), J. Phys. Oceanogr., 45, 657–677, 2015. a
Rascle, N. and Ardhuin, F.: A global wave parameter database for geophysical applications. Part 2: Model validation with improved source term parameterization, Ocean Model., 70, 174–188, 2013. a
Reichl, B. G., Ginis, I., Hara, T., Thomas, B., Kukulka, T., and Wang, D.: Impact of sea-state-dependent Langmuir turbulence on the ocean response to a tropical cyclone, Mon. Weather Rev., 144, 4569–4590, 2016a. a
Renault, L., Masson, S., Arsouze, T., Madec, G., and McWilliams, J. C.: Recipes for how to force oceanic model dynamics, J. Adv. Model. Earth Sy., 12, e2019MS001715, https://doi.org/10.1029/2019MS001715, 2020. a
Salih, A. A., Baraibar, M., Mwangi, K. K., and Artan, G.: Climate change and locust outbreak in East Africa, Nat. Clim. Change, 10, 584–585, 2020. a
Sauvage, C., Seo, H., Clayson, C. A., and Edson, J. B.: Impacts of waves and sea states on air-sea momentum flux in the Northwest Tropical Atlantic Ocean: parameterization and wave coupled climate modeling, ESS Open Archive [preprint], 30 pp., https://doi.org/10.1002/essoar.10512415.1, 2022. a
Saxby, J., Crook, J., Peatman, S., Birch, C., Schwendike, J., Valdivieso da Costa, M., Castillo Sanchez, J. M., Holloway, C., Klingaman, N. P., Mitra, A., and Lewis, H.: Simulations of Bay of Bengal tropical cyclones in a regional convection-permitting atmosphere–ocean coupled model, Weather Clim. Dynam. Discuss. [preprint], https://doi.org/10.5194/wcd-2021-46, 2021. a
Schultz, C., Doney, S. C., Zhang, W. G., Regan, H., Holland, P., Meredith, M., and Stammerjohn, S.: Modeling of the influence of sea ice cycle and Langmuir circulation on the upper ocean mixed layer depth and freshwater distribution at the West Antarctic Peninsula, J. Geophys. Res.-Oceans, 125, e2020JC016109, https://doi.org/10.1029/2020JC016109, 2020. a, b, c
Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Liu, Z., Berner, J., Wang, W., Powers, J. G., Duda, M. G., Barker, D. M., and Huang, X.-Y.: A description of the Advanced Research WRF Version 4, Tech. rep., NCAR Technical Note: NCAR/TN-556+STR, 145 pp., https://doi.org/10.5065/1dfh-6p97, 2019. a
Smith, S. D.: Coefficients for sea surface wind stress, heat flux, and wind profiles as a function of wind speed and temperature, J. Geophys. Res.-Oceans, 93, 15467–15472, 1988. a
Stramma, L., Cornillon, P., and Price, J. F.: Satellite observations of sea surface cooling by hurricanes, J. Geophys. Res.-Oceans, 91, 5031–5035, 1986. a
Sun, R., Subramanian, A. C., Miller, A. J., Mazloff, M. R., Hoteit, I., and Cornuelle, B. D.: SKRIPS v1.0: a regional coupled ocean–atmosphere modeling framework (MITgcm–WRF) using ESMF/NUOPC, description and preliminary results for the Red Sea, Geosci. Model Dev., 12, 4221–4244, https://doi.org/10.5194/gmd-12-4221-2019, 2019. a, b, c, d, e
Sun, R., Subramanian, A. C., Cornuelle, B. D., Mazloff, M. R., Miller, A. J., Ralph, F. M., Seo, H., and Hoteit, I.: The role of air–sea interactions in atmospheric rivers: Case studies using the SKRIPS regional coupled model, J. Geophys. Res.-Atmos., 126, e2020JD032885, https://doi.org/10.1029/2020JD032885, 2021. a, b
Sun, R., Villas Bôas, A. B., Subramanian, A. C., Cornuelle, B. D., Mazloff, M. R., Miller, A. J., Langodan, S., and Hoteit, I.: Focusing and defocusing of tropical cyclone generated waves by ocean current refraction, J. Geophys. Res.-Oceans, 127, e2021JC018112, https://doi.org/10.1029/2021JC018112, 2022. a, b, c
Tewari, M., Chen, F., Wang, W., Dudhia, J., LeMone, M., Mitchell, K., Ek, M., Gayno, G., Wegiel, J., and Cuenca, R.: Implementation and verification of the unified NOAH land surface model in the WRF model, in: 20th conference on weather analysis and forecasting/16th conference on numerical weather prediction, American Meteorological Society Seattle, WA, vol. 1115, 2165–2170, 2004. a
Thorpe, S.: Langmuir circulation, Annu. Rev. Fluid Mech., 36, 55–79, 2004. a
Tolman, H. L.: A third-generation model for wind waves on slowly varying, unsteady, and inhomogeneous depths and currents, J. Phys. Oceanogr., 21, 782–797, 1991. a
Tolman, H. L.: Subgrid modeling of moveable-bed bottom friction in wind wave models, Coast. Eng., 26, 57–75, 1995. a
Van Roekel, L., Fox-Kemper, B., Sullivan, P., Hamlington, P., and Haney, S.: The form and orientation of Langmuir cells for misaligned winds and waves, J. Geophys. Res.-Oceans, 117, https://doi.org/10.1029/2011JC007516, 2012. a, b, c, d
Warner, J. C., Sherwood, C. R., Signell, R. P., Harris, C. K., and Arango, H. G.: Development of a three-dimensional, regional, coupled wave, current, and sediment-transport model, Comput. Geosci., 34, 1284–1306, 2008. a
Weber, J. E. H., Broström, G., and Saetra, Ø.: Eulerian versus Lagrangian approaches to the wave-induced transport in the upper ocean, J. Phys. Oceanogr., 36, 2106–2118, 2006. a
Wu, L., Staneva, J., Breivik, Ø., Rutgersson, A., Nurser, A. G., Clementi, E., and Madec, G.: Wave effects on coastal upwelling and water level, Ocean Model., 140, 101405, https://doi.org/10.1016/j.ocemod.2019.101405, 2019b. a