<jupyter_start><jupyter_text>INTRODUCTION TO NUMERICAL ICE-SHEET MODELLING ***&copy; 2018, 2019 [<NAME>](mailto:<EMAIL>), Potsdam Institute for Climate Impact Research******- examples for simple diffusion and advection models (FTCS) as well as a simple ice-sheet model based on the Shallow Ice Approximation (SIA)- 2d ("flowline"): $x$ denotes horizontal extend of the ice sheet, $z$ vertical ice thickness***  Forward in Time, Central in Space (FTCS) schemes - based on finite differences (i.e. divide model domain in regularly-spaced nodes on which variables are calculated)- gradients/derivatives approximated by Taylor series expansion (omitt higher order derivatives)$$H(t+\Delta t)\approx H(t)+\Delta t\frac{\partial H}{\partial t}+\frac{(\Delta t)^2}{2}\frac{\partial^2 H}{\partial t^2}+\dots$$**Discretized Taylor approximations**- *temporal* derivatives forward in time:$$\frac{\partial H}{\partial t}\approx\frac{H_{t+\Delta t} - H_t}{\Delta t}$$- *spatial* derivatives centered in space:$$\frac{\partial H}{\partial x}\approx\frac{H_{i+1} - H_{i-1}}{2\Delta x}$$$$\frac{\partial^2 H}{\partial x^2}\approx\frac{H_{i+1} - 2H_i + H_{i-1}}{(\Delta x)^2}$$ The second derivative is determined by calculating local spatial gradients between consecutive grid points and taking gradient across these points, i.e., it is a symmetric derivative around $i$, involving only immediate neighbouring grid points $i+1$ and $i-1$<jupyter_code># load Python libraries and set up notebook

%matplotlib notebook

import numpy as np
import matplotlib.pyplot as plt
import time

plt.rcParams['figure.figsize']=10,4
plt.rcParams['lines.linewidth']=3
plt.rcParams['font.size']=12
plt.rcParams['axes.titleweight']='bold'
plt.close('all')<jupyter_output><empty_output><jupyter_text>Example 1: Diffusion - slow point-to-point transfer of material (steady creep)$$\frac{\partial H}{\partial t}=D \frac{\partial^2 H}{\partial x^2}$$$D$: diffusion coefficient (in the code we replace $D$ with $D\Delta t$)**FTCS scheme:**$$H_{i,t+1} = H_{i,t} + \frac{D\Delta t}{(\Delta x)^2} (H_{i+1,t} - 2H_{i,t} + H_{i-1,t})$$<jupyter_code>def diffusion(H,dx,m,D):
    
    Hn=np.array(H, copy=True) # copy existing H into new vector Hn
    
    for i in np.arange(1,m-1,1):
        Hn[i]=H[i]+D*(H[i+1]-2*H[i]+H[i-1])/dx**2
    H=np.array(Hn, copy=True) # copy new Hn back into old vector H
    
    return H<jupyter_output><empty_output><jupyter_text>Example 2: Advection - rapid transport of material with the motion of a fluid- magnitude of transport velocity determined by surface gradient of hillslope$$\frac{\partial H}{\partial t}=u \frac{\partial H}{\partial x}$$$u$: displacement speed**FTCS scheme:**$$H_{i,t+1} = H_{i,t} + \frac{u\Delta t}{(2\Delta x)} (H_{i+1,t} - H_{i-1,t})$$ Note that advection processes are inherently unstable. Solution: add a small amount of artificial diffusion in the code in order to stabilize the numerical integration ($D\ll u$).<jupyter_code>def advection(x,H,dx,m,D,u):
    
    Hn=np.array(H, copy=True) # copy existing H into new vector Hn

    for i in np.arange(1,m-1,1):
        if x[i]>=0:
            Hn[i]=H[i]+u*(H[i+1]-H[i-1])/(2*dx) + D*(H[i+1]-2*H[i]+H[i-1])/dx**2
        else:
            Hn[i]=H[i]-u*(H[i+1]-H[i-1])/(2*dx) + D*(H[i+1]-2*H[i]+H[i-1])/dx**2
    H=np.array(Hn, copy=True) # copy new Hn back into old vector H
    
    return H<jupyter_output><empty_output><jupyter_text>**Initial conditions**Slope at $t=0$:$$H_{i,t=0} = 1 - \left|x_i\right|$$**Boundary conditions***Dirichlet* boundary conditions (i.e. edge grid points remain constant over time):$$H_1=H_m=0$$<jupyter_code>### choose method
#method = 'diffusion'
method = 'advection'
%%time

fig,ax = plt.subplots(1,1)

plt.ion()
fig.show()
fig.canvas.draw()


m=101 # number of grid points
x=np.linspace(-1,1,m) # horizontal distance
dx=2/(m-1) # grid spacing
time_end=500. # duration (years)
H=1-np.abs(x) # initial slope
H0=np.array(H, copy=True)


# time stepping
for t in np.arange(0,time_end+1,1):
    
    ax.clear()
    
    # plot initial slope
    ax.plot(x,H0, c='black')
    
    ax.set_xlim([-1,1])
    ax.set_ylim([0,1])
    ax.set_xlabel('Scaled horizontal distance')
    ax.set_ylabel('Scaled elevation')
    ax.set_title(method)

    if method=='diffusion':
        H=diffusion(H,dx,m,D=2e-4)
    elif method=='advection':
        H=advection(x,H,dx,m,D=4e-6,u=4e-4) #D=1e-4,u=4e-3 #D=4e-6,u=4e-4
        # try D=0 :)
    
    # plot new thickness
    ax.plot(x,H, c='orangered')
    ax.text(0.93,0.95, 't='+str(t), va='top', ha='right')
    
    fig.canvas.draw()
    time.sleep(0.01)<jupyter_output><empty_output><jupyter_text>2d ice-sheet equation Vertical mean horizontal velocity (neglecting basal sliding) in flowline according to SIA:$$\bar{u}=-\frac{2A}{n+2}(\rho_i g)^n \left|\frac{\partial s}{\partial x}\right|^{n-1} \frac{\partial s}{\partial x} h^{n+1}$$$s=b+h$: surface elevation$A$: ice flow parameter (softness)$n$: Glen exponent$\rho_i$: ice density$g$: gravitation constantThickness evolution:$$\frac{\partial h}{\partial t}=-\frac{\partial (\bar{u}h)}{\partial x}+\dot{a}=-\frac{\partial q}{\partial x}+\dot{a}$$$\dot{a}$: accumulation rate$q$: flux<jupyter_code># Model setup
m=51 # number of grid points
L=1.5e6 # length of domain (m)
x=np.linspace(0,L,m) # horizontal distance
dx=L/(m-1) # grid spacing (m)
d=np.zeros(m) # diffusion coefficients
time_end=20000. # duration (years)

# Physical constants
A=1e-16 # Ice flow parameter
rho=910. # Ice density
grav=9.81 # Gravitation constant
n=3 # Glen index
a=0.3 # Surface mass balance<jupyter_output><empty_output><jupyter_text>Analytical solution (Vialov, 1958) Assumptions:- flat, rigid bed: $b=0$, i.e. $s=h$- no basal sliding- constant surface mass balance over whole  domain $\dot{a}$*Vialov profile* for steady state:$$h(x)=h_0\left[1-\left(\frac{x}{L}\right)^{(n+1)/n}\right]^{n/(2n+2)}$$with thickness at ice divide$$h_0=2^{n/(2n+2)}\left(\frac{\dot{a}}{A_0}\right)^{1/(2n+2)} \sqrt{L}$$and$$A_0=\frac{2A}{n+2}(\rho_i g)^n$$<jupyter_code>def vialov(L,m):
    
    Lv=L/2. # length of half-domain (m)
    xv=np.linspace(-Lv,Lv,m) # horizontal distance from ice divide (m)
    hv=np.zeros(m)
    
    A0=2./(n+2.)*A*(rho*grav)**n
    h0=2**(n/(2.*n+2.))*(a/A0)**(1./(2.*n+2))*Lv**0.5 # surface elevation at ice divide
    
    for j in np.arange(2,m-2,1):
        if xv[j]<=0:
            hv[j]=h0*(1-(np.abs(xv[j])/Lv)**((n+1.)/n))**(n/(2.*n+2))
        elif xv[j]>0:
            hv[j]=h0*(1-(xv[j]/Lv)**((n+1.)/n))**(n/(2.*n+2))
    
    return hv


# create plot
if True:
    plt.figure()
    
    hv=vialov(L,m)
    plt.plot(x,hv, c='steelblue')
        
    plt.xlim([0,L])
    plt.ylim([0,4000])
    plt.xlabel('Horizontal distance (m)')
    plt.ylabel('Elevation (m)')
    plt.title('Vialov profile')
    
    plt.show()<jupyter_output><empty_output><jupyter_text>Example 1: Explicit solution     diffusion problem    requires small time steps to remain stable    $\sim$ 20 min. for 20,000 years with $\Delta t$=2 years (on 2.8 GHz Intel Core i7 processor) **FTCS scheme:**$$h_{i,t+1}=h_{i,t}+\left[q_{i+\frac{1}{2},t}-q_{i-\frac{1}{2},t}\right]\frac{\Delta t}{\Delta x}+\dot{a}\Delta t$$with fluxes between grid points$$q_{i+\frac{1}{2},t}=\frac{2A}{n+2}(\rho_i g)^n \left(\frac{h_{i+1,t}+b_{i+1,t}-h_{i,t}-b_{i,t}}{\Delta x}\right)^n \left(\frac{h_{i,t}+h_{i+1,t}}{2}\right)^{n+2}$$$$q_{i-\frac{1}{2},t}=\frac{2A}{n+2}(\rho_i g)^n \left(\frac{h_{i,t}+b_{i,t}-h_{i-1,t}-b_{i-1,t}}{\Delta x}\right)^n \left(\frac{h_{i,t}+h_{i-1,t}}{2}\right)^{n+2}$$Pure advection equations can be numerically unstable in FTCS. In reality, the continuity equation is a diffusive equation due to the dependence of $q$ on $\partial h/\partial x$:$$\frac{\partial h}{\partial t}=-\frac{\partial (\bar{u}h)}{\partial x}+\dot{a}=\frac{\partial}{\partial x}\left(d\frac{\partial (h+b)}{\partial x}\right)+\dot{a}$$with diffusion coefficient$$d=-\bar{u}h\left(\frac{\partial (h+b)}{\partial x}\right)^{-1}=-\frac{2A}{n+2}(\rho_i g)^n \left|\frac{\partial (h+b)}{\partial x}\right|^{n-2} \frac{\partial (h+b)}{\partial x} h^{n+2}$$Diffusion terms (i.e. transformed flux terms) are evaluated between grid points ("staggered grid"):$$d_{i+\frac{1}{2},t}=\frac{2A}{n+2}(\rho_i g)^n \left|\frac{h_{i+1,t}+b_{i+1,t}-h_{i,t}-b_{i,t}}{\Delta x}\right|^{n-1} \left(\frac{h_{i,t}+h_{i+1,t}}{2}\right)^{n+2}$$$$d_{i-\frac{1}{2},t}=\frac{2A}{n+2}(\rho_i g)^n \left|\frac{h_{i,t}+b_{i,t}-h_{i-1,t}-b_{i-1,t}}{\Delta x}\right|^{n-1} \left(\frac{h_{i,t}+h_{i-1,t}}{2}\right)^{n+2}$$$$h_{i,t+1}=h_{i,t}+ \frac{\Delta t}{(\Delta x)^2} \left[d_{i+\frac{1}{2},t}(h_{i+1,t}+b_{i+1,t}-h_{i,t}-b_{i,t})-d_{i-\frac{1}{2},t}(h_{i,t}+b_{i,t}-h_{i-1,t}-b_{i-1,t})\right]+\dot{a}\Delta t$$ Figure: Staggered grid used in the flowline model.<jupyter_code>def explicit(b,h,dx,dt):
        
    slope=np.diff(b+h)/dx # slope on staggered grid
    
    for j in np.arange(1,m-1,1):
        hstag=(h[j+1]+h[j])/2. # ice thickness on staggered grid
        d[j]=2./(n+2.)*A*hstag**2*slope[j]**2*(rho*grav*hstag)**n
    
    h2=np.zeros(m)
    for j in np.arange(2,m-2,1):
        h2[j]=h[j]+(d[j]*slope[j]-d[j-1]*slope[j-1])*dt/dx+a*dt
    h=np.array(h2, copy=True)
    
    return h<jupyter_output><empty_output><jupyter_text>Example 2: Semi-implicit solution     greater numerical stability by evaluating right-hand side on timestep $t+1$ instead of $t$ ("implicit")    $\sim$ 10 min. for 20,000 years with $\Delta t$=2 years (on 2.8 GHz Intel Core i7 processor) **FTCS scheme:**$$h_{i,t+1}=h_{i,t}+\frac{\Delta t}{(\Delta x)^2}\left[d_{i+\frac{1}{2},t}(h_{i+1,t+1}+b_{i+1,t}-h_{i,t+1}-b_{i,t})-d_{i-\frac{1}{2},t}(h_{i,t+1}+b_{i,t}-h_{i-1,t+1}-b_{i-1,t})\right]+\dot{a}\Delta t$$Sort terms:$$h_{i,t+1}-\frac{\Delta t}{(\Delta x)^2}\left[d_{i+\frac{1}{2},t}(h_{i+1,t+1}-h_{i,t+1})-d_{i-\frac{1}{2},t}(h_{i,t+1}-h_{i-1,t+1})\right]=h_{i,t}+\frac{\Delta t}{(\Delta x)^2}\left[d_{i+\frac{1}{2},t}(b_{i+1,t}-b_{i,t})-d_{i-\frac{1}{2},t}(b_{i,t}-b_{i-1,t})\right]+\dot{a}\Delta t$$Define coefficients:$$\alpha_i=\frac{\Delta t}{(\Delta x)^2}d_{i-\frac{1}{2},t}$$$$\beta_i=1+\frac{\Delta t}{(\Delta x)^2}\left(d_{i-\frac{1}{2},t}+d_{i+\frac{1}{2},t}\right)$$$$\gamma_i=\frac{\Delta t}{(\Delta x)^2}d_{i+\frac{1}{2},t}$$$$\delta_i=\alpha_ib_{i-1,t}-(\beta_i-1)b_{i,t}+\gamma_ib_{i+1,t}+h_{i,t}+\dot{a}\Delta t$$Such that:$$-\alpha_ih_{i-1,t+1}+\beta_ih_{i,t+1}-\gamma_ih_{i+1,t+1}=\delta_i$$This system of equations can be written in matrix notation giving a tridiagonal matrix (all non-zero elements occur on the diagonal axis as well as the two adjacent axes). Such a system can easily be solved with the method below. Consider the following two equations that hold between two adjacent grid points:$$h_{i}=f_{i}h_{i+1}+g_{i}$$$$h_{i-1}=f_{i-1}h_{i}+g_{i-1}$$Multiply second equation by $\alpha_i$:$$\alpha_ih_{i-1}=\alpha_if_{i-1}h_{i}+\alpha_ig_{i-1}=\beta_ih_{i}-\gamma_ih_{i+1}-\delta_i$$Combine equations:$$h_i=\frac{\gamma_i}{\beta_i-\alpha_if_{i-1}}h_{i-1}+\frac{\delta_i+\alpha_ig_{i-1}}{\beta_i-\alpha_if_{i-1}}$$where$$f_i=\frac{\gamma_i}{\beta_i-\alpha_if_{i-1}}, \qquad g_i=\frac{\delta_i+\alpha_ig_{i-1}}{\beta_i-\alpha_if_{i-1}}$$<jupyter_code>def semi_implicit(b,h,dx,dt,dtdx,alpha,beta,gamma,delta,f,g):
    
    slope=np.diff(b+h)/dx # slope on staggered grid
    
    for j in np.arange(1,m-1,1):
        hstag=(h[j+1]+h[j])/2. # ice thickness on staggered grid
        d[j]=2./(n+2.)*A*hstag**2*slope[j]**2*(rho*grav*hstag)**n
    
    for j in np.arange(2,m-2,1):
        gamma[j]=d[j]*dtdx
        alpha[j]=d[j-1]*dtdx
        beta[j]=1.+gamma[j]+alpha[j]
        delta[j]=h[j]+a*dt+alpha[j]*b[j-1]-(beta[j]-1)*b[j]+gamma[j]*b[j+1]
    
    for j in np.arange(2,m-2,1):
        f[j]=gamma[j]/(beta[j]-alpha[j]*f[j-1])
        g[j]=(delta[j]+alpha[j]*g[j-1])/(beta[j]-alpha[j]*f[j-1])
    
    h2=np.zeros(m)
    for j in np.arange(m-2,0,-1):
        h2[j]=g[j]+f[j]*h2[j+1]
    h=np.array(h2, copy=True)
    
    return h<jupyter_output><empty_output><jupyter_text>**Initial conditions**No ice at $t=0$:$$H_{i,t=0} = 0$$**Boundary conditions***Dirichlet* boundary conditions (i.e. edge grid points remain constant over time):$$H_1=H_m=0$$Additionally, we have for the semi-implicit method:$$f_1=f_m=g_1=g_m=0$$  Case 1: Flat bed<jupyter_code>### choose method
#method = 'explicit'
method = 'semi-implicit'
%%time

fig,ax = plt.subplots(1,1)

plt.ion()
fig.show()
fig.canvas.draw()


# Initialize model
h=np.zeros(m) # initial thickness
b=np.zeros(m) # bedrock elevation

if method=='explicit':
    dt=2. # time step (years)
    
elif method=='semi-implicit':
    dt=10. # time step (years)
    dtdx=dt/dx**2
    
    # initialize vectors
    alpha=np.zeros(m)
    beta=np.zeros(m)
    gamma=np.zeros(m)
    delta=np.zeros(m)
    f=np.zeros(m)
    g=np.zeros(m)


# time stepping
for t in np.arange(0,time_end+1,dt):
    
    ax.clear()
    
    if method=='explicit':
        h=explicit(b,h,dx,dt)
    elif method=='semi-implicit':
        h=semi_implicit(b,h,dx,dt,dtdx,alpha,beta,gamma,delta,f,g)
        
    # plot bedrock
    ax.plot(x,b, c='black')
    
    # plot Vialov profile
    hv=vialov(L,m)
    ax.plot(x,hv, c='steelblue', label='analytical')
    
    # plot thickness profile
    ax.plot(x,h, c='orangered', label='numerical')
    ax.text(0.97,0.95, 't='+str(int(t)), transform=ax.transAxes, va='top', ha='right')
    
    ax.set_xlim([0,L])
    ax.set_ylim([0,4000])
    ax.set_xlabel('Horizontal distance (m)')
    ax.set_ylabel('Elevation (m)')
    ax.set_title(method+r' (dt='+str(int(dt))+')')
    ax.legend(loc="upper left")
    
    fig.canvas.draw()
    time.sleep(0.0001)<jupyter_output><empty_output><jupyter_text>Case 2: Sloped bed<jupyter_code>### choose bed shape
#bedshape = 'slope'
#bedshape = 'bump'
bedshape = 'trough'
%%time

fig,ax = plt.subplots(1,1)

plt.ion()
fig.show()
fig.canvas.draw()


# Initialize model
h=np.zeros(m) # initial thickness

# bed shape
angle=0.1 # bedrock slope in degrees
tan=np.tan(angle*np.pi/180)
if bedshape=='slope':
    b=-tan*x
elif bedshape=='bump':
    b=-tan*x+np.exp(-(x-L/2.)**2/(L/4.)**2)/(0.7*tan)
elif bedshape=='trough':
    b=-tan*x-np.exp(-(x-L/2.)**2/(L/4.)**2)/(0.7*tan)
else:
    b=np.zeros(m)

dt=5. # time step (years)
dtdx=dt/dx**2

# initialize vectors
alpha=np.zeros(m)
beta=np.zeros(m)
gamma=np.zeros(m)
delta=np.zeros(m)
f=np.zeros(m)
g=np.zeros(m)


# time stepping
for t in np.arange(0,time_end+1,dt):
    
    ax.clear()
    
    h=semi_implicit(b,h,dx,dt,dtdx,alpha,beta,gamma,delta,f,g)
        
    # plot bedrock
    ax.plot(x,b, c='black')
        
    # plot thickness profile
    ax.plot(x,b+h, c='orangered')
    ax.text(0.97,0.95, 't='+str(int(t)), transform=ax.transAxes, va='top', ha='right')
    
    ax.set_xlim([0,L])
    ax.set_ylim([np.min(b+h),3000])
    ax.set_xlabel('Horizontal distance (m)')
    ax.set_ylabel('Elevation (m)')
    ax.set_title('bed shape: '+bedshape+r' (dt='+str(int(dt))+')')
    
    fig.canvas.draw()
    time.sleep(0.0001)<jupyter_output><empty_output><jupyter_text>Excursus: "Human finite differences" Temperature diffusion according to$$\frac{\partial T}{\partial t}=k\frac{\partial^2 T}{\partial x^2}$$$k$: diffusion coefficient**FTCS scheme:**$$T_{i,t+1} = T_{i,t} + k\frac{\Delta t}{(\Delta x)^2} (T_{i+1,t} - 2T_{i,t} + T_{i-1,t})$$**Initial conditions**$T_{2,t=0}\dots T_{m-1,t=0}=10^\circ C$**Boundary conditions**$T_{1}=T_{m}=20^\circ C$ for all $t$<jupyter_code>m=10
k=0.5
Tmin=10.
Tmax=20.
T=np.ones(m)*Tmin
T[0]=T[m-1]=Tmax

plt.figure()

print('t=0: '+str(np.round(T,2)))

plt.plot(T)

# time stepping
for t in np.arange(1,11,1):
        
    Tn=np.array(T, copy=True) # correct
    #Tn=T # wrong
    
    for i in np.arange(1,m-1,1):
        Tn[i]=T[i]+k*(T[i+1]-2.*T[i]+T[i-1])
        
    T=np.array(Tn, copy=True) # correct
    #T=Tn # wrong
    
    print('t='+str(t)+': '+str(np.round(T,2)))

    plt.plot(T)
    
plt.xlim([0,m-1])
plt.ylim([Tmin-0.5,Tmax])
plt.xlabel('grid point $i$')
plt.ylabel('T ($^\circ C$)')
plt.title('temperature diffusion')

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.22728119758991241, 'fraction_numerical': 0.03230014286601652, 'mean_word_length': 4.4227012462108455, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 49, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '362706', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 6999, 'n_tokens_neox': 6727, 'n_words': 1177, 'path': 'introduction_to_numerical_ice_sheet_modelling.ipynb', 'repo_name': 'juliusgarbe/modelling_sandbox', 'size': '16125'}