- #1

- 11

- 1

## Main Question or Discussion Point

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.

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.

Python:

```
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)
```

Last edited: