Numerical modeling of a glacier's length -- Coding error(s)

[Mikkel]

Hey Physics Forum

I am currently doing my bachelor project in geophysics, with focus on the evolution of glaciers in Greenland. My project consists partly of programming, because I want to get better at it. I have, however, hit a wall. I can't seem to figure out what is wrong with my code and I have been trying different things for days. I would love any input/feedback.
The basic problem is describing how a glacier evolve over time, which is done solving the differential equation using Eulers Method.
I've tried to make a simple for loop, but I get negative values for the length of the glacier (L), which can not happen, so clearly something is wrong.

import numpy as np
import matplotlib.pyplot as plt
from math import exp # exponential function
from math import *

n = 5000    #Iterations (5000 years)
t = np.ones(n)
t[0] = 0    #Time(0)
dt = 1      #Time step

    ## Parameters ##
d0 = 200     #d0 (Start height)
s = 0.014   #Slope
lambd = 300     #lambda
sigma = 10000  
xs = 40000  #x-coordinate of the "bump"
ac = 0.0005 #Accumulation
alfaf = 0.7 #Constant
alfam = 2 #Constant
eps = 1     #Constant
delta = 1.127 #Constant
c = 2.4 #Climate sensitivity

    ## Grids ##
x = np.linspace(0,50000,5000) #My x grid
x[0] = 200              #Start position of glacier's head
L = np.ones(n)          #Glacier Length
B = np.zeros(n)
a = np.ones(n)*ac
Hf = np.ones(n)         #Defining vectors to loop through
F = np.ones(n)
dLdt = np.zeros(n)

for i in range(1, n):
    t[i] = i
    a[i] = a[i-1]+1*ac  #Accumulation (constant)
    x[i] = d0-s*x[i]+lambd*np.exp(-((x[i]-xs)/(sigma))**2) #Bed geometry
    B[i] = a[i-1]*L[i-1]    #Surface mass balance
    Hf[i] = np.max([alfaf*np.sqrt(L[i-1]), -eps*delta*x[i-1]]) #Ice thickness at glacier front
    F[i] = np.min([0, c*x[i-1]*Hf[i-1]])    #Flux of ice at the glacier front
    dLdt[i] = (2*(B[i-1]+F[i-1]))/(3*alfam*np.sqrt(L[i-1])) #Differential equation
    L[i] = L[i-1] + 1*dLdt[i-1]     #Updating L with Eulers Method

plt.plot(t,L)

 

[willem2]

Mod note: Fixed the array indexes. When you use i as an array index, in brackets, the browser interprets this as a BBCode for italics. Adding a space inside the brackets in front of 'i' prevents this from happening.

I don't really know what the computation for x[ i] "bed geometry" is. It does not depend on the other time-dependent variables. x[ i] starts out at 200 but goes to -389 in the end. This will make F[ i] go so negative that dLdt[ i] becomes negative and this will eventually make L negative around 2230, and you get square roots of negative numbers. L also becomes very big before it goes down again and becomes negative.

If the "bed geometry" is something like altitude or roughness of the surface under the glacier, it should probably depend on position, not on time and you should likely need to compute the solution of a partial differential equation, tracking ice thickness, accumulation and ice movement along the entire glacier.
I rather doubt that any 0-dimensional model of a glacier is of any use.

Accumulation ac is a bit strange also, this will simply go up linearly and not be constant.

 

[QuantumQuest]

Hey Physics Forum

I am currently doing my bachelor project in geophysics, with focus on the evolution of glaciers in Greenland. My project consists partly of programming, because I want to get better at it. I have, however, hit a wall. I can't seem to figure out what is wrong with my code and I have been trying different things for days. I would love any input/feedback.
The basic problem is describing how a glacier evolve over time, which is done solving the differential equation using Eulers Method.
I've tried to make a simple for loop, but I get negative values for the length of the glacier (L), which can not happen, so clearly something is wrong.

import numpy as np
import matplotlib.pyplot as plt
from math import exp # exponential function
from math import *

n = 5000    #Iterations (5000 years)
t = np.ones(n)
t[0] = 0    #Time(0)
dt = 1      #Time step

    ## Parameters ##
d0 = 200     #d0 (Start height)
s = 0.014   #Slope
lambd = 300     #lambda
sigma = 10000
xs = 40000  #x-coordinate of the "bump"
ac = 0.0005 #Accumulation
alfaf = 0.7 #Constant
alfam = 2 #Constant
eps = 1     #Constant
delta = 1.127 #Constant
c = 2.4 #Climate sensitivity

    ## Grids ##
x = np.linspace(0,50000,5000) #My x grid
x[0] = 200              #Start position of glacier's head
L = np.ones(n)          #Glacier Length
B = np.zeros(n)
a = np.ones(n)*ac
Hf = np.ones(n)         #Defining vectors to loop through
F = np.ones(n)
dLdt = np.zeros(n)

for i in range(1, n):
    t[i] = i
    a[i] = a[i-1]+1*ac  #Accumulation (constant)
    x[i] = d0-s*x[i]+lambd*np.exp(-((x[i]-xs)/(sigma))**2) #Bed geometry
    B[i] = a[i-1]*L[i-1]    #Surface mass balance
    Hf[i] = np.max([alfaf*np.sqrt(L[i-1]), -eps*delta*x[i-1]]) #Ice thickness at glacier front
    F[i] = np.min([0, c*x[i-1]*Hf[i-1]])    #Flux of ice at the glacier front
    dLdt[i] = (2*(B[i-1]+F[i-1]))/(3*alfam*np.sqrt(L[i-1])) #Differential equation
    L[i] = L[i-1] + 1*dLdt[i-1]     #Updating L with Eulers Method

plt.plot(t,L)

I would advise to first make sure that the steps you take together with the math involved, produce meaningful results from a physical perspective i.e. apply the model in an algorithmic approach using pen and paper. Then, you can go on and code it in Python, using the appropriate libraries and constructs. In other words, when the results we take from the execution of a program are wrong, there are one or more logical errors which are responsible for this.

 

[Mikkel]

Thank you both for your responses- really appreciate you taking your time!
You're probably right in, starting from scratch, and doing in by pen and paper, is a good idea. I will give this another try.

 

[Mark44]

And if you start from scratch, do a favor for yourself and anyone else who might have to read your code, by coming up with better, more descriptive variable names.

d0 #d0 (Start height) -- Why did you choose d0 to represent height?
lambd = 300 -- Why not go all the way and name this one lambda?
sigma = 10000 -- What's the purpose of this variable?
xs = 40000 #x-coordinate of the "bump" -- Why xs? x_bump or xBump would be more explanatory.
ac = 0.0005 #Accumulation -- If you called it, say accum or similar, it would be self-documenting.
alfaf -- What does this represent? It looks to be short for alfalfa, but there's probably none of that in a glacier.
alfam -- Ditto
eps -- Ditto
delta = 1.127 -- What is this significance of this constant?
c = 2.4 #Climate sensitivity -- You really should have a better name than c.
L -- Not too bad for length, but single-letter variables should be avoided except when used for loop counters such is i, j, and so on.
B = np.zeros(n) -- What does this array signify? There's not even a comment for it.
a -- Ditto
Hf = np.ones(n) #Defining vectors to loop through -- I have no idea what this comment means.
F -- What does this array signify? There's not even a comment for it.
dLdt -- This should have a comment. It's apparently the derivative of glacier length with respect to time.

You might think it takes too much time to create useful variable names, but having meaningful variable names goes a long way toward minimizing the number of bugs in a program. If you look back on the code in three months, it will be much easier to understand what the code is doing.

 

[pbuk]

You might think it takes too much time to create useful variable names, but having meaningful variable names goes a long way toward minimizing the number of bugs in a program. If you look back on the code in three months, it will be much easier to understand what the code is doing.

This is true in general, but in solving ODE Initial Value Problems where time is usually the independent variable there are a couple of conventions:

# Time
t = 0
# Step size (we could use delta_t or dt but h is a good one because it is what we usually
# use in the maths behind the code e.g. proofs of convergence, error bounds)
h = 0.1

There is quite a lot more wrong here than variable naming though I am afraid. Why don't you forget about Euler's method and use a Python library like scipy's solve_ivp? If you do this you will realize that solving an ODE IVP is all about defining a function that calculates first derivatives (which scipy's documentation unhelpfully calls fun; I prefer y_prime or dy_dt). This will also help you rethink how you can model a glacier as an IVP (initial value problem).

 

[Dr Transport]

python is forgiving, you can name/define the variables as you go along. That's what I have tended to do since I started playing around with it over the past year or so. When I am ready to pass the code/script along I clean up and drag all the definitions and the corresponding documentation up to the top for ease of reading. As I find I do not need a variable any longer, I'll do a delete block get rid of the overhead. As has been said previously, meaningful names are necessary for ease of debugging.