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.

1 Introduction

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 Methodology

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/\mathrm{4}\mathit{\pi }$ (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.

Figure 1 The schematic description of the SKRIPS regional coupled ocean–atmosphere model. The yellow block is the ESMF–NUOPC coupler; the white blocks are the ocean and atmosphere components; the red blocks are the implemented MITgcm–ESMF, WRF–ESMF, and WW3–ESMF interfaces.

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

The contribution of surface waves to the ocean momentum balance can be described by the wave-averaged momentum equations as follows (Suzuki and Fox-Kemper, 2016; Wu et al., 2019a):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{u}}{\partial t}}+(\mathit{u}\cdot \mathrm{\nabla })\mathit{u}& =-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}}}\mathrm{\nabla }p+{\mathbf{D}}^u+\mathit{b}-\mathit{f}\times \mathit{u}\\ & -\underset{\begin{array}{c}\text{Stokes}\\ \text{advection}\end{array}}{\underbrace{({\mathit{u}}^{\mathrm{S}}\cdot \mathrm{\nabla })\mathit{u}}}-\underset{\begin{array}{c}\text{Stokes}\\ \text{Coriolis}\end{array}}{\underbrace{\mathit{f}\times {\mathit{u}}^{\mathrm{S}}}}-\underset{\begin{array}{c}\text{Stokes}\\ \text{shear}\end{array}}{\underbrace{u_j^{\mathrm{L}}\mathrm{\nabla }{u}_j^{\mathrm{S}}}},\end{array}\end{array}$$

where t is time; $\mathrm{\nabla }=(\partial x,\partial y,\partial z)$; ρ_w is the density of water; p is the pressure; $\mathit{b}=-g\mathit{\rho }/{\mathit{\rho }}_{\mathrm{w}}\widehat{\mathit{z}}=b\widehat{\mathit{z}}$ is the buoyancy term; $\widehat{\mathit{z}}$ is the vertical unit vector; D^u is the diffusion; $\mathit{f}=f\widehat{\mathit{z}}$ is the Coriolis parameter; $-u_j^{\mathrm{L}}\mathrm{\nabla }{u}_j^{\mathrm{S}}$ is the Stokes shear force; $\mathit{u}=(u_{\mathrm{1}},u_{\mathrm{2}},u_{\mathrm{3}})=(u,v,w)$ is the wave-filtered Eulerian velocity (Eulerian velocity of the flow solved in MITgcm); u^S is the Stokes drift; and ${\mathit{u}}^{\mathrm{L}}=\mathit{u}+{\mathit{u}}^{\mathrm{S}}$ is the wave-filtered Lagrangian velocity. Here, the Einstein summation convention is used (e.g., $u_j^{\mathrm{L}}\mathrm{\nabla }{u}_j^{\mathrm{S}}=u_{\mathrm{1}}^{\mathrm{L}}\mathrm{\nabla }{u}_{\mathrm{1}}^{\mathrm{S}}+u_{\mathrm{2}}^{\mathrm{L}}\mathrm{\nabla }{u}_{\mathrm{2}}^{\mathrm{S}}+u_{\mathrm{3}}^{\mathrm{L}}\mathrm{\nabla }{u}_{\mathrm{3}}^{\mathrm{S}}$), although vector notation is used when it is unambiguous.

The tracer advection equation can be written as follows (Suzuki and Fox-Kemper, 2016; Wu et al., 2019a):

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial c}{\partial t}}+(\mathit{u}\cdot \mathrm{\nabla })c=-\underset{\begin{array}{c}\text{Stokes}\\ \text{advection}\end{array}}{\underbrace{({\mathit{u}}^{\mathrm{S}}\cdot \mathrm{\nabla })c}}+D^c,\end{array}$$

where c is a scalar quantity, such as potential temperature and salinity, and D^c 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 w_s is modified by multiplying with an enhancement factor ε as follows:

$$\begin{array}{}\text{(3)}& w_{\mathrm{s}}={\displaystyle \frac{\mathit{\kappa }{u}^{*}}{\mathit{\varphi }}}\mathit{\epsilon },\end{array}$$

where κ is the von Karman constant, u^* is the friction velocity, and ϕ is the stability function defined by Large et al. (1994). The enhancement factor is defined as follows:

$$\begin{array}{}\text{(4)}& \mathit{\epsilon }=\left|\cos \right(\mathit{\alpha }\left)\right|\sqrt{\mathrm{1}+(\mathrm{1.5}{\mathrm{La}}_{\mathrm{SLP}}{)}^{-\mathrm{2}}+(\mathrm{5.4}{\mathrm{La}}_{\mathrm{SLP}}{)}^{-\mathrm{2}}},\end{array}$$

where α is the angle between wind and Langmuir cells and La_SLP 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

$$\begin{array}{}\text{(5)}& {\mathit{\kappa }}_v=w_{\mathrm{s}}hG\left(\mathit{\sigma }\right),\end{array}$$

where h is the boundary layer depth, h is the KPP boundary layer depth, G(σ) is the shape function, and $\mathit{\sigma }=-z/h$ is the normalized depth. The boundary layer depth is determined based on the bulk Richardson number:

$$\begin{array}{}\text{(6)}& R{i}_{\mathrm{b}}\left(z\right)={\displaystyle \frac{(B_r-B(z\left)\right)\left|z\right|}{|{\mathit{u}}_r-\mathit{u}(z){|}^{\mathrm{2}}+V_t^{\mathrm{2}}(z)}},\end{array}$$

where B_r and u_r are the buoyancy and velocity averaged in the surface layer, B(z) and u(z) are the local buoyancy and local velocity, and V_t is the unresolved vertical shear proportional to w_s (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 Ri_cr=0.3.

Different from VR12-MA, LF17 parameterized the entrainment flux due to Langmuir turbulence by revising the bulk Richardson number:

$$\begin{array}{}\text{(7)}& R{i}_{\mathrm{b}}\left(z\right)={\displaystyle \frac{(B_r-B(z\left)\right)\left|z\right|}{|{\mathit{u}}_r-\mathit{u}(z){|}^{\mathrm{2}}+V_{\mathrm{tL}}^{\mathrm{2}}(z)}},\end{array}$$

the turbulence shear term $V_{\mathrm{tL}}^{\mathrm{2}}$ is defined as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{V}_{\mathrm{tL}}^{\mathrm{2}}\left(z\right)& ={\displaystyle \frac{C_{v}N\left(z\right)w_{\mathrm{s}}\left(z\right)\left|z\right|}{R{i}_{\mathrm{c}}}}\\ & {\left[{\displaystyle \frac{\mathrm{0.15}{w}_{*}^{\mathrm{3}}+\mathrm{0.17}{u}_{*}^{\mathrm{3}}(\mathrm{1}+\mathrm{0.49}{{\mathrm{La}}_{\mathrm{SL}}}^{-\mathrm{2}})}{w_{\mathrm{s}}(z{)}^{\mathrm{3}}}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}\end{array}$$

where dimensionless coefficient $C_v=\max (\mathrm{2.1}-\mathrm{200}\times \max (\mathrm{0},N),\mathrm{1.7})$, N(z) is the local buoyancy frequency, Ri_c is the critical Richardson number, $w_{*}=(-B_{\mathrm{0}}/h{)}^{\mathrm{1}/\mathrm{3}}$ is the convective velocity scale, and B₀ is the surface buoyancy flux. Here La_SL 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)

$$\begin{array}{}\text{(9)}& {\mathit{\tau }}_{\mathrm{oc}}={\mathit{\tau }}_{\mathrm{a}}-{\mathit{\tau }}_{\mathrm{aw}}+{\mathit{\tau }}_{\mathrm{ow}},\end{array}$$

where

$$\begin{array}{}\text{(10)}& {\mathit{\tau }}_{\mathrm{aw}}={\mathit{\rho }}_{\mathrm{w}}g\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathit{k}}{\mathit{\omega }}}{S}_{\mathrm{in}}\mathrm{d}\mathit{\omega }\mathrm{d}\mathit{\theta },\end{array}$$

and

$$\begin{array}{}\text{(11)}& {\mathit{\tau }}_{\mathrm{ow}}={\mathit{\rho }}_{\mathrm{w}}g\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathit{k}}{\mathit{\omega }}}{S}_{\mathrm{ds}}\mathrm{d}\mathit{\omega }\mathrm{d}\mathit{\theta }.\end{array}$$

Here, S_in and S_ds 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 z₀ 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 z₀ 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 Experimental design

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⁻¹, gusting to 200 km h⁻¹ (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.

1. CPL.AOW. Coupled ocean–wave–atmosphere (MITgcm–WW3–WRF) simulations.

2. 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.

3. 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 $\mathrm{1}/\mathrm{12}{}^{\circ }$ 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.

Table 1 The computational domain, WRF physics schemes, initial condition, boundary condition, and forcing terms used in the present simulations.

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).

4 Comparing coupled and uncoupled models

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⁻¹; IBTrACS: 254 km d⁻¹), 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.

Figure 2 Comparison between the tracks of Tropical Cyclone Mekunu obtained from IBTrACS data and the simulations. The thick solid lines indicate the tropical cyclone tracks obtained from averaging all ensemble members; the thin solid lines indicate the tropical cyclone tracks obtained from individual ensemble members. The text and markers highlight the location of the cyclone at specific times.

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⁻¹ for CPL.AOW, CPL.AO, and ATM.DYN (ensemble standard deviations are 3.13, 3.00, and 2.80 m s⁻¹). 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⁻¹, while in ATM.DYN the underestimation is as large as 14.1 m s⁻¹. 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.

Figure 3 The characteristics of Tropical Cyclone Mekunu obtained from the simulations plotted as functions of time. The solid lines indicate the ensemble-averaged simulation results; the shaded areas indicate the standard deviation of the results. Panel (a) shows the distance errors in comparison with IBTrACS data. Panel (b) shows the cyclone central pressure. Panel (c) shows the maximum wind speed. The simulations start on 20 May, but the results are not presented before the pressure starts to drop on 22 May.

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⁻¹) 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⁻²). 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.

Figure 4 The snapshots of the wind speed and latent heat fluxes at 00:00 UTC on 25 May obtained in the simulations. The ensemble-averaged fields are plotted, and we highlighted the region between 7 and 22^∘ N and between 46 and 62^∘ E. In panels (a)–(c) the 15 m s⁻¹ contour of wind speed is used to highlight the size of the tropical cyclone. The solid black lines indicate the ensemble-averaged track of the tropical cyclone in the simulations shown in Fig. 2. The red dots indicate the ensemble-averaged locations of the center of the tropical cyclone at the snapshot; the black dots indicate the ensemble-averaged locations of the center of the tropical cyclone each day at 00:00 UTC.

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⁻³ 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.

Figure 5 The evolution of SST during the tropical cyclone event. Panel (a) shows the SST when the simulation is initialized at 00:00 UTC on 20 May. Panel (b) illustrates the ensemble-averaged SST changes between 00:00 UTC on 20 May and 00:00 UTC on 27 May in CPL. Panel (c) shows the SST changes in the HYCOM analysis for the same period.

Figure 6 The evolution of MLD during the tropical cyclone event. Panel (a) shows the MLD at 00:00 UTC on 20 May when the simulation is initialized. Panel (b) shows ensemble-averaged MLD changes between 00:00 UTC on 20 May and 00:00 UTC on 27 May in CPL. Panel (c) shows the MLD changes in the HYCOM analysis for the same period.

4.3 Waves

Ocean surface waves are expected to affect air–sea interactions. The ensemble mean significant wave height (H_s) and the ensemble standard deviation obtained from the coupled simulation are shown in Fig. 7. The snapshots of the ensemble-averaged H_s are presented in Fig. 7a–c. On 26 May, the ensemble-averaged H_s 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 H_s 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 H_s from the ensembles is approximately 3 m, showing greater variance near the eye wall (Fig. 7d–f), while the spatial variability of H_s, with alternating high and low beams, is consistent throughout the ensemble. Although CPL.AOW captures the overall spatial variability of H_s, 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 H_s with altimeter data is shown in Appendix D.

Figure 7 Snapshots of the significant wave height H_s at 00:00 UTC on 22, 24, and 26 May 2018. Panels (a)–(c) show the ensemble-averaged H_s obtained from CPL.AOW. Panels (d)–(f) show the standard deviation of H_s of the ensembles from CPL.AOW. The 15 m s⁻¹ contour of wind speed is used to highlight the location of the tropical cyclone.

5 Sensitivity analysis to wave coupling

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).

Figure 8 The snapshots of the ensemble-averaged SST difference. Panels (a)–(c) show the SST difference between the simulations with Langmuir turbulence (CPL.LF17, CPL.VR12-MA, and CPL.LF17-ST) and without Langmuir turbulence (CPL.NoLT). Panels (d)–(f) show the same differences in the simulations with spectral nudging. The markers indicate the regions where the SST difference is significant (P<0.05).

Figure 9 The snapshots of the ensemble-averaged MLD difference. Panels (a)–(c) show the MLD difference between the simulations with Langmuir turbulence (CPL.LF17, CPL.VR12-MA, and CPL.LF17-ST) and without Langmuir turbulence (CPL.NoLT). Panels (d)–(f) show the same differences in the simulations with spectral nudging. The markers indicate the regions where the MLD difference is significant (P<0.05).

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 Ri_cr 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 $|{\mathit{u}}_r-\mathit{u}(z){|}^{\mathrm{2}}$ 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 u^L 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 V_tL(z) term is used in Eq. (7) for parameterizing the Richardson number. Although the velocity gradient $|{\mathit{u}}_r-\mathit{u}(z){|}^{\mathrm{2}}$ is also smaller, shown in Fig. 10d, the entrainment flux V_tl(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.

Figure 10 The snapshots of the vertical profiles at 00:00 UTC on 25 May obtained in the simulations. Panels (a)–(d) show the Richardson number, buoyancy term, density changes, and horizontal velocity. The solid lines indicate the vertical profiles obtained from CPL.NoLT, CPL.VR12-MA, and CPL.LF17. The horizontal dashed lines indicate the KPP boundary layer in MITgcm.

6 Discussion and summary

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.

Appendix A: Langmuir numbers

In this work, we used the Langmuir number summarized in Li et al. (2019). The Langmuir number La_LP is defined from the surface layer (a fraction of the mixed layer, which is the upper 20 % in their definition) averaged Stokes drift 〈u^S〉_SL and a reference Stokes drift $u_{\mathrm{ref}}^{\mathrm{S}}$ near the base of the mixed layer:

$$\begin{array}{}\text{(A1)}& {\mathrm{La}}_{\mathrm{SL}}={\left[u_{*}/(\langle u^{\mathrm{S}}{\rangle }_{\mathrm{SL}}-u_{\mathrm{ref}}^{\mathrm{S}})\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

The Langmuir number La_LP is used in the scaling of Langmuir-enhanced entrainment in LF17 in Eq. (8).

The projected Langmuir number La_SLP considers the misalignment between wind and waves:

$$\begin{array}{}\text{(A2)}& {\mathrm{La}}_{\mathrm{SLP}}={\left[{\displaystyle \frac{u_{*}\cos {\mathit{\theta }}_{\mathrm{wl}}}{\langle u^{\mathrm{S}}{\rangle }_{\mathrm{SL}}\cos ({\mathit{\theta }}_{\mathrm{ww}}-{\mathit{\theta }}_{\mathrm{wl}})}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

where θ_ww is the misalignment between wind and waves and θ_wl is the misalignment between wind and Langmuir cells. The Langmuir number La_SLP is used in VR12-MA, and the scaling of Langmuir-enhanced diffusivity is used in LF17 in Eq. (4).

Appendix B: Ocean roughness closures

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.

In Taylor and Yelland (2001), the ocean surface roughness is parameterized as follows:

$$\begin{array}{}\text{(B1)}& {\displaystyle \frac{z_{\mathrm{0}}}{H_{\mathrm{s}}}}=\mathrm{1200}(H_{\mathrm{s}}/L_{\mathrm{p}}{)}^{\mathrm{4.5}},\end{array}$$

where z₀ is the ocean roughness, L_p is the wavelength at the peak of the wave spectrum, and H_s is the significant wave height.

Drennan et al. (2003) proposed a wave-age-based formula to characterize the ocean roughness. The wind friction velocity is also considered in this formula:

$$\begin{array}{}\text{(B2)}& {\displaystyle \frac{z_{\mathrm{0}}}{H_{\mathrm{s}}}}=\mathrm{3.35}(u^{*}/C_{\mathrm{p}}{)}^{\mathrm{3.4}},\end{array}$$

where u^* is the wind friction velocity and C_p is the wave phase speed at the peak frequency.

Oost et al. (2002) also derived the following expression for the ocean roughness based on the experimental data based on wave age and wavelength:

$$\begin{array}{}\text{(B3)}& {\displaystyle \frac{z_{\mathrm{0}}}{L_{\mathrm{p}}}}={\displaystyle \frac{\mathrm{25}}{\mathit{\pi }}}(u^{*}/C_{\mathrm{p}}{)}^{\mathrm{4.5}}.\end{array}$$

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

Appendix C: Tropical cyclone characteristics using different setups

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).

Figure C1 The tracks of Tropical Cyclone Mekunu from the coupled simulations using different options to parameterize (a) Langmuir turbulence and (b) surface roughness. The thick solid lines indicate the locations of the center of the tropical cyclone obtained from averaging all ensemble members. The thin solid lines in the background denote the tracks of each ensemble member. The dashed lines denote the track of the tropical cyclone in IBTrACS data. The text and markers highlight the time and locations of the cyclone at specific times.

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.

Figure C2 The characteristics of Tropical Cyclone Mekunu obtained from the coupled simulations using different options to parameterize Langmuir turbulence. The solid lines indicate the ensemble-averaged simulation results; the shaded areas indicate the standard deviation of the results. Panel (a) shows the distance errors in comparison with IBTrACS data. Panel (b) shows the cyclone central pressure. Panel (c) shows the maximum wind speed.

Figure C3 The characteristics of Tropical Cyclone Mekunu obtained from the coupled simulations using different options to parameterize surface roughness. The solid lines indicate the ensemble-averaged simulation results; the shaded areas indicate the standard deviation of the results. Panel (a) shows the distance errors in comparison with IBTrACS data. Panel (b) shows the cyclone central pressure. Panel (c) shows the maximum wind speed.

Figure C4 The characteristics of Tropical Cyclone Mekunu obtained from the coupled simulations that do not consider the effect of Stokes forces (Stokes advection and Stokes Coriolis) and wind stress corrections. The solid lines indicate the ensemble-averaged simulation results; the shaded areas indicate the standard deviation of the results. Panel (a) shows the distance errors in comparison with IBTrACS data. Panel (b) shows the cyclone central pressure. Panel (c) shows the maximum wind speed.

Figure C5 The snapshots of the ensemble-averaged difference in 10 m wind speed and latent heat flux. Panels (a)–(c) show the 10 m wind speed difference between the simulations using different options to parameterize the surface roughness (TY2001, DGHQ, and OOST) compared with CPL.AOW (CHNK). Panels (d)–(f) show the latent heat differences for the same simulations. The markers indicate the regions where the differences are significant (P<0.05).

Figure C6 The snapshots of the ensemble-averaged SST difference when Stokes forces are considered implicitly in the simulations. Panels (a)–(c) show the SST difference between the simulations with Langmuir turbulence (CPL.LF17, CPL.VR12-MA, and CPL.LF17-ST) and without Langmuir turbulence (CPL.NoLT). Panels (d)–(f) show the same differences in the simulations with spectral nudging. The markers indicate the regions where the SST difference is significant (P<0.05).

Appendix D: Impact of velocity shear in KPP

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.

Figure D1 The snapshots of the ensemble-averaged SST and MLD difference. Panels (a)–(b) show the SST difference between the simulations with Langmuir turbulence (CPL.LF17, CPL.VR12-MA) and without Langmuir turbulence (CPL.NoLT). Panel (c) shows the SST difference using VR12-MA but does not enhance the diffusion coefficient for current velocity. Panels (d)–(f) show the differences in MLD. The markers indicate the regions where the SST difference is significant (P<0.05).

Appendix E: The significant wave height

To evaluate the simulation performance of surface waves, we compared the modeled H_s with along-track H_s measurements from the Jason‐3 and SARAL/AltiKa altimeters. We use quality‐controlled, unfiltered, and not resampled along‐track H_s 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 H_s 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 H_s are not completely consistent with the observational data.

Figure E1 The ensemble-averaged significant wave height in comparison with the altimeter data. Panels (a)–(c) show the comparison with Jason-3 data. Panels (d)–(f) show the comparison with SARAL data. The simulation results from Sun et al. (2022) are also presented. The shaded areas indicate the standard deviation of wave height in the ensemble simulations.

Appendix F: Validation against drifter 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.

Figure F1 Evolution of the SST during the tropical cyclone event in comparison with drifter data. Panel (a) shows the SST changes during the event from drifter data. Panels (b)–(d) show the ensemble-averaged SST changes from HYCOM, CPL.LF17, and CPL.NoLT, respectively. The red numbers indicate SST warming during the event; the blue numbers indicate SST cooling during the event.