unstable run with wave-enhanced bottom stress turn on

[Study_fish]

Dear ROMS Experts,

I am currently utilizing ROMS to simulate changes in water levels during the passage of a typhoon. In this simulation, I am using wave variable inputs as a forcing file rather than coupled with SWAN. My main purpose is to investigate the impact of wave-enhanced bottom stress on storm surge. Initially, I had planned to perform a nested simulation to obtain high-resolution results; however, the model frequently encountered instability during the initial step. Consequently, I have opted to begin with a simpler approach by simulating a single grid.

My CPP options are configured as follows:

#define ROMS_MODEL
#undef NESTING

#define UV_ADV
#define CURVGRID
#define UV_COR
#undef UV_LOGDRAG
#define LIMIT_BSTRESS
#define UV_VIS2

#define TS_U3HADVECTION
#define TS_C4VADVECTION
#define SALINITY
#define TS_DIF2
#define NONLIN_EOS

#define BULK_FLUXES
#define EMINUSP
#define ATM_PRESS
#define SOLAR_SOURCE

#define AVERAGES
#define DIAGNOSTICS_TS
#define DIAGNOSTICS_UV

#define SOLVE3D
#define SPLINES_VVISC
#define SPLINES_VDIFF
#define MASKING
#define SPHERICAL
#define DEBUGGING

#define MIX_GEO_UV

#define MIX_GEO_TS

#define LMD_MIXING
#define LMD_SKPP
#define LMD_BKPP
#define LMD_NONLOCAL
#define LMD_SHAPIRO
#define RI_SPLINES

#define WEC_VF
#define WET_DRY
#ifdef WEC_VF
# define WDISS_THORGUZA
# define WEC_BREAKING
# define ROLLER_RENIERS
# define WEC_ROLLER
# undef WEC_WCAP
#endif

#undef STREAMING
# ifdef STREAMING
# define SURFACE_STREAMING
# define BOTTOM_STREAMING
# define BOTTOM_STREAMING_YU
#endif

#define SSW_BBL
#ifdef SSW_BBL
# define SSW_CALC_ZNOT
# define SSW_CALC_UB
# undef ANA_SEDIMENT
# undef SSW_LOGINT
# undef SSW_ZORIP
#endif

#undef SSH_TIDES
#ifdef SSH_TIDES
# define UV_TIDES
# define ADD_FSOBC
# define ADD_M2OBC
#endif

#define HDF5

However, after running the model for 16 hours, it experienced a blow-up. I have several questions:

1. I have checked my initial and boundary files, and they appear to be normal. Since I am particularly concerned about flooding during storm surges, I enabled the wet_dry option and conducted a separate test without activating any other wave options, which ran without issues. Therefore, is the problem indeed related to the settings of ssw_bbl? I have defined #define SSW_CALC_UB because I am not coupling with the SWAN model and cannot provide bottom orbital velocity to ROMS. Additionally, SSW_BBL requires the user to specify the bottom wave period, and I have approximated this using the peak surface wave period. Could this be a potential issue?

2. Regarding the setting of the WEC_ALPHA parameter, I referred to the Sandy case study and set it to 0. Does this mean that it is unnecessary to define ROLLER_RENIERS and WEC_ROLLER? Could this be the cause of the blow-up?

3. I reviewed the run log and noticed something concerning: the side boundary condition for ubar_stokes appears as a string of garbled text. As I am not very familiar with Fortran, I find it difficult to understand the reason for this issue, even though I set the boundary condition for ubar_stokes to Gradient in the input file. I am worried that COAWST may not be correctly reading this boundary condition.

Any suggestions would be greatly appreciated. Thank you for your assistance.

[Study_fish]

Progress has been made in my running. Upon reevaluating the forcing files, I identified issues with my "Lwave" variable. Consequently, I performed iterative calculations of Lwave using the dispersion relation, employing the wave period at the spectral peak, which is publicly available data from Copernicus Marine Service. I created a new nested example; however, the model experienced a blow-up after running for one day.

I hypothesize that the issue may stem from my use of Pwave_top as a substitute for Pwave_bot. Given the presence of very deep topography within my model domain, this approximation may lead to significant inaccuracies in deep water regions. Could this be one of the factors contributing to the blow-up?

I reviewed the SWAN code for calculating Pwave_bot and noted that it requires the integration of the wave energy spectrum. Unfortunately, I do not have complete data on the wave energy density spectrum. I would appreciate any suggestions on how to accurately provide Pwave_bot.

Additionally, I examined the .rst file and found that the anomalous velocities occurred near the boundary of an island in the coastal region. As a rough corrective measure, would it be feasible to set the land mask for this island to represent ocean instead?

[lfparejar]

Just to confirm: you are using COAWST (e.g., https://github.com/DOI-USGS/COAWST) or ROMS (https://github.com/myroms)?

It's probably easier to make a SWAN domain based on your ROMS grid, figure out a SWAN input file, and run the models coupled - that way you would have all the wave variables you need. There are matlab tools to make these SWAN grid and bathymetry files.

What does your wave forcing file look like - what variables do you have and where do you get them from?

--
Regarding some of the wave variables:

-The bottom wave period "Pwave_bot" is ROMS' rewording of SWAN's "TMBOT" output variable.
-For your application, where you say there is deep topography, I wouldn't consider wave rollers (WEC_ROLLER and ROLLER_RENIERS) for now. I think of these as surfzone, very shallow water breaking processes. Same for WEC_ALPHA.
-I don't think the WEC_WCAP and WEC_BREAKING flags are actually coded right now, so no need.

I would also recommend activating wave options one by one to see what's causing the issue. Then it's easier to see if the issue is the vortex force formalism WEC_VF or the boundary layer scheme SSW_BBL (and associated flags).

Then again, making this a coupled model is probably better.

-
Fernando

[Study_fish]

Just to confirm: you are using COAWST (e.g., https://github.com/DOI-USGS/COAWST) or ROMS (https://github.com/myroms)?

Thanks to your reply! I'm using COAWST to conduct this running instead of ROMS, cuz VF approach is just implemented in COAWST.

It's probably easier to make a SWAN domain based on your ROMS grid, figure out a SWAN input file, and run the models coupled - that way you would have all the wave variables you need. There are matlab tools to make these SWAN grid and bathymetry files.

What does your wave forcing file look like - what variables do you have and where do you get them from?

I see that point. For certain reasons that cannot be elaborated upon here, I would like to run ROMS rather than the coupled model. My wave data is primarily obtained from https://data.marine.copernicus.eu/produ ... escription. For data that cannot be obtained (such as wave length and bottom wave period), I will use certain methods to approximate them. The wave length is calculated iteratively using the dispersion relationship of linear wave theory, and the bottom wave period is approximated using the peak wave period as a substitute.

Regarding some of the wave variables:

-The bottom wave period "Pwave_bot" is ROMS' rewording of SWAN's "TMBOT" output variable.
-For your application, where you say there is deep topography, I wouldn't consider wave rollers (WEC_ROLLER and ROLLER_RENIERS) for now. I think of these as surfzone, very shallow water breaking processes. Same for WEC_ALPHA.
-I don't think the WEC_WCAP and WEC_BREAKING flags are actually coded right now, so no need.

I would also recommend activating wave options one by one to see what's causing the issue. Then it's easier to see if the issue is the vortex force formalism WEC_VF or the boundary layer scheme SSW_BBL (and associated flags).

Then again, making this a coupled model is probably better.

These suggestions are very insightful. I realize that deploying wave rollers in deep terrain is not a good idea, and deep terrain can also lead to potential problems. Therefore, I modified my grid in subsequent trial simulations to avoid deep soundings as much as possible. Meanwhile, I also tried opening or closing different switches in combinations, with the closest I came to success being when the model exploded during a typhoon landfall. I will investigate the cause further and report back promptly.

[lfparejar]

In theory you can calculate Lwave based on the dispersion relation for surface gravity waves, but you would be missing a term of of k (dot) u for the Doppler effect because you have no currents. In SWAN, you may get bathymetric refraction (PROPTHETA) but not accurate frequency shift (PROPSIGMA), as currents are not received by SWAN.

I don't think the difference between Pwave_top and Pwave_bot is the biggest issue. For SSW_BBL, if I recall correctly, you can either define SSW_CALC_UB (which needs Pwave_bot) or not define it and provide Uwave_rms, which is the root-mean-square value (in m/s) of the orbital motion near the bottom (Ub). Ub and Uwave_rms are related like this: Ub = sqrt(2)*Uwave_rms. Ub here would be the bottom orbital velocity based on linear wave theory for a monochromatic wave.

Also keep in mind that using wave forcing instead of a dynamic wave model has some implications:

-Depth-dependent terms in the wave action density budget in SWAN (such as bottom friction, Sbf) will not see variations in sea surface elevation (zeta) computed by ROMS. SWAN will calculate these terms based on static 'h', not 'h + zeta'. If you have big tides or storm surges, these terms may not be accurate.

Without surface currents exported to SWAN, the wave energy will only be advected by the group velocity and not by the surface current.

[Study_fish]

I don't think the difference between Pwave_top and Pwave_bot is the biggest issue. For SSW_BBL, if I recall correctly, you can either define SSW_CALC_UB (which needs Pwave_bot) or not define it and provide Uwave_rms, which is the root-mean-square value (in m/s) of the orbital motion near the bottom (Ub). Ub and Uwave_rms are related like this: Ub = sqrt(2)*Uwave_rms. Ub here would be the bottom orbital velocity based on linear wave theory for a monochromatic wave.

useful information! Accroding to small-amplitude linear wave theory for a monochromatic wave, bottom orbital velocity is written as (Soulsby and Smallman, 1986):

Ub = wH/2sinh(kh)

so i can get the standard Uwave_rms for ROMS by using the relation Uwave_rms = Ub/sqrt(2), is it true? I will turn off the SSW_CALC_UB and try it.

Also keep in mind that using wave forcing instead of a dynamic wave model has some implications:

-Depth-dependent terms in the wave action density budget in SWAN (such as bottom friction, Sbf) will not see variations in sea surface elevation (zeta) computed by ROMS. SWAN will calculate these terms based on static 'h', not 'h + zeta'. If you have big tides or storm surges, these terms may not be accurate.

Without surface currents exported to SWAN, the wave energy will only be advected by the group velocity and not by the surface current.

That's the point. However, I'm not sure if the wave analysis products in CMEMS contain the surface currents.