To do idea 3: simulate ice flow to get H and u

[jkingslake]

Use a simple ice flow model to provide ice thickness (and hence basic hydraulic gradient psi) and sliding speed and grounding line speed.

A place to start could be Alex Robel's code here: https://github.com/aarobel/Ice-Stream-Flowline-Model-Thermal You could start by simulating a steady-state H, grounding line position and u, using this ice flow model and then using these as inputs to the hydrology model. (to get a non-oscillating solution to the ice flow model, @aarobel notes here that you will need to change some parameters.)

This would amount to one-way coupling between ice flow and hydrology. The next step would be to use the effective pressure to control sliding in the model, which will require some adjustment to the model. @aarobel, may have some ideas about how to do that :-)

[aarobel]

Hi @jkingslake! Sounds like a cool idea. If the goal is simply to have thickness, velocity, GL position from a flowline model of a marine-terminating glacier, then I would recommend using this model instead: https://github.com/aarobel/SSAsimpleM. Its simpler and better benchmarked against MISMIP and a good jumping off point for doing fancier things. However, if you'd like to do something with meltwater production and till, then the thermal model is definitely the way to do. Happy to provide more feedback or ideas if need be.

[jkingslake]

Thanks for the suggestion Alex! Yes that probably makes sense as a starting point. The overarching goal for this project is to couple hydrology and ice flow in as many ways as we can think of. You can see some of the other ideas in the other issues in this repo. So, ultimately it would be nice to have N effect the ice flow through the sliding. Would the simpler code you suggest be a good starting point for including a really simple N-dependent sliding law, do you think?

On Thu, Sep 30, 2021, 11:09 AM Alexander Robel @.***> wrote:

Hi @jkingslake https://github.com/jkingslake! Sounds like a cool idea. If the goal is simply to have thickness, velocity, GL position from a flowline model of a marine-terminating glacier, then I would recommend using this model instead: https://github.com/aarobel/SSAsimpleM. Its simpler and better benchmarked against MISMIP and a good jumping off point for doing fancier things. However, if you'd like to do something with meltwater production and till, then the thermal model is definitely the way to do. Happy to provide more feedback or ideas if need be.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/ldeo-glaciology/nye_eqns_bvp/issues/4#issuecomment-931410206, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALTXJ3OK6OLQNHZLWSWCU7LUER4RJANCNFSM5FCOSIXQ .

[aarobel]

Yes, I think SimpleSSA would be a good starting point as long as you don’t expect surging or things like that (in which case I would suggest the thermal flowline code which is built to accommodate large, rapid changes in velocity). The only thing you should need to change is the basal shear stress term in the momentum balance to something including N. Keep in mind that all model equations are scaled to the terminus position and are also non dimensionalized.

[glugeorge]

Hi! I've just gotten around to playing around with this code and have made a few changes here and there. Just for a quick sanity check, I'm taking that this sigma variable can be interpreted as like a non-dimensionalized distance between the origin (sigma=0) and the grounding line position (sigma=1)? I've attached a plot that I ran with the code and existing parameters here and I want to make sure I'm interpreting it correctly (here, I ran it in transient mode, and also changed the plot for thickness to a surface plot rather than contour):

[aarobel]

Hmm, something definitely seems off here. Looks like the glacier is growing and then running into some numerical issues just before year 3000. What exactly did you change?

I'm not sure if I've changed anything since I last updated the repo, but when I run the basic transient simulation which I have locally (initializing at steady-state with accumulation of 1 m/yr and then doing transient response to accumulation of 0.8 m/yr), the attached plot is what I get.

[glugeorge]

Oh ok, yeah that's what I got when I first ran it untouched - the relevant change I made that caused the blowup was just setting params.transient = 1 instead of 0 on line 33. So does that just initialize in a transient state rather than steady state? Thanks!

[aarobel]

Aha, yes, that explains that. The way that the code is currently setup is that I give it a pretty odd initial guess (just something thats convenient), and then I have it find a steady-state configuration that initial guess (lines 50-76). Then lines 79-100 involve running a transient simulation from that initial steady-state. So, if you are OK with that basic setup, I would change line 33 back to transient=0.

[jkingslake]

Thanks for the input @aarobel.

@glugeorge I would suggest trying to couple with the hydrology model without the lake initially, given that the lake causes complications when it over fills or empties. You could just ignore the lake equation and fix Q at s = 0 instead of fixing N.

Also, I would start with the steady state of the ice geometry because we will need to think carefully about how to scale the channel equations to take account of the moving grounding line.

Also, I guess you will be converting this to python?

[aarobel]

I would be very interested if you are planning to write a python version. I wrote a Julia version (also on my GitHub) that uses automatic differentiation for the solver (so it goes faster and is more accurate), but haven’t gotten around to writing a python version.

See the appendix of Schoof 2007 for more info on how to scale new prognostic equations. There are some typos in there, which I can always advise on.

[jkingslake]

OK thanks fo the point to Schoof 2007 for that. I will take a look.

Yes I saw that you had a Julia version. I wish I had the time to dedicate to just skipping python and going straight to learning Julia! :-)

[glugeorge]

Yeah I can give a shot at converting it to python! I'll first getting a version of the code without the lake and then make the conversion. I'll also look at the julia code since I feel it'd be useful to learn as well

[glugeorge]

Revisiting this, I've made some adjustments to the hydrology model to work with the sigma and tau coordinates in my latest commit (https://github.com/ldeo-glaciology/nye_eqns_bvp/blob/main/nf_bvp_python/nf_bvp.ipynb) , and for it to be in a steady state (without the lake). Next step is going through the ice sheet model and making the linkages

The existing code makes some assumptions about the glacier velocity and psi term as a function of sigma and tau, but I think connecting the two models will resolve it

[jkingslake]

Great @glugeorge !

It might be good to document the rescaling of the model equations on here somewhere to help me understand what you've done and considering you'll have to write it up one day anyway.

[glugeorge]

Sure thing, it's all documented on the overleaf document that I've been working on: https://www.overleaf.com/read/hzvvwvjkfwfj

But yeah these are the changed equations based off the rescaling from equations A1-A3 from Schoof (2007):

[jkingslake]

These all looks good to me. As you know, I have recently been going through exactly the same change of coordinates (with @rmskarbek 's help) for a firn compaction model (code here).

How have you scaled x_g? With the ice velocity scale?

Looking forward to seeing some of the results with a moving GL. I would start by running the hydro model to steady state with and then prescribing a constant GL migration and seeing how it effects N.

[glugeorge]

Yeah so x_g evolves non-dimensionally with the ice velocity over time, so I think it should be scaled with those. Here are the results with a stationary grounding line (no glacier velocity):

And then with a constant migration:

The only effect on N I see is that it kind of moves with the grounding line.

[aarobel]

A question about the above posts: in the original Schoof 2007 model, the grounding line position is an implicit function of bed topography and ice thickness at the grounding line (maintaining floatation at this grid point). The way the model equations are non-dimensionalized, there is a horizontal length scale [x], which is called xscale in the original code which scales the grounding line position. However, I'm not sure if you changed some of the scaling assumptions in your new coupled version of the model.

[glugeorge]

Good point, I think some of the scalings between SSA_simple and my current model vary a bit. I haven't coupled them yet, but I'll try to use consistent scalings when I go ahead to couple the two! Right now Im inclined to believe my existing scaling works since the prescribed constant evolution of the grounding line works out to 1km/y (which is the prescribed velocity) when I re-dimensionalize all the values

[aarobel]

I suppose it depends on what you would like the boundary conditions at the grounding line/terminus to be. If you are saying that there is no ice leaving the terminus, then it makes sense that dx_g/dt = u_g where x_g is the terminus position and u_g is the ice velocity at the terminus. An implicit boundary condition, like the floatation condition in SSA_simple, includes a grounding line ice flux removing ice from the domain in order to maintain floatation.

[jkingslake]

Yeah, I agree on these points Alex. But I think having the GL move at the ice velocity was just a simple assumption as a first try to check the hydrology model code. In general in the future George I'm guessing that you will be prescribing x_g from the ice flow model, in which case of course $\dot{x_g} \neq u$ in general.

On Mon, Jan 17, 2022, 11:26 AM Alexander Robel @.***> wrote:

I suppose it depends on what you would like the boundary conditions at the grounding line/terminus to be. If you are saying that there is no ice leaving the terminus, then it makes sense that dx_g/dt = u_g where x_g is the terminus position and u_g is the ice velocity at the terminus. An implicit boundary condition, like the floatation condition in SSA_simple, includes a grounding line ice flux removing ice from the domain in order to maintain floatation.

— Reply to this email directly, view it on GitHub https://github.com/ldeo-glaciology/nye_eqns_bvp/issues/4#issuecomment-1014712534, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALTXJ3JZYXMMGVL3Z5MIJVDUWQ7K3ANCNFSM5FCOSIXQ . You are receiving this because you were mentioned.Message ID: @.***>

[glugeorge]

Mostly for myself to keep track of progress. I've been trying to convert SSAsimple into python for easier communication between that script and mine, but am having difficulty with the solver. The matlab solver uses a Trust-Region-Dogleg method, which doesn't seem to have a straightforward implementation in python (I've tried all the scipy.optimize solvers and none seem to converge). The plan is to look into a python package that can calculate the jacobian and hessian and use that in scipy.optimize.minimize with the dogleg option, though I'm a bit skeptical of that working. I'm also going to look at the julia script because I understand it uses a different method of solving as well. If anyone has additional suggestions, would love to hear them!

[aarobel]

Yeah, I think you might be running into the same issue I was with getting Python solvers to converge. The Julia version uses an auto-differentiation package to calculate the Jacobian, and then a nonlinear solver package to calculated the solution. My informal testing indicates that the accuracy between the Julia and MATLAB versions is comparable, and after the first solve, the Julia solver runs about 4x faster (because it doesn't have to approximate the Jacobian). I think similar packages exist for Python, so you could try that.

[glugeorge]

Coupling process has been completed (to an extent)- new scripts and experiments at: https://github.com/glugeorge/coupled_ice_hydrology