model blowup due to very large velocity along topography

[sonaljit.m]

Hello,
I am facing an issue with my simulation using ROMS, and need some suggestions.

I am initializing my domain with HYCOM and NCODA analyzed fields from 2015 March, forced with atmospheric fluxes from MERRA 2. The bathymetry is obtained from ETOPO1.

https://www.dropbox.com/s/00jxl0wj8q67s ... p.png?dl=0

The figure above shows the bathymetry. The coastlines represent Puerto Rico and the Virgin Islands. Near the northern wall of my domain, there is a trench with depths exceeding 6km. Right along the border of that trench (as shown by the black n white spots near the northern wall), the velocity is exceeding to extremely large values, around 1 to 10 m/s which is obviously not reasonable. I'm showing only the zonal component, but the meridional component also shows the same thing at those locations. This is happening within the first few time steps no matter what time step I choose.

The figure below shows the same data, but focuses on one of the regions with the anomaly. You can see the velocity contour lines layered over the bathymetry.

https://www.dropbox.com/s/cuetkdglccwi5 ... d.png?dl=0

This sudden increase in the velocity is occuring along the entire depth profile at that location. The plot below shows the velocity profile at one of those locations, along the first 7 time steps. And then the simulation is blowing up (due to CFL criterion).

https://www.dropbox.com/s/000uwrfdi5l37 ... s.png?dl=0

Has anybody come across such a problem before? Can anyone suggest how to get over this anomaly?

Thanks,
Sonaljit.

[kate]

It looks like a barotropic instability mode, so I'd be looking at the barotropic timestep, not the baroclinic one. The back and forth nature of it looks like bottom drag. What are you doing for the bottom drag?

[sonaljit.m]

Thanks for the reply. Yes, turns out that the barotropic timestep was not small enough.

So, I increased NDTFAST, allowing more number of barotropic timesteps between DT, and those anomalies along the north topography disappeared. The simulation continues for a little while, and then again blows up on the 42nd timestep. Here is a snapshot at the 42 timestep, for the velocity u layered over the topography.

https://www.dropbox.com/s/00jxl0wj8q67s ... p.png?dl=0

In this simulation, the u is increasing to unrealistic values near the east and west walls, at regions where the topography is the shallowest. I set the minimum topography depth in my simulation at 10 m. The regions where u is unrealistic, have a topography depth within 10 to 20 m. (the actual topography at those regions is in the same range)

https://www.dropbox.com/s/dk92w8gm56cgo ... r.png?dl=0

This plot shows the horizontal advection for barotropic velocity. It is almost 0 everywhere, but attains large values at the regions where u is unrealistic (near west and east walls). The advection types I am using are as follows:

#define UV_C2ADVECTION (2nd-order centered)
#define TS_C2HADVECTION
#define TS_C2VADVECTION

For vertical mixing I am using GLS_MIXING, with CANUTO_A.

Unlike the back n forth nature in the earlier simulation, the u velocity gradually rises to such high values in this one. Regarding the bottom drag, I am using quadratic drag with coeffient 2.d-02. Do you have any suggestion for changing the drag type?

Also, is 10 m too low for the minimum topography depth?

Thanks.

[kate]

In this simulation, the u is increasing to unrealistic values near the east and west walls, at regions where the topography is the shallowest. I set the minimum topography depth in my simulation at 10 m. The regions where u is unrealistic, have a topography depth within 10 to 20 m. (the actual topography at those regions is in the same range)

The advection types I am using are as follows:

#define UV_C2ADVECTION (2nd-order centered)
#define TS_C2HADVECTION
#define TS_C2VADVECTION

Why those? I thought second-order was kept as an option just to show how bad it is compared to anything else.

For vertical mixing I am using GLS_MIXING, with CANUTO_A.

Probably not the culprit here.

Unlike the back n forth nature in the earlier simulation, the u velocity gradually rises to such high values in this one. Regarding the bottom drag, I am using quadratic drag with coeffient 2.d-02. Do you have any suggestion for changing the drag type?

Also, is 10 m too low for the minimum topography depth?

Quadratic drag is fine, especially with the LIMIT_BSTRESS option. I have used 10 m minimum depth, but it is asking a lot of your bottom smoothing in such a region. Did you do bottom smoothing?

Running for 42 steps before blowing up does not mean you have a short enough timestep yet - running for thousands of steps before blowing up means it's probably something else. Also, because it's blowing up near the boundary, what is your boundary treatment?

[sonaljit.m]

Okay, I made a few more changes based on your points. Advection is now 4th order central for all parameters.

#define UV_C4ADVECTION
#define TS_C4HADVECTION
#define TS_C4VADVECTION

Regarding bottom smooting, here's what I did: Obtained the bathy data h from ETOPO1, then lowpassed it with 10 km cutoff. After that I used the function smoothing_Positive_rx0(hgrd.mask_rho, h, rx0_max) (after masking), from bathy_smoother module. (If you have other suggestions about smoothing, please let me know)

Also, I am using an idealized wind stress 0.02 N/m2 (not varying spatially) throughout this domain.

On the boundary walls, I am using radiation boundary condition. I am providing climatology data processed from HYCOM/NCODA analysis (interpolated to ROMS grid).

Its still blowing up.

[kate]

On the boundary walls, I am using radiation boundary condition. I am providing climatology data processed from HYCOM/NCODA analysis (interpolated to ROMS grid).

So, you have radiation boundary conditions only? Those are not terribly well-behaved. What are you doing with the HYCOM/NCODA fields? Initial conditions only or climatology nudging of some sort? Everywhere or just around the edges? If you have it as a nudging at the edges, why not as part of the boundary condition too?

Does it run longer with a shorter timestep?

[sonaljit.m]

Here are the steps I'm taking for building the ROMS climat fields.

HYCOM/NCODA: I first design the ROMS grid, then interpolate the fields from the HYCOM/NCODA grid to the ROMS grid. During interpolation, I am first doing the horizontal 2D interpolation and then the vertical interpolation on the sigma surfaces.

Boundary condition: Radiation (Rad) (no nudging on edges). My climatology file and ROMS grid file both have the same grid-structure.

I went through a number of threads about model stability issues, and made a few more changes to my grid.
- Increased the vertical resolution near the bottom by modifying theta_b.
- Smoothing the bathymetry by lowpass filtering (the earlier grid was also lowpassed, but this time the lowpassing is stronger). (if there's a better suggestion for smoothing, please let me know)

After the 1st 4 timesteps (dt = 30 seconds), this is what it looks like:
https://www.dropbox.com/s/00jxl0wj8q67s ... p.png?dl=0

The gray contours show bathymetry, and colored lines are the u vel. Very large values near 2 m/s seem to occur in specific locations along the periphery of the white colored region (where bathymetry is a minimum set to 10 m).

Here's a closeup of the anomaly near east wall:
https://www.dropbox.com/s/cuetkdglccwi5 ... d.png?dl=0

The current run blowsup after 4 time-steps. However after reducing dt to 10 seconds (previous one was 30 sec), the run blowsup at the 6th timestep.