wetting and drying problem

[wpgong]

i am simulating an estuary dynamics using ROMS3.7 and activated wetting_and_drying. i tested different DCRIT depths from 0.1 to 1.0m. The model can run smoothly during flood tides. however, during ebb tides when some cells become dry, very large barotropic velocities occurred around these cells and the model blowed up.
i am wondering if anyone has experienced similar problems before and how to solve it?
thanks

Wenping 10/11

[kate]

Hernan told me he was working on a fix for WET_DRY - back in August. Any news, Hernan?

It is working for me in a Cook Inlet domain using the COAWST ROMS branch.

[wpgong]

The simulated domain in my situation is much complex, very irregular coastline, and the grid seems a little zigzag. does this affect the performance of wettting_and_drying process?
Hi, Kate, what are the values of Dcrit and time step in your Cook Inlet case?
thanks.

[lanerolle]

Have you tried using the LIMIT_BSTRESS CPP option when you have wetting-drying? I need it and without it, my runs are unstable and I use Dcrit=0.1 (the default value). Also, I use the wet/dry mask to mask out the net heat flux computed in bulk_flux.F.

I am having a problem where my wetting-drying simulation in Cook Inlet blows-up at isolated cells from time to time and the blow-up is only associated with T and S values (not water elevations or currents) and occurs during the summertime heating and not during the winter/spring/fall times! If I physically modify the bathymetry around the cell where the blow-up occurs, I can get the model running a little further but I need to keep doing this incessantly. When I look at where the unphysically large T, S values occur in the water column, it occurs in the middle where the grid stretching is most intense and the values appear to be reasonable near the top and bottom of the water column. So it appears that it might be some vertical advection issue or some issue like what Kate mentioned where T/S values at dry cells are getting updated somehow and we need to prevent this from happening. Any thoughts?

[kate]

I too use a DCRIT of 0.1. I tried 0.05 and it blew up. My timestep is 30 sec and I always run with LIMIT_BSTRESS.

For the T,S, I was getting blowups from the rotated mixing tensor. I hacked in a fix so that wet and dry cells don't diffuse with each other and that ran. Then I changed the bathymetry and I was back to having this blow up. I turned off explicit horizontal diffusion and now all is happy.

[wpgong]

thanks for you guys' reply.
i tried setting the CPP LIMIT_BSTRESS, and set the horizonatal viscosity and diffusivity as zero. The problem was still there. i will try more.

Wenping 10/13

[arango]

Hi Kate, I think that my big update to the wetting and drying will be released this week and the restart problem is now solved I am still running NOAA's Cook Inlet application but very slowly The grid is huge (1044x724x30) and I am running in my 8-CPUs desktop in case that I need to examine problem in the TotalView debugger. The problem is that this grid is so weird and it is very difficult to make sense inside the debugger when plotting or displaying tiled arrays. It is a nightmare and that's the reason that it is taking so long:



Anyway, I have been running this application for several days and restarting every 12 hours. Everything seems to go very smoothly and the model is not blowing-up with PERFECT_RESTART:

Code: Select all

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

      0     0 00:00:00  0.000000E+00  7.600860E+02  7.600860E+02  6.396795E+12
         (000,0000,00)  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00
      DEF_HIS   - creating history file, Grid 01: cook_inlet_his_0001.nc
      WRT_HIS   - wrote history  fields (Index=1,1) into time record = 0000001
      DEF_STATION - creating stations file, Grid 01: cook_inlet_sta.nc
      1     0 00:00:05  8.664375E-16  7.600860E+02  7.600860E+02  6.396795E+12
         (386,0113,01)  1.100680E-09  7.030073E-09  0.000000E+00  2.549995E-06
      2     0 00:00:10  2.402784E-10  7.600860E+02  7.600860E+02  6.396795E+12
         (181,0640,01)  5.215735E-04  5.364590E-05  2.222223E-11  2.209626E-02
      3     0 00:00:15  8.536992E-12  7.600860E+02  7.600860E+02  6.396795E+12
         (181,0640,30)  0.000000E+00  0.000000E+00  2.634209E+01  1.285078E-02
....

 NLM: GET_STATE - Read state initial conditions,             t =     0 12:00:00
                   (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1)

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

   8640     0 12:00:00  6.990734E-05  7.599744E+02  7.599745E+02  6.392954E+12
         (177,0823,30)  1.022442E-02  6.039011E-03  0.000000E+00  2.274181E-01
      DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc
   8641     0 12:00:05  6.995225E-05  7.599746E+02  7.599746E+02  6.392955E+12
         (178,0637,30)  0.000000E+00  0.000000E+00  9.772995E+00  2.270133E-01
   8642     0 12:00:10  6.999734E-05  7.599747E+02  7.599748E+02  6.392957E+12
         (228,0667,30)  2.714938E-06  0.000000E+00  1.091613E+01  2.266488E-01
   8643     0 12:00:15  7.004253E-05  7.599748E+02  7.599749E+02  6.392958E+12
...

 NLM: GET_STATE - Read state initial conditions,             t =     1 00:00:00
                   (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1)

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

  17280     1 00:00:00  1.628109E-03  7.595062E+02  7.595078E+02  6.385508E+12
         (172,0592,30)  3.321712E-02  1.269262E-02  0.000000E+00  7.566996E-01
      DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc
  17281     1 00:00:05  1.628419E-03  7.595064E+02  7.595080E+02  6.385506E+12
         (177,0636,30)  0.000000E+00  0.000000E+00  3.088772E+01  7.568594E-01
  17282     1 00:00:10  1.628729E-03  7.595065E+02  7.595081E+02  6.385503E+12
         (170,0619,30)  0.000000E+00  0.000000E+00  1.418715E+01  7.570190E-01
...

 NLM: GET_STATE - Read state initial conditions,             t =     1 12:00:00
                   (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1)

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

  25920     1 12:00:00  7.244958E-04  7.594712E+02  7.594719E+02  6.389436E+12
         (616,0511,30)  2.413032E-02  1.226487E-01  0.000000E+00  2.061889E+00
      DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc
  25921     1 12:00:05  7.252824E-04  7.594714E+02  7.594722E+02  6.389437E+12
         (179,0638,30)  1.161984E-03  3.170779E-06  8.484577E+00  2.060438E+00
  25922     1 12:00:10  7.260429E-04  7.594717E+02  7.594724E+02  6.389439E+12
         (207,0968,30)  6.414001E-03  0.000000E+00  9.312273E+00  2.059241E+00
...

 NLM: GET_STATE - Read state initial conditions,             t =     2 00:00:00
                   (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1)

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

  34560     2 00:00:00  5.806858E-03  7.586865E+02  7.586923E+02  6.375351E+12
         (617,0507,30)  6.565633E-03  1.592644E-01  0.000000E+00  2.073913E+00
      DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc
  34561     2 00:00:05  5.808851E-03  7.586864E+02  7.586922E+02  6.375339E+12
         (212,0974,30)  1.397617E-02  4.390969E-02  1.038334E+01  2.073022E+00
  34562     2 00:00:10  5.810832E-03  7.586862E+02  7.586920E+02  6.375327E+12
         (211,0975,30)  2.637443E-02  4.564990E-02  2.213804E+01  2.072133E+00
...
and so on

As you can see, the model is going stable restarting with PERFECT_RESTART. I think that the changes solved the problem that Lyon was having restating this application when WET_DRY was activated. I will release the update this week.

[kate]

Great! I'm clearly not having enough fun with grids.

[wpgong]

i checked the salinity and temperature data when the cells become dry, they are unreasonably large. At which place in the source code can we keep the salinity and temperature unchanged,i.e. as the values at the previous time step, when the cells become dry?
thanks.

[stvyng]

In Hernan's post,Cw occasionally appears to be rather large (i.e., > 1.0), while Cu/Cu are in the normal range. Would the large Cw cause any stability issue?

[nacholibre]

I am having a problem where my wetting-drying simulation in Cook Inlet blows-up at isolated cells from time to time and the blow-up is only associated with T and S values (not water elevations or currents) and occurs during the summertime heating and not during the winter/spring/fall times!

Did you try your simulation without atmospheric forcing? I have a case where I see unreasonably large T/S values even if I turn off radiation, solar source and rain. I believe the issue is occurring because of large slopes between the neighboring cells and it does not necessarily have to be in drying cells. When the water gets too shallow that the ratio of the depth in that cell with respect to its surrounding cells is greatly violating the gentle slope assumption (rx0 and more so rx1 ratios) tracers become unstable. For example, what happens if two neighboring cells have 2m and 3m depth initially, but during the simulation water level goes down to -1.9m in both and now it is 0.1m and 2.1m?

So far topographic smoothing, reducing number of layers, changing stretching of layers seem to be the easiest solutions. But they are not always a viable options as one might be diverging from the real physics of the problem. There is a relevant discussion going on under WET_DRY and Dcrit, although I think the issue is not only limited to drying cells.
Zafer

[kate]

i checked the salinity and temperature data when the cells become dry, they are unreasonably large. At which place in the source code can we keep the salinity and temperature unchanged,i.e. as the values at the previous time step, when the cells become dry?
thanks.

I tried doing this and it doesn't work. Things go bad when cells are barely wet, worse than if you let tracers evolve all the time.

[kate]

In Hernan's post,Cw occasionally appears to be rather large (i.e., > 1.0), while Cu/Cu are in the normal range. Would the large Cw cause any stability issue?

I get that too and I pretend all is well. That alone isn't enough to make ROMS blow up.

[lanerolle]

Hi Kate, Zafer and Wpgong,

Thanks Kate and Zafer for sharing your experiences with us about the extreme T/S values. As you know, I too am battling this problem currently and the most annoying thing this is that they occur at isolated (single cell/node) locations and then the whole simulation blows-up and crashes! According to what Zafer said, he found that (1) smoothing the bathymetry, (2) carefully choosing the vertical sigma-grid formulation and (3) the number of vertical levels had the greatest impact towards solving/alleviating this issue. I am pretty much settled with (2) and (3) and therefore, I only have (1) to play with in order to get my simulations to be robust and stable.

Kate and Zafer, could you please share some of your experiences/advice on how you stabilized your simulations? -e.g. what are the general rules of thumb you followed in fixing/smoothing the bathymetry? I will try out some of your techniques and see whether I too can stabilize my simulations.

On my part, I found that moving away from spline-based vertical advection gave me relatively greater numerical stability in my simulations. I now use C4 vertical advection for both momentum and tracers and U3 horizontal advection for both momentum and tracers.

Also, I am still puzzled as to why I get this T/S blow-up problem only in the warmer months. I begin my simulation on January 01 and there are no T/S related instabilities until ~May 15 and the tidal and sub-tidal forcing is pretty similar month-to-month. So the only other factor that I can think of is summer-time heating due to meteorological forcing. I have run several different grids in Cook Inlet, AK and in all cases, the instabilities occur in the summer-time. I use a constant Jerlov water type of 1 - do you also use this value or do you (as Kate mentioned to me in the past) use a spatially varying Jerlov number? Please let me know your thoughts on this.

Thanks,

Lyon.

PS. My model uses a DEM and hence has flooding over the mudflats in upper Cook Inlet. Do you think I need to do something special in my ROMS configuration/simulations to accommodate flooding?

[nacholibre]

My experience is that some sort of smoothing is usually required. You can tell it by looking at Beckmann & Haidvogel and Haney stiffness ratios (rx0 and rx1). The shallower the water, the more likely for the ratios to be exceeded more often. There are many options when it comes to smoothing, bin averaging, Shapiro filter, Mathieu's (Mathieu Dutour-Sikiric) smoothing toolbox for Matlab etc.

The reason that the blow up occurs in warmer months is probably that the extreme temperatures that you observe due to instability are enhanced when the heat input from the atmosphere is added on top. In other words the instability may be already there without the atmospheric forcing. Have you tried running without the atmospheric forcing to see if you still get unreasonably high tracer values?
Zafer

[kate]

Lyon - I did a lot of hand-tuning of my Cook Inlet bathymetry. My smoothing approach breaks down as depth approaches zero, so I told it to not smooth shallower than about 5 m. I then had to fix up areas where it blew up. I got something that ran with my old bathymetry and with the one I got from you, but blew up with an average of the two. A recent hack was to turn off explicit horizontal diffusion because the rotated mixing tensors were unstable at just one point.

[lanerolle]

Hi Zafer,

If I run my simulation without any surface/met forcing while also keeping T, S always constant (i.e. a constant density tidal simulation with or without sub-tidal forcing), the simulation is always stable and does not blow-up. I had to fix the bathymetry slightly to prevent blow-ups due to u-cfl, v-cfl violations as a result of excessively high velocities - but no blow-ups due to T/S values were seen anywhere. The T/S blow-ups came into existence after I ran the full synoptic hindcast simulation with surface/met forcing and during the warmer (May-September) months.

As the blow-ups in T/S occur in the shallows (< 5m), just like Kate, I too will aim to fix the bathymetry manually by hand.

Lyon.

[wilkin]

Is the wet/dry mask applied to the surface heat fluxes?

It would seem sensible that a cell designated "dry" is unable to heat or cool since there isn't really any water there. Surely disabling the heat flux would help in this situation?

[lanerolle]

Yes, I do apply the wet/dry mask to the heat flux in bulk_flux.F after it is computed. I assumed that would be the most likely place to apply the wet/dry mask although I have seen posts on the ROMS forum indicating another location in ROMS (step_t.F?) to apply this mask would be more appropriate.

Lyon.

[arango]

Applying the wetting and drying mask in step_t or to any tracer anywhere is really bad Lyon, this was actually the problem that you were having with the perfect restart in your Cook Inlet application. We were multiplying with such mask during I/O. I mentioned this before. If a particular cell is dry in the current timestep n (any tracer at that cell will be zero if multiplied by the wet/dry mask) and wet in the next timestep n+1 (any tracer there will not be zero), we will be in trouble when time stepping the solution and getting a huge gradient of values for tracers. Recall that with simple time stepping we need:

Code: Select all

    t(new) = t(old) + t_RHS

If you stand a the beach, you will notice that the temperature of the water indeed does not change too much every few seconds when the waves comes to our feet. In ROMS, the wetting and drying still have a minimum water thickness (Hz cannot be zero since with divide by Hz in the numerical kernel). We just set the velocities to zero in the dry cell but the tracers need to have similar values when dry or wet. We cannot have zero values because of multiplying by the wet/dry mask to the direct tracer value will yield huge time gradients (and instabilities) in the above equation.

[jivica]

Just wondering would it be possible to use one if statement in a way if cell is dry then use tracer old state (from the previous time step) without applying anything i.e. heatflux etc?
And keep doing that until it get wet again.
probably it is like that already in new version
Ivica

[wilkin]

I was thinking something doing something like this, probably at the end of set_vbc.F so that it catching all possible ways of setting the heat flux (whether by bulk fluxes, qcorrection, or prescribed).

Code: Select all

# ifdef MASKING
        DO j=JstrR,JendR
          DO i=IstrR,IendR
            stflx(i,j,itemp)=wetdry(i,j)*stflx(i,j,itemp)
# ifdef SALINITY
            stflx(i,j,isalt)=wetdry(i,j)*stflx(i,j,isalt)
# endif
          END DO
        END DO
# endif

[lanerolle]

I finally got my Cook Inlet ROMS model working in a stable fashion and the instabilities I saw were exactly as Zafer described. To stabilize my simulations, I had to do two things:

(1) hand-tune the bathy-topo to ensure that the vertical grid aspect ratios were limited in value during the drying phase and,
(2) re-distribute in time the river discharge of a couple of rivers in upper Knik Arm to ensure that discharge spikes (during Spring) over shallow mudflats did not cause the flow to excessively accelerate (this lead to both u-cfl/v-cfl violation + weird wetting-drying scenarios)

Thanks for all your help and suggestions - they worked!