Add tidal related external files

[minghangli-uni]

This PR complements #52. Additional related information is available in https://github.com/COSIMA/access-om3/issues/280

• Bottom-Xm-averaged stratification from WOA climatology
• Bottom roughness
• Tidal amplitude

[minghangli-uni]

@luweiyang, could you please review the script for generating the bottom-averaged frequency first? I will push the other scripts later.

[luweiyang]

The script looks great! For the barotropic (global 1-layer configuration) model I ran, I used depth-averaged N instead of bottom-500m-averaged N. I'd suggest using depth-averaged N if the model does not include tidal forcing, but this script will be really useful for generating N for a multilayer model with tidal forcing. I computed N from WOA18 data because it wasn’t available from the barotropic model. However, if the model has sufficient vertical resolution, N can also be computed from the model output. Thanks for putting this together!

[minghangli-uni]

Thanks @luweiyang Just to clarify, when you say “depth-averaged N,” do you mean averaging N over the entire water column, from surface to bottom? I’ve updated the code to allow both a bottom-Xm average and a full-column average, and the results look quite different. Could you confirm if this difference is expected?

[luweiyang]

Yes, @minghangli-uni the depth-averaged N is obtained by averaging from surface to seafloor. Both plots look reasonable, and the differences are expected — in the open ocean (deeper than 1000 m), the depth-averaged N is about an order of magnitude larger than the bottom-500m-averaged N.

The depth-averaged N (the larger N) is needed to estimate the energy in low-mode internal tides (large horizontal and vertical scales), whereas the bottom-500m-averaged N is more relevant for estimating the energy of high-mode internal tides (small horizontal and vertical scales). Low-mode internal tides are much more energetic than high-mode internal tides.

In a model with tidal forcing, low-mode internal tides will be resolved, and you can use bottom-averaged N to parameterize the dissipation of high-mode internal tides. However, if tidal forcing is absent, even if the resolution is high enough to resolve low-mode internal tides, they will not be present because the forcing was missing in the first place. So, this part of the internal tide energy will need to be included in the parameterization, and the depth-averaged N should be used.

Does it make sense?

[minghangli-uni]

Thanks, @luweiyang ! That makes perfect sense to me! Apologies for the delayed response—I was attending the Bluelink meeting over the past two days. I’ll update the other two scripts in this PR within a day or two.

[minghangli-uni]

The current approach to estimating bottom roughness is based on the polynomial fitting method proposed by Jayne, Steven R., and Louis C. St. Laurent. "Parameterizing tidal dissipation over rough topography." Geophysical Research Letters 28.5 (2001): 811-814. In this method,

Over eachgrid cell, a polynomial sloping surface is fit to the bottom topography (givenby H = a + bx+ cy+ dxy), and the residual heights are used to compute h2, the mean-square bottom roughness averaged over the grid cell.

An alternative approach using an FFT-based method, developed by @luweiyang, is also being considered. However, there are uncertainties in how this method handles coastal regions where resolution changes near the land-sea boundary. More sensitivity tests need to be done. So the use of Bottom-Xm-averaged stratification from WOA climatology may not be essential for now, but we will retain it for potential future improvements.

[minghangli-uni]

As discussed with @luweiyang and @aekiss, I tested three approaches for handling land cells when applying the polynomial fit to compute bottom roughness, using the high-resolution SYNBATH dataset (1/240 deg resolution):
1. All Depths: Retaining both positive (land) and negative (ocean) elevations as-is. While it preserves all data, it might risk introducing biases from land elevations (which are not meaningful for ocean bottom roughness) into the fit?
2. Land as 0 m: Setting all positive elevations to 0 m before fitting. But it may introduce sharp transitions between coastal and ocean regions. For example, the coastal areas could be overly smoothed because the fit enforces a hard zero rather than reflecting natural gradients?
3. Land as NaN: Masking positive elevations by setting them to NaN to exclude them from the fit. My understanding is this method may yield a relatively more accurate representation of the ocean floor because the polynomial fit is calculated solely using valid, ocean-relevant data, without the “affection” or influence of land cells?

---

Below shows the comparisons of the resulting HRMS (NB: the output is HRMS square) fields for each approach. For reference, I also include the bottom roughness from OM2, which lacks provenance. But I suspect it may have been derived from a 1/30 deg bathymetry dataset.

1. Four plots use a consistent color scale ranging from blue (0 m) to red (1000 m). The OM2 roughness field is capped at exactly 1000 m (this seems a bit suspicious to me and suggests that some postprocessing might have been applied), whereas the roughness from our polynomial fit reaches a maximum of 1382.10986328m.
   

2. The color ranges from blue (0m) to red (200m) for some more details.
   

3. Direct co
   
   mparisons between the three types.

[luweiyang]

Thanks @minghangli-uni ! So using both positive and negative elevations helps recover some hrms near the coast, but not a lot (~60m compared to 200m). Where do you see the extreme hrms values (>=1000m) in the OM2 file and your calculation? Do they appear in the same locations?

It would be nice to attach the figures showing the differences due to different source data (Smith and Sandwell 1997, GEBCO, and SYNBATH), and to better understand the differences caused by the sampling method - grid resolution (calculating on every grid cell vs. on a 1-deg grid and interpolating) and sampling window size (1/4-deg, 1-deg, etc), as we discussed.

[minghangli-uni]

Where do you see the extreme hrms values (>=1000m) in the OM2 file and your calculation? Do they appear in the same locations?

Thanks @luweiyang . Below the red crosses show the max HRMS values for each figure, while the yellow crosses mark the locations with values greater than 1000 m. This appears to confirm that the HRMS in OM2 did some forced postprocessing to cap the maximum value exactly at 1000 m. Not sure for the reason behind...

[minghangli-uni]

It would be nice to attach the figures showing the differences due to different source data (Smith and Sandwell 1997, GEBCO, and SYNBATH),

@luweiyang , below shows the results using GEBCO 2024 and SYNBATH and also their differences for each method. As expected, the differences arise from the inherent characteristics of the bathymetry sources. Similar to the above, red crosses show the max HRMS values for each figure, while the yellow crosses mark the locations with values greater than 1000 m.

[minghangli-uni]

and to better understand the differences caused by the sampling method - grid resolution (calculating on every grid cell vs. on a 1-deg grid and interpolating) and sampling window size (1/4-deg, 1-deg, etc), as we discussed.

@luweiyang The comparison above shows results from OM2, calculations on every grid cell (0.25deg), and calculations on 0.5deg and 1deg grids. The sampling window corresponds to each grid resolution respectively. Changing the resolution does have an impact, particularly near the coastlines. The OM2 results appear to be more similar to those from the 0.5deg resolution case?

Similar to previous plots, red crosses show the max HRMS values for each figure, while the yellow crosses mark the locations with values greater than 1000 m.

[minghangli-uni]

This commit cbfbcd3 estimates the deg-> meters conversion based on a geometric approximation. It assumes a spherical earth and uses a locally cartesian projection. Not sure about its accuracy near the poles or is there a more accurate method for this conversion? Appreciate any feedback on this.

    R = 6.371229e6  # earth radius in meters
    deg2rad = np.pi / 180
    lat0 = np.mean(y)  # a reference lat
    lon0 = np.mean(x)  # a reference lon

    xlon = (xv - lon0) * np.cos(yv * deg2rad) * deg2rad * R
    ylat = (yv - lat0) * deg2rad * R

[luweiyang]

Thanks @minghangli-uni for the updates.

When I used the plane fit method (as in JSL), I used longitude as x and latitude as y so both the original topography and fitted topography are functions of longitude and latitude. I then calculated the hrms using the original and fitted topography. My understanding is that your xv and yv represent the distance of a grid point from a reference grid point (the centre of the sampling domain). What is the reason for doing it this way? Also, do you plan to use the spectral method to calculate hrms at some point?

[minghangli-uni]

My understanding is that we are fitting a polynomial to bottom roughness (in meter), but our horizontal coordinates are still in degrees of lon/lat, so this means the fitted roughness comes out in meters per degree. To get a true roughness in meters per meter of horizontal distance, I think we should first convert our lon/lat offsets into linear distances (meters)?

And the other consideration is 1deg of latitude always spans roughly 111 km, but 1deg of lon spans only 111 km × cos(lat). That distortion grows as you move toward the poles, so using raw degrees would introduce lat-dependent errors in the roughness estimates. Correct me if I am wrong!

[minghangli-uni]

Below shows the comparisons of tidal speed amplitude between the current generation (left) and the one from OM2 (right).

The maximum values are 1.87796m/s and 0.99725m/s respectively. The unit displayed in the colorbar is incorrect. It should be m/s, not meters.

---

It takes around 36mins to generate tideamp using generate_tide_amplitude.py, with the majority of the time spent on the regridding step.

[minghangli-uni]

TPXO10 files and SYNBATH are now in ik11.

[luweiyang]

This commit cbfbcd3 estimates the deg-> meters conversion based on a geometric approximation. It assumes a spherical earth and uses a locally cartesian projection. Not sure about its accuracy near the poles or is there a more accurate method for this conversion? Appreciate any feedback on this.

    R = 6.371229e6  # earth radius in meters
    deg2rad = np.pi / 180
    lat0 = np.mean(y)  # a reference lat
    lon0 = np.mean(x)  # a reference lon

    xlon = (xv - lon0) * np.cos(yv * deg2rad) * deg2rad * R
    ylat = (yv - lat0) * deg2rad * R

It makes sense to me! xlon gets very small near the poles, but that’s expected. Maybe it would be nice to add a bit more explanation for future reference? xlon represents ∆x(=R*cos(lat)∆(lon)), and ylat represents ∆y(=R∆(lat)) away from the reference point (lon0, lat0), which is the grid point the hrms is calculated for and also the centre of the sampling box.

[luweiyang]

@minghangli-uni Could you add the reference here so others understand why you calculated the tidal velocity amplitude for OM3 this way?