Sea ice checkerboarding

[aekiss]

As mentioned in https://github.com/COSIMA/mom6-panan/issues/34#issuecomment-1487773967, the sea ice concentration displays checkerboarding, indicative of numerical instability. This occurs for DT = 300s in the 1/20th.

Are there SIS2 parameters we can tweak to fix this? It would be good to avoid reducing the timestep if we can.

[aekiss]

@adele157 what does SIC checkerboard look like with DT = 225?

[adele-morrison]

Don't know, it crashed. Is there some way to access the daily output from work of a crashed run?

[aekiss]

Not unless there are .nc files you can find in there. Might be enlightening to shorten the run, so it stops just before the crash, and save output from that final state

[adele-morrison]

I checked that there's no checkboarding in panan-01.

[adele-morrison]

It also occurs in the southern Weddell, but nowhere else. Maybe related to where the grid cells are small?

[adele-morrison]

Looking through the panan sea ice parameters, I notice we have: ICE_BULK_SALINITY = 0.0 ! [g/kg] default = 4.0 ! The fixed bulk salinity of sea ice.

Probably we want to change this to 4 or 5, so it matches our other ACCESS runs? Was there any reason for this @AndyHoggANU ?

Oh, actually I see then there is this: ICE_RELATIVE_SALINITY = 0.1 ! [nondim] default = 0.0 ! The initial salinity of sea ice as a fraction of the salinity of the seawater ! from which it formed.

which will do nearly the same thing. Though somewhat more complicated, because the sea ice salinity will vary in time and space.

[adele-morrison]

Ah, and this one that I didn't change going from 1/10 to 1/20: DT_ICE_DYNAMICS = 180.0 ! [seconds] default = -1.0 ! The time step used for the slow ice dynamics, including stepping the ! continuity equation and interactions between the ice mass field and ! velocities.

I had assumed it was the same as the ocean time step! Probably we should halve this going to 1/20th, yes? I'm a bit confused by the large number of different timesteps. I'm thinking they should all be multiples of one another, is that right @angus-g? Otherwise the model adjusts them so they are? There's only a quite small number of possibilities for different timesteps if they all need to be multiples of each other.

[AndyHoggANU]

Err. Yes. Good idea to check on the input files! I agree, it should be an integer fraction of the ocean timestep. It probably adjusts to be something like that if it's not, but we clearly need to document how many time steps to change!

[angus-g]

I'm thinking they should all be multiples of one another, is that right @angus-g? Otherwise the model adjusts them so they are?

I don't actually know what happens if they're not, but if it doesn't give a flat-out error it seems like it could manifest as some confusing behaviour. I think this is the difference between the "fast" and "slow" ice timesteps. Maybe it's the "fast" one that's locked to the ocean (dynamics) timestep? I might look for (or create) a diagram relating all the different timesteps we've got going on!

[adele-morrison]

There's also: DT_RHEOLOGY = 50.0 ! [seconds] default = -1.0 ! The sub-cycling time step for iterating the rheology and ice momentum ! equations.

and

NSTEPS_ADV = 1 ! default = 1 ! The number of advective iterations for each slow dynamics time step.

and

[adele-morrison]

I can confirm the checkerboarding does not go away after 19 days of DT = 225. @aekiss how long would you expect we need to run for with changed parameters to see a difference?

I also notice we are getting very unrealistic sea ice thickness at the Antarctic Peninsula (this is 19 June):

[adele-morrison]

Which sea ice time step should I try playing with?

1. DT_ICE_DYNAMICS = 180.0 "for the slow ice dynamics, including stepping the continuity equation and interactions between the ice mass field and velocities". This is equal to 1200 or 3600 in most of the GFDL (coarse-res) examples.
2. or DT_RHEOLOGY = 50.0 "The sub-cycling time step for iterating the rheology and ice momentum equations." This is equal to 50 or 100 in most GFDL examples.

Also, does anyone know of any high-res MOM6-SIS2 example configs I could compare against?

[adele-morrison]

Reducing the timestep gives a big improvement (good one on top, bad one below)! I changed both DT_ICE_DYNAMICS (halved from 180 to 90) and DT_RHEOLOGY (from 50 - which was not a factor of 180 before - to 30). Need to figure out which one solved the problem now.

How did we pick DT_ICE_DYNAMICS = 180.0 @AndyHoggANU ? This seems way lower than in any other examples I can find.

[AndyHoggANU]

Not clear to me how we picked DT_ICE_DYNAMICS but I would say it is one of those trial and error parameters. If this is the parameter that fixes the checkerboarding then we use that as our metric to set it!

[aekiss]

For comparison, CICE5 has 3 timesteps: thermodynamic, dynamic and elastic.

In ACCESS-OM2-01 the sea ice thermodynamic timestep matches the ocean baroclinic timestep, and the CICE5 dynamic timestep is shorter by a factor of 2, though we've also experimented with a factor of 3 when hitting CFL limits.

CICE5 uses an elasto-visco-plastic rheology, with an artificial elastic part which creates artificial elastic waves that are supposed to die away during 120 sub-timesteps within each dynamic timestep.

I guess SIS2 also uses EVP, and that may be what DT_RHEOLOGY is referring to? If so, your DT_ICE_DYNAMICS/DT_RHEOLOGY = 3 is a lot smaller than the ratio of 120 we use in ACCESS-OM2, and even smaller than the default of 240 in CICE6. This may be what is creating the waves in SIC visible in the first test you did (circled), with a ratio of 3.6:

[aekiss]

Can you look at your new SIC with a more detail-focused colormap to see if it has waves? Ideally using a snapshot or short time-average.

[aekiss]

What sort of time-averaging was used in the plots so far?

[adele-morrison]

They are daily outputs (I assume averages?).

On Fri, 31 Mar 2023 at 13:00, Andrew Kiss @.***> wrote:

What sort of time-averaging was used in the plots so far?

— Reply to this email directly, view it on GitHub https://github.com/COSIMA/mom6-panan/issues/35#issuecomment-1491182649, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACA44U7IXIGD3CAPLXTMZM3W6Y3CLANCNFSM6AAAAAAWLJJHJA . You are receiving this because you were mentioned.Message ID: @.***>

[aekiss]

Waves might be more conspicuous in SIS2 diagnostics that look at velocity gradients, eg divergence, strain rate etc.

SIS2 uses the EVP formulation of Bouillon et al 2009 and Bouillon et al 2013, which follows that in CICE.

Bouillon et al 2013 used a DT_ICE_DYNAMICS/DT_RHEOLOGY ratio of 300.

GFDL OM4p25 and OM4p5 use a rheology time step of 50 s and a dynamics time step of 1200s, ie ratio of 24 (Adcroft et al 2019).

So if we still have elastic waves maybe we should try (say) DT_RHEOLOGY=5?

Though oddly, Bouillon et al 2009 say checkerboarding can't occur on a C-grid...

[adele-morrison]

For the record, I'm now running the 1/20th with DT_ICE_DYNAMICS = 225.0, DT_RHEOLOGY = 9.0, based on scaling the GFDL values of these by the resolution factor. This is looking good (from very brief glances at concentration in ncview - this could definitely do with some more people looking in more detail and turning on more sea ice diagnostics!).

[schmidt-christina]

I was wondering what other variables to save to investigate this problem further when I run the 1/20th model with the new bathymetry in the Ross Sea.

But before, I just checked Adele's latest run /home/157/akm157/mom6/panan-005-zstar-ACCESSyr2/ and I still see the above pattern in the first year of her run and also emerging in the second year (not shown).

Do we need to reduce DT even further or save some more diagnostics first to see where it's coming from?

[adele-morrison]

I only changed the sea ice timesteps in the last run segment I did. I was waiting for the new bathymetry before running again.

On Thu, 6 Apr 2023 at 13:36, Christina Schmidt @.***> wrote:

I was wondering what other variables to save to investigate this problem further when I run the 1/20th model with the new bathymetry in the Ross Sea.

But before, I just checked Adele's latest run /home/157/akm157/mom6/panan-005-zstar-ACCESSyr2/ and I still see the above pattern in the first year of her run and also emerging in the second year (not shown).

Do we need to reduce DT even further or save some more diagnostics first to see where it's coming from?

[image: image] https://user-images.githubusercontent.com/38216156/230265629-c2c4f971-9291-414d-8d96-499c3d4d72c1.png

[image: image] https://user-images.githubusercontent.com/38216156/230265816-9a4d6987-1a55-49a0-a3dd-441b2309cc60.png

— Reply to this email directly, view it on GitHub https://github.com/COSIMA/mom6-panan/issues/35#issuecomment-1498445834, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACA44U4HINOYQGRXKT2FVFDW7Y24TANCNFSM6AAAAAAWLJJHJA . You are receiving this because you were mentioned.Message ID: @.***>

[schmidt-christina]

Ok, so this (last day of your run) is ok? There is still a bit of this pattern

[schmidt-christina]

I can then run the model again with this time step and new Ross Sea bathymetry. Any other sea ice diagnostics I should save?

[adele-morrison]

Yeah, looks like most of the checkerboarding is gone, but we still have waves. Is that a problem @aekiss? And how long would we need to run for with the smaller timesteps until these go away? That image above was only a 1 month run with shorter timesteps I think, starting from the bad restart with longer timesteps.

I reckon just start running @schmidt-christina. No doubt something else is going to break and we'll need to restart in any case. We can switch on more diagnostics as we go, but good to just get some more model years under our belt in the meantime, so we can learn a bit more and get the 1/40th going asap.

[adele-morrison]

@schmidt-christina can you confirm if the waves disappeared with the change in sea ice timesteps in your last 1/20th run? If so, let's close this.

[schmidt-christina]

Yes, there are no waves when running the 1/20th with DT_ICE_DYNAMICS = 225.0, DT_RHEOLOGY = 9.0

[aekiss]

I think we might need a bigger DT_ICE_DYNAMICS/DT_RHEOLOGY ratio than 225/9=25, i.e. a smaller DT_RHEOLOGY.

We could use the default DT_RHEOLOGY = -1.0 to avoid having to adjust DT_RHEOLOGY whenever we change DT_ICE_DYNAMICS. This will then use the ratio NSTEPS_DYN to determine the elastic timestep. NSTEPS_DYN defaults to 432, so will give a much smaller DT_RHEOLOGY of about 0.52s in our case with DT_ICE_DYNAMICS = 225.0.

https://github.com/NOAA-GFDL/SIS2/blob/c565ec390aba0346b26b8c0acedda99386008aad/src/SIS_dyn_cgrid.F90#L204-L218

  call get_param(param_file, mdl, "DT_RHEOLOGY", CS%dt_Rheo, &
                 "The sub-cycling time step for iterating the rheology "//&
                 "and ice momentum equations. If DT_RHEOLOGY is negative, "//&
                 "the time step is set via NSTEPS_DYN.", units="seconds", &
                 default=-1.0, scale=US%s_to_T)
  CS%evp_sub_steps = -1
  if (CS%dt_Rheo <= 0.0) &
    call get_param(param_file, mdl, "NSTEPS_DYN", CS%evp_sub_steps, &
                 "The number of iterations of the rheology and ice "//&
                 "momentum equations for each slow ice time step.", default=432)
  call get_param(param_file, mdl, "ICE_TDAMP_ELASTIC", CS%Tdamp, &
                 "The damping timescale associated with the elastic terms "//&
                 "in the sea-ice dynamics equations (if positive) or the "//&
                 "fraction of DT_ICE_DYNAMICS (if negative).", &
                 units="s or nondim", default=-0.2)

https://github.com/NOAA-GFDL/SIS2/blob/c565ec390aba0346b26b8c0acedda99386008aad/src/SIS_dyn_cgrid.F90#L776-L781

  if (CS%dt_Rheo > 0.0) then
    EVP_steps = max(CEILING(dt_slow/CS%dt_Rheo - 0.0001), 1)
  else
    EVP_steps = CS%evp_sub_steps
  endif
  dt = dt_slow/EVP_steps