Unexpected differences with Hybrid model compared to Elmer/Ice with Cuffey and Paterson A

[dlilien]

I am hoping I am just missing something silly here, but I am getting a substantial difference in velocity when using the Icepack hybrid model versus Elmer/Ice on test cases where I think there should be little or no difference. The SIA model in icepack matches Elmer.

The difference is caused by the rate factor--dividing the default value in icepack by a factor of ~2 ** 2.5 seems to get things almost identical. I have attached a minimal working example where the difference is seen; I am simply using built in rate_factor in icepack and I am using the same values in Elmer/Ice (or, rather, their suggested values but I have separately confirmed that those match Cuffey and Paterson). A reasonably recent Elmer/Ice is probably required for the viscosity formulation, but not bleeding edge and no need for any add-ins other than the elmerice.

For some reason I have been having problems with point evaluation in firedrake recently--for that reason, the last two plots in the notebook are not as clear as the could be, but I think it shows things pretty clearly anyway. Hopefully @danshapero can help point out if I am doing something wrong or if this is a bug in the hybrid solver. mwe.tar.gz

[danshapero]

One of the main tests for the hybrid solver is to take the vertical degree to be 0 and checking that it gives the same result as SSA. Given that this test passes, it at least confirms that the two models in icepack are consistent with each other, so if there's a bug in the hybrid model then there should also be a bug in SSA. Is it easy for you to repeat this experiment using the SSA in both and comparing the results? I ask because It looks like the Elmer simulation is using the Stokes model and not Blatter-Pattyn. This shouldn't be explainable by a scaling of 2^2.5 but it's still not completely an apples-to-apples comparison.

My hunch is that all this explainable by slightly different conventions about how you define the second invariant of a rank-2 tensor, some of which do or don't have some factors of 2 floating around. I can try to dig around in the Elmer source code to see what they do. I think that I followed pretty closely what's in both the Greve and Blatter or Cuffey and Paterson books but I won't rule out the possibility that I mis-transcribed an equation somewhere.

[RupertGladstone]

I heard from Thomas Zwinger that there is now a Blatter/Pattyn model available in Elmer/Ice, in case you wanted a direct comparison. Not sure if it is documented yet, but we can ask Thomas directly if you're interested.

Daniel's hunch sounds more likely to me, and getting Thomas' input for digging in source code could also be an option here. If one of these codes has a bugged implementation of viscosity I'm sure we all want to know!

BTW I just ran your test case for Elmer and get velocity norm 0.4 and max vel 1.0. Presume that is the same you get? This is with a very recent version of Elmer on my lapotp.

Note that you could also set a rate factor directly in Elmer instead of setting temperature; another way of making the comparison more direct, if your icepack rate factor is being set directly.

[dlilien]

My velocity norm and max match those numbers, @RupertGladstone.

I should have time to look at this in more detail next week. I think I can easily set up identical simulations with SSA in both Elmer and icepack. From the testing I did so far, though, I expect they will largely agree; using realistic geometries, in sliding dominated areas, the solutions from full Stokes in Elmer and Blatter/Pattyn in icepack were very similar.

[dlilien]

OK, I can confirm that icepack and Elmer match for SSA simulations. I still do not really understand why things look so different from Elmer full-Stokes and from the SIA when I use the hybrid solver with a really high value for C, but I think maybe I should just stop worrying about it. SSA comparison added here: mwe.tar.gz

[bryanvriel]

Hi @danshapero, I think the discrepancy may be due to the computation of the effective strain at this line:

def _effective_strain_rate(ε_x, ε_z, ε_min):
    return sqrt((inner(ε_x, ε_x) + trace(ε_x) ** 2 + inner(ε_z, ε_z) + ε_min**2) / 2)

I think inner(ε_z, ε_z) shouldn't be multiplied by 1/2 since ε_z already has a factor of 0.5 applied. Without the 1/2, I get an expanded effective strain rate that matches Equation 5.71 in Greve and Blatter.

In a flowband experiment with the hybrid model with a high C, if I pull inner(ε_z, ε_z) out of the parentheses as:

def _effective_strain_rate(ε_x, ε_z, ε_min):
    return sqrt((inner(ε_x, ε_x) + trace(ε_x) ** 2) / 2 + inner(ε_z, ε_z) + ε_min**2)

then I get a surface velocity that closely matches the theoretical SIA velocity. However, if I use the effective strain rate as currently implemented, I get a surface velocity that's 4x higher than SIA.

[dlilien]

Looking at Equation 14 in Pattyn, 2003, I think that @bryanvriel is right that there is a factor of 2 incorrect on ε_z in the effective strain rate. I agree with him that it probably makes sense to pull ε_min out of the parentheses as well so that the definition is intuitive. With Bryan's change, I get velocities that are within 50% of Elmer/Ice and icepack's SIA when C is large, as opposed to something like a factor of 5 before.

[danshapero]

Yes you're both exactly right, I can't believe I missed this -- it's the double counting because of the symmetry of the strain rate tensor. I'll write a test case and a fix today.

[danshapero]

I made a PR with the fix. If either or both of you can look at the test case and tell me whether or not you're convinced that would be great, or run your test cases with the patch applied. I'll merge it by end of day tomorrow either way.