Revisit bottom roughness generation

[minghangli-uni]

The existing bottom roughness calculation follows Jayne & St Laurent 2001,

The bottom roughness, h^2, is computed from the Smith and Sandwell [1997] 1/30° ocean topography database. Over each grid cell, a polynomial sloping surface is fit to the bottom topography (given by H = a + bx + cy + dxy), and the residual heights are used to compute h^2, the mean-square bottom roughness averaged over the grid cell", where the grid resolution is 0.5°.

but applying this approach directly on high-res ocean models is problematic. When roughness is computed at the model grid scale, the resulting h^2 becomes explicitly dependent on model resolution, which is unphysical. In the absence of explicit tidal forcing, resolved scale topography will not contribute to internal tide generation, so generation by resolved scales also needs to be parameterised. A clear summary of this issue is given by @aekiss in ACCESS-NRI/access-om3-configs#457

Below is the key part of that summary.

The key physical insight is that bottom roughness should instead be defined relative to the internal tide length scale, not the model grid spacing,

roughness should instead be calculated relative to the topography smoothed over a length scale of the lowest-mode internal tide wavelength, because roughness on finer scales than this is effective in generating internal tides. This is around 180km in the deep ocean, but depends on stratification and water depth, becoming shorter in shallow water. Importantly, this smoothing scale is independent of model grid resolution.

@CallumJShakespeare has provided a matlba implementation of this approach https://github.com/CallumJShakespeare/global_tide_inputs, that workflow computes roughness using:

• internal tide length scales derived from WOA18,
• MATLAB-specific routines (e.g. interp2).

To integrate this method into the ACCESS-OM3 workflow and update the stratification input to WOA23, this PR re-implements the method in Python with MPI parallelism.

Hence this script does following things,

• Compute the buoyancy frequency squared (N^2) from WOA23 T/S data
• From N2, computes the mode-1 internal tide wavelength (lambda1) at a given tidal frequency (e.g. M2)
• Smooth high-resolution bathymetry (SYNBATH) using a polar, Gaussian-weighted average with radius $\lambda1$, producing a physically meaningful mean depth.
• Compute depth variance relative to this smoothed depth, getting bottom roughness independent of model grid resolution.
• Interpolate the resulting roughness field onto model grids using esmf regridding.
• Complete metadata and provenance.

Notes:

• The computation involves interpolation over 670k ocean points, so a task-based MPI distribution is used. Work is distributed over valid ocean cells rather than latitude stripes. This avoids severe load imbalance due to land masks.
• The original matlab code uses interp2 for bathymetric sampling. But in this implementation bilinear interpolation is performed manually:
  • to minimise overhead in tight loops
  • to enforce strict NaN propagation (if any corner of the interpolation cell is nan, the result is nan), preventing land contamination of ocean roughness.

EDIT: 15 Jan 2026

• Regridding is performed separately by generate_bottom_roughness_regrid.py. This is intentionally split out due to a known xesmf + mpi4py issue in the analysis environment: Hanging with xesmf regridding in mpi on xp65 ACCESS-Analysis-Conda#207
• With 72 cpus, the full workflow completes in about 57 minutes of wall time. If the intermediate file already exists, the regridding step alone takes only about 1 minute.
• Once the intermediate bottom-roughness file is available, regridding to different model resolutions is very fast.

[minghangli-uni]

Overall the new \sqrt(h^2) (bottom left) show strong qualitative agreement in the large-scale features seen in both @CallumJShakespeare’s 100 km config (top left) and the OM2 configuration (top right). The peak magnitudes in the new field (bottom-left 25km config) are comparable to those in Callum’s matlab result, but are roughly twice those in OM2, which is capped at 1000m . The bottom-right figure (25km config) shows the roughness field currently used in OM3.

@aekiss @CallumJShakespeare would you mind taking a look and letting me know if you have any suggestions?

[minghangli-uni]

For completeness, the right panel shows the bottom roughness onto the OM3 100km grid, while the left shows Callum’s dataset.

Note:

• Callum's \sqrt(h^2) is based on WOA18, while the OM3 uses WOA23. The matlab bathymetric sampling uses the built-in interp2, whereas this implementation applies a manual interpolation approach to minimise overhead.

[minghangli-uni]

@aekiss thanks for raising this #105 (comment).

I've changed the workflow so that filling is no longer done purely on the woa ocean mask. Instead now the fill mask is obtained from the actual esmf regridding weights. After building the regridder, I've extracted which woa source cells contribute to at least one MOM6 wet target cell and then expand the fill mask to include both the original woa wet cells and these "needed" source cells. Then apply the laplacian fill over this combined mask before regridding. In this case any mom6 wet cells that draw from woa cells outside the original woa wet mask still receive valid values. So there''' be no nans caused by the mask mismatch and only extrapolating where the weights prove it is required.

Here's a notebook https://github.com/minghangli-uni/sandbox/blob/master/external_tidal/Update_bottom_roughness.ipynb showing the key steps and diagnostics. It also includes a comparison across different mapping methods - for the 25 km config, the resulting fields look qualitatively very similar though.

[aekiss]

OK thanks. That seems a very complicated way to go about it compared to the IC approach #105 (comment), but if it works I'm happy.
Revised patch seems to have problems, but we could use bilinear instead.

[minghangli-uni]

Revised patch seems to have problems, but we could use bilinear instead.

Correct. Here is what I see from the mask check.

# da1: conservative_normed regridding method (initial try)
# da2: conservative_normed regridding method (revised)
# da3: bilinear regridding method (revised)
# da4: patch regridding method (revised)

for i, da in enumerate([da1, da2, da3, da4], start=1):
    same_mask = np.array_equal(np.isfinite(da.values), mask.values)
    print(f"da{i} uses same mask:", same_mask)
da1 uses same mask: False
da2 uses same mask: True
da3 uses same mask: True
da4 uses same mask: False

So for lower resolution configs e.g. 100km configs, we will use the conservative_normed method. For moderate to high resolution configs, bilinear interpolation will be used.

If you're happy with these results, would you be able to have a final look and approve this @aekiss ?

[aekiss]

There are 2 unresolved comments - let me know when they're done or rejected

[aekiss]

Great work, thanks @minghangli-uni !

[minghangli-uni]

Really appreciate the careful review @aekiss Thanks!