DIVA is inaccurate when flwa varies strongly vertically

[billsacks]

From @matthewhoffman on October 4, 2016 20:33

The DIVA and BP velocity solvers calculate quite similar results for standard test cases (e.g. ISMIP-HOM), but recent testing has revealed that they generate quite different results when flwa varies significantly vertically.

Here is an example of two runs with identical setup except for which solver is used:

I isolated the issue - it only occurs when flwa varies vertically as in the above example. If flwa is made vertically constant, then the two solvers give similar answers:

It makes sense that the results of DIVA would be sensitive to the details of how flwa (or the effective viscosity) is integrated vertically. This may take some careful thought and testing to resolve in a way that allows DIVA to yield results comparable to BP. For example, flwa can vary vertically by two or more orders of magnitude, meaning a straight arithmetic average may be in appropriate. Similarly, consideration may be needed for the vertical arrangement of flwa values - presumably soft ice at the bed shuold affect the depth-integrated effective viscosity much more than soft ice near the surface.

For the record, I was using this flwa profile (with uniform vertical levels):

for i in range(nx):
   for j in range(ny):
     flwastag[0,:,j,i] = (
6.24821E-17,
5.98918E-17,
3.78746E-17,
1.86614E-17,
9.50214E-18,
7.66194E-18,
7.9094E-18,
1.06029E-17,
3.31514E-17,
#6.71515E-17)  #E=1
2.01e-16)  #E=3

which comes from the temperature profile in this paper: Ryser, C., M. P. Lüthi, L. C. Andrews, M. J. Hoffman, G. A. Catania, R. L. Hawley, T. A. Neumann, and S. S. Kristensen (2014), Sustained high basal motion of the Greenland ice sheet revealed by borehole deformation, J. Glaciol., 60(222), 647–660, doi:10.3189/2014JoG13J196. and using the updated Cuffey and Paterson flwa formula.

Copied from original issue: E3SM-Project/cism-piscees#61

[billsacks]

From @gunterl on October 4, 2016 20:46

Hi Matt,

Thank you for sharing these results, it's nice to see more work is done using more dycores. So far all the MIPS have been focused on obtaining more knowledge about the effect of basal friction law on ice sheet dynamics and ice has been considered to be isothermal. Just out of curiosity, where are you measuring the velocity that create these plots? (GL, ice divide, average).

Hope all is well.

Gunter

On Tue, Oct 4, 2016 at 2:34 PM, Matt Hoffman notifications@github.com wrote:

The DIVA and BP velocity solvers calculate quite similar results for standard test cases (e.g. ISMIP-HOM), but recent testing has revealed that they generate quite different results when flwa varies significantly vertically.

Here is an example of two runs with identical setup except for which solver is used: [image: screen shot 2016-10-04 at 11 05 27 am] https://cloud.githubusercontent.com/assets/4182034/19091091/61d8a0c4-8a3e-11e6-9233-eb580986187c.png

I isolated the issue - it only occurs when flwa varies vertically as in the above example. If flwa is made vertically constant, then the two solvers give similar answers: [image: screen shot 2016-10-04 at 10 59 45 am] https://cloud.githubusercontent.com/assets/4182034/19091144/960b1b4c-8a3e-11e6-8cf7-3f19626b4341.png

It makes sense that the results of DIVA would be sensitive to the details of how flwa (or the effective viscosity) is integrated vertically. This may take some careful thought and testing to resolve in a way that allows DIVA to yield results comparable to BP. For example, flwa can vary vertically by two or more orders of magnitude, meaning a straight arithmetic average may be in appropriate. Similarly, consideration may be needed for the vertical arrangement of flwa values - presumably soft ice at the bed shuold affect the depth-integrated effective viscosity much more than soft ice near the surface.

For the record, I was using this flwa profile (with uniform vertical levels):

for i in range(nx): for j in range(ny): flwastag[0,:,j,i] = ( 6.24821E-17, 5.98918E-17, 3.78746E-17, 1.86614E-17, 9.50214E-18, 7.66194E-18, 7.9094E-18, 1.06029E-17, 3.31514E-17,

6.71515E-17) #E=1

2.01e-16) #E=3

which comes from the temperature profile in this paper: Ryser, C., M. P. Lüthi, L. C. Andrews, M. J. Hoffman, G. A. Catania, R. L. Hawley, T. A. Neumann, and S. S. Kristensen (2014), Sustained high basal motion of the Greenland ice sheet revealed by borehole deformation, J. Glaciol., 60(222), 647–660, doi:10.3189/2014JoG13J196. and using the updated Cuffey and Paterson flwa formula.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/ACME-Climate/cism-piscees/issues/61, or mute the thread https://github.com/notifications/unsubscribe-auth/AH0K-MZHwOpXrpGu-hktW4Jsu21jmIqvks5qwrg4gaJpZM4KOIDR .

[billsacks]

From @matthewhoffman on October 4, 2016 23:4

@gunterl , thanks for the feedback. Yes, this was interesting to discover under a more realistic configuration and seems like an important detail for practical application.

My model domain is very simple. Here is a longitudinal profile of thickness (this is a "plastic" glacier shape where Tau_d=10^5 Pa everywhere):

At the upstream and downstream ends (far enough from the center to not matter), I apply Dirichlet b.c. of u=100 m/yr. Laterally, the domain is uniform and periodic. The results above are the surface speed from the cell in the center of the domain.

[billsacks]

From @gunterl on October 5, 2016 17:9

Hi Matt,

Thanks for continuing the conversation. My next question is if your test case is a marine terminated glacier? It doesn't look that way at a first glimpse. In which case have you tried to run SIA as well and see how it compares to DIVA? This would also be interesting to see.

On Tue, Oct 4, 2016 at 5:04 PM, Matt Hoffman notifications@github.com wrote:

@gunterl https://github.com/gunterl , thanks for the feedback. Yes, this was interesting to discover under a more realistic configuration and seems like an important detail for practical application.

My model domain is very simple. Here is a longitudinal profile of thickness (this is a "plastic" glacier shape where Tau_d=10^5 Pa everywhere): [image: screen shot 2016-10-04 at 4 59 20 pm] https://cloud.githubusercontent.com/assets/4182034/19095644/23d782a2-8a54-11e6-927c-e408a52cb72b.png At the upstream and downstream ends (far enough from the center to not matter), I apply Dirichlet b.c. of u=100 m/yr. Laterally, the domain is uniform and periodic. The results above are the surface speed from the cell in the center of the domain.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/ACME-Climate/cism-piscees/issues/61#issuecomment-251538782, or mute the thread https://github.com/notifications/unsubscribe-auth/AH0K-M8pxwGGyAIdIqzUmfjsHJyWQXURks5qwttsgaJpZM4KOIDR .

[billsacks]

From @gunterl on October 5, 2016 17:16

Or SSA for that matter?

On Tue, Oct 4, 2016 at 5:04 PM, Matt Hoffman notifications@github.com wrote:

@gunterl https://github.com/gunterl , thanks for the feedback. Yes, this was interesting to discover under a more realistic configuration and seems like an important detail for practical application.

My model domain is very simple. Here is a longitudinal profile of thickness (this is a "plastic" glacier shape where Tau_d=10^5 Pa everywhere): [image: screen shot 2016-10-04 at 4 59 20 pm] https://cloud.githubusercontent.com/assets/4182034/19095644/23d782a2-8a54-11e6-927c-e408a52cb72b.png At the upstream and downstream ends (far enough from the center to not matter), I apply Dirichlet b.c. of u=100 m/yr. Laterally, the domain is uniform and periodic. The results above are the surface speed from the cell in the center of the domain.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/ACME-Climate/cism-piscees/issues/61#issuecomment-251538782, or mute the thread https://github.com/notifications/unsubscribe-auth/AH0K-M8pxwGGyAIdIqzUmfjsHJyWQXURks5qwttsgaJpZM4KOIDR .

[billsacks]

From @matthewhoffman on October 5, 2016 20:14

No, this was totally grounded. I have not tried any other solvers on it. The specific setup was using the Schoof basal friction law with a prescribed time series of effective pressure. There was significant amounts of both sliding and deformation. It was kind of an ad hoc setup that I was using for something else, and I stumbled onto the problem unintentionally.

@DanFMartin and I have been discussing ideas about setting up a more standardized test that could be used to assess these differences. I was thinking a very simple geometry might be the "slab" test case I set up in CISM a couple years ago but never quite fully verified: https://github.com/ACME-Climate/cism-piscees/tree/develop/tests/higher-order/slab

It is described in sections 5.1-2 of: J.K. Dukowicz, 2012. Reformulating the full-Stokes ice sheet model for a more efficient computational solution. The Cryosphere, 6, 21-34. www.the-cryosphere.net/6/21/2012/

Then we would impose a temperature (or alternatively flwa) field that varies strongly vertically, but is horizontally uniform. I was using essentially the borehole temperature data from here: Lüthi, M. P., C. Ryser, L. C. Andrews, G. A. Catania, M. Funk, R. L. Hawley, M. J. Hoffman, and T. Neumann (2015), Heat sources within the Greenland Ice Sheet: dissipation, temperate paleo-firn and cryo-hydrologic warming, The Cryosphere, 9, 245–253, doi:10.5194/tc-9-245-2015. We could fit a curve to it to yield an easy to apply analytic expression that is PMP at surface and bed and cold in the interior (-10 or -15 C).

@gunterl , @DanFMartin , let me know if you have any other suggestions about how to set it up.

[billsacks]

From @gunterl on October 5, 2016 20:50

Hi Matt,

So I believe C (the powerlaw coefficient from Schoof's law) is spatially varying then, right? Or are you using something proportional to effective pressure somewhere else?

On Wed, Oct 5, 2016 at 2:14 PM, Matt Hoffman notifications@github.com wrote:

No, this was totally grounded. I have not tried any other solvers on it. The specific setup was using the Schoof basal friction law with a prescribed time series of effective pressure. There was significant amounts of both sliding and deformation. It was kind of an ad hoc setup that I was using for something else, and I stumbled onto the problem unintentionally.

@DanFMartin https://github.com/DanFMartin and I have been discussing ideas about setting up a more standardized test that could be used to assess these differences. I was thinking a very simple geometry might be the "slab" test case I set up in CISM a couple years ago but never quite fully verified: https://github.com/ACME-Climate/cism-piscees/tree/ develop/tests/higher-order/slab

It is described in sections 5.1-2 of: J.K. Dukowicz, 2012. Reformulating the full-Stokes ice sheet model for a more efficient computational solution. The Cryosphere, 6, 21-34. www.the-cryosphere.net/6/21/2012/

Then we would impose a temperature (or alternatively flwa) field that varies strongly vertically, but is horizontally uniform. I was using essentially the borehole temperature data from here: Lüthi, M. P., C. Ryser, L. C. Andrews, G. A. Catania, M. Funk, R. L. Hawley, M. J. Hoffman, and T. Neumann (2015), Heat sources within the Greenland Ice Sheet: dissipation, temperate paleo-firn and cryo-hydrologic warming, The Cryosphere, 9, 245–253, doi:10.5194/tc-9-245-2015. We could fit a curve to it to yield an easy to apply analytic expression that is PMP at surface and bed and cold in the interior (-10 or -15 C).

@gunterl https://github.com/gunterl , @DanFMartin https://github.com/DanFMartin , let me know if you have any other suggestions about how to set it up.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/ACME-Climate/cism-piscees/issues/61#issuecomment-251785947, or mute the thread https://github.com/notifications/unsubscribe-auth/AH0K-O3U4RNb7fx5pUxv6FR537TPrPc6ks5qxAUkgaJpZM4KOIDR .

[billsacks]

From @matthewhoffman on October 5, 2016 21:39

In my test, the Coulomb friction coefficent C is a constant, but effective pressure is a time varying 2d field that came from my subglacial hydro model.

I don't think any of that is relevant to producing the inaccurate DIVA velocities - I think that is all due to a depth-varying temperature field. So I'd suggest we start with something like no sliding or a moderate beta (say, 1000) and confirm we can reproduce this behavior in a simpler setup.

[billsacks]

From @gunterl on October 5, 2016 21:49

Yes I agree. I was just trying to get a better sense of the dynamics.

On Wed, Oct 5, 2016 at 3:39 PM, Matt Hoffman notifications@github.com wrote:

In my test, the Coulomb friction coefficent C is a constant, but effective pressure is a time varying 2d field that came from my subglacial hydro model.

I don't think any of that is relevant to producing the inaccurate DIVA velocities - I think that is all due to a depth-varying temperature field. So I'd suggest we start with something like no sliding or a moderate beta (say, 1000) and confirm we can reproduce this behavior in a simpler setup.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/ACME-Climate/cism-piscees/issues/61#issuecomment-251807501, or mute the thread https://github.com/notifications/unsubscribe-auth/AH0K-NkG-phn3kz-SVgPZz8EJvFsIGOnks5qxBkKgaJpZM4KOIDR .

[billsacks]

From @DanFMartin on October 5, 2016 23:46

I agree -- Let's start with a moderate beta and go from there. We can add complexity if needed later.

Sent from my iPhone

On Oct 5, 2016, at 2:39 PM, Matt Hoffman notifications@github.com wrote:

In my test, the Coulomb friction coefficent C is a constant, but effective pressure is a time varying 2d field that came from my subglacial hydro model.

I don't think any of that is relevant to producing the inaccurate DIVA velocities - I think that is all due to a depth-varying temperature field. So I'd suggest we start with something like no sliding or a moderate beta (say, 1000) and confirm we can reproduce this behavior in a simpler setup.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[billsacks]

From @matthewhoffman on October 25, 2016 22:31

I happened on this discussion in a Doug Brinkerhoff paper:

Second, depth-integrated models typically make the assumption that the vertical averages of strain rates (and by extension viscosity) are a sufficient approximation to the true vertical integral of the horizontal stress divergence. As demonstrated by Schoof and Hindmarsh [2010], this approximation holds in the high slip ratio context at all aspect ratios and in the low slip ratio case at shallow aspect ratios. It does not hold when aspect ratios are high and the slip ratio is low. Rather, our model uses numerical quadrature combined with a polynomial ansatz and can accommodate a vertically varying viscosity term in the membrane stresses. This also allows us to include a vertically varying temperature field with a reasonable amount of degradation in accuracy relative to the LMLa equations, a critical feature if this stress balance is tobe used for polythermal ice.

Brinkerhoff, D. J., and J. V Johnson (2015), Dynamics of thermally induced ice streams simulated with a higher-order flow model, J. Geophys. Res. Earth Surf., 120, doi:10.1002/2015JF003499.Received.

[billsacks]

From @whlipscomb on October 25, 2016 22:58

@matthewhoffman, Thanks for the lead. I'll re-read the paper and follow up with Doug or Jesse if needed.