Multi-layer shallow water model

[dhruvbhagtani]

Hi all,

I'd be very keen to implement a multi-layer version of the shallow water model. For example, here's the 2-layer version:

My ultimate ambition is to solve this set of equations that include a diapycnal term S_int:

where Q_net is a prescribed heat forcing at the surface. So, from what I understand and from what I discussed with @navidcy, if we have a multi-layer shallow water model like eqs. 56 above then I can implement the set of eqs with the diapycnal terms added forcing to each individual fluid layers. Does this sound right?

Is enhancing the ShallowWaterModel to have multiple layers feasible?

cc @francispoulin, @glwagner, @navidcy, @AndyHoggANU, @rmholmes

[dhruvbhagtani]

btw, forgot to mention that above: Uₙ = hₙ * uₙ and Vₙ = hₙ * vₙ. That is, the eqs are written in conservative form, similar to https://clima.github.io/OceananigansDocumentation/stable/physics/shallow_water_model/

[iuryt]

I don't have enough experience to say if this is feasible, although I have this impression that there is always this nothing-is-impossible atmosphere around Oceananigans. I think this is an awesome idea!!!

[glwagner]

@iuryt speaks truth... it's hard to say something isn't "feasible". However, I do think that polishing off an n-layer extension of the existing ShallowWaterModel would require some time. I think you want to write not only the continuous equations, but also the finite volume discretization (and there's some issues around discretization of the advection operator now see #1866).

There are some design questions regarding the user interface and the use of existing grids. One possibility is to use the existing grids (which all have a z-coordinate) to specify the number of layers (via the vertical size), but to ignore the vertical grid. Then users would also have to input the buoyancy of each layer (I think I'd prefer buoyancy to density, if that's possible...) But another design might be possible that "re-uses" whatever users provide for the vertical grid.

I just took a look at the code and was surprised to notice that all of the kernels and fields are already three-dimensional. That means that once the user interface is figured out and appropriate model modifications are made, the "only" thing left might be to generalize the kernels...

I'm not sure how you do immersed boundaries, or how you blend "explicit" bathymetry (added through the pressure gradient) and immersed boundaries. Does anyone know?

@apaloczy may also be interested in this discussion.

[francispoulin]

Great idea @dhruvbhagtani and I would be very happy to help with this.

It is definitely feasible and something that would be great to have, but @glwagner raises some good questions in terms of how to do this. If we can do 2-layer than n-layer is not very far away, and I think it would be good to build your model with the future in mind.

A few thoughts:

Geoff Vallis' has the n-layer equations in his textbook both with and without a free-surface. A free-surface model is what you mentioned in hour example, and what we currently have in the one-layer case, so probably best to focus on that.

Topography that can interact more than one layer is complicated because then we need wetting and drying, but it won't need immersed boundaries. At least i don't think so. Wetting and drying usually requires a positive preserving advection scheme. We have first-order upwinding (which I don't think anyone has every used), that is positive preserving and could do this. But it's not a great scheme. A limiter would be better. There are a lot of methods out thatk, but not the first thing to worry about. This if much longer term, I would say.

I agree with the idea of using the z coordiante for the layer number. That would be very convenient.

I think this is exciting and probably wouldn't be too long to code up after we decide how to set things up.

[glwagner]

Topography that can interact more than one layer is complicated because then we need wetting and drying, but it won't need immersed boundaries. At least i don't think so.

I think I see what you're saying. h just goes to 0 where a layer impinges on a boundary (including an island), we don't have to mask anything.

I agree with the idea of using the z coordiante for the layer number.

Seems like the main annoying thing is that Δz might be irrelevant (despite that it exists).

I've thought about this a bit in the context of eventually supporting generalized vertical coordinates. Some codes use the name "r" instead of "z", which they believe grants them license to interpret "r" as density (or buoyancy) rather than height. We can also "just" interpret "z" as buoyancy. Then every grid and model property is meaningful, which I like... but on the other hand, they may not have the meaning you expect, if the name "z" is hard coded in your brain.

[francispoulin]

Exactly! We have to allow for h getting close to zero, maybe even zero, but never anything negative.

Using this for the density seems like a good way to be efficient, but then we also have to make sure how it's used in the model Seems doable, one just has to be careful and worry about the details. I need to think more about Δz, but thanks for sharing your thoughts.

[glwagner]

Using this for the density

density buoyancy

[navidcy]

is buoyancy in each layer bₙ = - g ρₙ / ρ̄, where ρ̄ is the mean density?

[navidcy]

there is another issue similar to bathymetry: what happens when an interior layer outcrops

[glwagner]

I believe you'd end up ρ̄ being a reference density (which is mechanically irrelevant, but not thermodynamically irrelevant necessarily). But I could be wrong (I'm not sure why people would choose density if all else were equal, so maybe there's a good reason...)

[francispoulin]

In shallow water models we don't usually talk about buoyancy but you could I suppose. Since density is what appears in the momentum equation, that's what's typically used. I believe Greg was saying buoyancy since that's the variable we typically use in the other models, but I could be wrong.

Yes, you can have outcroppings at the surface. This can even happen in a one layer case. Imagine starting out with a one layer reduced gravity shallow water model that is in the shape of an inverted U. When perturbed, the interface will move and then you have to deal with height going to 0 and also becoming positive.

There are positive preserving schemes for WENO that we can code up and test with our current model to better understand how it works before moving to multiple layers. Again, very happy to talk about this too.

[francispoulin]

I believe you'd end up ρ̄ being a reference density (which is mechanically irrelevant, but not thermodynamically irrelevant necessarily). But I could be wrong (I'm not sure why people would choose density if all else were equal, so maybe there's a good reason...)

In shallow water we only have density appearing in front of the pressure gradient. There is no buoyancy term as we determine the hydrostatic pressure and plug that into the pressure.

[glwagner]

Buoyancy might not be worth the hassle...

I just know that density drops out of the Boussinesq equations; if we then convert to buoyancy coordinates and discretize in buoyancy space, we might obtain something identical or similar to the layered shallow water equations.

We've written Oceananigans to deal with buoyancy rather than density with great intention: a "reference density" does not always appear in every setup / numerical experiment (for example: using buoyancy itself as a tracer), and when it does appear it's only meaningful in a few specific places (equation of state, thermodynamic calculations at the surface). This not only clarifies the physics (eg what role does the reference density play?) but also helps reproducibility (we shouldn't have to report / set a reference density for experiments where its dynamically irrelevant). Yes, it's different from how we normally think about the problem, but there's always room for progress...

I'm not sure if the layered shallow water benefits from this philosophy too.

[apaloczy]

Also one thing that is important to get right in multi-layer shallow water models is layer-wise conservation of mass and other properties (see e.g., https://aronnax.readthedocs.io/en/latest/verification.html and http://www.nordet.net/beom.html).

Topography that can interact more than one layer is complicated because then we need wetting and drying, but it won't need immersed boundaries. At least i don't think so. Wetting and drying usually requires a positive preserving advection scheme. We have first-order upwinding (which I don't think anyone has every used), that is positive preserving and could do this. But it's not a great scheme. A limiter would be better. There are a lot of methods out thatk, but not the first thing to worry about. This if much longer term, I would say.

In addition to positive-preserving advection schemes, it seems one other way to deal with this is adding to the pressure gradients an extra term that prevents the thicknesses of outcropped/incropped layers from reaching zero following Salmon (2002). There is a deep discussion on the physics and practical implementation of this in different validation simulations in http://www.nordet.net/etc/doc_beom.pdf. Just wanted to put this out there, I don't have a good feel for which approach would be most desirable.

[francispoulin]

Buoyancy might not be worth the hassle...

I just know that density drops out of the Boussinesq equations; if we then convert to buoyancy coordinates and discretize in buoyancy space, we might obtain something identical or similar to the layered shallow water equations.

We've written Oceananigans to deal with buoyancy rather than density with great intention: a "reference density" does not always appear in every setup / numerical experiment (for example: using buoyancy itself as a tracer), and when it does appear it's only meaningful in a few specific places (equation of state, thermodynamic calculations at the surface). This not only clarifies the physics (eg what role does the reference density play?) but also helps reproducibility (we shouldn't have to report / set a reference density for experiments where its dynamically irrelevant). Yes, it's different from how we normally think about the problem, but there's always room for progress...

I'm not sure if the layered shallow water benefits from this philosophy too.

I don't know if using buoyancy instead of density affects the layered shallow water model much. One way to think of the model is that the density/buoyancy is piecewise constant, which is why one is just a scalar multiple of each other. We use hydrostatic balance to determine the pressure in terms of density/buoyancy and we can then use either variable.

If we are going to use buoyancy in the n-layer problem, then we might want to start with the one-layer version. I prefer to think of it as a reduced gravity shallow water model, since that is more general. In that context we define the reduced gravity as g' = Δρ g/ ρ₀, where Δρ is the change in density and ρ₀ is a reference density. It seems that the reduced gravity is proportional to the change in buoyancy between two layers. Not sure if that helps at all, but its an observation.

[francispoulin]

Also one thing that is important to get right in multi-layer shallow water models is layer-wise conservation of mass and other properties (see e.g., https://aronnax.readthedocs.io/en/latest/verification.html and http://www.nordet.net/beom.html).

Topography that can interact more than one layer is complicated because then we need wetting and drying, but it won't need immersed boundaries. At least i don't think so. Wetting and drying usually requires a positive preserving advection scheme. We have first-order upwinding (which I don't think anyone has every used), that is positive preserving and could do this. But it's not a great scheme. A limiter would be better. There are a lot of methods out thatk, but not the first thing to worry about. This if much longer term, I would say.

In addition to positive-preserving advection schemes, it seems one other way to deal with this is adding to the pressure gradients an extra term that prevents the thicknesses of outcropped/incropped layers from reaching zero following Salmon (2002). There is a deep discussion on the physics and practical implementation of this in different validation simulations in http://www.nordet.net/etc/doc_beom.pdf. Just wanted to put this out there, I don't have a good feel for which approach would be most desirable.

Thanks for the reference @apaloczy ! I know of Rick's work but didn't know whether people were using it very much. The reference seems like it will shed some light on the matter. This might very well be worth playing with to see what approach we prefer.

[kburns]

Just to add a reference, Kevlahan et al. 2015 doesn't discuss multiple layers, but seems to do a good job incorporating bathymetry and immersed-boundary continents that have proper wave reflection and don't increase wave speeds inside the continents. I'm trying this model out in Dedalus and maybe it's possible to extend this approach to multiple layers, but I'm not a shallow water expert.

[francispoulin]

Thanks @kburns for the reference.

Nicholas Kevlahan and I are collaborating on a problem now and we have used his wavelet based model. I agree that it does do a lot of things very well. It deals with topography and coastlines using Brinkman penalization, which is quite distinct from our immersed boundary method. I don't know how he preserves positivity in the height field but I will find out.

Update: their model uses a vertical coordiante in the vertical to avoid layer collapse. Also, in the multi-layer case they use a discretization of vertical diffusion, to also prevent layer collapse (depths going to zero).

[navidcy]

@glwagner & @dhruvbhagtani why don't we have a zoom and see how we can get things started and then continue the discussion either here or in a PR?

[francispoulin]

I am also happy to contribute if/when I can.

[dhruvbhagtani]

@iuryt you're right, the environment here is really amazing!

I will think and write down some FV discretisation equations by the time we zoom in. I'm not the best person to talk on how we should implement it in Oceananigans in terms of using z-coordinate as different levels, but it seems like a great and innovative idea.

Perhaps I'm asking a silly question but it's not clear to me that the equations I've written up (and I'll double check this with someone else too) seems to consider varying bathymetry, and I'm not sure if Oceananigans explicitly models bathymetry by adding a pressure gradient, and whether we would be double counting the effects of bathymetry using the equations I shared above.

I agree with @francispoulin, maybe it doesn't matter much whether we use density or buoyancy as our vertical coordinate. I do like the idea of using buoyancy instead of density, but that's just a personal preference.

[glwagner]

I'm not sure if Oceananigans explicitly models bathymetry by adding a pressure gradient, and whether we would be double counting the effects of bathymetry using the equations I shared above.

ShallowWaterModel is totally independent from the other models --- so this is up to how ShallowWaterModel is written. But in a shallow water system it seems (based on the discussion on this thread) that bathymetry is treated via pressure gradients (and also vanishing layers or other techniques to handle outcropping). So we won't use ImmersedBoundaryGrid, which is how complex domains are handled elsewhere in Oceananigans.

[navidcy]

Regarding vanishing layers: what MOM6 does (I think) is that layers never vanish but they are restricted to have height of 1e-10.

[francispoulin]

I believe that Rick Salmon's method essentially does the same thing. It never quite goes to zero but gets really small.

[glwagner]

The notes are crazy because 1 is surface and 2 is bottom.

[navidcy]

After #2522 me and @dhruvbhagtani will start working on the multi-layer. The plan is to extend the ShallowWaterModel to allow non-flat z dimensions with Nz the number of fluid layers.

It would be good if we add a regression test for the single-layer ShallowWaterModel to ensure we don't break things in the process.

[simone-silvestri]

The Bickley jet might be a good example to use as a regression test!

[navidcy]

Yeap, I agree.

So are you thinking a config similar to the examples (e.g., rectilinear), start with mean state + some prescribed perturbation and run for up to, e.g., t=20 and compare output? How does that sound?

[francispoulin]

Starting with a rectilinear grid makes sense, but it might also be fun to try a lat-lon grid.

[navidcy]

ok, now that #2522 is merged we can start thinking of multiple layers :) I put an overleaf doc with eqs at https://www.overleaf.com/read/mtyjxnnrjpqv

[francispoulin]

The notes look great!

One question though. In your definition for the reduced gravity, why not divide by $\rhoj$ instead of $\rho{j+1}$, since that's what naturally apprears in the pressure graident term? I know that because the densities are almsot the same, it won't matter much, but it is not clear why you make this approximation when you don't need to.

[navidcy]

That was probably just a typo :)

Fixed it!