`ShallowWaterModel.bathymetry`?

[navidcy]

Is this just a placeholder? Or does indeed ShallowWaterModel.bathymetry supports non-trivial bathymetry?

[francispoulin]

At the moment it is just a place holder but just created issue #1716 that will create a test for non-trival bathymetry.

Thanks @navidcy for asking about this.

[glwagner]

@francispoulin ?

[francispoulin]

This still needs to be done. I will raise it on list of priorities of my to-do list. It will be very easy when it's done.

[apaloczy]

Just found this issue while setting up a ShallowWaterModel experiment with nonzero bathymetry. Implementing this involves modifying src/Models/ShallowWaterModels/solution_and_tracer_tendencies.jl to include the terms associated with bathymetric gradients in x_pressure_gradient and y_pressure_gradient, correct? What else would be needed? I would be happy to help write a validation test case and/or an example for this.

[glwagner]

I think that's right --- it's essentially a matter of applying the transformation h -> h + b (where b is bathymetry` in the governing equations here:

https://clima.github.io/OceananigansDocumentation/stable/physics/shallow_water_model/

?

I think this might be achieved two ways:

1. Add new terms to the tendencies; or
2. Define the AbstractOperation

total_depth = solution.h + bathymetry

in the model constructor and add this property to ShallowWaterModel. Then, where appropriate in the tendency kernel functions, use total_depth rather than solution.h. I think the two appropriate places are the pressure gradient term and the mass flux term in the continuity equation. Method 2 might actually be numerically different than method 1 when using WENO advection for the mass flux (eg, we are interpolating a different function with different smoothness properties)... I'm not sure which is better (but something tells me that method 2 might be better?

As for how to handle the bathymetry I think something similar to what we do for immersed boundaries would be nice:

Oceananigans.jl/src/ImmersedBoundaries/grid_fitted_immersed_boundaries.jl

Lines 43 to 51 in 3322f18

	function ImmersedBoundaryGrid(grid, ib::GridFittedBottom)
	arch = grid.architecture
	bottom_field = Field{Center, Center, Nothing}(grid)
	set!(bottom_field, ib.bottom_height)
	fill_halo_regions!(bottom_field)
	offset_bottom_array = dropdims(bottom_field.data, dims=3)
	new_ib = GridFittedBottom(offset_bottom_array)
	return ImmersedBoundaryGrid(grid, new_ib)
	end

eg, in the constructor for ShallowWaterModel, something like

if !isnothing(bathymetry) # user has provided non-default bathymetry
    bathymetry_field = Field{Center, Center, Nothing}(grid)
    set!(bathymetry_field, bathymetry) # works for functions and arrays
    total_depth = solution.h + bathymetry_field
else # there's no bathymetry
    total_depth = solution.h
end

# Later, make sure to use `bathymetry_field` when instantiating `ShallowWaterModel`.

@apaloczy that would be pretty awesome if you want to tackle this! Happy to provide help along the way if you open a PR when you start down this road. I can also help with docs.

[francispoulin]

That's great that you want to do this @apaloczy ! I am also happy to help if I can.

I think that @glwagner explained it very well.

One question I will raise is, do we prefer to have to have height as a variable, or the free surface, which is more consistent with the hydrostatic model?

[apaloczy]

Thanks @glwagner and @francispoulin, I will give this a try and ping you both along the way. @francispoulin, I feel that having the free surface elevation rather than total height would be preferable too.

About option 2, I am missing how introducing the h -> h + b transformation would work for the pressure gradient term. The governing equations only apply for flat-bottom cases the way they are currently set up, right? I ask this because the total thickness h = H + 𝜂 - b is used instead of the free surface 𝜂 in the pressure gradient terms and -g∂ₓ𝜂 = -g∂ₓ(h + b - H) = -g∂ₓh only if b = 0. Here H, 𝜂 and b are the mean thickness, the free surface and the bathymetry, respectively.

So with bathymetry we would have -g∂ₓ(h + b - H) = -g∂ₓh -g∂ₓb, and in conservative form, -g∂ₓ(h²/2) -gh∂ₓb, right? so we would end up with this extra term proportional to the local bottom slope (as in e.g., FourierFlows/GeophysicalFlows.jl#163 (comment)).

[francispoulin]

Your are correct @apaloczy that we need to be careful. One thing I should point out is that we are solving the equations in conservation form, so it is a bit differente than the link you cited. If you look at equations 7 to 9 in this paper then you will see how this is written in conservation form. Note that the syntax is h is the total depth and h_B is the height of the topography.

Also, this paper simulates the instabilty of a geostrophically balanced Bickley jet over sloping topogrpahy. We could take the ShallowWaterModel example and modify the parameters slightly to reproduce the results in this paper, if we wanted some validation.

[glwagner]

I think the main reason to use total height is that it generalizes to a "stacked" shallow water model (where there are N layers rather than just 1). But perhaps the choice depends on whether shallow water model is valuable mostly as a stand-alone model for physics problems or whether it's intent is more as a testbed for develping numerical methods, etc (with perhaps the eventual possibility of generalizing to N layers).

[francispoulin]

It's probably easier to keep things in terms of h for the moment, and revisit this question later, if there is any interest in changing it.

[glwagner]

Right. Another reason to keep h is that someday we will hopefully have a nonlinear free surface with the hydrostatic model. In that case the physics overlap strongly, and so ShallowWaterModel will need to distinguish itself somehow so we can motivate continuing to maintain it.

[glwagner]

with bathymetry we would have -g∂ₓ(h + b - H) = -g∂ₓh -g∂ₓb, and in conservative form, -g∂ₓ(h²/2) -gh∂ₓb, right?

Oh, this seems right. Apologies for the confusion.

I'm still not sure we want to analytically integrate one of the pressure gradient terms h dx(h) = dx(h^2/2) but leave the other unchanged. If we treat them both similarly then we would write the pressure gradient term

-g h * ∂ₓ(h + b)

It might be worth testing both possibilities to see if one has favorable numerical properties?

Perhaps this is related to the whole issue re: well-balanced methods, etc... ?