I'm trying to model water currents in one natural bay of the Mideterranian Sea. The modelling area is relatively small (50x50 km approx). As the result of it, geometry is pretty ragged with several one-grid bays and peninsulas. Also, the bay has a form of boot with long nose. When wind forces are included in model, after approximately 4 hours of modeling time models blow up. Apparently, the surface elevation builds up in the 'nose' part of the bay.
What are the possible ways to avoid such behaviour? Also, what are the 'normal' values for kinetic and potential energies?
Thank you in advance
Artem Babayan
Model blows up in ragged region
The model should theoretically be able to handle one-grid bays. However, we don't generally allow them in our grids because the plotting software doesn't manage to plot those bays, making trouble there hard to diagnose. I have allowed one-grid passages between plottable bits into my grids.
Liz Dobbins had to increase the viscosity in one domain to damp out large tides that were causing trouble. Other things to try are bottom friction. Still, since we don't have experience with one-grid bays, there could be a possible ROMS bug involved. You'd think a large surface gradient would push the wter back out. How shallow are you getting it? Or is it trying to get to a negative water thickness? We clamp the water depth to something like 10 m in Cook Inlet even though the reality involves wetting/drying.
Liz Dobbins had to increase the viscosity in one domain to damp out large tides that were causing trouble. Other things to try are bottom friction. Still, since we don't have experience with one-grid bays, there could be a possible ROMS bug involved. You'd think a large surface gradient would push the wter back out. How shallow are you getting it? Or is it trying to get to a negative water thickness? We clamp the water depth to something like 10 m in Cook Inlet even though the reality involves wetting/drying.
I have several model domains with one-grid point wide bays. At 1 km resolution for the New Jersey shore I have the narrow back bay of Barnegat Bay and Little Egg Harbor, both with narrow entrances that get reasonably large tidal velocities and wind set-up within the bay.
At 10 km resolution the Chesapeake Bay part of a US East Coast model has a lot of bays that are 1 grid point wide. Some have rivers entering and others don't - no stability problems there.
So I don't believe there is any fundamental problem with narrow channels in ROMS. If it's critical you achieve correct dynamics in a narrow channel then you might want to re-evalaute your choice of resolution.
It the Bay of Fundy I found it necessary to make the head of the bay deeper than I originally had chosen because the large tide amplitudes could take the depth plus sea level (h+zeta) close to zero on the low tide, and also drive excessively large velocities. But this is not related to the bay width.
With regard to Kate's comment about plotting, this is why I use Matlab to diagnose these sorts of problems - you can look at individual cell values and find where things are blowing up if you do a few restarts close to the time of instability.
John.
At 10 km resolution the Chesapeake Bay part of a US East Coast model has a lot of bays that are 1 grid point wide. Some have rivers entering and others don't - no stability problems there.
So I don't believe there is any fundamental problem with narrow channels in ROMS. If it's critical you achieve correct dynamics in a narrow channel then you might want to re-evalaute your choice of resolution.
It the Bay of Fundy I found it necessary to make the head of the bay deeper than I originally had chosen because the large tide amplitudes could take the depth plus sea level (h+zeta) close to zero on the low tide, and also drive excessively large velocities. But this is not related to the bay width.
With regard to Kate's comment about plotting, this is why I use Matlab to diagnose these sorts of problems - you can look at individual cell values and find where things are blowing up if you do a few restarts close to the time of instability.
John.
John Wilkin: DMCS Rutgers University
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu
-
artem_vb
Thank you all for your replies. It seems like the part of the problem arises from ragged coastline -- the water is blown into bays and stays there, which causes numerical instablity. After I artificially smoothed the coastline the model was running for 16 hours of instead of 4, as before, but still blows up eventually. Note, that without wind, when tidal forces only are present the models runs smoothly for several days.
The reply by John Wilkin (thank You again) assumes that the roots of problems lie somewhere else.
I suspect the main sources of problem may be:
a) Relatively shallow simulated area (max. depth ~150m, but average is around 40-45m).
b) Coarse resolution (500m) for such an area (50x50km) -- at some points the depth in adjacent cells is 14m, 3m and 0m.
So the questions are:
1) Am I right and can it really be the reasons for model to blow up (as you can probably guess I'm not an expert in this field)?
2) If I use finer grid (100m) will ROMS be able to cope with it (that is will the pprimitive-equations model still valid under such conditions)?
3) Are the ways to smooth the bathymetry somehow?
The reply by John Wilkin (thank You again) assumes that the roots of problems lie somewhere else.
I suspect the main sources of problem may be:
a) Relatively shallow simulated area (max. depth ~150m, but average is around 40-45m).
b) Coarse resolution (500m) for such an area (50x50km) -- at some points the depth in adjacent cells is 14m, 3m and 0m.
So the questions are:
1) Am I right and can it really be the reasons for model to blow up (as you can probably guess I'm not an expert in this field)?
2) If I use finer grid (100m) will ROMS be able to cope with it (that is will the pprimitive-equations model still valid under such conditions)?
3) Are the ways to smooth the bathymetry somehow?
None of the parameter ranges you use:
Can we assume your C-grid velocity and rho points masks are consistent with your coastline, i.e. how do you build your land/sea mask arrays?
ROMS solves the hydrostatic primitive equations. If you believe the hydrostatic assumption is still valid (depth scale/horizontal scale << 1) then using a 100 m grid size can be valid.
John.
are particularly different from our typical coastal applications. Water should not accumulate in the coastal cells.a) Relatively shallow simulated area (max. depth ~150m, but average is around 40-45m). b) Coarse resolution (500m) for such an area (50x50km) -- at some points the depth in adjacent cells is 14m, 3m and 0m.
Can we assume your C-grid velocity and rho points masks are consistent with your coastline, i.e. how do you build your land/sea mask arrays?
ROMS solves the hydrostatic primitive equations. If you believe the hydrostatic assumption is still valid (depth scale/horizontal scale << 1) then using a 100 m grid size can be valid.
John.
John Wilkin: DMCS Rutgers University
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu
-
artem_vb
Thank you for your time. I really appreciate it.
Er, do you mean, whether am I use curvilinear grid? Then the answer is 'no'. The area is divided with ordinary rectangular grid. I had a feeling, that it can cause problems as well, because quite a large chunk of the cells are land cells, so this is not very effective way of doing things... Can it be helped somehow?Can we assume your C-grid velocity and rho points masks are consistent with your coastline, i.e. how do you build your land/sea mask arrays?
The cell depth of 0 m could very well be the problem. I run a program that clamps the depth to some fixed minimum depth, say 30 m or 10 m, depending. If the water depth goes negative, you will be in trouble.artem_vb wrote: a) Relatively shallow simulated area (max. depth ~150m, but average is around 40-45m).
b) Coarse resolution (500m) for such an area (50x50km) -- at some points the depth in adjacent cells is 14m, 3m and 0m.
So the questions are:
1) Am I right and can it really be the reasons for model to blow up (as you can probably guess I'm not an expert in this field)?
2) If I use finer grid (100m) will ROMS be able to cope with it (that is will the pprimitive-equations model still valid under such conditions)?
3) Are the ways to smooth the bathymetry somehow?
We do smooth the bathymetry. There is a metric for the change in depth that you can reasonably get away with:
rvalue = |h1-h2| / (h1+h2)
Getting this larger than 0.4 is considered unsafe. For instance, 14 m adjacent to 3 m gives 0.65. Smoothing - it's an art, not a science. I've got some Fortran codes that John Wilkin wrote for doing this.
When I come across very rugged and highly varying bathymetry, I usually filter it
using a simple Laplacian filter and even with a little bit of filtering, I find a massive
improvement in the model performance - numerical stability and accuracy of results.
In Matlab, you can filter very easily by using the conv2(H,f) command (H is your
bathymetry and f is the filter stencil eg. 3x3, 5x5, etc.). I vary the weights on the
stencil using the discrete function f=1/(A*(i*i +j*j) where i, j=-1, 0, 1 (for a 3x3 stencil)
and A=2^(K-1) and K=6 is sufficient usually. Speciallly for small bodies of water, it is
important to ensure that you maintain the volume of water when filtering the bathymetry
and the Laplaciian filter seems to do a good enough job.
using a simple Laplacian filter and even with a little bit of filtering, I find a massive
improvement in the model performance - numerical stability and accuracy of results.
In Matlab, you can filter very easily by using the conv2(H,f) command (H is your
bathymetry and f is the filter stencil eg. 3x3, 5x5, etc.). I vary the weights on the
stencil using the discrete function f=1/(A*(i*i +j*j) where i, j=-1, 0, 1 (for a 3x3 stencil)
and A=2^(K-1) and K=6 is sufficient usually. Speciallly for small bodies of water, it is
important to ensure that you maintain the volume of water when filtering the bathymetry
and the Laplaciian filter seems to do a good enough job.
-
artem_vb
Thank you all very much for your help. I'm now trying to implement several options in order to improve stability.
May I ask another question, may be not so directly connected with this topic -- I want to save the surface wind stresses in order to ensure, that I calculate wind influence right. How can I do it? It seems, like ROMS does not save it by default. Are there any switches in the cpp definition files? I could not locate any which is directly responsible for this.
May I ask another question, may be not so directly connected with this topic -- I want to save the surface wind stresses in order to ensure, that I calculate wind influence right. How can I do it? It seems, like ROMS does not save it by default. Are there any switches in the cpp definition files? I could not locate any which is directly responsible for this.
Artem,
Please read all the way through the comments in the ocean.in file. All the input/output switches you need to select output of the wind stresses are documented there, as are the options that control how often history, averages and diagnostics files etc are created and how often data is saved.
Though there is no formal manual for ROMS yet, most of the basic information for getting started is described in ocean.in and cppdefs.h.
John.
Please read all the way through the comments in the ocean.in file. All the input/output switches you need to select output of the wind stresses are documented there, as are the options that control how often history, averages and diagnostics files etc are created and how often data is saved.
Though there is no formal manual for ROMS yet, most of the basic information for getting started is described in ocean.in and cppdefs.h.
John.
John Wilkin: DMCS Rutgers University
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu