Hydrostatic consistency and the rx0 rx1 factors

[rduran]

It is well known that the bathymetry should be smoothed to avoid fictitious pressure gradients and that is where the conditions

3 < rx1 (Haney) < 7
0 < rx0 (Beckman and Haidvogel) < 0.4

come from.

in Robert Miller's "Numerical modeling of Ocean Circulation" 2007 Cambridge University Press, one can read regarding the pressure gradient error in sigma coordinates (page 136 after some analysis) that (please use LaTeX to read below expression):

We therefore obtain the expected behavior of the error decreasing as $\delta \sigma$ decreases as long as $\left| \frac {\sigma_j D_x} {D} \right| \Delta x \leq \Delta \sigma $ . This condition is often referred to as hydrostatic consistency

Also in page 141 after discussing Haney's and Beckman and Haidvogel's work, Miller states specifically about Beckman and Haidvogel's r factor:

If we recast (hydrostatic consistency expression) as $ r \leq 2 \Delta \sigma / \sigma_j$ we can see that this choice will result in violation of hydrostatic consistency unless r is chosen << 0.2 or the resolution in $\sigma$ is very coarse

My point is that pressure gradient errors depend also on the vertical resolution. I would like to test a new vertical discretization for hydrostatic consistency before I interpolate a year worth of daily boundary conditions only to find I have a fictitious pressure gradient spicing up my dynamics.

In conclusion I have a couple of questions:

1) Why is only smoothing of the bathymetry (to reduce D_x) mentioned (in the forum or wiki roms) as a solution to avoid this problem when the vertical resolution is also a factor and the same smoothed bathymetry may work for one vertical discretization and fail for another one? ... Am I missing something?

2) I wrote a few lines of matlab code to check for hydrostatic consistency given a roms_grd.nc but I am not sure I am understanding hydrostatic consistency criterion correctly, it seems the way I calculated it might be a little to restrictive or perhaps I should be using Cs_w: can someone take a look and comment? Ideally if this ends up being a useful check then others can benefit from this code:

Code: Select all

function count2=check_hydro_consistency(grd)
% count=check_hydro_consistency(grd)
% Uses roms_get_grid.m
% grd=roms_get_grid(roms_grd.nc); 
% $Id: roms_get_grid.m 393 2010-07-20 15:11:52Z wilkin $
h=grd.h;
dx=1./grd.pm;
dy=1./grd.pn;
Cs_r=grd.Cs_r;
dsigma=diff(grd.z_r,1); 
% I use 1 in diff command because in my case depth is the first dimension
% in the z_r array, if depth is the third dimension in your array change 1
% for 3; Notice this also would affect dsigma(:,k,kk) below

[dhdx dhdy]=gradient(h,dx(:,1),dy(:,1)); %this is D_x and D_y
%in my case I only need a column because spacing in Xi and Eta directions
%depends only on latitude [i.e. Eta direction - note the dimensions for my
%2D arrays are  (dim(Eta),dim(Xi)) ]
        
% check for hydrostatic consistency from Miller 2007:
% I think that this would be
% max(abs(Cs_r*dhdx(k,kk)/h(k,kk))*dx(k,kk)) <= min(dsigma(:,k,kk))
% so if  
% max(abs(Cs_r*dhdx(k,kk)/h(k,kk))*dx(k,kk)) > min(dsigma(:,k,kk)) 
% the test fails

count=0;
count2=zeros(size(h));

disp('Starting checkout -----------------')
for k=1:size(h,1)
    for kk=1:size(h,2)

a= max(abs(Cs_r*dhdx(k,kk)/h(k,kk))*dx(k,kk));
b= min(dsigma(:,k,kk));
c= max(abs(Cs_r*dhdy(k,kk)/h(k,kk))*dy(k,kk));

%test for failure in Xi direction
if a > b 
count=count+1;
count2(k,kk)=1;
end

%test for failure in Eta direction
if c > b 
count=count+1;
count2(k,kk)=1;
end

    end
end
disp(['Consistency not met ', num2str(count),' times'])
disp('Finished checkout -----------------')
figure
pcolor(count2), shading flat
colorbar
title('One values are where consistency is not met','Fontsize',14)
end

[staalstrom]

Hi

I study internal tides in a sill fjord (Oslofjord, Norway).
The topography is complicated, and I use a horizontal resolution of 75m.
In this application I get currents of up to 9 cm/s when I run a test case with a strong pycnocline at 20 m depth (same as the sill depth) and no horizontal gradients in the stratification. I managed to reduce this error to about 1 cm/s, which is of the same order as the accuracy of a typical ADCP.

I programmed a more complicated transform that is dependent on the water depth.
In roms the vertical transform is independent of depth. The depth dependency is handled by the stretching function. I have set the stretching to C=s (traditional sigma coordinates). I constructed the transform so I always have high vertical resolution around 20 m depth. This reduce the spurious pressure gradient error current a lot. But I get problems in the depth outside the sill where I have stratification below sill depth and coarse vertical resolution. In areas where the water depth is less than 20 m I get very thin layers and the haney number (rx1) is more than 100, I see no problems.

I have not understood what is the problem with the hydrostatic inconsistency (but I am still young). My experience is that high vertical resolution helps a lot. Coarse vertical resolution on the other hand is not good idea at all. So in my opinion the haney (rx1) number is not very useful.

Best regards
Andre Staalstrom (ans@niva.no)

[mathieu]

Yes the number given for rx0, rx1 correspond to the other one given in this forum.
I have not checked the formula but as far as I remember the hydrostatic consistency means that two adjacent vertical cells share one level and it is equivalent to rx1 <= 1.
Yes the erroneous pressure gradient depend on the vertical discretization, but a priori this error is taken into account by the rx1 which depends on the vertical discretization as opposed to the rx0 which does not depend on it.
Having hydrostatic consistency is very hard. It is equivalent to having rx1<=1 so you will have to make compromise. And yes you might have to run a model only to discover some problem later on. That is the nature of this business.
For computing hydrostatic consistency you need Sc_w, Cs_w, the Vtransform, Vstretching in order to get the Z_w. See the ROMS code for the formula.
Of course rx1 is an approximation to the true behavior. The real measure of error is the error itself and GETM now has a vertical dynamic adaptation mechanism that allow them a 10fold decrease in the horizontal pressure gradient.

[rduran]

Andre, Mathieu,

thanks for the helpful replies.

[wilkin]

For a comprehensive and up-to-date review of the pressure gradient truncation error issue, various historical approaches to reducing this error in sigma coordinate models, and how it manifests in modern versions of ROMS, I highly recommend reading section 5 of Sasha Shchepetkin's article in the Handbook of Numerical Analysis:

Alexander F. Shchepetkin, James C. McWilliams "Computational Kernel Algorithms for Fine-Scale, Multiprocess, Longtime Oceanic Simulations," Handbook of Numerical Analysis, Volume 14, 2009, Pages 121-183, P. G. Ciarlet (Ed.), Elsevier, doi:10.1016/S1570-8659(08)01202-0

The section concludes "Although PGF cannot be eliminated entirely, it can be verified that for flat stratification, the cancellation of terms ... is fourth-order accurate, and the new scheme is robustly tolerant of “hydrostatically inconsistent” grids with (delta-x/delta-z) · ∂z/∂x|s > 1 (Haney [1991])."