Simplify problem to explain the C grid/remapping checkerboard instability

[JFLemieux73]

The goal here is to use numerical experiments to guide us for the mathematical explanation of the checkerboard instability.

Here are comments sent by email by Mats and Elizabeth:

'About demonstrating the checkerboard pattern mathematically: If linearised shallow water equations were applicable, I imagine investigating the dispersion relation of a C-grid discretisation with traditional incremental remapping would reveal stationary waves (= grid scale noise) due to the velocity averaging to obtain the corner velocities. Not necessarily trivial because of the geometric nature of incremental remapping, but I think simplifications could be done that still keep the essence. Since CICE don’t solve shallow water equations, I wonder if there is a thickness/concentration or strain rate feedback in the divergence of the ice stress tensor that, in a linearised sense, would trigger something similar.'

Mats

JF, I think this can be shown in 1D. I suspect it's the sequential spatial derivatives (of the velocity components and stresses and the ice pressure) that cause the problem. The pressure term of course brings in the ice thickness and area. The location of those quantities on the mesh can enhance the null space solution or damp it, because of the averaging as Mats says.
e

[JFLemieux73]

./cice.setup --case Eastblock -g gbox80 -m ppp6 -e intel -p 1x1 -s boxforcee,boxclosed,buildincremental

I changed things in ice_in:
npt=15 (days)
ndte = 1200
deltaminEVP = 2e-9
capping_method = 'sum'
kridge = 1
ktransport = 1
ice_data_type = 'eastblock'
ice_data_conc = 'p8'
Pstar = 1e4
grid_ice = 'C'

[JFLemieux73]

Ok I am able to reproduce what we had in issue #792. Here is aice after 15 days:

[JFLemieux73]

Note for me...The figs are here:

/home/jfl001/Lemieux_et_al_CICE_Cgrid/Checkerboard/Eastblock/FIGS

The problem is already a bit simplified as Coriolis is zero (Exp1).

Exp2: let's make the ocean-stress linear. In stepu_C and stepv_C I change

vrel = aiX(i,j)rhowCw(i,j)*sqrt((uocn(i,j) - uold)**2 + &
(vocn(i,j) - vold)**2) ! m/s

by vrel = aiX(i,j)rhowCw(i,j)*0.05d0

Here is aice after 15 days:

[JFLemieux73]

Exp3: the same term can be simpler. I set aiX(i,j)=1...vrel = rhowCw(i,j)*0.05d0. Here is aice after 15 days:

The checkerboard is still there...

[JFLemieux73]

Exp4: The default value for the ellipse ratio is 2 (e_yieldcurve = e_plasticpot = 2).

If I set the ellipse ratio to a very large value, shear stress should disappear (i.e. eta = 0).

Here is aice after 15 days for e_yieldcurve = e_plasticpot = 2000:

Different but still there...

[JFLemieux73]

Exp5: Let's go back to e_yieldcurve = e_plasticpot = 2. In stepu_C and stepv_C, vold and vvel are forced to be zero. Only u components are non-zero.

Note that I changed the min value (0.8 instead of 0.9). So the problem is still there. Is it caused by the du/dy and or dP/dy?

[JFLemieux73]

Exp6: Do we still see the problem with set ice_data_conc = 'c1'?

Here is aice after 15 days:

Still there...

[JFLemieux73]

ASIDE: I ran the case above with kstrength=1 and didnt see the noise...Maybe it is just because the ice is more rigid with kstrength=1. I went back to the first exp above (with ice_data_conc = 'p8') and the checkerboard pattern is present with kstrength=1 (Cf was set to 6) as shown below:

So it is not a question of kstrength=0 or 1 but maybe kstrength=0 is more prone to lead to the checkerboard.

[JFLemieux73]

@phil-blain helped me to continue this investigation. The idea is to simplify the problem before doing the stability analysis. We ran the eastblock test case with C-grid and l_fixed_area = T (no checkerboard) and stored the velocity field (uvel,vvel,uvelE and vvelN) after 15 days. We then used these constant (in time but spatially varying fields) for the remapping. We ran the eastblock test case with l_fixed_area = T or l_fixed_area = F. The results look the same and there is no checkerboard observed in both simulations.

we used this branch for the test:

Checkerboard_Lemieux_et_al_2024

The test cases were ran here:

/home/jfl001/Lemieux_et_al_CICE_Cgrid/CHECKERBOARD/Investigate_checkerboard

[JFLemieux73]

The next step is to run with l_fixed_area = T and store the velocity fields for all the time levels. These velocity fields will then be used for a simulation with l_fixed_area = F. If the checkerboard does not appear with this test case would mean that it is really an interaction between the dynamics and the transport (more complicated to understand...).

[JFLemieux73]

Ok Philippe and I just worked on this. As mentioned above, a run with l_fixed_area = F used the velocity fields (for transport) at all time levels from a run with l_fixed_area = T. The checkerboard disappears...So clearly the checkerboard is due to an interaction between the dynamics and remapping. This is more complicated to explain mathematically...

We used this branch for the test:

Checkerboard_B_Lemieux_et_al_2024

The test cases were ran here:

/home/jfl001/Lemieux_et_al_CICE_Cgrid/CHECKERBOARD/Investigate_checkerboard_B

[JFLemieux73]

Ok going back to experiments with both transport and dynamics.

here: /home/jfl001/Lemieux_et_al_CICE_Cgrid/CHECKERBOARD/VISCOUS/VISCOUS_EB

v=0
e=1 (eta=zeta)
viscous (in visc_replpress):

af=1d-03
zetax2 = afstrength/(DminArea)
rep_prs = strength
etax2 = epp2izetax2

We still see the checkerboard with this simplified exp:

[eclare108213]

We still see the checkerboard with this simplified exp

This is with l_fixed_area = T? and using VP? My initial thought was that elastic waves were kicking off the checkerboard in the high-concentration region, but if this is VP then there's something about the VP discretization (or VP itself) that the advection can't fix. It's entirely possible that the VP discretization creates a checkerboard all by itself, or that it still produces a null space solution when run with the advection. I went through LOTS of EVP discretizations that checkerboarded on the B grid with advection turned off. Does this checkerboard appear with advection turned off? Maybe not, if ice needs to converge at the wall for it to become apparent.

That said, the other possibility is that the VP solution in this test case is not well-defined, e.g. see appendix B in https://www.sciencedirect.com/science/article/pii/S0021999101967105. Two-dimensional solutions usually overcome this problem with VP, but in this test case the solution might not be "strong enough" (whatever that means) in 2D to overcome the tendency to checkerboard. That's pretty hand-wavy, but I'm quite sure we can come up with a simple 1D, theoretical explanation of the checkerboarding in all of these cases.

[JFLemieux73]

Hi Elizabeth,

This is the same good old checkerboard we have seen before (see issues #791, #792 and #811). This is C-grid (EVP because VP does not exist yet for the C-grid) with advection turned on. All the experiments above that have the checkerboard use l_fixed_area = F.

We know how to fix the problem by using l_fixed_area = T.

Intuitively we understand why this fixes the problem: l_fixed_area = T ensures that the divergence implied by the remapping is consistent with the divergence implied by the EVP solver.

The goal here is to simplify the problem so that we can do a stability analysis for the paper...i.e. come up with a more solid explanation instead of intuition.

[eclare108213]

I thought you had found a new case that was still checkerboarding even with the new fixes. Thanks for the clarification!

[JFLemieux73]

yes sorry my notes are not really clear!

[JFLemieux73]

Ok I can further simplify the problem by setting eta=0 (no shear stress). The problem still appears with l_fixed_area =F (no checkerboard with l_fixed_area = T). In conclusion we could have for the stability analysis:

v=0
no water stress
ice_data_conc = 'c1'
f=0
af=1d-03
zetax2 = af*strength/(DminArea)
rep_prs = strength
etax2 = 0 (or e tends toward infinity)

[JFLemieux73]

note: as the noise appears at the wall, maybe it is simpler to assume it is plastic. I have inscreased ndte to 2400, reduced elasticDamp to 0.12 and set deltaminEVP=2-12 (very small...). The checkerboard is still present with l_fixed_area = F.

[JFLemieux73]

Problem has been simplified for modal analysis. Checkerboard is associated with the spatial averaging of C-grid velocities at the U point. See Lemieux et al., in prep for details.

[JFLemieux73]

Special experiment to answer to reviewer 1:

/cice.setup --case Eastblock -g gbox80 -m ppp6 -e intel -p 1x1 -s boxforcee,boxclosed,buildincremental

I changed things in ice_in:
npt=1 (days)
ndte = 2400
deltaminEVP = 2e-11
elasticDamp = 0.12d0
capping_method = 'sum'
kridge = 1
ktransport = 1
ice_data_type = 'eastblock'
ice_data_conc = 'c1'
Pstar = 1e4
grid_ice = 'C'

I want to look at the beginning of the checkerboard formation. So that it arrives 'quickly' I set ice_data_conc = 'c1'. This means that the ice is strong right from the start (rheology term is important).

[JFLemieux73]

It takes a few hours to see the checkerboard. This is aice after 13h:

Note that to see the checkerboard, the range is between 0.996 and 1.0 (for aice)

[JFLemieux73]

Notice the checkerboard that starts to develop close to the wall. I calculated delta to see if it is plastic or viscous. I converted the values to s-1. It is the beginning of the simulation and the ice is still pilling up against the eastern wall. It is clearly plastic close to the wall and in the region of the checkerboard. The figure of log(delta) shows values between 10^-7 s-1 to 10^-5 s-1in the region of the checkerboard (recall that deltamin=2e-11 s-1).

[JFLemieux73]

The divergence (s-1) is also consistent with ice converging at the wall and in the plastic regime. This is also after 13h. The displayed range is -1e-7 to 0. The zone of convergence) 'contains' the region of checkerboard.

[JFLemieux73]

Note that during the period between 0h and 13h when the initial checkerboard shown above was created, the same conclusions apply for delta (plastic) and the divergence (strong convergence in the region close to the wall).