C-grid crash with velocity instabilities

anton-seaice commented 2 weeks ago

This may well be related to #941 but Ill write it up seperately for clarity.

I ran CICE6 c-grid, coupled to MOM6 and using data atmosphere forcing. Its a nominal 0.25 degree resolution, we've been running with a CICE6 / MOM6 Baroclinic / Coupling timestep of 1350s (possibly too-long) and a short MOM6 barotropic timestep (675s).

After ~15 months, CICE crashes with negative ice area. This is the crash location, on the Greenland coast.

The last instantaneous output shows clear patterning in sea ice concentration at the crash location:

And instabilities in the C grid velocities. Note the main instability shown is in uvelE and vvelN

It possible ofcourse this is not an issue with C-grid advection - our ocean model configuration is still undergoing changes. However we are confident that this configuration wouldn't crash with CICE B-grid. I am happy to do re-runs, provide extra plots etc to investigate further.

Crash log:

New mass < 0, i, j =          50          31
 Old mass =  1.584993725340139E-004
 New mass = -5.727752649206490E-006
 Net transport = -1.642271251832203E-004
 istep1, my_task, iblk, cat =       28799          80           4           2
 (print_state) negative area (ice)�
 (print_state) istep1, my_task, i, j, iblk:       28799          80          50
          31           4
 (print_state) Global block:         450
 (print_state) Global i and j:        1069        1002
 (print_state) Lat, Lon (degrees):   81.6865092532964     
  -14.8643879243489     

 aice    1.00000000000000     
 aice0  1.110223024625157E-016

 n =           1
 aicen  0.000000000000000E+000
 vicen  0.000000000000000E+000
 vsnon  0.000000000000000E+000
 Tsfcn  -1.79028623325174     

 n =           2
 aicen  1.584993725340139E-004
 vicen  1.541077117507753E-004
 vsnon  2.501025765687183E-009
 hin  0.972292251300263     
 hsn  1.577940483739432E-005
 Tsfcn  -20.1079585645912     

 n =           3
 aicen  2.710599464180833E-004
 vicen  5.327318905188697E-004
 vsnon  4.277164629733349E-009
 hin   1.96536558631641     
 hsn  1.577940483739432E-005
 Tsfcn  -22.7879628146431     

 n =           4
 aicen  1.058556202804489E-002
 vicen  3.713120284791096E-002
 vsnon  4.344231959587402E-003
 hin   3.50772143694754     
 hsn  0.410392187781622     
 Tsfcn  -26.5975737336460     

 n =           5
 aicen  0.988984878653003     
 vicen   8.07225859281393     
 vsnon  0.575890443169084     
 hin   8.16216583999580     
 hsn  0.582304598987849     
 Tsfcn  -26.7320323153733     

 qice, cat            1  layer            1  0.000000000000000E+000
 qice, cat            1  layer            2  0.000000000000000E+000
 qice, cat            1  layer            3  0.000000000000000E+000
 qice, cat            1  layer            4  0.000000000000000E+000

 qice, cat            2  layer            1  -307755649.562915     
 qi/rhoi  -335611.395379406     
 qice, cat            2  layer            2  -305955239.880947     
 qi/rhoi  -333648.026042472     
 qice, cat            2  layer            3  -305976569.446496     
 qi/rhoi  -333671.286201195     
 qice, cat            2  layer            4  -305655112.966541     
 qi/rhoi  -333320.733878453     

 qice, cat            3  layer            1  -306509750.424470     
 qi/rhoi  -334252.726744242     
 qice, cat            3  layer            2  -305988097.504735     
 qi/rhoi  -333683.857693277     
 qice, cat            3  layer            3  -305997732.789325     
 qi/rhoi  -333694.365091957     
 qice, cat            3  layer            4  -305915232.195830     
 qi/rhoi  -333604.397160120     

 qice, cat            4  layer            1  -330320531.213201     
 qi/rhoi  -360218.681802837     
 qice, cat            4  layer            2  -321686608.131194     
 qi/rhoi  -350803.280404791     
 qice, cat            4  layer            3  -312374812.336550     
 qi/rhoi  -340648.650312486     
 qice, cat            4  layer            4  -298970272.633928     
 qi/rhoi  -326030.831661863     

 qice, cat            5  layer            1  -335702127.033382     
 qi/rhoi  -366087.379534768     
 qice, cat            5  layer            2  -324703839.344245     
 qi/rhoi  -354093.608881401     
 qice, cat            5  layer            3  -313196657.706246     
 qi/rhoi  -341544.882994815     
 qice, cat            5  layer            4  -297369967.905740     
 qi/rhoi  -324285.679286521     

 qice(i,j)  -4984078201.07574     

 qsnow, cat            2  layer            1  -110121000.000000     
 qs/rhos  -333700.000000000     
 Tsnow  0.000000000000000E+000

 qsnow, cat            3  layer            1  -110121000.000000     
 qs/rhos  -333700.000000000     
 Tsnow  0.000000000000000E+000

 qsnow, cat            4  layer            1  -125310280.183219     
 qs/rhos  -379728.121767329     
 Tsnow  -21.7393727617777     

 qsnow, cat            5  layer            1  -126251426.824129     
 qs/rhos  -382580.081285240     
 Tsnow  -23.0863712635802     

 qsnow(i,j)  -471803707.007348     

 sice, cat            1  layer            1  0.000000000000000E+000
 sice, cat            1  layer            2  0.000000000000000E+000
 sice, cat            1  layer            3  0.000000000000000E+000
 sice, cat            1  layer            4  0.000000000000000E+000
 sice, cat            2  layer            1  0.649201856400313     
 sice, cat            2  layer            2   2.35458111192744     
 sice, cat            2  layer            3   3.03109206979846     
 sice, cat            2  layer            4   3.18929777221700     
 sice, cat            3  layer            1  0.649201856400314     
 sice, cat            3  layer            2   2.35458111192744     
 sice, cat            3  layer            3   3.03109206979846     
 sice, cat            3  layer            4   3.18929777221700     
 sice, cat            4  layer            1  0.649201856400311     
 sice, cat            4  layer            2   2.35458111192744     
 sice, cat            4  layer            3   3.03109206979848     
 sice, cat            4  layer            4   3.18929777221700     
 sice, cat            5  layer            1  0.649201856400316     
 sice, cat            5  layer            2   2.35458111192745     
 sice, cat            5  layer            3   3.03109206979847     
 sice, cat            5  layer            4   3.18929777221700     

 uvel(i,j)  0.421109861482645     
 vvel(i,j)  0.626738511376921     
 uvelE(i,j)   2.01817027718648     
 uvelN(i,j)  0.109513019152795     

 atm states and fluxes
             uatm    =    11.4565684617028     
             vatm    =   -5.35149824880375     
             potT    =    247.552100916471     
             Tair    =    247.552100916471     
             Qa      =   2.939509003865902E-004
             rhoa    =    1.39554159484570     
             swvdr   =   0.000000000000000E+000
             swvdf   =   0.000000000000000E+000
             swidr   =   0.000000000000000E+000
             swidf   =   0.000000000000000E+000
             flw     =    173.710396449191     
             frain   =   6.817549163993367E-009
             fsnow   =   3.857187849140834E-006

 ocn states and fluxes
             frzmlt  =    1130.89758885969     
             sst     =   -1.78566706771846     
             sss     =    32.9096734053628     
             Tf      =   -1.79028623325174     
             uocn    =   7.000659626299037E-002
             vocn    =   4.925604188204023E-002
             strtltxU=  -7.371310560342537E-003
             strtltyU=   3.605963438896476E-002

 srf states and fluxes
             Tref    =    247.355015644017     
             Qref    =   3.016600190419180E-004
             Uref    =    12.5024800827170     
             fsens   =    33.0991790281445     
             flat    =   -3.38962097696834     
             evap    =  -1.175569424330643E-006
             flwout  =   -207.304897334251     

 (abort_ice)ABORTED: 
 (abort_ice) error = (diagnostic_abort)ERROR: negative area (ice)
 ERROR: (abort_ice)(diagnostic_abort)ERROR: negative area (ice)

(Built cice using https://github.com/CICE-Consortium/CICE/commit/12dd204349090058a66715163932ae3243f9632c , I can't see any commits since then they would impact this).

phil-blain commented 2 weeks ago

Hi Anton!

@JFLemieux73 is not at the office this week, but I think his first question would be: does it also crash with upwind advection ?

anton-seaice commented 2 weeks ago

Thanks Phil - I will try that, and I will also make some figures of spatial plots of uvelE/vvelN similar to https://github.com/CICE-Consortium/CICE/issues/976#issuecomment-2383594287

anton-seaice commented 2 weeks ago

Both re-running with upwind advection or re-running with a 900s timestep didn't crash the model. Neither are a workable solution for us, although I could experiment with timesteps some more if there is nothing else strange going on.

JFLemieux73 commented 1 week ago

Hi Anton, I have the impression this is an instability related to the coupling. Can you restart just before the oscillations develop and show plots of uvel,vvel as above with a smaller time step? I am adding @dupontf to the conversation. Fred any idea?

dupontf commented 1 week ago

Hi Anton, what is your ice-ocean drag? and when you said above

re-running with a 900s timestep didn't crash

I suppose you mean that both the coupling (CICE) and MOM time-steps were 900s?

I am trying to eliminate a coupling instability: which time-level is passed to CICE from MOM?

anton-seaice commented 1 week ago

Thanks both

Here is the plot from the 900s timestep, ( the coupling , CICE and MOM all run at this 900s timestep ).

You can see the instability but it resolves - it doesn't show up in the sea-ice concentration. Maybe this is normal ?

dupontf commented 1 week ago

Thanks Anton, this still could correspond to a coupling instability to me as 900s corresponds to a reduction of the coupling timestep, so again, which value of drag and which time-level are exchanged between CICE and MOM?

I have found in the past for instance in NEMO and LIM2 that having the ocean current at time-level n given to CICE for calculating ice velocity at time-level n+1 and the associated stress given back to NEMO for computing current at time-level n+1 (and with Leapfrog time-stepping) yields an unstable situation akin to have $$u^{n+1}=u^{n-1}+Cd (u{ice}-u^n)$$

dupontf commented 1 week ago

During the CICE virtual meeting, it was determined that MOM is using a RK2 time-steppîng, which should not be as sensitive as Leapfrog to the time level exchanged. So less likely a coupling issue.
Elizabeth raised the possibility of a Coriolis mode. This would require more investigation.
D. Bailey clarified that the momentum exchange is done on T-point through the coupler.
A. Roberts suggested to plot a phase lag plot between ice and ocean.

JFLemieux73 commented 4 days ago

Hi @anton-seaice,

It would be good to have an idea if this comes from the remapping. Could you please restart just before the oscillations develop and use upwind instead? Use your usual time steps: coupling timestep of 1350s and MOM6 barotropic timestep of 675s. Do we still see the oscillations in the velocity (please show plots as above)?

JFLemieux73 commented 4 days ago

About Coriolis...I also coded a semi-implicit approach:

https://github.com/JFLemieux73/CICE/tree/semi_imp_Coriolis

This is not up to date with the current code. If you want to try it you might have to cherry pick some parts in ice_dyn_evp and ice_dyn_shared.

eclare108213 commented 4 days ago

A few thoughts from me:

What we know (anything else?):

Configuration ~0.25 deg resolution 1350 s CICE6 / MOM6 Baroclinic / Coupling timestep 675 s MOM6 barotropic timestep
C-grid fails, B-grid does not Oscillations in space and time lead to negative area and mass ~15 months into the simulation Location on the northern Greenland coast
Using upwind instead of remapping does not fail, but it’s very diffusive so that does not tell us much.
Using a 900 s time step for CICE, MOM and their coupling does not fail but still shows some instability.

Some suggestions:

To debug, it’s helpful to have a restart file immediately before the oscillations develop.
In the failing models, is the momentum exchange done at cell centers or at the C-grid points? If it’s at the C-grid points, does it still fail when coupling at cell centers?
Could this be coming from residual elastic waves? Crank up EVP subcycling to ~1200 or more.
To further understand time step dependence of coupling versus component model time steps, run CICE and MOM at 900 s (or smaller, to eliminate potential instabilities coming from them) and couple less often, e.g. 1350 s. Andrew's coupling diagnostics could also help here.
Try JF’s semi-implicit Coriolis approach.

eclare108213 commented 4 days ago

Considering the location, what happens if you turn on landfast ice? Maybe the increased tensile stress would help?

If none of those things provide any insight, start simplifying your configuration as much as possible while maintaining the failure mode. E.g. Turn off the thermodynamics. Turn off the ice-ocean coupling, or individual terms associated with the dynamics. Turn off advection. Try to set up a similar configuration in a box model (grid size, thickness, velocity magnitude and angle with the coast, etc)....

anton-seaice commented 3 days ago

Thankyou all for the suggestions!

Re: upwind advection

It would be good to have an idea if this comes from the remapping. Could you please restart just before the oscillations develop and use upwind instead? Use your usual time steps: coupling timestep of 1350s and MOM6 barotropic timestep of 675s. Do we still see the oscillations in the velocity (please show plots as above)?

These are the plots from the experiment with upwind advection:

This appears to rule out the advection scheme somewhat as the instability is still there with upwind advection.

JFLemieux73 commented 3 days ago

Good! This is very useful.

CICE-Consortium / CICE

C-grid crash with velocity instabilities #992