CROCO blowing up at 500m resolution

Hello,
I am facing a problem with my CROCO setup, I was hoping if anyone familiar with this problem might have a possible solution. I am trying to simulate an idealized problem of an island within the open ocean. The island is elliptical shaped, as shown here in green landmask:

The min depth is 10 m, and steeps down to 1500 m towards the open ocean. The grid details are given below.

  • Horizontal resolution: 500 m
  • Min depth 10 m, max depth 1500 m
  • Haidvogel number: 0.203, Haney number: 2.96
  • theta_s = 5, theta_b = 0.0, Tcline = 10.0
    Based on what I’ve read, the Haidvogel no should be < 0.25 and Haney no should be < 3.0, in order to give stable results.
    I initialized it with climatology from HYCOM 2018 and forcing from ECMWF.

I started with a time step dt of 10 sec and NDTFAST 120, but it is blowing up within 5 steps. Then I reduced the time step to 2 sec keeping NDTFAST the same, but it is still blowing up after 66 sec. Here is a timelapse video showing the zonal velocity transect and the SSH at each timestep till it blows up:
video timelapse attached

I noticed that around the 54th time step, the velocity forms zig-zag patterns near the landmask. These patterns amplify to unrealistically large values. Also the SSH forms such patterns in the same area, and blows up eventually.

I further reduced the time step to 1 sec and it seems to run fine without blowup (1 sec timelapse). But running the model at such small time steps is quite unrealistic, so I am wondering what am I doing wrong here?

Here is a list of preprocessors I am using for my setup:
TIDES
OBC_EAST
OBC_WEST
OBC_NORTH
OBC_SOUTH
CURVGRID
SPHERICAL
MASKING
WET_DRY
NEW_S_COORD
SOLVE3D
UV_COR
UV_ADV
SALINITY
NONLIN_EOS
UV_HADV_UP3
UV_VIS2
UV_VIS_SMAGO
UV_VADV_SPLINES
TS_HADV_RSUP3
TS_DIF4
TS_VADV_AKIMA
SPONGE
GLS_MIXING
GLS_KEPSILON
CANUTO_A
QCORRECTION
SFLX_CORR
CLIMATOLOGY
ZCLIMATOLOGY
M2CLIMATOLOGY
M3CLIMATOLOGY
TCLIMATOLOGY
ZNUDGING
M2NUDGING
M3NUDGING
TNUDGING
ANA_BSFLUX
ANA_BTFLUX
SSH_TIDES
UV_TIDES
POT_TIDES
TIDERAMP
OBC_M2CHARACT
OBC_M3ORLANSKI
OBC_TORLANSKI
SPLIT_EOS
TS_HADV_C4
MPI_COMM_WORLD
M2FILTER_POWER
HZR
VAR_RHO_2D
UV_MIX_S
VIS_COEF_3D
DIF_COEF_3D
TS_MIX_ISO
TS_MIX_IMP
TS_MIX_ISO_FILT
NTRA_T3DMIX
SPONGE_GRID
SPONGE_DIF2
SPONGE_VIS2
NF_CLOBBER

Thanks,
Sonaljit.

Follow up - based on some reading from the documentation, I made the following changes:

  • used 5th order WENO advection schemes for lateral and vertical advection.
  • for subgrid turbulence, used the 3D mixing schemes GLS_MIX2017_3D and UV_VIS_SMAGO_3D
  • Used NBQ, with DIAG_CFL for active diagnosis of the CFL parameters

With these changes, I used a time step dt 10 sec and NDTFAST = 480. The setup ran for over 1300 time steps without any errors. I am gonna run it for more time steps to see what happens. Hopefully, this is will go smoothly.

Anyone familiar with running fine resolution ROMS (~500 m) with variable terrain? Can you confirm if I am using the right preprocessors?

Thanks.
Sonaljit.