Instabilities with SWAN_COUPLING and WET_DRY together

[drews]

I am applying wind stress at 20 m/s to Lake Erie at an angle of 20° north of due east (from the WSW, aligned with the axis of the lake). I get a blowup error at time 01:34:45 (hours), and upon further investigation I see unrealistic instabilities arising along the north shore, where the waves impact the cliffs (42.6 North, 80.7 West). The lake level is 174.1 meters above sea level, and these cliffs rise to 200 meters asl. The ROMS log shows the Courant number Cw rising above 5, which of course is not good.

ROMS/TOMS version 3.5
SVN Revision : 563M

Here is the final time step:
5685 0 01:34:45 3.279316E-02 6.458590E+06 6.458590E+06 4.933649E+11

(498,148,1) 5.903344E-04 1.435162E-03 5.985260E+00 2.005215E+01


Here is a snapshot of the free-surface zeta at the blowup time:

Figure 1.

Here is the lake level in a couple of cells along the northern shore near the cliffs:

Figure 2.

ErieUnstable02.jpg (242.49 KiB) Viewed 9581 times

My C++ #defines are:
#undef UV_ADV
#define UV_COR
#define UV_QDRAG
#define SOLVE3D
#define SPLINES
#define AVERAGES
#define ATM_PRESS
#define ANA_STFLUX
#define ANA_BTFLUX
#define WET_DRY
#define MASKING

#define GLS_MIXING
#define TKE_WAVEDISS

#define SWAN_COUPLING
#define NEARSHORE_MELLOR08

#ifdef SWAN_COUPLING
# define MCT_LIB
# define SVENDSEN_ROLLER
#endif

Other values for ROMS are:
1.000 dt Timestep size (s) for 3-D equations.
10 ndtfast Number of timesteps for 2-D equations between
Lm == 719 ! Number of I-direction INTERIOR RHO-points

Mm == 239 ! Number of J-direction INTERIOR RHO-points

N == 3 ! Number of vertical levels

DCRIT == 0.10d0 ! m


Minimum X-grid spacing, DXmin = 6.88508034E-01 km
Maximum X-grid spacing, DXmax = 6.88508034E-01 km
Minimum Y-grid spacing, DYmin = 9.26478945E-01 km
Maximum Y-grid spacing, DYmax = 9.26478945E-01 km
Minimum Z-grid spacing, DZmin = -3.33333001E+05 m
Maximum Z-grid spacing, DZmax = -3.69999999E+01 m

Minimum barotropic Courant Number = 5.97131259E-03
Maximum barotropic Courant Number = 8.96145073E-03
Maximum Coriolis Courant Number = 9.75860076E-05

Maximum grid stiffness ratios: rx0 = 2.271973E-01 (Beckmann and Haidvogel)
rx1 = 1.135987E+00 (Haney)

The model coupling is:
Ocean Export: bath:SSH:Ubar:Vbar:ZO
Waves Export: Wamp:Wlen:Wdir:Wptop:Wpbot:Wbrk:Wdiss:Wubot

Some SWAN parameters are:

SET level 174.1 &
depmin 0.1 &
rho 1000.0 &
inrhog 1
MODE DYNAMIC
CGRID REG 0 0 0 495726 222355 720 240 &
CIRcle 36 0.0418 1.0 24
INPgrid BOTtom REG 0 0 0 720 240 688.5 926.5 EXC 9999
INPgrid WInd REG 0 0 0 720 240 688.5 926.5 &

NONSTATIONARY 20061201.000000 1 HR 20061203.000000

READinp BOTtom -1 'bathymetry.bot' 3 2 FREE
READinp WInd 1 'winds.wnd' 3 2 1 FREE
BREaking
FRICtion
FRAME 'LakeErie' 0 0 0 495726 222355 720 240
OUTPut OPTIons BLOCK 8 1000
BLOCK 'LakeErie' NOHEADer 'ErieBlockWaves-hsign.txt' LAY 3 HSIGN &
OUTPUT 20061201.000000 1 MIN
COMPute (NONSTat 20061201.000000 30 SEC 20061203.000000)
STOP


The outputs from SWAN to ROMS are stable. The wave height hsign smoothly increases up to 2.1 meters before the blowup. Wave heights calculated by SWAN are routinely greater than the water depth in cells along the coast, which puzzles me (see Cutoff).

That's the Baseline case. Here is what I have tried so far:

NoWetDry: Turn off the ROMS WET_DRY option. Adjust water level to a sea-level datum. Mask off all water shallower than 2.1 meters to avoid cells going dry. This configuration results in a blowup error at time 07:11:24 hours in cell (352,198) when that cell went dry, so the model ran much longer. Along the northern shoreline the values for zeta are stable but not realistic; I get surge heights greater than 1 meter, with a level reduction of 2 meters in the adjacent cell. The basin also loses about 4% of its total volume. The vbar currents are about 2 meters/second up to the northern shore. Something really weird is happening when the waves meet the coast there. The wave ratio (height/water depth) climbs smoothly to 0.8 and 1.0 and beyond along that northern coast.

WET_DRY is very useful for studying storm surge. Waves contribute to storm surge. It would be valuable to get the combination working.

TinyWaves: Force SWAN with winds that are 1/1000 the strength of the baseline winds. Force ROMS with normal baseline winds. Result: The model runs with no instabilities up to about 12:50 hours, at which point it quits with a "File too large" error. The problem is that the SWAN output files exceed 2 Gigabytes. Nevertheless, I can conclude that the model will run further.

Cutoff: Permit no wave height larger than 0.8 of the water depth. Reduce the wave length by the same fraction. I accomplished this by adding custom code to ocean_coupler.f90 after retrieving the SWAN values. Result: Blowup error at time 01:35:11 hours in cell (495,147) at a sea cliff on the southern shore. The same large wavelike instabilities are present along the northern shore and at a few points along the southern shore. This result really surprised me, because I thought wave heights greater than the water depth were causing the problem.

Switching from UV_QDRAG to UV_LOGDRAG.
Result: Blowup error at time 01:00:14 hours in cell (393,195) at the wave impact zone along the north shore.

Thanks to John Warner for his help so far. Any ideas about what to try next?

[drews]

It looks like my combined model has a direction problem. I had set up SWAN for Cartesian coordinates, but ROMS requires the Nautical convention for wave direction, as John pointed out in this recent thread:

viewtopic.php?f=14&t=2301
For SWAN use:
SET INRHOG 1 NAUTICAL

The WikiROMS also says:
https://www.myroms.org/wiki/index.php/Variables#D
"Dwave: wind-induced wave direction. Direction the waves are coming from; measured clockwise from geographic North. (nautical convention)."

In my case the wave direction Dwave should be 250 degrees. (By the way, my coupling interval between ROMS and SWAN is 1 minute.) I'll make the directional revision to the SWAN forcing winds and run it again.

[drews]

Converting SWAN to use the Nautical convention for wind was not easy. I am forcing SWAN from a file of wind components (all constant for now). The wind is blowing from the West-Southwest along the axis of the lake (20 degrees north of due east), from Toledo to Buffalo. This wind direction at 20 m/s requires the following components in the SWAN forcing file:
x = 18.7939, y = 6.8404

SWAN behaves as if the "first guess" uses the Cartesian convention for wind direction, then all following calculations use the Nautical convention. So the initial waves are in the wrong direction. Over the first hour of simulation they swing around to come from Nautical 250 degrees (WSW). This would be a pretty big bug to have been overlooked by the SWAN developers, so I think something in my SWAN configuration is wrong instead. Nevertheless, the work-around is to suppress the initial guess in the SWAN configuration file:

INITial ZERO
GEN3 KOMEN AGROW

Now the waves start in the correct direction and don't take long to build. They are mostly directed at 250 degrees, and they curl outward along the north and south shores as they should:

Figure 3.

ErieUnstable03.jpg (73.62 KiB) Viewed 9495 times

Figure 3. SWAN-calculated waves on Lake Erie.

Back to the original problem: The shoreline instabilities are still present. The free-surface zeta looks similar to Figure 1, although the oscillations are slightly bigger. The blowup occurs at 01:27:18 hours in cell (495,147). Here is the plot of zeta for that cell (along the southern shore):

Figure 4.

ErieUnstable04.jpg (254.13 KiB) Viewed 9495 times

Figure 4. Instability in zeta along the coast.

The jump in the water level in this plot is small, only 0.07 meters (from 217.1 to 217.17 meters). But this cell is 43 meters above the lake surface (at 174 m asl). Why is it getting wet? Surely the water is not really splashing up that high.

[drews]

I tried running the Lake Erie model with DCRIT = 0.5 meters and time step DT = 0.5 seconds. The maximum Courant value Cw reported in the ROMS log file is 0.939, which removes the possibility of a CFL violation. The model terminated at 01:53:56 hours with a blowup error in cell (399,192). That cell has oscillations in zeta on the order of 2 meters in amplitude.

Then John Warner set me up with the COAWST version of ROMS and SWAN. When I ran the model, all the former problems went away! There are no more wavelike instabilities coming off the northern shore. A few grid cells go dry in Sandusky Bay. The storm surge piles up nicely at Buffalo like it should. John, you rock!

SWAN eventually got stuck at 06:25:59 hours when one of the coupling files exceeded 2 Gigabytes, but I know how to fix that. Hwave stabilized at about 4 meters in the center of the lake.

There must be some procedure for writing COAWST modeling improvements back into the ROMS code base. How do I request this?

[drews]

Here is what I did to reduce the size of the SWAN coupling/output files to keep them smaller than 2 Gigabytes:

1. Increase the coupling interval from 1 to 6 minutes. I can't remember exactly where I changed this run parameter.

2. The block output from SWAN has to be changed, too. That command in erie.swn has to reflect the 6-minute coupling interval as well:
BLOCK 'LakeErie' NOHEADer 'ErieBlockWaves-ubot.txt' LAY 3 UBOT &
OUTPUT 20061201.000000 6 MIN

[AlexisEspinosa]

Dear drews,

Do you know if ROMS incorporated the stability procedure from COAWST?
I'm curious about it.

I'm suffering some blow-ups myself with ROMS and I do not know if I should ask Prof. John Warner access to his ROMS-COAWST.

a) Do you know what was the main difference that made COAWST ran smoothly compared to the main ROMS version?

On the other hand, I've been trying to find something about vertical courant number (Cw) restrictions but there is not much in the formum (except your comment above). I'm also using wet-dry and, as you posted, to keep Cw<1 implies the use of very small time steps.

b) Does Cw should always be Cw<1?

I'm asking this because I have obtained good results compared to field measurements even when Cw~40 (Cu~0.05 Cv~0.05). But on the other hand, I recently made some changes to my bathymetry and now the model is blowing-up. New bathymetry is not that different from the previous one, although it incorporates new data, so it is more realistic. Then, the question is, if Cw>1 is the cause of the recent blow-ups, then:

c) why the model was not blowing up previously if the maximum Cw's are pretty much the same?


Many thanks,
Alexis Espinosa

[jcwarner]

i try my best to keep a version of roms in coawst that is consistent with the rutgers web site. however, we do make modifications in coawst that may take some time to get back to rutgers. the same is for wrf and swan - we make changes that work for the coupled system and sometimes these features make their way back to those root versions.

if you using swan coupling, then you should use coawst.
if you are using wet-dry, we have a few new options that seem to work better for very shallow water. let me know if you want access to coawst.
thanks,
john

[AlexisEspinosa]

Yes, of course. I'm very interested.

[CharlesJames]

Converting SWAN to use the Nautical convention for wind was not easy.

I was just browsing the Forum to see if anyone had dealt with this before although it looks like it hasn't been discussed for a while. I was just faced with the same problem after having developed a SWAN model with a large number of Cartesian convention forcing files and realizing that when coupling to COAWST (or ROMS) the hydrodynamic model is expecting Nautical convention Dwave from SWAN.

Before attempting this I had a look to see how often Dwave is used in COAWST and it turns out that there are only 12 files which use the Dwave variable and all of them actually transform it back into the Cartesian convention. I created a small cppdef wrapper called DWAVE_CART around each one that removed the '1.5_r8*pi-' transformation whenever it was implemented and it seems to have worked fine (and saves two floating point operations per modification!). Now my Stokes drifts are finally going in the right direction and I'm able to use the original SWAN configuration.

These are the files I modified although I've only tested wec_vf.F:

Nonlinear/mb_bbl.h
Nonlinear/sg_bbl.h
Nonlinear/ssw_bbl.h
Nonlinear/Sediment/sed_bedload.F
Nonlinear/Wec/radiation_stress.F_mellor05
Nonlinear/Wec/radiation_stress.F_mellor08
Nonlinear/Wec/radiation_stress.F_mellor11
Nonlinear/Wec/radiation_stress.F
Nonlinear/Wec/wec_roller.F
Nonlinear/Wec/wec_stokes.F
Nonlinear/Wec/wec_streaming.F
Nonlinear/Wec/wec_vf.F