Investigate thick ice in C-grid simulations

[JFLemieux73]

gx3 simulations show thicker ice with C-grid (compared to B-grid). See issue #702. Many idealized test cases show comparable results between C and B grid. Why is it not the case for gx3? Here are a few things to investigate:

-Coriolis
-advection
-deformations (used for redistribution)
-thermo
-gx 3 restart

I will use commit bffe4e8. I will use a combination of test cases to investigate this.

[JFLemieux73]

boxsyme with:
ndte=1200
capping=0
visc_method = 'avg_zeta'
deltaminEVP =2e-9
Pstar=1e4
elasticDamp=0.36 (default)

I finally set coriolis='constant'. Here is the uvel field after one day for the B grid:

and the uvelE field (C grid)

[JFLemieux73]

The vvel and vvelN fields are respectively:

[JFLemieux73]

The patterns are similar but not exactly the same. Note that the 4 figures above are steady-state solutions (after one day). Same fields after one time step (1 hour) are also very similar (not shown).

[JFLemieux73]

I now want to test advection:
./cice.setup -c boxwallp5 -m daley -e intel -s debug,boxwallp5,gridc -g gbox80
ndte=1200
capping=0
visc_method = 'avg_zeta'
deltaminEVP =2e-9
Pstar=1e4
elasticDamp=0.1 (I decreased it...)
I have set close_boundaries=.true. and ktransport=1.

@eclare108213 does it make sense to have ktransport = 1 if kridge = -1 ?

[eclare108213]

Generally transport on with ridging off will not work - convergence will make aice > 1 and the code is likely to crash. However if you have a setup without convergence, e.g. a uniform velocity that just advects a blob of ice through the domain, or maybe a divergent flow field with open boundaries where the ice could just disappear from the domain, then it could work.

[JFLemieux73]

Thanks. Ok I will use kridge=1. Here is the B grid hi after 7 days (hist_avg=.true., daily outputs)

and the C grid hi after 7 days:

[JFLemieux73]

It looks like there is some noise in both south-east north-east corners in the C grid solution. I change the range of values: 1.5 to 1.80.

[JFLemieux73]

So we see the noise in both corners. This might be the problem we are seeing in the gx3 simulations (there are a lot of corners...). Is it due to advection or ridging?

[JFLemieux73]

aice does exceed 1.0 during the first three days. I want to separate the effect of advection and ridging. Below is hi after 3 days with ktransport=1 but kridge=-1. I set the range to 1.0 to 1.7 m. For the B grid:

and for the C grid:

[JFLemieux73]

The patterns are the same. For both the thick band on the wall is 1.6626 m thick. Note however that the thickest ice is ~1.68 m for B and ~1.72 for C (is that an issue?). My impression though is that the problem is in the ridging...probably in the deformation_T subroutine. Note that we dont calc DeltaT in deformationT the same way we calc DeltaT for the viscosities (Bouillon approach). This is a legacy of the CD grid...I guess this is what we need to fix.

[JFLemieux73]

Hum...but the ridging should not change hi (mean thickness). Most of the tests done for the C grid are with uniform conditions. Is it possible there is a problem when aice, hi (and strength) vary spatially?

[eclare108213]

In the gx3 cases, are you running with the Hibler strength or the complicated ITD formula?

[JFLemieux73]

It is kstrength=1 (complicated formula) in gx3 while it is kstrength=0 (Hibler) here. I could test both.

[dabail10]

I'm looking here at deformations_T. Have you updated to the latest code where Tony has refactored these? I think it is doing the right thing for the C grid. It calls strain_rates_T which is now a module procedure that includes two versions of this. The one that should be called for the C grid is strain_rates_Tdtsd. This computes divT, tensionT, shearT, and DeltaT.

[JFLemieux73]

Yes but it computes shearT as:

   shearT(i,j) = (dxT(i,j)**2)*(uvelN(i,j)/dxN(i,j) - uvelN(i,j-1)/dxN(i,j-1)) &
                 + (dyT(i,j)**2)*(vvelE(i,j)/dyE(i,j) - vvelE(i-1,j)/dyE(i-1,j))

That's not the Bouillon approach that we use for DeltaT for the viscosities at the T point.

[JFLemieux73]

In fact we don't even need to recalculate the strain rates as they are already stored. I will try that.

[JFLemieux73]

I now use a smaller boxwallp5 test case (12x12) for debugging. Here is a snapshot of hi after 10 days for B grid:

and for the C grid

[JFLemieux73]

Not sure if this helping...

[JFLemieux73]

I also ran it with upwind advection.

B grid:

C grid:

[JFLemieux73]

So it is thicker with the C grid. I am going back to the first few time levels to see what is going on at the beginning.

[JFLemieux73]

I have coded a new deformations_T subroutine for which DeltaT is calc using shearTsqr (the avg of the shearU**2). As expected this does not change the volume. I ran the same test as above with upwind and the hi field is the same. This could however change the results when thermo and /or kstrength=1 are used. This will need to be tested in gx3 experiments for example.

[JFLemieux73]

Here is dvidtd (change of volume due to dynamic) for the B grid after one time step (1 hour):

C grid:

Note that this is with upwind. We already see a difference...

[JFLemieux73]

I guess it makes sense as uvel (for B grid) is interpolated to the E point (uvel + 0)/2 for the upwind advection...

[JFLemieux73]

But I guess with remap it should be very similar:
B grid:

C grid:

[JFLemieux73]

Going back to the 80x80 boxwallp5 test case. I want to see the impact over longer sims. Here is hi after 29 days with upwind advection.
B grid:

C grid:

[JFLemieux73]

Here are also the pressure (sigP) fields after 29 days:
B grid:

C grid:

[JFLemieux73]

They look good. The pressure increases toward the coast to balance the wind stress. I am going to do the same thing with the remapping scheme.

[JFLemieux73]

Here is hi (C-grid) after 29 days with remapping:

[JFLemieux73]

It's a bit less smooth near the corners but it's hard to tell if something is wrong. I think I should focus on gx3 simulations. I will use issue #702 for this.

[JFLemieux73]

I also tested a northeast forcing. Here is the B grid hi after one month:

and C grid hi:

They are very similar.

[JFLemieux73]

Idealized test cases do not show this thickening. This should be investigated in gx3 and gx1 simulations (see issue #702).