Better way to smooth bathymetry than "repeated interpolation"?

[glwagner]

Right now we rely on a sort of crude and uncontrolled method to smooth bathymetry, eg from the bathymetry generation example

h_rough = regrid_bathymetry(grid; interpolation_passes = 3)
h_smooth = regrid_bathymetry(grid; interpolation_passes = 40)

It doesn't really seem possible to predict what the effect of the "interpolation passes" is. Obviously more of them smooths, but how much do they smooth? How much do we need? Is there any way to understand what the net outcome is, regardless of the grid we are on?

I think this method was implemented not because it is a good method to use, but because it was convenient to simply "keep interpolating". Convenience isn't always the best motivation...

I think diffusion or a moving average might be a better way to smooth bathymetry. This would preserve the total ocean volume. Also, rather than specifying "number of passes" I feel like it would be more useful to be able to specify something like the "maximum gradient" or "maximum curvature" and iterate until that criterion is reached.

[NoraLoose]

I think diffusion or a moving average might be a better way to smooth bathymetry. This would preserve the total ocean volume. Also, rather than specifying "number of passes" I feel like it would be more useful to be able to specify something like the "maximum gradient" or "maximum curvature" and iterate until that criterion is reached.

That's basically what we do when we prepare the bathymetry for a ROMS simulation. We perform two consecutive steps, first a domain-wide smoothing (to avoid grid-scale instabilities) and then an iterative local smoothing until the maximum slope criterion is satisfied. For more details, see the ROMS-Tools documentation.

[NoraLoose]

The domain-wide smoothing step might be more relevant for a terrain-following coordinate (as in ROMS) than for a z-coordinate.

[navidcy]

While working on #475 and in an attempt to understand what the regrid_bathymetry method does, we were chatting with @glwagner in the morning trying to understand what the effect of the step-by-step interpolations do (eg interpolate into intermediate resolutions before arriving at your final resolution versus interpolating straight from high-resolution to final coarse resolution.) The figure at https://clima.github.io/ClimaOceanDocumentation/stable/literated/generate_bathymetry/ suggests that intermediate interpolations lead to horizontal smoothing...

We both think that that might be more reasonable is to apply the smoothing to the native_grid of the topography (at the highest resolution) and then interpolate to the user-desired resolution.

#477 is a first attempt to allow for new smoothing operations

[glwagner]

Curious @NoraLoose do you smooth the source bathymetry, or only the "regridded" field (eg the bathymetry that you can after interpolation / regridding)? I'm also curious whether we should naively interpolate, or use a proper conservative remapping, which conserve the average depth.

[NoraLoose]

We first regrid, then smooth the regridded field. I think this approach has a couple of advantages. Smoothing after regridding—using a fixed kernel size (say, 4–8 times the grid scale)—serves a few purposes:

• It reduces aliasing effects when regridding from very high resolution to a coarser grid.
• It ensures the final bathymetry is smooth at the grid scale, which is especially important for terrain-following coordinates (used in ROMS), but maybe less so for z-coordinates.

If you smooth the source bathymetry before regridding, you have less control over smoothness at the target grid scale, and you’re still at risk of aliasing—unless you’re very careful about how the initial smoothing is done, and how it’s tuned to the target resolution.

Overall, I think conservative regridding followed by conservative smoothing is a great strategy.

[NoraLoose]

For the sake of completeness, I will just say that we are actually performing three steps:

1. Regridding the source bathymetry to the target grid.
2. Domain-wide smoothing with a fixed kernel size.
3. Iterative local smoothing to reduce steep slopes until a steepness criterion is met. See also Section 4.2.1 in this paper.

Not sure how critical step 3 is for z-coordinate models—it’s essential for terrain-following coordinates to reduce pressure gradient errors, but maybe less so here.