Near-bottom numerical instability over the slope

[apaloczy]

Hello all,

I'm having problems with pressure gradient errors over an idealized continental slope application. My initial field is just flat stratification, with an N2 profile derived from a climatological dataset. I have baroclinic (barotropic) velocity errors of 10 cm/s (4 cm/s) after a steady-state is reached, when a clear along-slope jet develops (attached).

The horizontal resolution is ~1.65 km in both xi and eta directions. The number of vertical levels is 50, and the coarsest vertical resolution is 80 m, near the bottom (attached). Vtransform = 2, Vstretching = 4, Tcline = theta_s = theta_b = 0. Sponge layers that ramp the harmonic viscosity and diffusivity coefficients up to 5 times their inner values are set along all 3 open boundaries.

Boundary conditions are: Chapman explicit for the free-surface, Shchepetkin for the normal 2D velocities, Gradient for the TKE field and Radiation + Nudging for the tangential 2D velocities, 3D velocities and tracers. All nudging is done towards the initial conditions at the boundaries, that is, flat T and S profiles, and zero everywhere for all velocities and the free-surface, reflecting the initial state of rest.

The Beckmann and Haidvogel number (rx0) is 0.075, and the Haney number (rx1) is 3.7, which I understand are conservative values. The current values for the harmonic viscosity and diffusivity coefficients is 200 m2/s. I experimented with 50 m2/s, and the error is only about 1 cm/s larger. With a biharmonic viscosity of 4.5e+07, the error is about the same. With biharmonic viscosity and biharmonic diffusivity though, the error oscillates violently between 10 cm/s and 100 cm/s.

The CPP options I'm using are:

UV_ADV, UV_COR, UV_QDRAG, DJ_GRADPS, SOLVE3D
UV_U3HADVECTION, TS_U3HADVECTION
UV_C4VADVECTION, TS_C4VADVECTION
UV_VIS2, TS_DIF2
CURVGRID, SPHERICAL, MASKING
MIX_GEO_UV, MIX_GEO_TS, VISC_GRID, DIFF_GRID, MY25_MIXING, SPONGE, RADIATION_2D, ANA_SMFLUX, ANA_STFLUX, ANA_BTFLUX, T_PASSIVE, ANA_BPFLUX, ANA_SPFLUX

Talking to the author of this post (viewtopic.php?f=17&t=3323), we believe this might be a pressure gradient noise issue similar to the one he reported there. I implemented many of his useful suggestions, which brought improvements to the solution, but the pressure gradient errors are still much larger than what I would expect for my conservative values of rx0 (0.075) and rx1 (3.7). In fact, this paper (http://www.sciencedirect.com/science/ar ... 0302000033) compares different pressure gradient algorithms in a steep seamount test case, and shows baroclinic velocity errors of O(0.1 cm/s), with rx0 = 0.07 and rx1 = 2.7 and using the same pressure gradient algorithm I have been using (the parabolic splines density jacobian, CPP option DJ_GRADPS) and a similar maximum bottom slope. Also, their vertical resolution is 225 m, much coarser than my 80 m.

The things I have already tried so far are: Playing with the vertical stretching curve, increasing the vertical resolution, increasing the horizontal resolution, changing the mixing coefficients between harmonic and biharmonic (and their values), decreasing the actual slope of the continental slope (currently the maximum slope is 0.03), capping the topography at 2000 m (it was 4000 m before), using time-dependent horizontal mixing coefficients (UV_SMAGORINSKY and TS_SMAGORINSKY options), changing the rotation of the viscosity/diffusivity tensors between along S-surfaces and along geopotential surfaces (MIX_S_TS, MIX_GEO_TS, MIX_S_UV and MIX_GEO_UV options), changing the boundary conditions, putting the sponge on and changing the pressure gradient algorithm (which showed that the DJ_GRADPS is the algorithm that produces the smallest error for my application. This agrees with the result obtained for the study I mentioned above for the seamount test case, which has similar values of rx0 and rx1). Most of these changes have had small beneficial effects.

I have also checked for errors in the initial field itself, but it is pretty flat, as it should be, only with very small-scale sawed features that always appear in the areas of maximum bottom slope due to interpolation errors (attached).

Any suggestions on ways to decrease the errors on the 3D velocity field to an acceptable level (I think, lower than perhaps 1.0 cm/s) is greatly appreciated.

Thank you all for your time,

[lanerolle]

Have you tried another sigma grid formulation which refines the grid both at the top (surface) and bottom? What about something like V_transform=1, V_stretching=1, theta_s=4.5, theta_b=0.9, Tcline=10? V_transform=1, V_stretching=1 is more robust (but may not necessary be more accurate) than V_transform=2, V_stretching=4 (the recommended one) which is sensitive to the rx1 bathymetry smoothness parameter and can lead to numerical instabilities - at least, this is my own experience.

[apaloczy]

Hi Ianerolle, thank you for your reply.

With the parameters you suggested, the vertical stretching curve becomes the red line in the attached figure. The maximum Haney number rx1 increases from 3.7 to 14.1:

(ORIGINAL) N, theta_b, theta_s, tcline: 50.0, 0.0, 0.0, 0.0
(MODIFIED) N, theta_b, theta_s, tcline: 50.0, 0.9, 4.5, 10.0

Haney number for original case:
Maximum r-value = 3.69139749326
Haney number for modified case:
Maximum r-value = 14.1338503033
Coarsest dz for original case: 79 m
Coarsest dz for modified case: 84 m

However, I have been suggested to prioritize finer vertical resolution (dz, something like 20-30 m) over the hydrostatic consistency of the vertical grid (low rx1). What do you think of this approach? Playing with the stretching curve, I found it difficult to get both low rx1 and fine dz. For example, if I increase N from 50 to 120 with my original configuraton (V_transform=2, V_stretching=4, theta_s = theta_b = Tcline = 0), my coarsest dz goes from 79 m to 33 m, but rx1 increases from 3.7 to 8.9.

I tried lowering rx1 as much as possible in previous tests, by decreasing the number of vertical levels and coarsening the vertical resolution near the bottom, but it did not seem to decrease the errors.

Also, why do you say Vtransform = 1, Vstretching = 1 is a more robust configuration? I mean, do you get less noise with it, relative to the recommended Vtransform = 2, Vstretching = 4 configuration, even though rx1 gets big because of the near-bottom refinement?

Thank you for your time.

[rsoutelino]

Hi Andre,

This post has additional info on grid stiffness: viewtopic.php?f=14&t=2841&p=10698&hilit ... rx1#p10698

[saulo]

The rx1 is a bit tricky as it seems application dependent. A value considered high may work on a given setup while a "low" value may not (a setup that depends largely on the topographical slope, which looks real steep in your case and the sharpness of gradients but also on the algorithms your choosing to solve it with may have an impact). Like Ivica would say thgough, its all in the grid. In a high resolution application I had to smooth the bathymetry significantly more than I was expecting, I had a reasonable rx1 to start with, something like 7, but had to bring it down to 5 before it started working with "no" PGE. The bad rx1 were on the steep mid slope part near the bottom layers. With the high horizontal resolution I was using, it required even more vertical resolution in those layers than I could afford, otherwise the changes in thickness were still too strong it seems.

[apaloczy]

Thanks everyone for your replies.

Following Ianerolle's suggestion, I did two test runs changing only the vertical stretching curve as depicted on my previous post. It seems that, for my application, the configuration with Vtransform = 2 and Vstretching = 4 produces less Horizontal Pressure Gradient Error (HPGE). I've attached plots of mid-water velocity and kinetic energy/maximum velocity histories.

ct). Like Ivica would say thgough, its all in the grid. In a high resolution application I had to smooth the bathymetry significantly more than I was expecting, I had a reasonable rx1 to start with, something like 7, but had to bring it down to 5 before it started working with "no" PGE. The bad rx1 were on the steep mid slope part near the bottom layers. With the high horizontal resolution I was using, it required even more vertical resolution in those layers than I could afford, otherwise the changes in thickness were still too strong it seems.

I see, saulo. That seems to be where the HPGE errors are largest in my application too. But I'm afraid I didn't follow everything you said. You said that you had to increase vertical resolution? But that increases rx1, right? So that's something that has been confusing me: To decrease HPGE, do you need to have low rx0 or fine dz? Or both? Because the way I understand it, if you refine dz, you increase rx1. Was the decrease in HPGE in your case only because you decreased rx1 or you tried other things too?

I have experimented with the bottom slope too. the maximum slope was ~0.08, and now is down to ~0.03.

Would it be a reasonable mitigation approach to subtract this flat solution from the solutions with forcing? Evidently, this would not correct any errors due to nonlinear interactions between the flat HPGE flow and the physical flow, I don't know. Has this been tried or looked at?

Thanks,

[saulo]

The error decreased because I decreased rx1 alone.
Actually my case was I had high vertical resolution, a large number of vertical levels and I would had to increase the horizontal resolution more than I could aford to bring rx1 down (so the reverse of what I said, sorry about that, it was a while ago). rx1 is sort of how well your grid points are aligned, so if you have lots of points high dz you need more dx,dy to maintain things leveled. Your case seems that you have high dx,dy but the slope may be still too steep. Not sure there would be a number of vertical levels that would mitigate it significantly. rx1 is easy to calculate offline, you can play around with different dx, dz to get the best possible value but it is still no guarantee it will work. rx1 is not everything in terms of hpge and some slopes may be just too steep. Once with POM I think had to use a 6th order scheme to minimize it enough. Boa sorte

[apaloczy]

Got it, saulo. I will fiddle with the grid a little more to decrease rx1. So based on the information I found so far and considering that I have tried pretty much everything, it looks like my slopes are just too steep. At least the largest errors are not located around the shelfbreak, which is the area I am most interested in analyzing.

Does anyone have any other suggestions? They will be very useful.

[jivica]

viewtopic.php?f=17&t=1080#p12021

see at the bottom -> my last response, hope it helps.
cheers
I.

[apaloczy]

Thanks jivica, this procedure seems very useful. I will follow that in another application where I'm seeking a realistic topography.

But in this application with idealized topography, I don't understand why I'm getting this much HPGE with this configuration. I'm using safe rx0 and rx1 values (0.075 and 3.7), and the bottom slope is not overly steep, I guess (0.03).

Also, I have already done the equivalent of the procedure you describe by decreasing the bottom slope (i.e., smoothing) between the 200 m and 2000 m isobaths, which is where the largest HPGE velocities are located. As you can see in the attached figures, the maximum slope was originally 0.12, then I tried 0.08, then 0.03, and still have HPGEs of O(10 cm/s).

Thanks for looking at this.

[mastrorocco]

Hi Andre,

I suspect there might an inconsistency between the s-coordinate formulation used to create your initial condition and the one used in ROMS. This could lead to horizontal density gradients in your computations.

I don't know why you are using Tcline = 0 and I am not sure if this is fine (perhaps someone else can comment on that!). Noticed that ROMS uses the following to control the thickness of the stretching hc = min(hmin,Tcline). See the documentation https://www.myroms.org/wiki/index.php/V ... coordinate. Make sure the Python/Matlab script you are using to create your initial condition is consistent with what ROMS is doing. I would also recommend closing all the boundaries when assessing the pressure gradient errors.

Hope it helps. Good luck!

Gustavo

[hugobastos]

Hi there,
seems that you had a lot of fun with this...and i think that gustavo have a good point. Are you really sure about your initial conditions and their configuration in the sigma coordinates? Because your rx values are already too conservative and even in other grids where the rx are higher its suppose to not raise such errors after a week or so.

Regarding the processing of the initial and boundary files when you change the Vfunctions for the lanerolle sugestion, did you remake your initial and bry conditions? If not, you need too...it's a basic thing, but chance of this to occurs when your are doing a lot of tests with vertical functions is very high!

Also note that your maximum hpge will be close to the north since there is a large slope there...it will be nice to see your velocity horizontal distribution close to the bottom ( or the layer that have the maximum value...).

I think a good point to make sure any errors regarding the treatment of the initial conditions is to provide a simpler T profile. I usually just use the ANA_INITIAL from another case ( say SEAMOUNT), close the boundaries (clamped) and run for 10days ( your results seems to be "stable" after 20 days, so between 1week and 20 days its a reasonable time window). Note that the time here is important since you dont want your basin to be contaminated with the reflection over the boundaries.

You can also test this is another grid (say your grid is a problematic one, or maybe its easier to check in a light one). Say the WC13 example. Pickup the grid, change to an analytical initial condition (again SEAMOUNT or any other, just remember about the cppdefs and the linear equation of state and its parameters - they change across the .in files) and run the test with the boundary closed. It's suppose to raise a low hpge error since the bathymetry is very smooth (I've done this in the past). You can also try your initial condition function on this grid to check for consistency and the vertical stable profile in sigma coordinates.

Again your boundaries are open, and this can pose another problem if they were created with wrong functions or if you change the vertical transform (assuming the same tools are used for ini/bry).

dont give up hehehe

Cheers

[apaloczy]

Thanks mastrorocco and hugobastos for bringing this vertical coordinate issue to my attention. It seems like a very good source of trouble.

don't know why you are using Tcline = 0 and I am not sure if this is fine (perhaps someone else can comment on that!).

I set Tcline to 0 in order to minimize bottom refinement, and hence rx1. But looking at the two transformation equations (eqs. 1 and 2) between sigma and S coordinates in the wiki, it seems that using Tcline=0 (hence hc=0) would not cause any problems. Just in case, I am now using non-zero Tcline. Thanks for pointing this out.

Noticed that ROMS uses the following to control the thickness of the stretching hc = min(hmin,Tcline). See the documentation https://www.myroms.org/wiki/index.php/V ... coordinate. Make sure the Python/Matlab script you are using to create your initial condition is consistent with what ROMS is doing. I would also recommend closing all the boundaries when assessing the pressure gradient errors.

From the wiki, hc = min(hmin,Tcline) is used only if Vtransform=1. For Vtransform=2 (my current case), it says that Tcline can be any positive value (so I guess this also means trouble if Tcline=0).

I am preparing the initial conditions files by taking a single T,S profile and linearly interpolating it to the model vertical grid, point-wise in (xi,eta). To do that, I need to get the depths of all rho points in each (xi,eta) point. So I checked the code I'm using to do that. I'm using the pyroms package, the pyroms.vgrid.s_coordinate* group of functions (https://github.com/kshedstrom/pyroms/bl ... s/vgrid.py). Comparing this code against the stretching functions and transformation equations in the wiki, it seems to be correct.

I kept the boundaries open to try to isolate HPGEs from boundary reflection noises. I did a lot of tests with the open boundary conditions and the sponge layer, and it seems to be working fine now. The noise in the free-surface field at day 100 is smaller than 0.5 cm, for example.

Regarding the processing of the initial and boundary files when you change the Vfunctions for the lanerolle sugestion, did you remake your initial and bry conditions? If not, you need too...it's a basic thing, but chance of this to occurs when your are doing a lot of tests with vertical functions is very high!

Yes, I did. Double-checked that now.

Also note that your maximum hpge will be close to the north since there is a large slope there...it will be nice to see your velocity horizontal distribution close to the bottom ( or the layer that have the maximum value...).

Yes, I have been looking at the horizontal velocity in the layer where its maximum is found (which is almost always in the first or second S-level near the bottom), see my first post. Regarding the maximum slope in the north part of the domain, it is gone now, I switched to a near-even slope all along the continental slope (see the "second decreased bottom slope" figure on my last post). I simply linearly interpolated between the 200 m and 2000 m isobaths. After that switch, the maximum HPGE decreased by ~4 cm/s and spreaded all along the continental slope. But the maximum HPGE indeed was in that steepest slope in the north prior to this switch.

Again your boundaries are open, and this can pose another problem if they were created with wrong functions or if you change the vertical transform (assuming the same tools are used for ini/bry).

I see your point. But to make the *_bry files I just slice the boundaries of the *_ini files and copy the vertical stretching variables, so that should be ok.

I think a good point to make sure any errors regarding the treatment of the initial conditions is to provide a simpler T profile. I usually just use the ANA_INITIAL from another case ( say SEAMOUNT), close the boundaries (clamped) and run for 10days ( your results seems to be "stable" after 20 days, so between 1week and 20 days its a reasonable time window). Note that the time here is important since you dont want your basin to be contaminated with the reflection over the boundaries.

You can also test this is another grid (say your grid is a problematic one, or maybe its easier to check in a light one). Say the WC13 example. Pickup the grid, change to an analytical initial condition (again SEAMOUNT or any other, just remember about the cppdefs and the linear equation of state and its parameters - they change across the .in files) and run the test with the boundary closed. It's suppose to raise a low hpge error since the bathymetry is very smooth (I've done this in the past). You can also try your initial condition function on this grid to check for consistency and the vertical stable profile in sigma coordinates.

Thanks for this suggestion, I will try it out to make sure that it isn't the initial conditions and try out different grids.

[saulo]

Its best to close the boundaries to evaluate the error for sure.
Looking back at some of your plots there seems to be a lot more activity near the boundaries than expected.
For an application with no forcing and no horizontal variability in T/S you really don't have to worry about boundary reflections. They are irrelevant, there are no real dynamics in the model anyway, the error is very likely not going to grow exponentially because of that. But bad BCs can def cause weirdness and to a lesser extend very unbalanced initial conditions too. Since you smoothed things and it didn't make a lot of difference, I think Gustavo is right and you should check your initial/bdry profiles.

[apaloczy]

Updates:

I capped the topography at 1200 m instead of at 2000 m, keeping everything else the same, and the HPGE decreased from 10-30 cm/s to ~3 cm/s. Also, the sudden acceleration to ~40 cm/s in the first 5 days is gone. The HPGE now increases more smoothly (attached).

After day ~50, weird banded structures just like the ones described in this post (viewtopic.php?f=17&t=3323) start to develop (see example attached), and increase the noise on the continental slope to ~4-6 cm/s. Since I am not too interested in integrating beyond day ~50, I should be fine with those. But following the suggestion of that post's author, I tried switching from MIX_GEO_UV to MIX_S_UV, and there was a decrease in ~1.5 cm/s in this banded-like noise, and its development was delayed in ~10 days (not shown).

So, apparently, there is something weird happening on the slope between 1200 m and 2000 m. Capping at 1200 m caused the bulk of the noise to dissapear, but there is still a 3 cm/s noise.

Following the mastrorocco's and saulo's suggestions, I tried analytical initial conditions with ANA_INITIAL and closed all boundaries. I started with the stratification of the SEAMOUNT case, which is a very gentle exponential profile of T with constant S. The HPGE in this case grows linearly up to 1.2 cm/s on day 50 (attached).

After that, I tried using a more steep (larger N2) linear T profile with T(z=-1200 m) = 4 C and T(z=0 m) = 26.5 C and constant S = 35.5, in an attempt to mimic my original profile. The HPGE also grew linearly until day 50 but this time it reached 3 cm/s (attached), the same level of HPGE I am getting with my non-analytical T,S profiles. So it seems that the HPGE is sensitive to N2.

After that, I calculated the depths of RHO-points, z_r(i,j,k), with both the matlab module (function set_depth.m) and the pyroms.vgrid submodule. The maximum absolute difference (point-wise) between the two arrays was O(1e-13 m), while the finest dz I have is 2e-1 m. So it seems we can rule the stratification out. Besides, from the tests I mentioned above, I'm getting the same amount of HPGE both with z_r being computed internally (ANA_INITIAL) and being specified externally (calculated with pyroms).

Any other ideas? Do you think 3 cm/s is the lowest HPGE I can get?

[hugobastos]

Hello again andré,

all right, so your function/initial distribution is right and there is no error associated with the initial conditions. And yes, the hpge is sensible to the profile you impose, since the pressure at two z point would be different in steeper sigma coordinates. The thing with the gentle seamount case would just to put a "standard" value on your hpge error and check for your function. I think you in the right way,imposing the local profile to watch for magnitude of the errors, but i have something in mind regarding your plot of velocities:

This plot of maximum points is a single point, or is the maximum in every timestep? I'm thinking that you probably used the maximum per timestep, which indicates that the location of the point can be varying right?

My tip is that after the ~50day, the point where you find the maximum value of velocity should have change close to the boundary. So maybe your conclusions about the interior error being the double that you see after sometime it's just a dynamic response of a signal that reach a cross solid boundary and need to flow along. Albeit the magnitude of the signal generated is low, they indeed propagate around your basin , closing the "gyre". In the boundary corners of your grid they should generate some higher velocities (for example in the shelf break of the north/south boundary).

Sometimes this maximum velocity may not be the real maximum velocity of your hpge, since the interior values should be less than that. If you decrease the time of integration you will not be able to watch the signal arrive at the boundary so your maximum velocities will be smaller.

If this is happening, maybe you need to define a fix region for evaluating your maximum velocities, think about the dynamics of that region and decide about the value based on the mean generated locally. The manifestation of the hpge in the simulations is kinda of a local feature around the topography (and you usually spot this way before this signal in some way can propagate through the domain/boundary), so your estimations of the error in the global maximum velocity can be biased.


Cheers

[apaloczy]

Hi hugobastos, thanks for your reply.

This plot of maximum points is a single point, or is the maximum in every timestep? I'm thinking that you probably used the maximum per timestep, which indicates that the location of the point can be varying right?

Yes.

My tip is that after the ~50day, the point where you find the maximum value of velocity should have change close to the boundary. So maybe your conclusions about the interior error being the double that you see after sometime it's just a dynamic response of a signal that reach a cross solid boundary and need to flow along. Albeit the magnitude of the signal generated is low, they indeed propagate around your basin , closing the "gyre". In the boundary corners of your grid they should generate some higher velocities (for example in the shelf break of the north/south boundary).

Sometimes this maximum velocity may not be the real maximum velocity of your hpge, since the interior values should be less than that. If you decrease the time of integration you will not be able to watch the signal arrive at the boundary so your maximum velocities will be smaller.

If this is happening, maybe you need to define a fix region for evaluating your maximum velocities, think about the dynamics of that region and decide about the value based on the mean generated locally. The manifestation of the hpge in the simulations is kinda of a local feature around the topography (and you usually spot this way before this signal in some way can propagate through the domain/boundary), so your estimations of the error in the global maximum velocity can be biased.

That seems to be what is happening here, I think. I realized that by examining the whole velocity field, but now the errors due to momentum getting trapped in the northern boundary (which are not HPGE, only an indirect effect of it) and propagating around it ciclonically vanished with some changes in the OBCs and the sponge layers (see attached). Now the white dot that marks the cell with the maximum velocity doesn't jump around anymore, it stays near the bottom, and within the region of the domain with the maximum bottom slope, where one would expect the HPGE to be maximum. The maximum HPGE is down to ~1.35 cm/s at day 100 (attached) anyway, so that is good enough, I guess.

Thanks again,