sea_level roughness depends on timestep

[aekiss]

I thought I should make a note of this issue that Josué discovered and investigated earlier this year.

The top panel shows the sum of the squared sea_level grid cell differences in i and j, calculated daily over 1958-2019 for the three 0.1 deg IAF cycles. The 3 other panels show the timestep. As you can see, the SSH roughness goes up whenever the timestep is reduced (open image in a new tab and zoom in for a better view).

This timestep-dependence is rather disconcerting, as it affects geostrophic kinetic energy calculations.

The fact that the roughness increases with decreasing timestep surprises me - if this was due to numerical instability in the barotropic solver I would have expected it to go the other way around.

I have checked that the correct files are being loaded, so this isn't related to https://github.com/COSIMA/cosima-cookbook/issues/252.

[aekiss]

Some things that could be good to know:

• [ ] does this signal appear at lower resolution (would be much easier to investigate)?
• [x] does it appear in monthly averages? (if not, signal in sea surface gradient may be oscillatory on timescales between 1 day and 1 month)
• [ ] how does it show up in histograms? (ie is signal widespread, or a few extreme values?)
• [ ] what geographical regions dominate the signal?
• [ ] does it show up in wavenumber spectra? or frequency spectra?
• [ ] does it appear in any other fields?
• [ ] does this signal go away if barotropic_split is adjusted to keep the barotropic timestep the same?

[aekiss]

@StephenGriffies have you seen anything like this before?

[StephenGriffies]

I must admit not knowing just what I am seeing.

What is the math expression that we expect to be diagnosing here?

What is the timestep that is being changed?

[aekiss]

The timeseries is

(sea_level_{i+1} - sea_level_{i})^2 + (sea_level_{j+1} - sea_level_{j})^2

summed over all grid cells. It's not intended to be physically meaningful, but is just a simple diagnostic with similar properties to the surface geostrophic KE that Josué was looking at.

By "timestep" I mean baroclinic timestep, which equals the tracer timestep (baroclinic_split=1). barotropic_split is fixed at 80.

The barotropic dynamics use a predictor-corrector method, with the recommended method (barotropic_time_stepping_a=true) and default dissipation parameter pred_corr_gamma=0.2 and smoothing smooth_eta_t_laplacian=true, vel_micom_lap=0.05m/s.

[StephenGriffies]

This is odd. I do not have prior experience. Will need to think on this.

[aekiss]

Yes, it's puzzling.

Note that the timestep is reduced to deal with numerical instability, so it is actually correlated with model state. However, with the normal timestep the model would typically go unstable (usually by a CFL error in CICE) only on a particular day or two within a 3-month run (e.g. due to a storm near a tripole), whereas the timeseries show enhanced slopes throughout the full 3 months in which the timestep was reduced. So I don't think the signal could be due to the particular instability that motivated the reduction in timestep.

[aekiss]

@josuemtzmo do you have anything to add?

[josuemtzmo]

I haven't looked into this issue since I raised it. As far as I remember, the model surface velocities didn't show steps (I'm checking it atm, I will update you once I get more info), it only seems to affect the sea_level and potentially eta_t, which I believed is used to compute energy budgets. I looked in regions and hemispheres, and the step seems to occur everywhere spatially:

Difference between 29-06-2010 and 02-07-2010 Eddy Kinetic Energy computed from sea_level.

Time series of KE for 2010.

[aekiss]

thanks @josuemtzmo - our old slack discussions are still available here https://arccss.slack.com/archives/C6PP0GU9Y/p1614235878023200 and here (PM) https://climatefluidphysics.slack.com/archives/D0119RUBT5G/p1614729405001300

[StephenGriffies]

Is the output saved every time step or daily time averages. If time averages, then perhaps something is amiss with the diag manager time averaging functionality.

[aidanheerdegen]

@aekiss is there a jupyter notebook available with the code for reproducing this?

[aekiss]

Yes, just uploaded here: https://github.com/aekiss/notebooks/blob/master/ssh_gradient_timeseries_steps/ssh_gradient_timeseries_steps.ipynb the daily timeseries files take a long time to calculate, so it would be best to grab mine from /g/data/v45/aek156/notebooks/github/aekiss/notebooks/ssh_gradient_timeseries_steps/01deg_jra55v140_iaf*_sea_level_diffsq_sum_????-01-01_????-01-01.nc

[aekiss]

Thanks @StephenGriffies, they are daily time averages. We do have have some monthly sea_level snapshots in cycle 1, so I'll see if the steps show up there as well.

[aekiss]

Unfortunately monthly sea_level snapshots weren't saved over the dates in which the timestep changed. I will look around for other experiments instead.

[aekiss]

Here's what the timeseries look like when calculated from monthly mean data. The roughness is still dependent on timestep, suggesting it isn't due to intra-monthly variability.

[aekiss]

Here are timeseries of the horizontal sum of monthly mean u and v squared at 30m depth (chosen to be below the Ekman layer, so the velocity is presumably close to the surface geostrophic current due to sea_level gradients).

There is no evidence of the timestep dependence we would expect if the sea_level roughness timestep dependence was genuine. So it seems more likely that this is some sort of bug in the diagnostic manager, as @StephenGriffies suggested.

[aekiss]

Unfortunately I don't see this timestep-dependence in SSH roughness in a 1 degree experiment, so this may require higher resolution to investigate properly.

FYI this run is in /home/156/aek156/payu/1deg_jra55_iaf_iss243

[aekiss]

@russfiedler suggested this might be a symptom of the checkerboard null mode that affects the barotropic equations on a B-grid. We might want to see if smoothing parameters have any effect (see above for the parameters we currently use at all resolutions).

[russfiedler]

I don't fully understand your python code but am I right in saying that you are computing forward (backward?) differences on the T points rather than centered differences on the U points?

[aekiss]

They are simple 1-sided backward differences on the t-points (using diff), so would highlight the checkerboard mode

            sea_level_xdiff = sea_level.diff('xt_ocean')
            sea_level_ydiff = sea_level.diff('yt_ocean')
            sea_level_diffsq = sea_level_xdiff.isel(yt_ocean=slice(0, -1))**2 \
                             + sea_level_ydiff.isel(xt_ocean=slice(0, -1))**2
            sea_level_diffsq_sum = sea_level_diffsq.sum('xt_ocean').sum('yt_ocean')

[aekiss]

The barotropic mass convergence conv_rho_ud_t shows checkerboarding in some marginal seas - e.g. North Sea, Baltic, Mediterranean, Red Sea, Persian Gulf, Black Sea (click to zoom in). The x transects are from the Baltic. These plots are from restart375, ending in 1991-01-01 in cycle 2 with dt= 360s /g/data/ik11/restarts/access-om2-01/01deg_jra55v140_iaf_cycle2/restart375/ocean/ocean_barotropic.res.nc

Checkerboarding is also present (and visually very similar) in /g/data/ik11/restarts/access-om2-01/01deg_jra55v140_iaf_cycle2/restart379/ocean/ocean_barotropic.res.nc (ending in 1992-01-01) which used dt=540s.

So it's not clear to me that checkerboarding could explain the dependence on timestep. Also the fact that the roughness increases with decreasing timestep surprises me - if this was due to numerical instability in the barotropic solver I would have expected it to go the other way around.

[aekiss]

Also worth noting that some 0.25° configs also show checkerboarding: https://github.com/COSIMA/access-om2/issues/274#issuecomment-1544942926