Waves in SKRIPS: WaveWatch

. In this work, we integrated the WaveWatch (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58) WAVEWATCH III model into the regional coupled model SKRIPS (Scripps–KAUST Regional Integrated Prediction System). The WaveWatch (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58) WAVEWATCH (cid:58) III model is implemented with flexibility, meaning the coupled system can run with or without the wave component. (cid:58)(cid:58) In (cid:58)(cid:58)(cid:58)


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 and increase socio-economic implications for coastal communities in that region (Henderson-Sellers et al., 1998;Murakami et al., 2017;Bhatia et al., 2018).Cyclone Mekunu ::: We :::::::::: investigated :::::::: Cyclone ::::::: Mekunu ::::::: because :: it : was the strongest tropical cyclone in the north Indian Ocean in 2018 (Government of India, 2018).It had a clear signature in SST cooling and MLD deepening, which can be used for testing the model.We investigate :::::::::: 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 Section 2. An overview of cyclone Mekunu, the design of the experiments, and the validation data are presented in Section 3. Section 4 details the numerical simulation results and Section 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.

Model Description
In this work, the version 6.07 (Tolman, 1991;WW3DG, 2019) of the WAVE-height, WATer depth and Current Hindcasting third generation wave model (WaveWatch ::::::::::::: 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).
Following our previous implementations (Sun et al., 2019) and the WW3-ESMF interface (WW3DG, 2019) for the prototype case, we implemented the WW3-ESMF interface in the coupled modeling system.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, the surface boundary fields are exchanged via subroutine calls during the execution of the simulations by the WW3-ESMF interface.In addition, WW3 grid information is provided to the coupler in the initialization subroutine for online re-gridding of the exchanged boundary fields, although online regridding is not needed in the present work.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.

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): du dt ∂u ∂t ::: where t is time; ∇ = (∂x, ∂y, ∂z); ρ w is the density of water; p is the pressure; b = −gρ/ρ w ẑ = bẑ is the buoyancy term; ẑ is the vertical unit vector; D u ::: , although vector notation is used when it is unambiguous.The tracer advection equation can be written as (Suzuki and Fox-Kemper, 2016;Wu et al., 2019a): where c is a scalar quantity, such as potential temperature and salinity; D c is the diffusion.

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 0 : 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 discussed in ::: that :::::  2018).In addition to the extremely heavy rainfall in Oman and Yemen, 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).
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.The wave ::::::::::::::: ocean-atmosphere : model is not coupled to the ocean or the atmosphere in the simulations :::: wave ::::: model, aiming to demonstrate the impact of the wave model on the simulation results.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 assimilated HYCOM /NCODA 1/12 • daily global analysis ::::::: HYCOM : data as initial and boundary conditions for ocean temperature, salinity, and horizontal velocities (Chassignet et al., 2007).The :: At :::: each :::: time ::::: step, ::: the boundary conditions for the ocean are updated by linearly interpolating between the daily HYCOM /NCODA analysis 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 days, 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)  The time step of the ocean model is 120 seconds.The horizontal sub-grid mixing is parameterized using nonlinear Smagorinsky viscosities, and the K-profile parameterization (KPP) 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).

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(Knapp et al., , 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 SST is 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)(Lumpkin and Pazos, 2007).The simulated SST fields are also vali-  speed by about 5.9 m/s and 7.3 m/s, while in ATM.DYN the underestimation is as large as 13.5 m/s.The ATM.DYN does not capture the intensification of the TC between May 24 and 26 that is present in the IBTrACS observations and in the coupled simulations.We hypothesize that this is due to the SST cooling in ATM.DYN being greater than in the coupled simulations (see Section 4.2).Despite CPL.AOW better simulates the minimum pressure and maximum wind speed compared with the IBTrACS data than CPL.AO, the overall performance of CPL.AOW is not better but somewhat worse than CPL.AO, although both coupled runs better simulate the tropical cyclone than ATM.DYN.
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 further :::::::::: investigated ::: in :::::: Section :: 5. :

SST and Mixed Layer Depth
The :: 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 is presented in Fig. 5to highlight the impact of the tropical cyclone on the ocean.The results obtained from CPL.AO are not presented here, but in Section 5 we investigate the effect of wave coupling on the coupled system on SST and MLD .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).The SST cooling in HYCOM shown in : It

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 snapshot :::::::: snapshots : of the ensembleaveraged H s at 00 UTC May 26 is :: are : presented in Fig. 7(a).The :::: a-c).::: On :::: May ::: 26, ::: the : ensemble averaged H s is as high as Near the eye wall of the tropical cyclone, the standard deviation of H s from the ensembles is approximately 3 m, showing that the wave height is sensitive to the wind speed ::::: greater :::::::: variance ::: near ::: the :::: eye ::: wall : (Fig. 7b) :::: (d-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 the appendix :::::::: Appendix :: D.
Although the simulations with spectral nudging have smaller internal variability, they underestimate the intensity of the tropical cyclone in the simulations.

The Vertical Profiles
To investigate the SST warming and cooling in the wake zone , we have plotted :::   To test the coupled ocean-wave-atmosphere model, we performed a series of simulations of 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.Then we 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 the 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.

Figure 1 .
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.

Figure 2 .
Figure 2. Comparison between the tracks of 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 time and locations of the cyclone at specific times.

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

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

Figure 8 .
Figure 8. 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.

Figure 10 .
Figure 10.The snapshot ::::::: 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).

Figure 11 .
Figure 11.The snapshots of the vertical profiles at 00 UTC May 24 obtained in the simulations.Panel (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 mixed :::: KPP ::::::: boundary layer depth in MITgcmto illustrate the mixing layer.

Figure A1 .
Figure A1.The tracks of 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.
:::::::::::: WAVEWATCH : III model driven by ERA5 wind only; WAV.CUR indicates the simulation using stand-alone WaveWatch ::::::::::::: WAVEWATCH III model driven by ERA5 wind and HYCOM currents.It can be seen from Fig. A7 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 A2 .
Figure A2.The characteristics of 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 A3 .
Figure A3.The characteristics of 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 A7 .
Figure A7.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.
. To initialize the wave model, we allowed the wave field to spin-up for 19 days from May 01, 2018 and then we analyze the period from May 20, 2018.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 • E to 78 • E. The horizontal grid has 448×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 sigma 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 On ::: the :::: other ::::: hand, ::: we ::: did ::: not :::: spin ::

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