Modules Python utilisés dans ce notebook :
In [1]:
import numpy as np
from scipy.sparse import dia_matrix, csc_matrix
from scipy.sparse.linalg import spsolve
from scipy import stats
from matplotlib import pyplot as plt
from matplotlib import animation
import seaborn as sns
from JSAnimation import IPython_display
from IPython.display import YouTubeVideo
%matplotlib inline
sns.set_context('notebook')
sns.set_palette('winter')
Dans la plupart des cas, on veut modéliser l'évolution de la topographie, $z(x,y,t)$, au cours du temps sur des échelles allant de $10^2$ à $>10^6$ années. Cette évolution est réprésentée sous sa forme la plus générique par :
$$ \frac{\partial z}{\partial t} = U - E $$où $z$ (L) est l'élévation de la surface topographique, $U$ (L T$^{-1}$) est le taux de variation de l'élevation de la surface topographique sous l'action des processus tectoniques (> 0), et $E$ (L T$^{-1}$) est le taux de variation de l'élévation de la surface topographique sous l'action des processus érosifs (> 0), de dépôt (< 0).
Parfois, on cherche à modéliser l'épaisseur de sol (régolithe), l'épaisseur de la couche active dans les chenaux, etc... L'évolution de ces variables au cours du temps ont un impact direct sur l'évolution de l'élévation. Certains modèles intégrent aussi les composantes latérales des processus tectoniques et/ou érosifs.
Pour la majorité des modèles, $U$ est considéré comme une "condition à la limite", imposée au système (forçage), alors que $E$ est paramétré notamment en fonction de variables dérivées - localement ou non - de la surface topographique (pente, courbure, aire drainée...).
Outre leur paramétrisation des processus érosifs, les modèles se distinguent par le nombre de dimensions spatiales prises en compte :
Quelques liens intéressants:
Un bon article général sur cette thématique: Tucker & Hancock (2010)
Quelques illustrations:
In [2]:
YouTubeVideo('T0BWWjSvK30', width=600, height=450)
Out[2]:
In [3]:
YouTubeVideo('U7IC-L2fq2o', width=600, height=450)
Out[3]:
D'autres animations (modèles CASCADE, DOUAR - Jean Braun).
Comme beaucoup d'autres systèmes physiques, au départ d'un état initial ou d'une perturbation (d'origine climatique ou tectonique par exemple), la surface topographique peut en théorie évoluer jusqu'à atteindre un certain équilibre dynamique (steady-state), lors duquel la surface n'évolue plus en apparence mais est toujours sous l'influence des processus actifs d'érosion et de soulèvement (parfaitement balancés en tout point de l'espace) :
$$ \begin{align} \frac{\partial z}{\partial t} & = 0 \\ U & = E \end{align} $$Cet état permet donc de s'affranchir du terme $\frac{\partial z}{\partial t}$ et d'établir une relation fonctionnelle claire - et invariante dans le temps - entre les différentes variables du système (par exemple, topographie vs. paramètres du modèle). C'est pourquoi la plupart des études théoriques s'attachent à décrire cet état en premier lieu. Une telle relation fonctionelle peut aussi permettre une confrontation bien plus aisée du modèle aux données de terrain (en supposant toutefois que cet état est atteint dans le milieu naturel étudié).
L'érosion des versants est le plus souvent modélisée comme étant un processus purement diffusif (diffusion linéaire et homogène). Considérant le profil à une dimension d'un versant, l'équation de diffusion s'écrit :
$$ \frac{\partial z}{\partial t} = U + K \frac{\partial^2 z}{\partial x^2} $$où $x$ (L) est la position le long du profil, et $K$ (L$^2$ T$^{-1}$) est une constante (diffusivité).
Postulat de base et développement
L'usage de l'équation de diffusion repose sur une idée vieille de plus d'un siècle (Davis, 1892; Gilbert, 1909) selon laquelle la quantité de sédiments transportés le long des versants est linéairement proportionelle à la pente locale :
$$ q_s = - K \frac{\partial z}{\partial x} $$où $q_s$ (L$^2$ T$^{-1}$) est le flux volumique sédimentaire par unité de largeur (> 0 si la pente est "descendante"), et $\frac{\partial z}{\partial x}$ est le gradient topographique (variation locale de l'élévation, ou pente locale, < 0 de haute vers basse élévation). A noter que cette simple relation est aussi utilisée dans d'autres domaines, par exemple pour le transport de chaleur vs. température (loi de Fourier), pour le transport d'un fluide à travers un milieu poreux vs. charge (loi de Darcy), pour le courant électrique vs. potentiel électrique (loi d'Ohm) ou de manière générale en théorie de la diffusion (première loi de Fick). En géomorphologie, on parle souvent de Geomorphic Transport Laws - GTLs (Dietrich et al., 2003).
En appliquant la conservation de la masse (voir figure ci-dessous), on a
$$ \rho \frac{\partial z}{\partial t} = \rho U - \rho \frac{\partial q_s}{\partial x} $$où $\frac{\partial q_s}{\partial x}$ représente la variation locale du flux sédimentaire (notion de divergence) et $\rho$ (M L$^{-3}$) représente invariablement la masse volumique du substratum et de la couche sédimentaire en surface (simplification). Après développement, on retrouve l'équation de diffusion écrite ci-dessus :
$$ \begin{align} \frac{\partial z}{\partial t} & = U - \frac{\partial}{\partial x} \left(- K \frac{\partial z}{\partial x}\right) \\ & = U + K \frac{\partial^2 z}{\partial x^2} \end{align} $$A noter que $\frac{\partial^2 z}{\partial x^2}$ est la variation locale de la pente, c'est à dire la courbure locale du terrain ($> 0 \rightarrow$ concave, $< 0 \rightarrow$ convexe).
La figure ci-dessous schématise l'évolution d'un profil de versant résultant de l'action combinée du soulèvement tectonique $U$ et du transport sédimentaire $q_s$ (conservation de la masse appliquée à un élement 1D de versant).
Observations
Gilbert (1909) avait fait l'hypothèse d'un transport linéarement dépendant à la pente locale pour expliquer l'origine de la forme convexe de la partie supérieure des versants. Cette hypothèse va de pair avec celle d'un équilibre dynamique atteint sur ces parties de versants (voir exercices ci-dessous).
Certaines observations, bien que peu nombreuses, soutiennent ce modèle. Par exemple, McKean et al. (1993) ont calculé des flux sédimentaires à partir de concentrations en isotopes $^{10}$Be mesurées sur les parties convexes de versants situés en Californie. Considérant une épaisseur de sol uniforme sur le versant et un équilibre entre apport (production in-situ) et perte (transport sédimentaire) de $^{10}$Be, leurs résultats montrent une dépendence linéaire du flux sédimentaire par rapport à la pente locale :
Small et al. (1999) (figure ci-dessous) et Heimsath et al. (1999, 2000) ont également effectué des mesures en isotopes cosmogéniques qui supportent le modèle (du moins, en partie).
Par souci de simplification, on ne considère ici aucun soulèvement tectonique ($U = 0$). L'équation de diffusion à une dimension devient alors
$$ \frac{\partial z}{\partial t} = K \frac{\partial^2 z}{\partial x^2} $$c'est à dire une équation aux dérivées partielles linéaire et homogène (pas de terme source).
Comme d'autres équations aux dérivées partielles ou différentielles, cette équation peut avoir une solution analytique - appelée fonction de Green - pour autant que l'on spécifie des conditions initiales et aux bords du domaine particulières. Si on considère que le profil initial du versant est donné par la distribution de Dirac, $z(x,0) = \delta(x)$, alors on peut montrer qu'au temps $t$, la distribution de l'élévation $z$ pour $x \in [-\infty, +\infty]$ est donnée par :
$$ z(x,t) = G(x,t) = \frac{1}{\sqrt{4 \pi K t}} \exp \left(- \frac{x^2}{4 K t} \right) $$Cette fonction gaussienne $G(x,t)$ est la fonction de Green propre à l'équation de diffusion sans conditions aux bords particulières et sans terme source. Elle est aussi appelée solution fondamentale de l'équation de diffusion.
Exercice a: calculer et illustrer sur un graphique la solution fondamentale de l'équation de diffusion, pour différentes valeurs de $t$.
Exercice b: reproduire le graphique en faisant varier le paramètre $K$. Comment $K$ influence-t-il la solution ?
Exercice c: on a vu que $K$ détermine la magnitude du flux sédimentaire linéairement proportionel à la pente, ce paramètre étant considéré constant dans le temps et l'espace. Au vu de l'influence de ce paramètre sur la solution de l'équation de diffusion et sur base d'une bonne connaissance de la géomorphologie des versants, de quels facteurs $K$ depend-il ? Ces facteurs sont-ils uniquement liés au transport sédimentaire ? Sont-ils constants dans le temps et l'espace ?
In [4]:
# diffusivity
K = 0.01
# --------------
# ** solution **
# set array of x values
x = np.linspace(-20, 20, 200)
# init figure
ax = plt.subplot(111)
# plot z(x,t) for different values of t (log steps)
for t in np.logspace(2, 4, 5):
z_t = (1. / np.sqrt(4 * np.pi * K * t)) * np.exp(- x**2 / (4 * K * t))
ax.plot(x, z_t)
Solution exercice b: on voit que la diffusivité $K$ influence l'efficacité globale avec laquelle le versant est érodé. Plus la valeur de $K$ est élevée, plus l'érosion est intense et donc plus la courbe de Gauss est élargie pour les mêmes valeurs de $t$.
Solution exercice c: fonction de la pente, le transport sur les versants sera aussi plus ou moins efficace en fonction de l'activité biologique, du climat (cycles gel/dégel ou humide/sec), du type de sol ou de la nature des sédiments transportés. $K$ mesure tous ces facteurs. $K$ est exprimé en L$^2$ T$^{-1}$ et intègre aussi la profondeur sous la surface topographique, le long de laquelle le transport s'effectue. Enfin, l'effacité érosive ne dépend pas uniquement de l'efficacité avec laquelle les sédiments sont transportés sur le versant, cela dépend aussi de l'efficacité avec laquelle la roche mère est dégradée en sédiment mobile. $K$ mesure aussi cette efficacité, qui dépend du climat, du type de roche, etc... Le paramètre $K$ dépend donc de tous ces facteurs. Considérer $K$ constant dans le temps et l'espace revient donc à considérer tous ces facteurs comme constants.
Le profil de versant initial défini pour calculer la solution fondamentale ci-dessus n'est évidemment pas du tout réaliste, il signifie que toute la matière sous le versant (substratum + sol) est concentrée en un seul point de l'espace $x=0$ ! Mais ce cas de figure est tout de même utile. En effet, on peut montrer que l'on peut calculer analytiquement l'érosion par diffusion à partir de conditions initiales quelconques $z(x, 0) = g(x)$, la solution étant donnée par la convolution de $g(x)$ avec la solution fondamentale :
$$ z(x, t) = G(x, t) \ast g(x) $$$\ast$ étant l'opérateur de convolution.
Dans l'exemple ci-dessous (d'après Avouac, 1993), on applique cette solution convoluée à partir d'un escarpement de faille défini de la manière suivante:
$$ g(x) = a H(x) + b x $$où $a$ est la moitié de la hauteur de l'escarpement, $b$ est la pente du versant environnant et $H(x)$ est donné par
$$ \begin{align} H(x) & = - \frac{1}{2}, & \text{si} \quad x < 0, \\ & = \frac{1}{2}, & \text{si} \quad x > 0, \end{align} $$Le code ci-dessous génère une figure représentant le profil initial de l'escarpement (rouge), le profil après érosion par diffusion (vert) ainsi que les profils intermédiaires (bleu).
Exercice: lancer la cellule de code ci-dessous et faire varier les paramètres en entête.
In [5]:
# half scarp offset [m]
a = 20.
# regional slope [#]
b = 0.2
# diffusivity [m^2 yr^-1]
K = 0.01
# erosion duration [yr]
t = 1e4
# --------------
# set x array
# (set L >> domain of interest
# to avoid convolution boundary effects)
L = 200.
x = np.linspace(-L, L, 2*L)
# set initial elevation
H = np.ones_like(x) * 0.5
H[np.nonzero(x < 0)] = -0.5
z0 = 2. * a * H + b * x
# plot initial elevation (in red)
ax = plt.subplot(111)
ax.plot(x, z0, c='r')
# diffusion kernel (for convolution)
def diffusion_kernel(x, t):
return ((1. / np.sqrt(4 * np.pi * K * t))
* np.exp(- x**2 / (4 * K * t)))
# calculate and plot the solution for t and intermediate times
t_array = np.linspace(0, t, 4)[1:]
for ti in t_array:
k = diffusion_kernel(x, ti)
z_ti = np.convolve(z0, k, mode='same')
ax.plot(x, z_ti)
# plot final solution in green
ax.plot(x, z_ti, c='g')
# plot axis labels and automatic limits
z0_subset_ind = np.nonzero(
np.logical_and(x > -L/3, x < L/3)
)
z0_subset = z0[z0_subset_ind]
_= plt.setp(ax,
xlim=[-L/4, L/4],
ylim=[z0_subset.min(), z0_subset.max()],
xlabel='x [m]',
ylabel='relative elevation [m]')
Comme c'est souvent le cas avec les équations aux dérivées partielles, la résolution analytique de l'équation de diffusion peut très vite se compliquer - voire devenir impossible - lorsqu'on considère des cas de figures un peu plus élaborés que ceux présentés ci-dessus (conditions aux bords, termes sources, variation spatiale et/ou temporelle de certains paramètres...). On utilise alors des méthodes numériques pour résoudre l'équation différentielle (par ex. les différences finies - voir prochaine séance de TP, les volumes finis ou les éléments finis).
Le code ci-dessous implémente un schéma numérique aux différences finies (Euler implicite, centré spatialement d'ordre 2) pour l'équation de diffusion incluant le terme $U$. Aux pieds du profil du versant modélisé, l'élévation est fixée et reste constante dans le temps (conditions aux bords dites de Dirichlet). Ces aspects de la résolution numérique de l'équation de diffusion seront vus en détails lors de la prochaine séance de TP.
Plusieurs conditions initiales, z_init, sont possibles: un profil plat d'élévation nulle ('flat'), une fonction sinus ('sin'), deux pentes rectilignes ('tri') ou le résultat de la simulation précédente ('prev').
Deux figures sont produites: le profil du versant (élévation) et la courbure (négative) du terrain. Le profil initial est dessiné en rouge, la solution finale en vert, et les solutions intermédiaires en bleu.
Exercice a: Lancer une simulation avec z_init = 'flat', L = 300, Dx = 5, K = 0.01 et U = 0.0001. Le système évolue-t-il vers un état d'équilibre ? Si oui, après combien de temps et comment expliquer la forme du profil du versant ?
Exercice b: Faire varier les paramètres K et U. Comment influencent-ils le système (forme du versant, temps requis pour atteindre l'équilibre...) ? Montrer qu'il est possible de prédire exactement le profil du versant à l'équilibre, si $K$ et $U$ sont connus.
Exercice c: Choisir différentes conditions initiales ('tri', 'sin'). Quel impact sur la solution finale ?
Exercice d: Lancer une permière simulation avec les paramètres de l'exercice a. Lancer une deuxième simulation avec z_init = 'prev' et changeant la valeur de K et/ou U. Comment évolue le système ?
Exercice e: Lancer une série de simulations avec z_init = 'tri', L = 300, H = 60, dist_tb = 30 K = 0.01 et U = 0, et en faisant varier Dx. Est-ce que cela à un impact sur la solution ? Quelle conclusion importante en tirer ?
In [6]:
# -- hillslope parameters
# hillslope length [m]
L = 300.
# initial elevation {'flat', 'sin', 'tri', 'prev'}
z_init = 'flat'
# hillslope height [m] (for z_init = 'sin' or 'tri')
H = 60.
# distance between hillslope top-base [m] (for z_init = tri')
dist_tb = 30.
# fixed node-spacing [m]
Dx = 2.
# -- model parameters
# diffusivity [m^2 yr^-1]
K = 0.01
# uplift rate [m yr^-1]
U = 0.0001
# -- time parameters
# fixed time step duration [yr]
Dt = 10000.
# simulation duration [yr]
T = 2e6
# -- output parameters
# nb. of plotted intermediate solutions
n_plot = 5
# nb. of saved solutions (for animation)
n_save = 50
# --------------
def neg_curvature(z, Dx):
"""
Return negative curvature given elevations `z`
and node-spacing `Dx` (boundary effects!)
"""
c = np.zeros_like(z)
c[1:-1] = (z[2:] - 2 * z[1:-1] + z[:-2]) / (Dx**2)
return c
# create mesh (add two boundary nodes as
# the first and last elements of the `x` vector)
x = np.arange(-Dx, L + 2 * Dx, Dx)
# set initial conditions
if z_init == 'flat':
z = np.zeros_like(x)
elif z_init == 'sin':
z = (-np.cos(x * 2 * np.pi / L) + 1.) * H / 2.0
elif z_init == 'tri':
x_top = L / 2.
x_left, x_right = x_top - dist_tb, x_top + dist_tb
s_left = H / (x_top - x_left)
d_left = s_left * x + H - s_left * x_top
s_right = -H / (x_right - x_top)
d_right = s_right * x - s_right * x_right
z = np.minimum(d_left, d_right)
z[np.nonzero(z < 0)] = 0.
elif z_init == 'prev':
try:
z = z # 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))
c = neg_curvature(z, Dx)
# set boundary conditions and working arrays
z[0] = z[-1] = 0.
z0 = z.copy()
z1 = z.copy()
b = np.zeros_like(z[1:-1])
b[0] = z[0]
b[-1] = z[-1]
c0 = c.copy()
# calculate the CFL
CFL = K * Dt / Dx**2
print "CFL: {}".format(CFL)
# build the following tridiagonal matrix (N-2, N-2)
# N is the number of nodes
# ( 1+2CFL -CFL 0 ... 0 )
# ( -D 1+2CFL -CFL ... 0 )
# ( 0 -CFL 1+2CFL ... 0 )
# ( ... ... ... ... ... )
# ( ... ... -CFL 1+2CFL -CFL )
# ( 0 0 0 ... -CFL 1+2CFL)
diag_12CFL = np.ones_like(x) * (1 + 2 * CFL)
diag_CFL = np.ones_like(x) * -CFL
A = dia_matrix(
(np.array([diag_CFL, diag_12CFL, diag_CFL]),
np.array([-1, 0, 1])),
shape=(x.size - 2, x.size - 2)
)
A_csc = csc_matrix(A)
# init and plot initial condition in red
ax = plt.subplot(211)
ax2 = plt.subplot(212)
ax.plot(x, z, c='r')
ax2.plot(x[2:-2], c[2:-2], c='r')
# time init, main loop and saving/plot 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,]
c_out = [c,]
while t < T:
# right side of linear system for implicit scheme
r = z[1:-1] + CFL * b
# solve linear system (sparse matrix)
z1[1:-1] = spsolve(A_csc, r)
# add uplift
z1[1:-1] += U * Dt
# update elevations (Z)
z = z1.copy()
# compute negative curvature
c = neg_curvature(z, Dx)
# plot/save intermediate solution
if plot_inc > plot_max_inc:
clr_index = int(t / T * 255.)
clr = plt.cm.winter(clr_index)
ax.plot(x, z, color=clr)
ax2.plot(x[2:-2], c[2:-2], color=clr)
plot_inc = 0
if save_inc > save_max_inc:
t_out.append(t)
z_out.append(z)
c_out.append(c)
save_inc = 0
# increment time step
t += Dt
plot_inc += 1
save_inc += 1
# save/plot final solution in green
ax.plot(x, z, c='g')
ax2.plot(x[2:-2], c[2:-2], c='g')
t_out.append(t)
z_out.append(z)
c_out.append(c)
# print elevation at hill top
# at the start and end of the run
print "hilltop (start): {} m".format(z0.max())
print "hilltop (end): {} m".format(z.max())
# plot axis labels
_= plt.setp(ax,
xlim=[0, L],
xlabel='x [m]',
ylabel='elevation [m]')
_= plt.setp(ax2,
xlim=[0, L],
xlabel='x [m]',
ylabel='negative curvature [m$^{-1}$]')
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 [7]:
# 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)
c_arr = np.array(c_out)
t_arr = np.array(t_out) / 1e3
fig, ax = plt.subplots(nrows=2, sharex=True)
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='x [m]',
ylabel='elevation [m]')
plt.setp(ax[1],
xlim=[0, L],
ylim=[c_arr[:, 2:-2].min(), c_arr[:, 2:-2].max()],
xlabel='x [m]',
ylabel='negative curvature [m$^{-1}$]')
def init():
line1.set_data(x, z_arr[0])
line2.set_data(x[2:-2], c_arr[0][2:-2])
tlabel.set_text('time: {} kyr'.format(t_arr[0]))
return line1, line2, tlabel
def animate(i):
line1.set_data(x, z_arr[i])
line2.set_data(x[2:-2], c_arr[i][2:-2])
tlabel.set_text('time: {} 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[7]:
Solution exercice a: Le système évolue bien vers un état d'équilibre et le profil du versant tend vers une forme convexe. La convexité constante s'explique en développant l'équation de diffusion à l'état d'équilibre ($\frac{\partial z}{\partial t} = 0$)
Solution exercice b: Les paramètres $U$ et $K$ ont un impact sur le temps mis pour atteindre l'équilibre et sur le degré de convexité de la forme du profil à l'équilibre. Plus précisément, la convexité (ou courbure locale négative), constante à l'équilibre, est donnée par le rapport $\frac{U}{K}$.
Solution exercice c: Les conditions initiales différentes n'ont pas d'impact sur la forme du profil de versant à l'équilibre, qui ne dépend que de $U$ et $K$. Par contre, le temps mis pour atteindre l'équilibre diffère.
Solution exercice d: Le système évolue vers un nouvel état d'équilibre, défini par un nouveau rapport $\frac{U}{K}$.
Solution exercice e: Dans la simulation définie, modifier Dx a un impact sur l'aspect du profil initial du versant. Par exemple, de Dx = 5 m à Dx = 50 m, les pentes du versant initial sont plus faibles, pour la même fonction donnée $z(x)$. Par conséquent, l'érosion - linéarement dépendante à la pente - est moins efficace et le temps mis pour éroder tout le versant est plus long. Ceci montre que l'efficacité de l'érosion dépend de l'échelle à laquelle on discrétise la topographie ou avec laquelle on mesure la pente. La calibration de la diffusivité $K$ est donc dépendante de cette échelle !
Le géomorphologue n'est, a-priori, sans doute pas pleinement satisfait de l'utilisation de l'équation de diffusion pour expliquer les formes convexes des versants. Il est attaché à la description des processus de transport sédimentaire, et l'explication du transport linéairement proportionel à la pente locale ne suffit pas, cela reste un modèle empirique vu comme une "boite noire".
Il est sans doute plus parlant pour un géomorphologue de décrire ce qu'il se passe au niveau "micro", c'est à dire au niveau des déplacements individuels des particules de sol mobile.
Comme modèle alternatif, on pourrait ainsi postuler que le transport sédimentaire sur les versants aux échelles de temps "géomorphologiques" est le résultat d'un grand nombre de déplacements successifs - à occurence régulière - de chacune des particules sédimentaires disposées sur le versant, la longueur de chaque déplacement étant aléatoire, indépendante des précédents déplacements (et des déplacements des autres particules) et réalisée à partir d'une certaine densité de probabilité donnée. Ce processus stochastique est bien connu en physique statistique, il est appelé "marche aléatoire" (random walk).
Considérant certains processus de transport sur les versants comme l'effet splash ou la bioturbation, on peut en effet raisonnablement penser que chaque déplacement individuel d'une particule est aléatoire - bien que dirigé dans le sens de la pente -, indépendant et restreint à son environnement immédiat (quelques centimètres).
NOTE: Ce modèle de marche aléatoire adopte en fait une approche lagrangienne, dans le sens où on cherche à décrire l'évolution du système à travers l'évolution des particules en mouvement qui forment ce système. L'approche utilisée ci-dessus pour l'équation de diffusion est, au contraire, une approche eulérienne car on cherche à décrire l'évolution du système à travers l'évolution de ses états (flux sédimentaire, élévation) en chaque point de l'espace.
Exercice a: compléter le code ci-dessous afin d'implementer une marche aléatoire à une dimension pour $P$ particules. Toutes les particules sont initialement situées en un seul et même point de l'espace $x = 0$. Le code est sensé créer une figure montrant la distribution de la position finale des particules, après un grand nombre $N$ de déplacements indépendants. Chaque déplacement a une longueur $\Delta x$ générée aléatoirement à partir de la loi normale $Z(\mu, \sigma)$, de moyenne $\mu = 0$ et d'écart-type $\sigma$. Suggestion : utiliser le module numpy.random pour générer des valeurs indépendantes à partir de la distribution $Z(\mu, \sigma)$.
Exercice b: Faire varier N. Quelle forme prend la distribution de la position des particules après $N$ déplacements ? Le montrer si possible en utilisant un test statistique disponible dans le module scipy.stats, importé sous le nom stats dans ce notebook.
Exercice c: Utiliser d'autres distributions disponible dans le module scipy.stats et ayant une variance finie (par exemple Gumbel, exponentielle...). Quel impact ont-elles sur la distribution obtenue (pour une grande valeur de N) ?
Exercice d: Si on compare la distribution obtenue avec la solution fondamentale de l'équation de diffusion, que peut-on en conclure ?
In [30]:
# variance of the probability density function
# for particle jump lengths [m]
sigma = 10e-2
# number of particles
P = 10000
# number of jumps
N = 1000
# --------------
# ** solution **
# generate P*N independent realizations of the normal law
# for each jump of each particle
jumps = np.random.normal(scale=sigma, size=P * N)
jumps = jumps.reshape(P, N)
# calculate the final x position after N jumps
# for each particle (considering initial x position
# equals 0 for all particles)
x_final = jumps.sum(axis=1)
# plot the distribution of `x_final`
ax = plt.subplot(111)
sns.distplot(x_final, color='b', ax=ax)
_ = plt.setp(ax,
yticks=[],
xlabel='x [m]')
# perform the Anderson-Darling test for
# `x_final` coming from a normal distribution (see doc)
A2, crit_vals, conf_levels = stats.anderson(x_final, dist='norm')
print "Anderson-Darling test"
print "A2: {}".format(A2)
print "critical values: {}".format(crit_vals)
print "confidence levels: {}".format(conf_levels)
Solution exercice b: La distribution prend la forme d'une distribution normale lorsque N est grand (on ne peut pas rejeter l'hypothèse de normalité).
Solution exercice c: Quelle que soit la distribution de probabilité (de variance finie) utilisée pour les déplacement des particules, la distribution obtenue est une distribution normale pour grand N.
Solution exercice d: On remarque que la distribution normale obtenue est semblable à la solution fondamentale de l'équation de diffusion. Or le problème posé ici est similaire au problème posé pour calculer la solution fondamentale (toutes les particules concentrées initialement en $x = 0$ $\sim$ distribution de Dirac). La distribution obtenue représente aussi un échantillon de la densité de probabilité de la position d'une particule après $N$ déplacements. En fait, on peut montrer que pour toute série de déplacements indépendants et aléatoires de variance finie, les deux modèles sont équivalents à la limite (c'est à dire pour $N \rightarrow \infty$). L'un (marche aléatoire) décrit le système au niveau "micro", c'est à dire les déplacements aléatoires des particules, répétés sur de courtes échelles de temps. L'autre (diffusion) décrit le système au niveau "macro", c'est à dire l'évolution du versant à long terme.
Question: En quoi le modèle présenté dans la section 3 peut être discutable, concernant le transport sur les versants en général ?
Réponse: Le modèle comporte pas mal de simplifications, comme par exemple la non prise en compte de la taille ou de l'enfouissement des particules conditionant leur déplacement. En outre, d'autres processus de transports peuvent produire, lors d'un seul événement, un déplacement allant généralement bien au délà de l'environnement immédiat de la particule (par exemple, glissements de terrain, coulées boueuses, détachements de la couche active, etc...).
Si certaines observations, comme celles présentées précédemment, soutiennent le modèle de diffusion linéaire simple, d'autres observations s'en écartent.
Heimsath et al. (2005) ont notamment montré qu'une loi de transport incluant l'épaisseur de sol, $h$, explique mieux les observations collectées dans différents sites d'études (Californie, Australie), sur des étendues qui ne sont plus seulement limitées aux parties convexes sommitales des versants. La loi de transport proposée par les auteurs est la suivante :
$$ q_s = - K h \frac{\partial z}{\partial x} $$Furbish et al. (2009) relient cette relation avec un modèle particulaire (niveau "micro") bien plus élaboré que celui décrit dans la section 3, qui prend notamment en compte la profondeur à laquelle se trouvent les particules en mouvement.
Par ailleurs, Roering et al. (1999) ont montré que les pentes rectilines observées au sein de certains versants ne peuvent pas être expliquées par la diffusion linéaire. Ils proposent que le flux sédimentaire varie plutôt non-linéairement avec la pente locale :
$$ q_s = - \frac{K \frac{\partial z}{\partial x}}{1 - (\frac{\partial z}{\partial x} / S_c)^2} $$$S_c$ étant la pente critique à partir les laquelle de flux sédimentaire est infini. Cette loi de transport est mieux adaptée aux flux sédimentaires calculés indépendament sur base de données topographiques qu'ils ont collectées (figure ci-dessous).
Foufoula-Georgiou et al. (2010) ont eux montré qu'un modèle plus général pouvait également expliquer ces mêmes observations. Le modèle est ici basé sur un transport non-local, c'est à dire un transport qui ne dépend plus seulement de la pente locale mais également des pentes situées en amont du versant. Cela se traduit par une "super-diffusion" :
$$ \frac{\partial z}{\partial t} = U + K \frac{\partial^\alpha z}{\partial x^\alpha} $$où $1 < \alpha < 2$ et $\frac{\partial^\alpha z}{\partial x^\alpha}$ est une dérivée fractionelle (opérateur non-local). Le cas $\alpha = 2$ correspond à l'équation de diffusion.
Ce modèle a l'avantage d'être théoriquement lié à l'approche de marche aléatoire présentée dans la section 3 (équivalence "macro" / "micro"), mais où la distribution de la longueur des déplacements des particules n'a ici plus de variance finie lorsque $1 < \alpha < 2$ (la variance diverge lorsque le nombre de déplacements augmente). Ce cas de figure a du sens pour des processus rapides comme les glissements de terrain.
Braun et al. 2001 ont enfin montré que plusieurs lois de transports combinées pouvaient expliquer aussi bien les formes convexes/divergentes en haut de versant que les parties concaves/convergentes en bas de versant.
Tucker & Bradley (2010) présentent également une critique de la diffusion, sur base d'un modèle particulaire simple.
In [8]: