Problem with ROMS wetting/drying and turbulence model(s)

[lanerolle]

I have been struggling really hard trying to get my ROMS application for Cook Inlet, AK working. This application needs wetting/drying and is quite challenging because the tidal water elevation range is ~10m (~33ft) - so there is heavy and rapid wetting/drying going on.

Irrespective of what I do, the model blows-up even with a baroclinic time step of 1 s and a barotropic time step of 1/30 s (NDTFAST=30) and for this configuration, the CFL conditions will not be violated. The run I have been trying to do is a simpler version with M2 tides only and having a constant density (T = 15 C, S = 35 PSU always). I always use Dcrit of 0.1 m although I have tried both 0.05 m and also 0.2 m.

I find that if I do not use any vertical turbulence model (MY2.5, GLS, KPP, etc.) but instead using #define ANA_VMIX and employ a constant value of 0.002 m2/s for the vertical eddy-viscosity (as in the ESTUARY test case), the model runs successfully and also does the wetting/drying well without any problems!!!

I tried both the Quadratic and Logarithmic bottom stress formulations and the outcomes are the same.

Could someone please let me know whether I need to do anything special with the vertical turbulence models when there is heavy and rapid wetting/drying and/or whether they become invalid in certain situations? Are some of the ROMS turbulence models better suited for wetting/drying than others?

Thanks.

PS. I have applied ROMS successfully with wetting/drying in Chesapeake Bay (where the tidal range is ~2m) and everything works well both with the MY2.5 and GLS turbulence models. The Cook Inlet application configuration is not significantly different to this Chesapeake Bay configuration but of course is much more challenging in terms of wetting/drying.

[lanerolle]

Sorry, I should have elaborated on the nature of the model blow-up. The model blows-up due to very large velocities which develop at isolated horizontal locations on sigma level #1 (i.e. the bottom most sigma-level on the rho-grid in the vertical). The blow-up also manifests itself in the turbulent kinetic energy (tke) but not in the general length scale (gls). Furthermore, the large velocities are confined to the bottom most ~5 sigma levels and their magnitudes diminish when moving from level #1 to say level #5 and by the time, we reach level N(ng) (top) there is no sign of high velocities.

So it appears that when there is heavy and rapid wetting/drying, the turbulence models do something strange at (or near) sigma level #1.

[nganju]

what does your vertical stretching look like?

[jcwarner]

I had an issue with this exact same problem, and fixed it by not allowing the computation of the bottom stress to exceed a value that would make the velocity change sign. This is all related to the barotropic/clinic mode splitting and the fact that the bottom stress needs to be computed before barotropic mode. This similar issue was recently discussed by Sasha in Hawaii, but not for the wet/dry. I would like to give you a modification to set_vbc.F and see if it works for you. Can you email to me your version of ROMS/Nonlinear/set_vbc.F? I will make one modification, and then you try it and see if this helps.
-john
jcwarner@usgs.gov

[lanerolle]

I tried your modified code and the model runs beautifully without encountering any numerical stability issues (or blow-ups)!! The wetting-drying also appears to be happening nicely. Thanks very much for your help. I have a few questions and could you please let us know what you think?

(a) Why does your code modification work - in particular with the MY2.5 mixing scheme? Are you addressing a known shortcoming associated with vertical turbulence closures? Do you anticipate its need for the GLS models also?

(b) What is the physical reasoning behind the code modification - in particular the use of the SIGN function and the u*Hz/dt term? - Is this from boundary layer theory?

(c) Will your code modification be included in the future official releases of ROMS?

Thanks again.

[jcwarner]

well.... lets start the bidding at ..
ooops. I cant do that.
Thanks for your compliments. It is nice when we work hard on something and it helps people out. This was an issue identified by a few other people here at the usgs (Maitane, and Neil) and we worked hard to figure this out. It is a short coming of the way the bottom stress is imposed in the model (at least this is my understanding of it). I think someone like Sasha would offer a better explanation, and he touched on this in Hawaii during his presentation, but for reasons other than wetdry. I dont remember all the gory details, but what happens is that when a cell first becomes wet, especially if the bathy is very irregular with strong gradients, there can be a large local pressure gradient. The bottom stress term becomes very large. So strong that it can change the magnitude of the veolcity in the predictor step. This causes an instability in step3d_uv and is not recoverable in the corrector step. It is rare, but of course when it does happen and one little cell somewhere gets bad, the whole simulation stops. So my fix was to limit the bottom stress magnitude so that it can not change the sign of the veolcity in the predictor step. So in the predictor step:

... FC(i,0)=dt(ng)*bustr(i,j)
...
cff1=u(i,j,k,nstp)*0.5_r8*(Hz(i,j,k)+Hz(i-1,j,k))
cff2=FC(i,k)-FC(i,k-1)
u(i,j,k,nnew)=cff1+cff2

so this is basically u_new = u_old*dz + bustr*dt (m2/s)
I dont want to allow the bustr term to change the sign of u_new. This is unphysical. How can the extraction of energy make the current go the other way? imho, it can not. It can only slow it down to zero. So i took the liberty to limit the magnitude of bustr*dt so that it would be smaller than u_old*dz, and therefore not allow u_new to change sign:
bustr*dt < uold*dz
bustr < u_old*dz/dt
The few lines of code added to set_vbc just limit the bustr, for example:
# ifdef WET_DRY
cff2=0.75_r8*0.5_r8*(Hz(i-1,j,1)+Hz(i,j,1))
bustr(i,j)=SIGN(1.0_r8,bustr(i,j))*MIN(ABS(bustr(i,j)), &
& ABS(u(i,j,1,nrhs))*cff2/dt(ng))
# endif
The 0.75 is my arbitrary assigment limitation. Similar stuff for bvstr, and this is coded for all 3 bottom stress options (log, quad, lin).

I am not sure how this exactly affects the performance of the turbulence closure. If there are problems with bottom stress, that is used as a BC for tke so the closures would be affected.

As far as adding new capabilities to Rutgers version - this is always an issue. I can only provide strong guidance, but have no direct control over how things are distributed.

[kate]

I don't think there's any reason to limit this to the case of WET_DRY. I've had to put in similar code without WET_DRY for tides in Bristol Bay (in which my tides exceed those of the real world). I got help from Sasha on this after hearing his talk.

[arango]

Yes, Sasha brought this to our attention for applications with bottom boundary layers (BBL) where the bottom drag is an important term in the governing equations. If you look his presentation, you will see his statement that having an implicit bottom drag requires a substantial redesign of the kernel in our version of ROMS. Sasha gave us a nice lecture about ROMS algorithms and we allocated a lot of time for his presentation

The ROMS Workshops are very useful and we encourage users to attend when possible. The quality of the scientific presentation is impressive

Adding this remedy is in my to do list and will be available soon as a CPP option.

[rduran]

I dont want to allow the bustr term to change the sign of u_new. This is unphysical. How can the extraction of energy make the current go the other way?

Greg Holloway would argue that:

We consider eddies interacting with bottom topography. It is seen that traditional ‘‘eddy viscosity’’ and/or ‘‘topographic drag’’, which
would reduce large scale flows toward rest, are wrong. The second law of thermodynamics is violated; the ‘‘arrow of time’’ is running backwards! From statistical dynamics, approximate corrections are obtained, yielding a practical improvement to the fidelity of ocean models.

and further:

Use of the term ‘‘form drag’’ has been largely replaced by ‘‘form stress’’ or ‘‘topographic stress’’ to recognize that this force may not simply retard the mean flow but may also force the mean flow. Numerical experiments confirmed that, with mean wind forcing applied to U of sign U > 0, the response tended toward U < 0.

reference:
From Classical To Statistical Ocean Dynamics
http://www.springerlink.com/content/t31252tt68g9784l/