To do ideas 1: replace u in u.dS/dx with a sliding law

[jkingslake]

Replace the ice flow speed u in the $u.dS/dx term in the channel evolution equation with a sliding law.

First step would be to use something very simple like u = C tau_b / N, where C is a constant and tau_b is the basal shear stress and N is the effective pressure. This will require you to use a non-zero boundary condition on N at the terminus (otherwise u will --> inf), and to define tau_b. Start by assuming tau_b is uniform and equal to the yield stress of ice, 100 kPa.

A different approach that would allow you to keep N = 0 and terminus wold be to use a sliding law that does not give u = inf when N = 0, e.g., Schoof, 2005.

[glugeorge]

I opted to use a basic sliding law as suggested in the form u = C tau_b / N. To get a rough sense of how this would work, I used tau_b = 100kPa, and then looked at my previous results for N (approx 600kPa) and an assumed u of 100 m/year to get an estimate for C = 6 m/year. My non-dimensionalization strategy is attached below:

Next I just used a very small N_t as the terminus boundary condition 0.0001 - I noticed if I used rho_ig*H (H=ice thickness, which isnt referred to in the rest of the equations if my memory serves me correct), the lake just drains after 1 cycle (also interestingly enough S at the terminus shrunk back too), and I wanted to keep the other parameters the same.

First, I wanted to check to see if the nonzero boundary condition would cause S at the terminus to decrease without a glacier velocity adjustment. The output, shown below, looks more or less the same as the simulation with Nt=0 and no velocity adjustment.

Afterwards, with the velocity adjustment described above, we get: Here we can see that the terminus S eventually reaches a stable cycle as well.

To do: I want to output the velocities to make sure things make sense, and I would love to take suggestions on determining C in a more rigorous manner + steps to looking for a non-uniform tau

[jkingslake]

Looks great! Just a clarification, the top set of plots is without the new term, but with N_t = 0.0001, and the bottom one is with the new term and N_t = 0.0001, correct?

If so, then this looks really promising, the only issue is that when you plot out the sliding speed as a function of x (i.e. u(x)) you will find that u get unrealistically large at the terminus, where N_t* is so tiny. The next step will be to use a sliding law that doesn't have this unrealistic blow-up in u at small N.

But for now it would be interesting is to see all four variables as functions of x and t, as color maps, i.e.

imagesc(x,t,S)
imagesc(x,t,N)

etc.

[jkingslake]

also, have you updated the repo with this new version of the code?

[jkingslake]

and @agstub can weigh in on getting realistic values for C. He recently used values for a sliding parameter from an inversion by Arthern et al., 2015. Although be careful that the definitions of the parameters align exactly with the ones you want (i.e. the sliding law is actually the same as the one you are using).

For now I would simply use a constant value from the literature. For example, in this 2013 paper a adopted a value from Hewitt and Fowler, 2008, although again, my sliding law is slightly different than the one you've used so far (there is an exponent on the shear stress (my eqn 10 in the 2013 paper).

[glugeorge]

Just a clarification, the top set of plots is without the new term, but with N_t = 0.0001, and the bottom one is with the new term and N_t = 0.0001, correct?

Yes, that is correct. I will try to add the color maps with my next commit.

also, have you updated the repo with this new version of the code?

I just pushed it, will also try adding C values from literature in my next commit.

[glugeorge]

After fixing some bugs and using c values from literature and a more realistic boundary condition, here are some plots:

Evidently, even after converting to logarithmic colorbar for velocity, it's still not good. If I filter values greater than 1km/y and less than -100m/y, the figure ends up looking like:

I find this one quite interesting we can see that the glacier at the lake outlet speeds up and slows down/backs up with each flood pulse, but then the terminus is just consistently fast (presumably due to boundary conditions). I'm going to move on to other to-do ideas, but if I were to revisit this one, my next steps would be incorporating a non-constant tau and using a sliding law that allows for the N_t = 0 boundary condition - though this may be eventually resolved when I couple it with an ice flow model

[jkingslake]

Ok, looks good. One unrealistic thing that is happening is the lake level is getting much larger than the flotation depth (~90m?) and therefore N goes negative. This causes the blow up in u during lake filling and the region of negative N and negative u at the lake-end of the domain while the lake is at it greatest depth.

The model simply doesnt apply when N = 0, for the reasons we discussed. I would simply adjust psi to avoid this happening, as it is a regime that our model cant be applied to.

Good idea to move on to the other ideas though.

Also, if. you have updated the repo after having made these new plots its a good idea to link to the commit and mention the scripts which produce the figures, in case others want to reproduce them.

[glugeorge]

Added the adjustment to avoid velocity blowing up at lake terminus. Closing issue for now to pivot to issue #2

[jkingslake]

Great, good work @glugeorge! Can you stick in a figure, so we can show @agstub ?

[glugeorge]

For @agstub to look at: Plots of lake depth, channel flow, and channel cross sectional area (I realized theres no legend but the orange lines are for channel exit, whereas the blue ones are the lake exit. Some colormaps of effective pressure, channel flow, channel cross sectional area, and glacier basal velocity:

[agstub]

This looks awesome!