No slip boundary conditions not being implemented

[suyashbire1]

We are running a zonally re-entrant setup on a spherical grid forced by restoring buoyancy at the surface. It seems that the model is not conserving zonal angular momemtum. We have NOSLIP = True in the MOM_input file. But when we look at the output, the model seems to have free slip turned on instead. This is crucial because any numerically generated spurious angular momentum is expected to be dissipated by friction at the boundaries. We believe this is not happening because of free slip condition.

The documentation does say that implementation of no slip condition is not as clean as free slip. So I'm wondering if the implementation of no slip conditions still not complete. Or is there something we might be missing from our setup? Let me know if you need to take a closer look at our setup.

[pittwolfe]

A few more details: the domain is zonally periodic, 22 degrees wide and 60 degrees N/S extent. The forcing is relaxation to a zonally symmetric buoyancy profile equivalent to a uniform temperature of 20ºC south of 20ºN, 6ºC north of 40ºN and a linear ramp between the two. Even though the initial conditions are horizontally uniform and at rest, the model develops a rather hefty zonal current of transporting around 120 Sv. (See attached figure) Since the zonal momentum is initially zero and there are no sources of zonal momentum, one would expect a purely baroclinic current or, if there's enough nonlinearity, banded jets with a net zero zonal transport over the whole domain. (Note that the final flow is steady and zonally symmetric, so there are no eddies or physical instabilities.)

Obviously, MOM is not conserving zonal momentum. This is not itself too worrying, since it's hard to conserve momentum in spherical coordinates. One usually hopes that dissipation would destroy the excess momentum generated by the numerics and you'd end up with a somewhat physical final state. In this geometry, the only way to dissipate momentum are drag against the sides and bottom. We are using quadratic bottom drag, but the vertically sheared flow has arranged itself so the bottom velocities are zero. This leaves side drag. We have NO_SLIP = True, but inspection of the zonal velocity (also on attached figure) shows that the boundary velocity is not zero, though its normal derivative is. Thus, the model seems to be using free slip boundary conditions even though no slip conditions are requested. Comments in the configuration file recommend against using no-slip because free-slip conditions are more "natural" on a C grid, but in this relatively simple test case, the combination of free-slip and lack of momentum conservation leads to a completely unphysical solution. No slip conditions are compatible with C grid finite volume codes (even for biharmonic friction) and are implemented in, e.g., the MITgcm.

Are there plans to implement no-slip in MOM6? Am I just misinterpreting these results? If no-slip conditions are really implemented in MOM6, how can they be turned on (since setting NO_SLIP = True seems not to work)?

Thanks, Christopher

[adcroft]

@suyash309 Can you make the configuration directory available so we can look at the setup? Or post MOM_parameter_doc.short ...

[suyashbire1]

I'm trying to upload it all on Github but am not able to. I'll figure it out. For the time being i'm pasting the MOM_parameter_doc.short here.

UPDATE: I figured out what was going wrong with my git.

Anyway, here is the configuration directory.

COORD_CONFIG = "USER" : This user defined initialization is here.

THICKNESS_CONFIG = "USER_witheta" : This user defined initialization is here.

SPONGE_CONFIG = "USER_southonly" : Sponges are defined here.

BUOY_CONFIG = "USER" : Buoyancy forcing is here.

You can see some of the figures from this run after about 20 years here.

[MJHarrison-GFDL]

@suyash309 Can you also provide the sponge/forcing and initial condition files? You should be able to create your own fork on Github then push your changes, otherwise just send the files as an attachment.

[suyashbire1]

@MJHarrison-GFDL, @adcroft I edited my previous post. There are links to everything there. Thank you for your help.

[adcroft]

I'm wondering about BT_STRONG_DRAG = True. This is an option within the barotropic-baroclinic split which has a lot of options due to the extended development of the algorithm. It could be overestimating the bottom drag and causing the weirdness. Can you see if setting it to False helps?

[suyashbire1]

I've changed it to BT_STRONG_DRAG = False. I'll update here sometime tomorrow when I get some results.

UPDATE: Hi @adcroft unfortunately, that does not seem to have changed anything. Here are the results after about 11 years of the run.

[MJHarrison-GFDL]

Did you consider moving the sponge away from the boundary?

On Mon, Oct 5, 2015 at 10:46 PM, Suyash Bire notifications@github.com wrote:

I've changed it to BT_STRONG_DRAG = False. I'll update here sometime tomorrow when I get some results.

— Reply to this email directly or view it on GitHub https://github.com/NOAA-GFDL/MOM6/issues/230#issuecomment-145728147.

[suyashbire1]

You mean extending it meridionally to, say, 10N instead of the current 4N?

[MJHarrison-GFDL]

Sorry, what I meant to say was that the jet should be away from the boundary.

I don't yet fully understand the problem here. If the desired solution is no flow at the boundary, this seems inconsistent with having a jet only a few grid points away.

On Wed, Oct 7, 2015 at 10:54 AM, Suyash Bire notifications@github.com wrote:

You mean extending it meridionally to, say, 10N instead of the current 4N?

— Reply to this email directly or view it on GitHub https://github.com/NOAA-GFDL/MOM6/issues/230#issuecomment-146219815.

[pittwolfe]

The problem is the jet. We only specify buoyancy forcing at the surface and southern sponge and forcing and sponges are zonally uniform, so there should be no source of zonal momentum. Therefore the jet is not there because of the forcing and there is no incompatibility between the forcing and boundary condition. Even if the jet were formed by some turbulent inverse cascade, the tangential flow should go to zero at the boundary due to no slip. What we see, however, is that the normal derivative of the tangential velocity goes to zero, suggesting that the model is giving us free-slip boundary conditions even though we asked for no-slip boundary conditions.

[Hallberg-NOAA]

Hi Christopher and Suyash,

  Yes, you are absolutely right that the coastal lateral stresses
probably are being zeroed out in your idealized test case, depending
on how the layer thicknesses are being set in the halo points on
land.  

  This has been a good discussion, that others might find helpful. 
I originally encoded the no-slip boundary conditions, so perhaps I
can comment on a number of issues that have arisen in this
conversation thread.

  1) MOM6 is written with the intention that it conserves angular
momentum pretty well. For example, conservation of angular momentum
dictates the structure of the lateral stress tensor calculations in
MOM6. However, a quadratic or linear bottom drag term or flow around
bathymetric features can be effective ways to extract angular
momentum from the solid earth if there is any near-bottom flow (as
is usually the case). However I do not think that we have carefully
examined this in all possible cases, so if there is clear evidence
that this is not the case, we would like to have that drawn to our
intention.  (The "pretty well" has to do with the deliberate use of
different definitions of the layer thickness at velocity points
depending on the flow direction and speed and relative thicknesses
of adjacent points, but that is a very long a detailed technical
discussion, much of which will appear in the MOM6 documentation
paper that we are writing.

  2) The no-slip boundary conditions appear in the code in two
places - in the estimate of relative vorticity at boundaries, and in
the estimate of the shears at the boundaries (as determined by a 2-d
land-mask). On a C-grid they do not dictate that the tangential
velocity is actually 0 at any explicitly modeled location.

  3) The ocean has no vertical sides. Density surfaces end where
they intersect the (usually sloping) bottom or the free surface. To
conserve angular momentum, we write the horizontal viscosity term on
a C-grid as 
  du / dt = 1/h_u Div([h_q K_q, h_h K_h] Grad u)
where u is the velocity, h_u is an estimate of the layer thickness
at velocity points, h_q is an estimate of the thickness at vorticity
points, h_h is the thickness at thickness points, and [K_h, K_q] are
the Laplacian viscosities at the thickness and vorticity points. 
Everything except for h_q and h_u are well defined from 2nd order
staggered discretizations or the C-grid variables themselves. 
Conservation of angular momentum requires that the same h_q appear
in the u- and v- equations, and we have chosen to use simple
arithmetic mean thicknesses for h_u in this term. Ensuring a finite
acceleration when one or more of these thicknesses at velocity
points go to 0 requires that the h_q also go to 0 in these
circumstances. Therefore, we use the harmonic mean of the harmonic
mean of the 2 values of h_u and harmonic mean of the 2 values of h_v
(see line ~466 of MOM_hor_visc.F90, which starts hq = ...).  The
harmonic mean of A and B is 2 A B / (A + B), and it has the nice
property that it is a second order accurate average of A and B, but
goes to 2 A when A << B.  This works well, but we do not
currently apply land masks, and setting h over land to 0 has the
effect of setting h_q at coastlines to 0, and disabling the fluxes
with a no-slip boundary condition.

 4) For the most part, in MOM6 we implement strong "lateral" drag
via a vertical boundary condition with the CHANNEL_DRAG scheme,
which calculates the bottom drag on a cell that partially intersects
a sloping bottom.  This seems like the right approach for realistic
simulations, but it obviously is not doing what we would like for
idealized vertical-walled oceans.

 5) Properly implementing a no-slip boundary condition will require
changing the definition of h_q at coastal points using the (2-d)
land-sea masks to determine when to replace the thicknesses of
masked velocity points with values that are temporarily extrapolated
from unmasked neighbors. This should work well enough for straight
or convex coastlines, but will introduce some (probably acceptable)
cross-contamination of u- and v- thicknesses at concave corners in
coastlines.  Once we have written and tested this code modification,
we will be closing out this issue.

  Hopefully all this is helpful!

  - Bob Hallberg

p.s., The BT_STRONG_DRAG only pertains to bottom drag was a red
herring in this case.
  
On 10/07/2015 11:59 AM, Christopher L. Pitt Wolfe wrote:

  The problem is the jet. We only specify buoyancy
    forcing at the surface and southern sponge and forcing and
    sponges are zonally uniform, so there should be no source of
    zonal momentum. Therefore the jet is not there because of the
    forcing and there is no incompatibility between the forcing and
    boundary condition. Even if the jet were formed by some
    turbulent inverse cascade, the tangential flow should go to zero
    at the boundary due to no slip. What we see, however, is that
    the normal derivative of the tangential velocity goes to zero,
    suggesting that the model is giving us free-slip boundary
    conditions even though we asked for no-slip boundary conditions.
  —
    Reply to this email directly or view
      it on GitHub.

[StephenGriffies]

Hi,

I am only vaguely following the stencil details. But the discussion of spurious momentum sources makes me ask the following, perhaps off-the-mark questions, but wishing to be sure we are passing the "laugh test".

With flat isopycnals and no momentum or buoyancy forcing at the ocean surface, does MOM6 keep the flow static if there is a vertical wall? What about sloping walls? Do answers to these questions depend on whether free-slip, no-slip, bottom drag, etc?

Seems that the physically correct behaviour is static flow with horizontal isopycnals will stay static, regardless the frictional boundary conditions. Capisce?

Steve

On Thu, Oct 8, 2015 at 5:30 PM, Robert Hallberg notifications@github.com wrote:

Hi Christopher and Suyash,

Yes, you are absolutely right that the coastal lateral stresses probably are being zeroed out in your idealized test case, depending on how the layer thicknesses are being set in the halo points on land.

This has been a good discussion, that others might find helpful. I originally encoded the no-slip boundary conditions, so perhaps I can comment on a number of issues that have arisen in this conversation thread.

1) MOM6 is written with the intention that it conserves angular momentum pretty well. For example, conservation of angular momentum dictates the structure of the lateral stress tensor calculations in MOM6. However, a quadratic or linear bottom drag term or flow around bathymetric features can be effective ways to extract angular momentum from the solid earth if there is any near-bottom flow (as is usually the case). However I do not think that we have carefully examined this in all possible cases, so if there is clear evidence that this is not the case, we would like to have that drawn to our intention. (The "pretty well" has to do with the deliberate use of different definitions of the layer thickness at velocity points depending on the flow direction and speed and relative thicknesses of adjacent points, but that is a very long a detailed technical discussion, much of which will appear in the MOM6 documentation paper that we are writing.

2) The no-slip boundary conditions appear in the code in two places - in the estimate of relative vorticity at boundaries, and in the estimate of the shears at the boundaries (as determined by a 2-d land-mask). On a C-grid they do not dictate that the tangential velocity is actually 0 at any explicitly modeled location.

3) The ocean has no vertical sides. Density surfaces end where they intersect the (usually sloping) bottom or the free surface. To conserve angular momentum, we write the horizontal viscosity term on a C-grid as du / dt = 1/h_u Div([h_q K_q, h_h K_h] Grad u) where u is the velocity, h_u is an estimate of the layer thickness at velocity points, h_q is an estimate of the thickness at vorticity points, h_h is the thickness at thickness points, and [K_h, K_q] are the Laplacian viscosities at the thickness and vorticity points. Everything except for h_q and h_u are well defined from 2nd order staggered discretizations or the C-grid variables themselves. Conservation of angular momentum requires that the same h_q appear in the u- and v- equations, and we have chosen to use simple arithmetic mean thicknesses for h_u in this term. Ensuring a finite acceleration when one or more of these thicknesses at velocity points go to 0 requires that the h_q also go to 0 in these circumstances. Therefore, we use the harmonic mean of the harmonic mean of the 2 values of h_u and harmonic mean of the 2 values of h_v (see line ~466 of MOM_hor_visc.F90, which starts hq = ...). The harmonic mean of A and B is 2 A B / (A + B), and it has the nice property that it is a second order accurate average of A and B, but goes to 2 A when A << B. This works well, but we do not currently apply land masks, and setting h over land to 0 has the effect of setting h_q at coastlines to 0, and disabling the fluxes with a no-slip boundary condition.

4) For the most part, in MOM6 we implement strong "lateral" drag via a vertical boundary condition with the CHANNEL_DRAG scheme, which calculates the bottom drag on a cell that partially intersects a sloping bottom. This seems like the right approach for realistic simulations, but it obviously is not doing what we would like for idealized vertical-walled oceans.

5) Properly implementing a no-slip boundary condition will require changing the definition of h_q at coastal points using the (2-d) land-sea masks to determine when to replace the thicknesses of masked velocity points with values that are temporarily extrapolated from unmasked neighbors. This should work well enough for straight or convex coastlines, but will introduce some (probably acceptable) cross-contamination of u- and v- thicknesses at concave corners in coastlines. Once we have written and tested this code modification, we will be closing out this issue.

Hopefully all this is helpful!

• Bob Hallberg

p.s., The BT_STRONG_DRAG only pertains to bottom drag was a red herring in this case.

On 10/07/2015 11:59 AM, Christopher L. Pitt Wolfe wrote:

The problem is the jet. We only specify buoyancy forcing at the surface and southern sponge and forcing and sponges are zonally uniform, so there should be no source of zonal momentum. Therefore the jet is not there because of the forcing and there is no incompatibility between the forcing and boundary condition. Even if the jet were formed by some turbulent inverse cascade, the tangential flow should go to zero at the boundary due to no slip. What we see, however, is that the normal derivative of the tangential velocity goes to zero, suggesting that the model is giving us free-slip boundary conditions even though we asked for no-slip boundary conditions. — Reply to this email directly or view it on GitHub.

— Reply to this email directly or view it on GitHub https://github.com/NOAA-GFDL/MOM6/issues/230#issuecomment-146691422.

Dr. Stephen M. Griffies NOAA Geophysical Fluid Dynamics Lab 201 Forrestal Road Princeton, NJ 08542 USA

[Hallberg-NOAA]

  There are no spurious momentum (or angular momentum) sources in
MOM6.

  We routinely verify this with the "resting" test case.

  - Bob Hallberg

On 10/08/2015 09:55 PM, Stephen Griffies wrote:
Hi,

  I am only vaguely following the stencil details. But the
  discussion of
  spurious momentum sources makes me ask the following, perhaps
  off-the-mark
  questions, but wishing to be sure we are passing the "laugh test".

  With flat isopycnals and no momentum or buoyancy forcing at the
  ocean
  surface, does MOM6 keep the flow static if there is a vertical
  wall? What
  about sloping walls? Do answers to these questions depend on
  whether
  free-slip, no-slip, bottom drag, etc?

  Seems that the physically correct behaviour is static flow with
  horizontal
  isopycnals will stay static, regardless the frictional boundary
  conditions. Capisce?

  Steve

  On Thu, Oct 8, 2015 at 5:30 PM, Robert Hallberg
  <notifications@github.com>
  wrote:

  >
  >
  >
  >
  >
  > Hi Christopher and Suyash,
  >
  > Yes, you are absolutely right that the coastal lateral
  stresses
  > probably are being zeroed out in your idealized test case,
  depending
  > on how the layer thicknesses are being set in the halo points
  on
  > land.
  >
  > This has been a good discussion, that others might find
  helpful.
  > I originally encoded the no-slip boundary conditions, so
  perhaps I
  > can comment on a number of issues that have arisen in this
  > conversation thread.
  >
  > 1) MOM6 is written with the intention that it conserves
  angular
  > momentum pretty well. For example, conservation of angular
  momentum
  > dictates the structure of the lateral stress tensor
  calculations in
  > MOM6. However, a quadratic or linear bottom drag term or flow
  around
  > bathymetric features can be effective ways to extract angular
  > momentum from the solid earth if there is any near-bottom
  flow (as
  > is usually the case). However I do not think that we have
  carefully
  > examined this in all possible cases, so if there is clear
  evidence
  > that this is not the case, we would like to have that drawn
  to our
  > intention. (The "pretty well" has to do with the deliberate
  use of
  > different definitions of the layer thickness at velocity
  points
  > depending on the flow direction and speed and relative
  thicknesses
  > of adjacent points, but that is a very long a detailed
  technical
  > discussion, much of which will appear in the MOM6
  documentation
  > paper that we are writing.
  >
  > 2) The no-slip boundary conditions appear in the code in two
  > places - in the estimate of relative vorticity at boundaries,
  and in
  > the estimate of the shears at the boundaries (as determined
  by a 2-d
  > land-mask). On a C-grid they do not dictate that the
  tangential
  > velocity is actually 0 at any explicitly modeled location.
  >
  > 3) The ocean has no vertical sides. Density surfaces end
  where
  > they intersect the (usually sloping) bottom or the free
  surface. To
  > conserve angular momentum, we write the horizontal viscosity
  term on
  > a C-grid as
  > du / dt = 1/h_u Div([h_q K_q, h_h K_h] Grad u)
  > where u is the velocity, h_u is an estimate of the layer
  thickness
  > at velocity points, h_q is an estimate of the thickness at
  vorticity
  > points, h_h is the thickness at thickness points, and [K_h,
  K_q] are
  > the Laplacian viscosities at the thickness and vorticity
  points.
  > Everything except for h_q and h_u are well defined from 2nd
  order
  > staggered discretizations or the C-grid variables themselves.
  > Conservation of angular momentum requires that the same h_q
  appear
  > in the u- and v- equations, and we have chosen to use simple
  > arithmetic mean thicknesses for h_u in this term. Ensuring a
  finite
  > acceleration when one or more of these thicknesses at
  velocity
  > points go to 0 requires that the h_q also go to 0 in these
  > circumstances. Therefore, we use the harmonic mean of the
  harmonic
  > mean of the 2 values of h_u and harmonic mean of the 2 values
  of h_v
  > (see line ~466 of MOM_hor_visc.F90, which starts hq = ...).
  The
  > harmonic mean of A and B is 2 A B / (A + B), and it has the
  nice
  > property that it is a second order accurate average of A and
  B, but
  > goes to 2 A when A << B. This works well, but we do not
  > currently apply land masks, and setting h over land to 0 has
  the
  > effect of setting h_q at coastlines to 0, and disabling the
  fluxes
  > with a no-slip boundary condition.
  >
  > 4) For the most part, in MOM6 we implement strong "lateral"
  drag
  > via a vertical boundary condition with the CHANNEL_DRAG
  scheme,
  > which calculates the bottom drag on a cell that partially
  intersects
  > a sloping bottom. This seems like the right approach for
  realistic
  > simulations, but it obviously is not doing what we would like
  for
  > idealized vertical-walled oceans.
  >
  > 5) Properly implementing a no-slip boundary condition will
  require
  > changing the definition of h_q at coastal points using the
  (2-d)
  > land-sea masks to determine when to replace the thicknesses
  of
  > masked velocity points with values that are temporarily
  extrapolated
  > from unmasked neighbors. This should work well enough for
  straight
  > or convex coastlines, but will introduce some (probably
  acceptable)
  > cross-contamination of u- and v- thicknesses at concave
  corners in
  > coastlines. Once we have written and tested this code
  modification,
  > we will be closing out this issue.
  >
  > Hopefully all this is helpful!
  >
  > - Bob Hallberg
  >
  > p.s., The BT_STRONG_DRAG only pertains to bottom drag was a
  red
  > herring in this case.
  >
  >
  > On 10/07/2015 11:59 AM, Christopher L. Pitt Wolfe wrote:
  >
  > The problem is the jet. We only specify buoyancy
  > forcing at the surface and southern sponge and forcing and
  > sponges are zonally uniform, so there should be no source of
  > zonal momentum. Therefore the jet is not there because of the
  > forcing and there is no incompatibility between the forcing
  and
  > boundary condition. Even if the jet were formed by some
  > turbulent inverse cascade, the tangential flow should go to
  zero
  > at the boundary due to no slip. What we see, however, is that
  > the normal derivative of the tangential velocity goes to
  zero,
  > suggesting that the model is giving us free-slip boundary
  > conditions even though we asked for no-slip boundary
  conditions.
  > —
  > Reply to this email directly or view
  > it on GitHub.
  >
  > —
  > Reply to this email directly or view it on GitHub
  >

https://github.com/NOAA-GFDL/MOM6/issues/230#issuecomment-146691422.

  -- 
  Dr. Stephen M. Griffies
  NOAA Geophysical Fluid Dynamics Lab
  201 Forrestal Road
  Princeton, NJ 08542
  USA
  —
    Reply to this email directly or view
      it on GitHub.

[adcroft]

Good point - I forget about these small tests too often. The 'seamount' example also tests this as well as demonstrates the absence of spontaneous motion in several coordinate systems over topography.

A forced Taylor cap over a seamount might be a good candidate for exercising sponges? I vaguely recall someone using that test in their thesis...