Modules Python utilisés dans ce notebook
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
import seaborn as sns
from JSAnimation import IPython_display
%matplotlib inline
sns.set_context('notebook')
sns.set_palette('winter')
Modèle stream-power
In [2]:
# -- channel parameters
# channel length [m]
L = 30000
# channel initialization distance [m]
x_crit = 300
# Hacks's law parameters (drainage area -> x)
ka = 6.69
h = 1.67
# number of nodes
n_nodes = 100
# -- model parameters
# m/n ratio
m_n = 0.5
# n (slope exponent)
n = 1.0
# uplift rate [m yr^-1]
U = 0.004
# bedrock erodibility
K = 0.00005
# -- initial / boundary conditions
# boundary conditions at baselevel {'fixed', 'lowering', 'fixed_grad'}
b_level = 'fixed'
# baselevel lowering rate for b_level = 'lowering'
b_lower = 0.001
# initial channel profile {'ramp', 'flat', 'steady', 'prev'}
z_init = 'steady'
# initial channel slope for z_init = 'ramp'
islope = 0.01
# initial plateau elevation for z_init = 'flat'
plat_z = 1.0
# initial uplift rate for z_init = 'steady'
Ui = 0.002
# initial bedrock erodibility for z_init = 'steady'
Ki = 0.00005
# -- time parameters
# simulation duration [yr]
T = 100000
# stability parameter (CFL)
stabil = 1
# -- output parameters
# nb. of plotted intermediate solutions
n_plot = 5
# nb. of saved solutions (for animation)
n_save = 50
# ---------------
# This code has been adapted from the Matlab script
# "detach_ftfs.m" originally written by K.X. Whipple.
# check if the given boundary condition is valid
if b_level not in ['fixed', 'lowering', 'fixed_grad']:
raise ValueError("unknown boundary condition '{}'"
.format(b_level))
# functions to calculate some arrays
def log_slope(z, dx):
"""
Compute log slope from elevations.
"""
return np.log10((z[0:-1] - z[1:]) / dx)
def steady_profile(x, U, K):
"""
Compute the steady-state channel profile
(elevations) given U and K.
"""
x_temp = x**(1 - (h * m_n))
z_coef = (U / K)**(1. / n)
z_coef *= ka**(-m_n) * (1. - h * m_n)**-1
z = ((-1.0 * x_temp) + L**(1.0 - h * m_n)) * z_coef
return z
# determine values of derived constants and
# determine dt for standard stability condition
dx = 1.0 * L / n_nodes
m = m_n * n
dt = dx / (K * (ka**m) * (L**(h*m))) / stabil
# Set the x array and adjust n_nodes if needed
n_nodes = int(n_nodes - round(x_crit / dx) + 1)
x = np.linspace(x_crit, L, n_nodes)
# set the initial channel profile
if z_init == 'ramp':
z = (L - x) * islope
elif z_init == 'flat':
z = np.ones_like(x) * plat_z
z[-1] = 0
elif z_init == 'steady':
z = steady_profile(x, Ui, Ki)
elif z_init == 'prev':
try:
z[-1] = 0 # try taking the array already in memory
except NameError:
raise ValueError("can't use the previous solution, "
"no simulation run yet.")
else:
raise ValueError("unknown initial conditions '{}'"
.format(z_init))
z0 = z.copy()
z1 = np.empty_like(z)
# compute constants in the equation to be solved
# outside of the time loop to maximize efficiency.
dU = U * dt
Kx = (K * (ka**m) * dt / ((dx)**n)) * (x**(h*m))
area = ka * x**h
# compute log(slope) and log(area)
log_s = log_slope(z0, dx)
log_a = np.log10((area[0:-1] + area[1:]) / 2.)
# plot initial condition (in red)
fig, ax = plt.subplots(nrows=2)
ax[0].plot(x, z0, c='r')
ax[1].plot(log_a, log_s, c='r')
# init time increments and outputs
t = 0.
plot_inc = 0
plot_max_inc = int((T / dt) / n_plot)
save_inc = 0
save_max_inc = int((T / dt) / n_save)
if z_init != 'prev':
t_out = [t,]
z_out = [z,]
log_s_out = [log_s,]
# Main time loop
while t < T + dt:
# calculate channel profile after one time step
z1[:-1] = z[:-1] + dU
z1[:-1] -= Kx[:-1] * np.abs(z[1:] - z[:-1])**n
# stability => no increase in elevation
# along the profile
z1[:-1] = np.maximum(z1[:-1], z[1:])
# boundary conditions
if b_level == 'fixed':
z1[-1] = z[-1]
elif b_level == 'lowering':
z1[-1] = z[-1] - (b_lower * dt)
elif b_level == 'fixed_grad':
z1[-1] = z1[-2] - (z[-2] - z[-1])
# update z
z = z1.copy()
# log(slope)
log_s = log_slope(z, dx)
# plot intermediate solutions
if plot_inc > plot_max_inc:
clr_index = int(255. * t / T)
clr = plt.cm.winter(clr_index)
ax[0].plot(x, z, c=clr)
ax[1].plot(log_a, log_s, c=clr)
plot_inc = 0
if save_inc > save_max_inc:
t_out.append(t)
z_out.append(z)
log_s_out.append(log_s)
save_inc = 0
# increment time step
t += dt
plot_inc += 1
save_inc += 1
# save / plot the final solution (in green)
log_s = log_slope(z, dx)
ax[0].plot(x, z, c='g')
ax[1].plot(log_a, log_s, c='g')
t_out.append(t)
z_out.append(z)
log_s_out.append(log_s)
# calculate and plot the analytic solution at steady-state
z_steady = z[-1] + steady_profile(x, U, K)
log_s = log_slope(z_steady, dx)
ax[1].plot(log_a, log_s, c='k', alpha=0.7, ls='--')
ax[0].plot(x, z_steady, c='k', alpha=0.7, ls='--')
# plot refinements
fig.set_size_inches(9, 9)
plt.setp(ax[0],
xlabel='distance from divide [m]',
ylabel='elevation [m]')
plt.setp(ax[1],
xlabel='log drainage area [$\log$(m$^2$)]',
ylabel='log channel gradient [#]')
plt.tight_layout()
Le code dans la cellule ci-dessous crée une animation à partir du (des) résultat(s) de la (des) dernière(s) simulation(s) effectuée(s) ci-dessus. L'exécution peut prendre du temps !
In [3]:
# create an animation from a matplotlib figure
# using JSAnimation: https://github.com/jakevdp/JSAnimation
# may take a while to complete
z_arr = np.array(z_out)
log_s_arr = np.array(log_s_out)
t_arr = np.array(t_out) / 1e3
fig, ax = plt.subplots(nrows=2)
fig.set_size_inches(9, 9)
line1, = ax[0].plot([], [], c='b', lw=2)
line2, = ax[1].plot([], [], c='b', lw=2)
tlabel = plt.text(0.8, 0.9, '', transform=ax[0].transAxes)
plt.setp(ax[0],
xlim=[0, L],
ylim=[z_arr.min(), z_arr.max()],
xlabel='distance from divide [m]',
ylabel='elevation [m]')
plt.setp(ax[1],
xlim=[log_a.min(), log_a.max()],
ylim=[log_s_arr.min(), log_s_arr.max()],
xlabel='log drainage area [$\log$(m$^2$)]',
ylabel='log channel gradient [#]')
plt.tight_layout()
def init():
line1.set_data(x, z_arr[0])
line2.set_data(log_a, log_s_arr[0])
tlabel.set_text('time: {:.3f} kyr'.format(t_arr[0]))
return line1, line2, tlabel
def animate(i):
line1.set_data(x, z_arr[i])
line2.set_data(log_a, log_s_arr[i])
tlabel.set_text('time: {:.3f} kyr'.format(t_arr[i]))
return line1, line2, tlabel
animation.FuncAnimation(fig, animate, init_func=init,
frames=z_arr.shape[0]-1, interval=50,
blit=True)
Out[3]:
In [3]: