a possible improvement of the spline vertical viscosity

[mizuta]

Hi:
I have a question about the spline interpolation used to get the stress by the vertical viscosity in step3d_uv.F. I guess that derivatives at the boundary (or alternatives) should be prescribed in the spline interpolation. What kind of the boundary condition is assumed in ROMS? It seems to me that the wind stress is missing in the boundary condition, while the viscous stress should coincide with the wind stress at the surface. If this is the case, I wonder why such things can happen,

Just as a test of a possible candidate of modification, I added the following four lines (numbered 1,2,3,4) to the step3d_uv.F,

Code: Select all

            DC(i,k)=cff*(u(i,j,k+1,nnew)-u(i,j,k,nnew)-                 &
     &                   FC(i,k)*DC(i,k-1))
          END DO
        END DO

1       DO i=IstrU,Iend
2         DC(i,N(ng)-1)=DC(i,N(ng)-1)                                   &
3    &                  -CF(i,N(ng)-1)*sustr(i,j)/AK(i,N(ng)-1)
4       END DO

!
!  Backward substitution.
!
        DO i=IstrU,Iend
          DC(i,N(ng))=0.0_r8
        END DO
        DO k=N(ng)-1,1,-1
          DO i=IstrU,Iend
            DC(i,k)=DC(i,k)-CF(i,k)*DC(i,k+1)
...

and similar lines for the meridional velocity.
Then, I performed an idealized run of wind-driven Ekman flows (which I made by myself). The thickness of the Ekman layer, sqrt(akv/f), was resolved by three grid points with the total grid numbers, N=30, in the vertical direction. The Ekman spiral was reasonably reproduced, and
compared to the theory (Pedlosky 1987) the velocity at the shallowest grid point was;

(1) 0.900 times of the theory, for the case without modification
(2) 0.999 times of the theory, for the case with the above noted modification
(3) 0.995 times of the theory, for the case "undef SPLINE" in cppdefs.h

The results was improved by the modification. Actually, the result in example (1) implies that the energy input by winds is underestimated by 10 percent in this example. This is worse than example )3), in which a simpler scheme is employed.
Thanks,

Genta

[shchepet]

Code: Select all

     DC(i,N(ng)-1)=DC(i,N(ng)-1)-CF(i,N(ng)-1)*sustr(i,j)/AK(i,N(ng)-1)

Yes, this is potentially a strong move, but it is also dangerous and can backfire.

What you are proposing is essentially to set the uppermost grid-box-side derivative
of the velocity field to du/dz = sustr/Akv at the free surface. This is the only natural
physical boundary condition, provided that
(1) Akv remains finite at surface; and
(2) the vertical profile of u (Ekman spiral or whatever) is adequately resolved.

The problem is that most vertical mixing parameterization yield vanishing Akv
at free surface, Akv --> 0, as z--> zeta, (this is often known as law of wall), or
if some regularization associated with roughness of free surface is assumed,
then Akv goes to something small, thought >0. This smallness may cause
disastrous oscillations of vertical spline within the upper grid box, and
a completely non-physical solution.

So instead of the natural physical boundary condition, several artificial
variants were tried. The one we found work the best (on the argument of
compromise of accuracy vs. robustness) is the assumption that the vertical
profile within the uppermost grid box is linear.

This is a very simple ad hoc condition, and it can be possibly improved.
Note that the blue word adequately above is not mathematically strictly
defined. There must be a some kind of slope-limiting algorithm: at first try
to set du/dz = sustr/Akv at free surface, but restrict it to a smaller value of
the same sign, if some limit based on local curvature of u-profile is exceeded.
However, at this time I do not know how to specifically implement this.


...It occurs to me that despite being used in ROMS for about 10 years, the
splined version of vertical viscosity/diffusion operator was never openly
described in any document. Some long time ago I made a write-up, mostly
for myself, and, accordingly, it is far from being complete,
http://www.atmos.ucla.edu/~alex/ROMS/vert_visc.pdf
but I believe it is still useful to explain what is going on in the code,
and how the tri-diagonal splined implicit operator is formed.

[mizuta]

Hi:
Thank you for your valuable comments and document.

of the velocity field to du/dz = sustr/Akv at the free surface. This is the only natural physical boundary condition, provided that
(1) Akv remains finite at surface; and
(2) the vertical profile of u (Ekman spiral or whatever) is adequately resolved.
The problem is that most vertical mixing parametrization yield vanishing Akv at free surface, Akv --> 0, as z--> zeta, (this is often known as law of wall), or

That is true that I did not consider the law of wall. It will be so difficult to fit such velocity as du/dz=1/(z-zeta) by polynomials. According to your document, there are two methods for calulating the vertical viscosity. One is with the spline interpolation and the other is without the interpolation, and the latter has no singularity at the surface. So, I got the impression that this latter method may not be so bad, if for example,

(1b) Akv=c*(zeta-z) and
(2) the vertical profile of u (Ekman spiral or whatever) is adequately resolved.

In what parameter ranges (the vertical resolution etc) the method with the spline interpolation is recommended to be used? Is it necessary to worry about the small Akv in such parameter ranges?
Thanks,
Genta