Add free-slip, partial-slip and no-slip momentum BCs

[dengwirda]

This PR updates the way horizontal momentum boundary conditions are computed --- highlighting that the current MPAS discretisation applies no-slip conditions, introducing new free-slip and partial-slip conditions, and cleaning-up the way cell-vertex remapping is handled across partially masked cells/duals in general.

No-slip and free-slip BCs can be expressed in terms of $\nabla \times \mathbf{u}$ on the boundaries --- no-slip conditions (the current implementation) correspond to a 'full' relative vorticity computed on boundary vertices, while free-slip conditions are imposed by setting $\nabla \times \mathbf{u} = 0$ on boundary vertices (corresponding to no induced torque when a free-stream velocity flows tangential to a wall). I've also introduced so-called 'partial-slip' conditions (see e.g. NEMO) as an interpolation between the no-slip and free-slip cases, leading to a reduced vorticity on boundaries. See e.g. Roullet and Gaillard for confirmation of how momentum BCs should be applied in a vector-invariant formulation.

A new lateral_walls namelist section for BC settings is introduced, where config_wall_slip_factor controls the amount of slip at horizontal walls (both coastlines at the surface and masked boundaries against bathymetry in general). config_wall_slip_factor = 0.0 is a no-slip condition, config_wall_slip_factor = 1.0 is a free-slip condition, and values between 0.0 and 1.0 correspond to partial-slip cases.

No-slip (left) vs free-slip (right) behaviour can be seen in the channel case below (developed with @hetland) in which the initially uniform zonal flow generates strong vorticity sheets at the north/south boundaries in the no-slip case, while the flow remains quiescent in the free-slip case.

This PR also updates the way vorticity and layer-thickness is computed and remapped at partially masked cells and duals at boundaries, ensuring that the numerical stencils reflect the masking. In the case of computing potential vorticity $\frac{1}{h}(\nabla \times \mathbf{u} + f)$ it's necessary to remap layer thickness from cells to duals (triangles) using the intersecting cell-dual 'kite' areas. Currently, I feel we are not accounting for masked cells/duals correctly. Consider a case at a boundary in which one of the three cells surrounding a boundary vertex is masked, and where a spatially uniform flow with layer-thickness $h$ exists inside the domain.

In the current implementation the thicknesses from the two valid cells is remapped along with a 0 value for the masked cell. The result is divided by the full area of the dual (triangle) resulting in a vertex layer-thickness of only $\frac{2}{3} h$, inconsistent with the spatially uniform interior thickness. An alternative approach implemented here is to accumulate only the unmasked area of the dual (triangle) and use this in the remapping, which will results the expected vertex thickness of $h$.

$$h_{v_{old}} = \frac{A_1 h_{c_1} + A_2 h_{c_2} + 0}{A_{tria}} = \frac{A_{tria}}{A_{tria}} (\frac{1}{3}h_{c_1} + \frac{1}{3}h_{c_2}) = \frac{2}{3}h$$

$$h_{v_{new}} = \frac{A_1 h_{c_1} + A_2 h_{c_2}}{A_1 + A_2} = \frac{A_{tria}}{A_{tria}} (\frac{\frac{1}{3}h_{c_1} + \frac{1}{3}h_{c_2}}{\frac{1}{3} + \frac{1}{3}}) = h$$

$$A_{1,2,3} = \frac{1}{3}A_{tria}$$

I've tested both in channel and global standalone spin-ups. Results are non-BFB, with somewhat increased kinetic energy developing in global spin-ups with free-slip BCs due to the reduced dissipation compared to no-slip lateral boundaries.

• Roullet, G. and Gaillard, T., 2022. A fast monotone discretization of the rotating shallow water equations. Journal of Advances in Modeling Earth Systems, 14(2), p.e2021MS002663.

• https://www.nemo-ocean.eu/doc/node58.html

[xylar]

@dengwirda, same questions as for #48: Is this ready to review? If so, who would you like to review it before it goes to E3SM?

[dengwirda]

@xylar, let me know how you'd like the 'discussion' aspect of this to proceed --- this branch is ready to run, but I'd say that additional testing is required before review.

[xylar]

@dengwirda, typically these PRs on E3SM-Ocean-Discussion only proceed if you ping specific people to indicate that you want their feedback. You can do that by adding them as a reviewer or just by @ing them and saying what kind of feedback you're looking for. Neither of those typically requires the code to be tested, though it could make things a bit clearer if you just ping people to initiate discussion and then flag them for a review when things are tested. The review here is typically cursory -- are things far enough along that we don't think we'll have a lengthy discussion (both in time and number of posts) on E3SM?

[dengwirda]

@xylar please feel free to add extra reviewers and/or review yourself here.
The slip BCs are hopefully not controversial, though the implications of using them in different config.'s is something to discuss.
The updates to layerThicknessVertex at boundaries will/may change behaviour in general.

[xylar]

Is boundaryVertex also needed in the present() arguments?

[xylar]

@dengwirda, your list of reviewers looks good to me.

[xylar]

@dengwirda, I suggest you make one PR to E3SM to fix dualArea as you have done here. I think we just want to do that once and for all, and it will be a non-BFB PR for sure. Then, we can do a second PR with the boundary conditions. That one can leave no-slip as the default for now but we can test the other options and consider switching to them before the v3 code freeze.

[mark-petersen]

@dengwirda this looks like a good addition. When you compare the current code with this PR using config_wall_slip_factor = 0.0, are the differences reasonably small? It that case, it looks like order-of-operations differences, so should be <1e-10 for the first few time steps. If that is the case, I think we could proceed to a PR to E3SM-Project with the default set to config_wall_slip_factor = 0.0.

[xylar]

@mark-petersen, I could be wrong but my understanding was that, because the area for vertices was being computed with full triangles rather than from the sum of the kites that actually exist, the differences will not be small and we will need to treat this as a (hopefully) non-climate-changing PR. That is why I was suggesting that we make the bug fix for the vertex area first and then the rest.

[mark-petersen]

@xylar I see what you are saying. This is the change where invAreaTri1 includes all three kites, while areaDual only includes the kites from active cells. This would be answer changing, and the changes in the loop above could go into a single PR.

[mark-petersen]

The changes in this loop are what is needed to include wall_slip_factor. For this section, it is true that wall_slip_factor = 0.0 should only change the solution due to order of operations. For a second PR for just the addition of wall_slip_factor and these lines, it will be easier to show that only minor changes occurred.

[sbrus89]

@dengwirda, is there a reason to exchange loop order here? It seems like the order of operations issue could be avoided by applying the (1.0_RKIND - wall_slip_factor * boundaryVertex(k, iVertex)) factor in a separate loop k loop following the original do i = 1, vertexDegree loop.

[mark-petersen]

I think the invAreaTri1 here needs a correction. The logic is similar to the areaDual below, but opposite. Currently, next to a full land boundary, areaTriangle(iVertex) does not include the non-existent cells. But beside a stair-step boundary below water where the neighbor exists, areaTriangle(iVertex) is the full area. So here invAreaTri1 is treating coastal boundaries versus underwater boundaries in a different way. Thinking about it, it seems each active edge should accumulate an associated area with it (the sub-triangle about that edge), and we should divide by that accumulated area.

[xylar]

I think you're suggesting exactly the same correction as what happens below with areaDual and that one is needed for the same reason -- handling topography.

[xylar]

I wonder, though, if the desired behavior is that the non-existent neighboring cells contribute zero vorticity, rather than that they don't contribut. And thus we would want to divide by the full triangle area, which is unfortunately not available for coastal vertices after culling.

[xylar]

Do you think it's accurate to think of the actual boundary as the one connecting the middle of each edge where normalVelocity = 0 (except when the edge is parallel to the boundary)?

I think the boundary really follows the boundary edges so it connects vertices, not edge centers.

[xylar]

If so, then we would only want relativeVorticity = 0 at the vertices with all boundary edges and relativeVorticity != 0 for vertices with one non-boundary edge because those vertices "stick out" into the flow.

To me, the vertices with 2 boundary edges and one sticking out into the flow are precisely the ones that implement the choice of boundary conditions. Whether that free edge is allowed to contribute to the curl or not determines the boundary condition.

[xylar]

I'm wondering if we should be aiming for an implementation that would provide the same shear in the case where a straight boundary has a sawtooth shape (2 sides of hexagons) and one where the boundary has a snaking shape (3 sides of hexagons).

I think that's a really good question, and also whether either of these can be expected to approach the straight-wall solution at very high resolution.

[cbegeman]

If the boundary really follows the vertices, it seems unlikely that we could get both wall geometries to approach the no-slip solution. We could test the full dual area and the partial dual area (@xylar's suggestion for area lying in the active mesh) and see if either leads to better agreement. Given that, as @dengwirda mentions, most ocean models use free slip, we should leave the no-slip development there and move on to testing the free-slip case in realistic configs.

[xylar]

I totally agree that focusing too much on getting no-slip right is silly if free-slip is the way to go. Free slip seems much easier to get right.

[dengwirda]

@mark-petersen, @xylar, yes I think you're both right here --- the wall_slip_factor change by itself should hopefully be non-climate-changing (though perhaps not BFB due to order-of-op type differences) but the boundary vertex thickness issue will change the way vorticity is computed in general --- at all coastlines for the surface layer, and at any cells/duals where layers are maksed against bathymetry in the deep ocean.

One reason I've grouped this work together here is that implicitly both changes relate to the way momentum boundary conditions are imposed. The wall_slip_factor should allow users to interpolate between no-slip and free-slip conditions explicitly, but the vertex thickness change will alter the way PV is computed on boundaries in all cases. In the example I have above, it looks like we may currently be computing too small a boundary vertex thickness (due to the masked kites), which would give too large a boundary PV value. Considering that free-slip has PV=0 on boundaries, and that no-slip uses a full PV value, this may mean that some kind of more-than-no-slip condition is being used currently. What do you think?

[xylar]

@dengwirda, I see the logic in including them together.

I think it will still be simpler in terms of how to frame these changes to the broader project if we first introduce a bug fix to the area weighting of layer thickness on vertices, knowing that the changes will not be small but that it's a bug. Then, we can still point to that PR when we introduce the other changes here but we would point out that the change in behavior after the bug fix is small as long as you are using the no-slip boundary condition. We would further point out that the "stealth" option to have free-slip or partial slip boundary conditions would be expected to introduce changes of a similar magnitude a and in a similar direction to the previous bug fix. If testing is sufficient at that point, we could also advocate for making something other than no-slip the default at that point.

But I think it really will be necessary to introduce these changes in those 2 stages, rather than in one PR.

[mark-petersen]

I'm happy for these changes to move to E3SM-Project. Thanks.

[xylar]

After over 2 years, I have had a chance to return to thinking about this.

In the PR description, Darren states:

In the current implementation the thicknesses from the two valid cells is remapped along with a 0 value for the masked cell. The result is divided by the full area of the dual (triangle) resulting in a vertex layer-thickness of only $2/3 h$, inconsistent with the spatially uniform interior thickness. An alternative approach implemented here is to accumulate only the unmasked area of the dual (triangle) and use this in the remapping, which will results the expected vertex thickness of $h$.

I do not think this is correct. Instead, I think the poorly named areaTriangle is actually the sum of the areas of the valid kite areas. Here is an example from the oQU240 mesh:

I therefore believe that the changes to the relative vorticity calculation in this PR are not needed. We could make a note that areaTrinagle is not, in fact, the full area of the triangle on a vertex but rather the sum of the valid kites. But this need not be recomputed.

[cbegeman]

@xylar Thank you!

I just wanted to document that @mark-petersen, myself, and Alaina Hogan are testing this branch and will also test a version without the dual mesh area changes. There is no need for anyone to make changes to this PR right now.

[dengwirda]

@xylar and all re: the dual-area issue, I'd say that it's worth considering what happens throughout the water column (where land cells are present in the mesh but masked) rather than just at the surface (where cells are fully culled out).

It looks like I have both the wall_slip_factor and layerThicknessVertex in the same PR here, though they can be thought about separately. Both relate to how vorticity is generated at boundaries though, and the consequences for energy + enstrophy balance.

[xylar]

@xylar and all re: the dual-area issue, I'd say that it's worth considering what happens throughout the water column (where land cells are present in the mesh but masked) rather than just at the surface (where cells are fully culled out).

@dengwirda, ah, of course. That's a very good point.

[cbegeman]

My preference for a path forward would be to move the areaDual computation to mpas_ocn_mesh.F and save it as a variable called something like areaVertex with dimensions nVertLevels, nVertices. To be super clear, areaVertex would be the sum of only the active cell areas.

We would use areaVertex for relativeVorticity = circulation / areaVertex and layerThicknessVertex = layerThicknessVertex / areaVertex

How does that sound?

I also want to confirm that we should be determining active cells on a vertex using cellMask. circulation is summed over minLevelVertexTop : maxLevelVertexBot so I think cellMask would reflect that (maxLevelVertexBot being the one that is inclusive of inactive neighbors).

[xylar]

@cbegeman, sorry for waffling on this. I made a note about about why I think the area to use in the circulation is not the same as areaVertex. If you agree with my reasoning, I am not sure it makes sense to store areaDual after all, since it will only be used to interpolate layer thickness to vertices and not in the denominator of the relative vorticity calculation.

[dengwirda]

All --- a few extra thoughts:

• In my model these days I've implemented things similarly to what @cbegeman mentions above --- with the hexagonal mesh being 'the mesh' at each layer, and control volumes associated with vertices and edges clipped to the boundary of the active hexagons. No triangles (or edge CVs) are allowed to stick out into the land essentially.

• Differently to this PR (and along the lines of the discussion in the other tab here) these days I use these clipped CVs (and their fractional areas, etc) for everything --- both computing curl as well as remapping variables between cells, vertices, edges, etc. Based on the clipping above, the set of CVs for all entities should always occupy the same overall area, and use of the fractional areas means that all of the transfers from one entity to another is an integral remapping that satisfies conservation.

• There is some work of ours looking at mesh convergence with different BCs on some standard test cases with shear layers and whatnot. That paper isn't out yet though. I'd say the main impact is the free- vs no-slip aspect (the wall_slip_factor here). Understanding the correct boundary CV treatment is something I'm interested in and happy to discuss though.

• A quick word about why I started this originally --- most GCMs use free-slip BCs, whereas MPAS has always been no-slip. With no-slip, there's a shear layer generated against bathymetry that tends to spool-up a lot of vorticity / eddies and tends to dissipate energy. The width of the shear layer is based on viscosity and unless there's enough resolution to resolve it (pretty unlikely in the case of a global model) no-slip seems like the wrong choice to me. Of course, whether the momentum BCs have a big impact on global biases is hard to say, but it's another thing I'd be interested to understand better.

[xylar]

Thanks for staying engaged with this, @dengwirda ! We'll look forward to your paper. I think we're stuck with our jagged boundaries for the time being but what you describe sounds really promising for a future improvement!

[cbegeman]

That sounds like an interesting approach, @dengwirda. Please do share your paper when it is available!

[cbegeman]

Suggested change

	!$acc private(invAreaTri1, i, iEdge, k, r_tmp, wall_slip_factor)
	!$acc private(invAreaTri1, i, iEdge, k, r_tmp)

This change was made prior to testing. At least on some compilers, including wall_slip_factor in private results in it being assigned to 0 in the loop.

[dengwirda]

Fixes a bug of mine.

[cbegeman]

Idealized testing

Barotropic gyre:

• NCC wall slip = 0. No slip case: differences in streamfunction between master and wallslip = 0 are O(1e-9) after 2 years
• CC wall slip = 1. Free slip case: error near the boundary is reduced from master to wallslip = 1 but error away from the boundary is increased. This may not be a cause for concern because the analytic solution is only accurate in the limit of $\nu_2 / L_x \beta = \infty$
  • Master

• • Wallslip = 1

Barotropic channel: with z-level and trough geometry to test free-slip when deeper levels terminate against bathymetry

• NCC wall slip = 0. Circulation and thus relative vorticity remains non-zero at boundary vertices.
  • Master

• • Wallslip = 0 zTop

• • Wallslip = 0 zBot

• CC wall slip = 1.
• • Wallslip = 1 zTop

• • Wallslip = 1 zBot

[cbegeman]

E3SM testing

v3.LR.piControl, 50 years, climo 1-50
https://e3sm.atlassian.net/wiki/spaces/pd/pages/5545886223/v3.LR.piControl-wallslip1
Increased throughflow in narrow passages:

No significant diffs in any other fields:

v3.SORRM G-case, 20 years, climo 11-20 (free-slip conditions on left with this branch, control with alfred3 branch on right)
https://e3sm.atlassian.net/wiki/spaces/pd/pages/5545885697/20251002.v3.SORRME3r3.GMPAS-JRA1p5-DIB-PISMF.wallslip1
No significant diffs in SST:

Significant diffs in seafloor velocity magnitude statistics (see upper right for min, mean, max):

Too early to say if AMOC diffs are significant:

Throughflow in narrow passages increased:

No significant sea ice diffs:
https://portal.nersc.gov/project/m2136/bin/iice/mpas-a.cgi?url1=https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.cbegeman/BlueTip/20251002.v3.SORRME3r3.GMPAS-JRA1p5-DIB-PISMF.wallslip1/mpas_analysis/ts_0001-0020_climo_0011-0020&label1=wallslip1&url2=https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.dcomeau/E3SMv3_dev/20250918.v3.SORRME3r3.GMPAS-JRA1p5-DIB-PISMF.alfred3/mpas_analysis/ts_0001-0020_climo_0011-0020&label2=alfred3&category=sea_ice&component=sea_ice&plot=&diff=0&dconfig=cryo&

[cbegeman]

I'd like to migrate this over to an E3SM PR. If anyone has any objections or thinks further testing is needed before that review stage, let me know before Wed, Nov 5.

[xylar]

Sorry for the slow response. Yes, go for it!