C-grid crash with velocity instabilities #992
Comments
|
Hi Anton! @JFLemieux73 is not at the office this week, but I think his first question would be: does it also crash with |
|
Thanks Phil - I will try that, and I will also make some figures of spatial plots of uvelE/vvelN similar to #976 (comment) |
|
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. |
|
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? |
|
Hi Anton, what is your ice-ocean drag? and when you said above
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? |
|
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 ? |
|
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 |
|
|
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)? |
|
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. |
|
A few thoughts from me: What we know (anything else?):
Some suggestions:
|
|
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).... |
|
Thankyou all for the suggestions! Re: upwind advection
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. |
|
Good! This is very useful. |
|
I made a new test after updating to CICE 6.6.0. All exchange variables are on T grid. CICE on B-Grid runs a 384 hour forecast to completion. CICE on C-Grid crashes after 56 hours with error below. I noticed the print state is for the open water, not ice. The location is in Hudson Bay. Ice is melting, as the date is July 13, 2021. New mass < 0, i, j = 6 3032 aice 5.974142390434661E-004 A few notes:
I'd like to start making the ice phase plots mentioned above. Could I get help on what data to look at? Where/when I should be writing what ice information? Thanks for the help! |
|
So you are coupling every 16 ice timesteps? I think most other coupled systems couple the ice model every timestep because of averaging challenges when fields vary in area. Could the time averaging operation you're using for coupling be generating errors in the C-grid coupling that don't exist when running on the B-grid? |
|
@proteanplanet - Did you have a paper for how to examine the phase / coupling diagnostics ? |
|
Hi Tony. Sorry if I am mistaken or misunderstand, but I am not aware of anything like time-averaging in our setup. At the coupling frequency (one hour), instantaneous data are exchanged. that data is used for the next forecast hour until data is exchanged again. I have not tried with a higher coupling frequency, it would be interesting to see the impact on runtime if we coupled more often than one hour. |
|
Thanks @daveh150. So you are coupling instantaneous data and not conserving, that's clearer now. I guess you are resolving the diurnal cycle and inertial period pretty well, coupling once per hour. Maybe the timestepping/coupling is not playing a role. |
Hello @anton-seaice and @daveh150. I sent an email to Anton suggesting we meet online to go over your coupling configuration, since there are a number of things that aren't clear to me from the configuration described here. It would be easier to go over this in person so you can explain how your model is configured. I invited @eclare108213, and @apcraig may also like to listen in. |
|
For anyone playing along at home, daveh's setup and ours (ACCESS) are different. Yes we both use CICE6, but our coupling is different (and not sure about the ocean)
As noted the B-grid runs. This is the B-grid velocity from the same time that the c-grid version crashes, the sharp spike might be odd ? 
These are the plots in #992 (comment)
The exchange is done at the cell-centers. Its not something which is easy to change at this point.
I ran with
It would take some work do to this exactly as described. Although I think I achieved something similar: I set the thermodynamic timesteps in MOM and CICE and the coupling timestep to 2700, and then the MOM baroclinic timestep to 1350 and 
This is possible, we are currently using Ktens = 0 and the basal stress parameterisation is off.
There are no instabilities if run with advection turned off ( |
|
Hi @anton-seaice. There's a bunch of things potentially going on here that aren't related to anything mentioned on this thread. I can talk now if you are available. I've emailed a WebEx link. |
|
Hi guys! back from Europe and just to add to the noise, I have spoken to the SI3 guys in Paris and they have had the explicit Coriolis treatment for the C-grid implementation for about 10 years (i.e., the LIM3 period) with no adverse report. |
|
Sorry if I missed the discussion on coupling @proteanplanet suggested. I certainly would be interested in learning from y'all about your coupling setup and trying to correct/improve ours if another opportunity arises. |
|
@anton-seaice Thank you for posting your results. Just to clarify, did your run with ndtd=2 run to completion on C-Grid? From what I can tell that halves the dynamic timestep. Did you also have ndte=1200 in that test? Also If I understood the model did not crash with ndte=1200, though you still saw spikes in velocity? Thanks! |
|
Thanks @proteanplanet for helping. I would guess the problem is related to the ice-ocean stress. @anton-seaice I would like to suggest an additional test to verify if rheology is involved or not. Could you please use your default setup (and time steps) and set Pstar=0 (if kstrength=0) or Cf=0 (kstrength=1)? Doing this sets the rheology term to zero. Please restart just before the instability and show plot the velocities as above. |
|
I see two potential issues related to the ice-ocean stress:
For point 2), @daveh150, @dabail10 and @anton-seaice could you please confirm in you coupled setup that aice on the left is equal to aice on the right? In other words, are they at the same time level? |
|
@JFLemieux73 sorry I need some hand-holding here... can you point out where that is in the code? Please forgive my inexperience in finding where the momentum balance is. I am missing where the stresses are multiplied by aice. Thanks! |
Elizabeth suggested lengthening the coupling timestep, however in our set-up the coupling & model timesteps are linked. So instead I doubled the normal (thermodynamic) timesteps and then set ndtd=2 (and similar in the ocean). This was with the default ndte (120).
Yes correct - ice velocities were unstable / oscillating but recovered (probably due to the atmosphere wind-stress forcing changing) |
|
Thanks again @JFLemieux73
I am very suspicious of this but can't turn it into anything certain This is uvelE just before the instability. The cell which fails is to the adjacent (to the left) to a cell with 0 velocity. I guess the cell with 0 velocity (grey) is sensible enough, as the motion (positive to the right) is pushing the ice onto a land cell (white).
My feeling is this is worth fixing if possible even if unrelated to this instability.
I don't quite follow but before the crash aice is close to 1.
We use relative stresses, so I assume both tau_air and tau_io are calculated from wind/ocean/ice velocities from the previous timestep. We use a data atmosphere, so i would have to double check for the wind.
Results are completely smooth with Cf=0
|
|
Thanks @anton-seaice. The ice is compact so I don't think my points above (ice edge velocity and aice in water and air stresses) are good candidates. Can you try with your default setup (and time steps) and use kstrength=0? (don't set Pstar to zero). |
|
@eclare108213 I am wondering if it is not related to the stability issue described in Lipscomb, W. H., E. C. Hunke, W. Maslowski, and J. Jakacki (2007), Ridging, strength, and stability in high-resolution sea Maybe the Cgrid is more sensitive to this problem... |
|
@JFLemieux73 that's an excellent thought. It would make sense, too, that @anton-seaice's runs with 2 dynamics iterations per thermo time step would be more stable. Using the Hibler ice strength formulation should also be a good indicator, as you suggest. If this is the issue, then we can focus on the ridging parameters and if necessary, the participation / redistribution functions. |
|
Thanks @eclare108213. @anton-seaice if this is the problem it is due to the fact that the ice strength is explicit when solving for the momentum equation. It should be less sensitive with the Hibler formulation. Due to the time difference I would like to ask right away for a few more experiments. Could we see time series of the velocity (as above) for kstrength=0 with
Thanks! |
|
@anton-seaice, it is in region with some of the smallest grid cells? What are roughly the dimensions of the grid cells (in km) where it crashes? |
|
Thankyou for all the help on this @JFLemieux73 - the cells are small but not smallest. X dimension is ~7km and 7 dimension is ~11km. The instability grows to the left / negative X direction which is the small dimension in the grid in this area. The crash point shows excessive ridging prior to the crash in the original run, and the internal ice-stress history output also looks erroneous (as always - hard to know what is cause vs effect). Which I think is consistent with the suggestion of an instability related to ridging. From original run, blue is atmosphere-ice stress, orange is ocean-ice stress, green is internal stress, all in x direction at the c-grid u point. Note the internal stress field is larger magnitude than the surface stresses:
With the hibler ice strength formulation, the model runs however the velocities are greatly reduced:
|
|
Thanks a lot @anton-seaice. It really makes sense with what is described in the paper cited above (Lipscomb et al. 2007). I also have a paper about this is you are interested: Lemieux, J.-F., D. A. Knoll, M. Losch, and C. Girard (2014), A second-order accurate in time IMplicit-EXplicit (IMEX) integration scheme for sea ice dynamics, J. Comput. Phys., 263, 375–392, doi:10.1016/j.jcp.2014.01.010. No worries if the velocities are small at this location with the Hibler formulation. With Pstar=2.75e4 the ice is probably strong enough for not failing in compression (given the current ice and wind conditions). The fact that the blue and green curves are the same is because aice~1. As @eclare108213 wrote we might have to adjust the ridging parameters to fix the problem with kstrength=1. For the moment I would suggest to use kstrength=0. Note that this is what we use at ECCC in all our operational systems. With kstrength=1, we get unrealistic thicknesses in tidally active regions (it might be the same problem that is described here...). Thicknesses are more realistic with kstrength=0. We use Pstar=2.25e4. |
|
Thank you @anton-seaice for the clarification. Another data point: Before the holiday last week I ran a test with ndtd=2. I kept all our other parameters the same, including ndte=120 and our standard dt=225. This test ran successfully for 16 days. This seems to agree with @anton-seaice results. I'll make another test with ndtd=1 and kstrength=0 to see if that works too. Thank you also @JFLemieux73 for the info on using Pstar=2.25e4 with kstrengh=0 and the issue with thick ice in tidally active regions. I haven't looked at that in our system but will do so now! |
|
Hi @daveh150 - is your crash occurring at the coast with high winds (>10m/s) leading to high ridging / internal ice strength values ? |
|
Apologies for the very-long-post. The TLDR is approximately, the ridging when using C-grid kicks off the instability, which can be worsened by coupling with the ocean. Prior to the instability, the internal ice stress shows more sensitivity to high wind speeds when using the B-grid compared to the ice-velocities being more sensitive when using the C-grid. Thanks to @proteanplanet , @eclare108213 and @apcraig for help in identifying we have a two-timestep lag in the ocn-ice coupling. This is related to us running to ice&ocn models concurrently. I ran some more tests to try and determine the impact of this.
For all cases, there are a few timesteps where the ice internal stress term oscilates, however it doesn't grow unstably. There is minor patterning in the
These are two figures from the ice-ocn serial run:
These results are a bit counter-intuitive. Both reducing the timestep and running without feedbacks from the ocean model give results which shows some instability which recovers. The initial perturbation for the instability comes from ridging term in CICE. I went back and looked at the ridging term in the original B-grid runs. This configuration uses a data-atmosphere, and all the coupling is on the A-grid/center/tracer points (i.e.
However the ridging rate in the C-grid run is much much larger. In the c-grid results, prior to instability, the velocity is higher and the internal stress is lower prior to the ridging occuring and then the ridging term appears to cause the instability.
In the B-grid results, before the time of the instability (i.e before 16:00) it looks like internal stress term is more sensitive to changes in strair and the velocity is lower. Then ridging occurs in a more plausible/stable way.
|
|
Hi @anton-seaice, just to clarify...are these results above all with kstrength=1 (Rothrock)? |
|
Yes - with Rothrock and default Cf value (17) |
|
Hi @anton-seaice, Following our discussion yesterday, I think you should have a look at section 4 in Lipscomp et al. 2007. You could set up a similar idealized experiment. Let me know if you want to proceed and I could help for setting up the experiment. That would be a nice way to see if the C-grid is more sensitive than the B-grid for this instability issue. |
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:
(Built cice using 12dd204 , I can't see any commits since then they would impact this).
The text was updated successfully, but these errors were encountered: