Blowing up problem using GLS mixing

[lcbernardo]

Dear ROMS users,

I'm encountering a problem using the GLS mixing option (k-epsilon in my case). In most ROMS configurations I've seen in the forums and other sites, if GLS mixing is used, the following combination of options is usually set in the header file:

# if defined GLS_MIXING
#define KANTHA_CLAYSON
#define N2S2_HORAVG
# endif

However, whenever I try to activate N2S2_HORAVG, my model always blows up near the start. Even if I change KANTHA_CLAYSON to CANUTO_A or CANUTO_B, the same situation occurs. However, if I don't activate N2S2_HORAVG, my model runs fine. I have used the KANTHA_CLAYSON and N2S2HORAVG combination successfully using MY25_MIXING, so I'm confused as to what the problem might be. I already tried using a smaller time-step, but the blow-up still occurs. Might anyone have any thoughts on what the problem might be? I would greatly appreciate any help.

My options in the header file are as follows:

#define UV_ADV
#define UV_COR
#define UV_VIS2
#define UV_SMAGORINSKY
#define TS_DIF2
#define TS_SMAGORINSKY

#define TS_U3HADVECTION
#define TS_C4VADVECTION

#define MIX_S_UV
/*#define MIX_GEO_UV*/
#define MIX_S_TS
/*#define MIX_GEO_TS*/
/*#define MIX_ISO_TS*/

#undef ANA_GRID
#define MASKING

#define NONLIN_EOS
#define SALINITY
#define ATM_PRESS

/*#define PJ_GRADPQ4 */
#define DJ_GRADPS
#define SOLVE3D

#define WET_DRY

/*** Option for Boundary condition ***/
#define RADIATION_2D

/*** Option for tidal forcing ***/
#define SSH_TIDES
#define UV_TIDES
/*#define RAMP_TIDES*/
#define ADD_FSOBC
#define ADD_M2OBC

#define SOLAR_SOURCE

#define BULK_FLUXES
#ifdef BULK_FLUXES
# define LONGWAVE
# define EMINUSP
# define COOL_SKIN
# define ANA_SRFLUX
# define ALBEDO
#else
# define ANA_SMFLUX
# define ANA_STFLUX
#endif

/* define only one of the following 5 */
#define UV_LOGDRAG
/*#define UV_QDRAG*/
/*#define MB_BBL*/
/*#define SG_BBL*/
/*#define SSW_BBL*/

#ifdef MB_BBL
/*# define MB_CALC_ZNOT*/
#endif

#ifdef SSW_BBL
/*# define SSW_CALC_ZNOT*/
# define SSW_LOGINT
#endif
#define LIMIT_BSTRESS

#define ANA_NUDGCOEF

#ifdef SOLVE3D

/* define only one of the following 4 */
/*# define BVF_MIXING*/
# define GLS_MIXING
/*# define MY25_MIXING*/
/*# define LMD_MIXING*/

# if defined GLS_MIXING
/*# define CRAIG_BANNER*/
# define KANTHA_CLAYSON
/*# define CANUTO_A*/
/*# define CANUTO_B*/
/*# define K_C4ADVECTION*/
/*# define K_C2ADVECTION*/
# define N2S2_HORAVG
/*# define ZOS_HSIG*/
/*# define TKE_WAVEDISS*/
# endif

# if defined MY25_MIXING
# define KANTHA_CLAYSON
# undef CANUTO_A
/*# define K_C4ADVECTION*/
# define N2S2_HORAVG
# endif

# ifdef LMD_MIXING
# define LMD_RIMIX
# define LMD_CONVEC
# define LMD_SKPP
# define LMD_BKPP
# define LMD_NONLOCAL
# endif

# define ANA_BPFLUX
# define ANA_BTFLUX
# define ANA_BSFLUX
# define ANA_SPFLUX

#endif

[jcwarner]

this is a situation where you will need to dig in a little deeper. You say with n2s2 avg turned on, it blows up near the beginning. Suggest you see how many steps it takes before it blows up, lets say it is a 100. then save the model output in a his file every 5 times (or so). that will give you 20 time steps. then see what is happening to the region of instability. there should also be a restart file that has the final time step in it. that should show where the problem is. At that location is it have strong lateral gradients in bathy - temp -salt ?? maybe the initial condition is strange at that location (interp temp from a global model but this is near a land point that was not in the global model)??? dig in a little deeper. of course you can run the model without that option, but it will probably eventually develop other issues, so i suggest you figure this out first.

-john

[wilkin]

As John Warner says, you should save output nearing the time of blow-up and try to catch where and for what reason the model is going awry. I suggest capturing the same time for the case that does not blow-up, and compare.

Perhaps repeat the Richardson number calculation offline to visualize the difference that the smoothing makes.

Generally speaking, I've found k-kl within GLS to be more robust than k-epsilon, with no real difference in the physics solution. You might explore using k-kl settings unless you have a particular reason to expect k-epsilon to be a more appropriate closure for your problem.

I also notice you have:

Code: Select all

#define MIX_S_TS
/*#define MIX_GEO_TS*/
/*#define MIX_ISO_TS*/

If your blow-up is near steep bathymetry the false diapycnal mixing due to horizontal along S-surfaces could be a culprit, and your solution in the no N2S2 smoothing case, while stable, may be poor.

Consider activating

Code: Select all

#define MIX_GEO_TS

[lcbernardo]

Dear Dr. Warner and Dr. Wilkin,

Thank you for sharing your thoughts and suggestions regarding how to possibly resolve this issue. Indeed my modeling domain includes some relatively steep bathymetric features, and I'll try to investigate the blowup more closely using the restart file as basis, and saving data at finer intervals nearing the blowup.

As for the proper GLS setting, I have been testing various configurations to see which one gives the best reproduction of some of the field data we have. So far I've tested the model with MY25_MIXING and LMD_MIXING without experiencing blowups, and have just started testing the GLS_MIXING options. So I'll be giving the k-kl some consideration as well. I will also try activating MIX_GEO_TS, though I guess my main task if to figure out if the steep bathymetric features are leading to some unrealistic solutions. If I am able to resolve this, I'll post about it in this topic.

Thanks,
Lawrence

[lcbernardo]

I'd like to post an update on this problem, after trying out some of the suggestions. First, after checking the log file, I found that the blowup occurs at a grid cell which is at the westernmost tip of an island, and at the surface level (NaN values appear). The depths in the immediate area are shown in the following figure, where the problem grid cell is indicated:

depth_values.png (37.9 KiB) Viewed 6506 times

I have tried a bit of smoothing in this area, but the blowup still occurs, and at the same time step and location. I also tried the k-omega and gen settings for GLS, as well as switching on MIX_GEO_TS instead of MIX_S_TS, and the same problem occurs. However, when I set it up to use k-kl, the model is able to continue running. For reference, here are excerpts from the respective log files:

k-epsilon:
8588161 5964 00:01:00 2.419131E-02 1.586564E+04 1.586567E+04 2.739042E+14
(056,119,05) 1.563521E-02 2.946967E-02 1.069506E-01 1.195170E+00
8588162 5964 00:02:00 2.413290E-02 1.586560E+04 1.586562E+04 2.739043E+14
(055,119,15) 1.849584E-02 3.710627E-02 1.981509E-01 1.192609E+00
8588163 5964 00:03:00 2.408643E-02 1.586557E+04 1.586560E+04 2.739043E+14
(022,062,15) 1.074742E-03 4.076535E-03 3.344242E-01 1.264081E+00
8588164 5964 00:04:00 2.404642E-02 1.586558E+04 1.586560E+04 2.739044E+14
(112,094,15) 1.545827E-03 1.696484E-03 2.903132E-01 1.250261E+00
8588165 5964 00:05:00 2.401343E-02 1.586559E+04 1.586561E+04 2.739044E+14
(159,085,15) 1.857205E-03 6.168238E-03 2.637254E-01 1.189504E+00
8588166 5964 00:06:00 2.398749E-02 NaN NaN 2.739043E+14
(159,085,15) 7.452478E-04 6.328074E-03 5.408304E-01 1.188985E+00

Blowing-up: Saving latest model state into RESTART file

WRT_RST - wrote re-start fields (Index=1,1) in record = 0000002

k-kl:
8588161 5964 00:01:00 2.419131E-02 1.586564E+04 1.586567E+04 2.739042E+14
(056,119,05) 1.563521E-02 2.946967E-02 1.069506E-01 1.195170E+00
8588162 5964 00:02:00 2.413315E-02 1.586560E+04 1.586562E+04 2.739043E+14
(055,119,15) 1.849586E-02 3.710625E-02 1.981509E-01 1.192609E+00
8588163 5964 00:03:00 2.408686E-02 1.586557E+04 1.586560E+04 2.739043E+14
(022,062,15) 1.074746E-03 4.076525E-03 3.344241E-01 1.264080E+00
8588164 5964 00:04:00 2.404696E-02 1.586558E+04 1.586560E+04 2.739044E+14
(112,094,15) 1.545869E-03 1.696480E-03 2.903131E-01 1.250259E+00
8588165 5964 00:05:00 2.401404E-02 1.586559E+04 1.586561E+04 2.739044E+14
(159,085,15) 1.853111E-03 6.169464E-03 2.637265E-01 1.189504E+00
8588166 5964 00:06:00 2.398816E-02 1.586559E+04 1.586562E+04 2.739043E+14
(159,085,15) 7.391761E-04 6.329990E-03 5.408346E-01 1.188984E+00
8588167 5964 00:07:00 2.396657E-02 1.586561E+04 1.586563E+04 2.739043E+14
(181,081,15) 9.293810E-04 4.238111E-04 5.872347E-01 1.192991E+00
8588168 5964 00:08:00 2.394613E-02 1.586563E+04 1.586565E+04 2.739043E+14
(181,081,15) 1.853887E-03 1.785394E-03 4.559670E-01 1.225200E+00
8588169 5964 00:09:00 2.392591E-02 1.586565E+04 1.586567E+04 2.739044E+14
(183,082,15) 6.380530E-03 3.828488E-03 3.301767E-01 1.257228E+00
8588170 5964 00:10:00 2.390668E-02 1.586566E+04 1.586568E+04 2.739043E+14
(183,082,15) 7.763664E-03 6.004263E-03 5.780641E-01 1.232996E+00
8588171 5964 00:11:00 2.388978E-02 1.586566E+04 1.586568E+04 2.739042E+14
(159,085,15) 2.890900E-03 2.703220E-03 5.309864E-01 1.208348E+00
(and so on...)

I guess for now I'll be trying to see what more I could do to smooth the bathymetric features and see what works.

Lawrence

[jcwarner]

if the closure is making it unstable, then look at temp, salt, u, v, and AKv and AKt. perhaps the options of

#define LIMIT_VDIFF --> use to impose an upper limit on vertical diffusion **
#define LIMIT_VVISC --> use to impose an upper limit on vertical viscosity **
might help ???

Also, when the depth is -100, are you using wet/dry? are those cells land masked?
I see you have wet/dry on. that option was not really designed to function in a setting where the bathy goes from 60 to -100. it is more for a nearshore tidal setting. So if you are using wet/dry to just remove points that are -100 (above sea level) i would suggest you use land/sea masking instead and just make the -100 values to be 5 (or something) annd just mask them out.

[lcbernardo]

Thank you for your reply Dr. Warner. I have wetting and drying turned on, but the cells with a -100 are part of the permanent land mask. Even so, will it be useful to still change the value to something like 5 or so? But anyway, I did try to turn off WET_DRY, but the model still blew up, and at a different point and sigma layer this time. And this point is not adjacent to a land cell, which makes me think that I could still do more to maybe smooth the bathymetric features. Might you have any suggestions on some tools to help facilitate this? So far, I have tried the Shapiro filter matlab script included in the ROMS matlab folder and even the QUICKIN module of the Delft3D modelling system, but I feel they are not working well since I still could mostly not get the recommended rx0 and rx1 whenever I run the model. Or if I do, it would be for a bathymetry smoothed in such a way that some of the real features are no longer recognizable.

[kate]

I feel they are not working well since I still could mostly not get the recommended rx0 and rx1 whenever I run the model. Or if I do, it would be for a bathymetry smoothed in such a way that some of the real features are no longer recognizable.

Indeed, you show neighboring points with rx0 of almost 0.58, which is larger than I like to go. To go from 60 m to 250 m in that distance, you need another grid point in between to keep the r-values down.

I can't recommend the tools I use since I can barely keep them working as systems change. When they really break for good I'll have to focus on the Python stuff. Still, the algorithms of smoothing more where the r-values are large should work in any language. But yes, that does change the nature of the shelf-break and you have to live with some compromises, possibly including more resolution.

[lcbernardo]

Thanks for your thoughts on this Kate. I don't think I mentioned in on this post, but the grid cell size is 1.5 km on a side. But indeed, I guess I'll have to consider using more resolution. I've also come across your bathymetry tools on GitHub, which I would like to look into more, even if I'm not a Python user. Might be worth my while to learn the language.