fix metric terms in sea-ice stress divergence

[mjlosch]

What changes does this PR introduce?

Bug fix

What is the current behaviour?

When I implemented the stress divergence of sea ice, I thought that treating the individual components of the stress tensor like scalars would suffice to take care of all metric terms (see documentation and Losch et al 2010), but little did I know about tensor calculus: you also need to differentiate the unit vectors and this leads to terms that cannot be expressed in flux form. These terms are missing.

What is the new behaviour?

Two extra metric terms of the stress divergence each are added to the u- and v-equation:

u: k_{2}\sigma_{12} - k_{1}\sigma_{22}
v: k_{1}\sigma_{12} - k_{2}\sigma_{11}

For all solvers these terms are added to the RHS of the equation. For EVP and the implicit Krylov/JNFK solvers this is simple, but for the LSR solver this means that actually some part of the u/v-terms are not treated implicitly. This may seem inconsistent (as other metric terms are treating implicitly), but it is sooo much easier this way.

Does this PR introduce a breaking change?

All results on curvilinear grids change (a little). Fortunately most sea-ice simulations are carried out on cubed sphere or llc-grid where the Arctic is on some sort of 'near-cartesian' grid and the metric terms are small. As an example I attach results from a test on a 36km grid (carved out from a cubed sphere, basically coarse grained from the 18km grid that was used in Nguyen et al 2011, etc).

I ran 3 12-year-simulations: "default" (with the extra metric terms missing), "all metric terms", "all k1/2=0"; the latter means that the metric terms in the strain rates and the extra metric terms in the stress divergence are set to zero, but the metric terms that are covered by the finite volume discretisation are retained (I cannot change that without large incisions, which I do not want to make). From the timeseries of standard deviation of differences it's clear that the differences do not increase much after 2-3 years, so that it is probably OK to compare April 2012 (the 12th year).

In general the differences I see are small (order 1-10 centimetres, although some larger values appear every now and then). A little surprising: the order of magnitude of differences does not really change between the different comparisons, I had expected that the effect of the metric terms in the strain rates would be larger. My interpretation is that in these runs there are not really any systematic differences, just different realisations of "noise" that lead to a similar orders of magnitude of differences.

Other information:

We need to decide, what the defaults should be. There's a new runtime flag SEAICEuseExtraMetricTerms, that is currently .True., which can be set to .FALSE. to recover old results.

Suggested addition to tag-index

o pkg/seaice: add missing metric terms to stress divergence (bug fix)

• changes all results on curvilinear coordinates
• SEAICEuseExtraMetricTerms = F turns off new terms to recover old results

[mjlosch]

I can easily recover old (and wrong) results by setting SEAICEuseExtraMetricTerms = .FALSE., in

verification/global_ocean.cs32x15/input.icedyn/data.seaice
verification/global_ocean.cs32x15/input.seaice/data.seaice
verification/lab_sea/input.hb87/data.seaice
verification/lab_sea/input.salt_plume/data.seaice
verification/lab_sea/input/data.seaice
verification/seaice_itd/input.thermo/data.seaice
verification/seaice_obcs/input.regDenom/data.seaice
verification/seaice_obcs/input/data.seaice

and for the AD-expeirments

verification/global_ocean.cs32x15/input_ad.seaice/data.seaice
verification/global_ocean.cs32x15/input_ad.seaice_dynmix/data.seaice
verification/lab_sea/input_ad/data.seaice
verification/lab_sea/input_tap.noecco/data.seaice