# Exploring the diffusion equation with Python

Ever since I became interested in science, I started to have a vague idea that calculus, matrix algebra, partial differential equations, and numerical methods are all fundamental to the physical sciences and engineering and they are linked in some way to each other. The emphasis here is on the word vague; I have to admit that I had no clear, detailed understanding of how these links actually work. It seems like my formal education both in math and physics stopped just short of where everything would have nicely come together. Papers that are really important in geomorphology, sedimentology or stratigraphy seemed impossible to read as soon as they started assuming that I knew quite a bit about convective acceleration, numerical schemes, boundary conditions, and Cholesky factorization. Because I didn’t.

So I have decided a few months ago that I had to do something about this. This blog post documents the initial – and admittedly difficult – steps of my learning; the purpose is to go through the process of discretizing a partial differential equation, setting up a numerical scheme, and solving the resulting system of equations in Python and IPython notebook. I am learning this as I am doing it, so it may seem pedestrian and slow-moving to a lot of people but I am sure there are others who will find it useful. Most of what follows, except the Python code and the bit on fault scarps, is based on and inspired by Slingerland and Kump (2011): Mathematical Modeling of Earth’s Dynamical Systems (strongly recommended). You can view and download the IPython Notebook version of this post from Github.

Estimating the derivatives in the diffusion equation using the Taylor expansion

This is the one-dimensional diffusion equation:

$\frac{\partial T}{\partial t} - D\frac{\partial^2 T}{\partial x^2} = 0$

The Taylor expansion of value of a function u at a point $\Delta x$ ahead of the point x where the function is known can be written as:

$u(x+\Delta x) = u(x) + \Delta x \frac{\partial u}{\partial x} + \frac{\Delta x^2}{2} \frac{\partial^2 u}{\partial x^2} + \frac{\Delta x^3}{6} \frac{\partial^3 u}{\partial x^3} + O(\Delta x^4)$

Taylor expansion of value of the function u at a point one space step behind:

$u(x-\Delta x) = u(x) - \Delta x \frac{\partial u}{\partial x} + \frac{\Delta x^2}{2} \frac{\partial^2 u}{\partial x^2} - \frac{\Delta x^3}{6} \frac{\partial^3 u}{\partial x^3} + O(\Delta x^4)$

Solving the first Taylor expansion above for $\frac{\partial u}{\partial x}$ and dropping all higher-order terms yields the forward difference operator:

$\frac{\partial u}{\partial x} = \frac{u(x+\Delta x)-u(x)}{\Delta x} + O(\Delta x)$

Similarly, the second equation yields the backward difference operator:

$\frac{\partial u}{\partial x} = \frac{u(x)-u(x-\Delta x)}{\Delta x} + O(\Delta x)$

Subtracting the second equation from the first one gives the centered difference operator:

$\frac{\partial u}{\partial x} = \frac{u(x+\Delta x)-u(x-\Delta x)}{2\Delta x} + O(\Delta x^2)$

The centered difference operator is more accurate than the other two.

Finally, if the two Taylor expansions are added, we get an estimate of the second order partial derivative:

$\frac{\partial^2 u}{\partial x^2} = \frac{u(x+\Delta x)-2u(x)+u(x-\Delta x)}{\Delta x^2} + O(\Delta x^2)$

Next we use the forward difference operator to estimate the first term in the diffusion equation:

$\frac{\partial T}{\partial t} = \frac{T(t+\Delta t)-T(t)}{\Delta t}$

The second term is expressed using the estimation of the second order partial derivative:

$\frac{\partial^2 T}{\partial x^2} = \frac{T(x+\Delta x)-2T(x)+T(x-\Delta x)}{\Delta x^2}$

Now the diffusion equation can be written as

$\frac{T(t+\Delta t)-T(t)}{\Delta t} - D \frac{T(x+\Delta x)-2T(x)+T(x-\Delta x)}{\Delta x^2} = 0$

This is equivalent to:

$T(t+\Delta t) - T(t) - \frac{D\Delta t}{\Delta x^2}(T(x+\Delta x)-2T(x)+T(x-\Delta x)) = 0$

The expression $D\frac{\Delta t}{\Delta x^2}$ is called the diffusion number, denoted here with s:

$s = D\frac{\Delta t}{\Delta x^2}$

FTCS explicit scheme and analytic solution

If we use n to refer to indices in time and j to refer to indices in space, the above equation can be written as

$T[n+1,j] = T[n,j] + s(T[n,j+1]-2T[n,j]+T[n,j-1])$

This is called a forward-in-time, centered-in-space (FTCS) scheme. Its ‘footprint’ looks like this:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_format = 'svg' # display plots in SVG format

plt.figure(figsize=(6,3))
plt.plot([0,2],[0,0],'k')
plt.plot([1,1],[0,1],'k')
plt.plot([0,1,2,1],[0,0,0,1],'ko',markersize=10)
plt.text(1.1,0.1,'T[n,j]')
plt.text(0.1,0.1,'T[n,j-1]')
plt.text(1.1,1.1,'T[n+1,j]')
plt.text(2.1,0.1,'T[n,j+1]')
plt.xlabel('space')
plt.ylabel('time')
plt.axis('equal')
plt.yticks([0.0,1.0],[])
plt.xticks([0.0,1.0],[])
plt.title('FTCS explicit scheme',fontsize=12)
plt.axis([-0.5,2.5,-0.5,1.5]);

Now we are ready to write the code that is the solution for exercise 2 in Chapter 2 of Slingerland and Kump (2011). This is an example where the one-dimensional diffusion equation is applied to viscous flow of a Newtonian fluid adjacent to a solid wall. If the wall starts moving with a velocity of 10 m/s, and the flow is assumed to be laminar, the velocity profile of the fluid is described by the equation

$\frac{\partial V}{\partial t} - \nu \frac{\partial^2 V}{\partial y^2} = 0$

where $\nu$ is the kinematic viscosity of the fluid. We want to figure out how the velocity will change through time as a function of distance from the wall. [Note that I have changed the original 40 m/s to 10 m/s — the former seems like an unnaturally large velocity to me].

We can compare the numerical results with the analytic solution, which is known for this problem:

$V = V_0 \Big\{ \sum\limits_{n=0}^\infty erfc\big(2n\eta_1+\eta\big) - \sum\limits_{n=0}^\infty erfc\big(2(n+1)\eta_1+\eta\big) \Big\}$

where

$\eta_1 = \frac{h}{2\sqrt{\nu t}}$

and

$\eta = \frac{y}{2\sqrt{\nu t}}$

dt = 0.0005 # grid size for time (s)
dy = 0.0005 # grid size for space (m)
viscosity = 2*10**(-4) # kinematic viscosity of oil (m2/s)
y_max = 0.04 # in m
t_max = 1 # total time in s
V0 = 10 # velocity in m/s

# function to calculate velocity profiles based on a
# finite difference approximation to the 1D diffusion
# equation and the FTCS scheme:
def diffusion_FTCS(dt,dy,t_max,y_max,viscosity,V0):
# diffusion number (has to be less than 0.5 for the
# solution to be stable):
s = viscosity*dt/dy**2
y = np.arange(0,y_max+dy,dy)
t = np.arange(0,t_max+dt,dt)
r = len(t)
c = len(y)
V = np.zeros([r,c])
V[:,0] = V0
for n in range(0,r-1): # time
for j in range(1,c-1): # space
V[n+1,j] = V[n,j] + s*(V[n,j-1] -
2*V[n,j] + V[n,j+1])
return y,V,r,s
# note that this can be written without the for-loop
# in space, but it is easier to read it this way

from scipy.special import erfc

# function to calculate velocity profiles
# using the analytic solution:
def diffusion_analytic(t,h,V0,dy,viscosity):
y = np.arange(0,h+dy,dy)
eta1 = h/(2*(t*viscosity)**0.5)
eta = y/(2*(t*viscosity)**0.5)
sum1 = 0
sum2 = 0
for n in range(0,1000):
sum1 = sum1 + erfc(2*n*eta1+eta)
sum2 = sum2 + erfc(2*(n+1)*eta1-eta)
V_analytic = V0*(sum1-sum2)
return V_analytic

y,V,r,s = diffusion_FTCS(dt,dy,t_max,y_max,viscosity,V0)

# plotting:
plt.figure(figsize=(7,5))
plot_times = np.arange(0.2,1.0,0.1)
for t in plot_times:
plt.plot(y,V[t/dt,:],'Gray',label='numerical')
V_analytic = diffusion_analytic(t,0.04,40,dy,viscosity)
plt.plot(y,V_analytic,'ok',label='analytic',
markersize=3)
if t==0.2:
plt.legend(fontsize=12)
plt.xlabel('distance from wall (m)',fontsize=12)
plt.ylabel('velocity (m/s)',fontsize=12)
plt.axis([0,y_max,0,V0])
plt.title('comparison between explicit numerical
\n(FTCS scheme) and analytic solutions');


The dots (analytic solution) overlap pretty well with the lines (numerical solution). However, this would not be the case if we changed the discretization so that the diffusion number was larger. Let’s look at the stability of the FTCS numerical scheme, by computing the solution with different diffusion numbers. It turns out that the diffusion number s has to be less than 0.5 for the FTCS scheme to remain stable. What follows is a reproduction of Figure 2.7 in Slingerland and Kump (2011):

dt = 0.0005 # grid size for time (m)
dy = 0.0005 # grid size for space (s)
y,V,r,s = diffusion_FTCS(dt,dy,t_max,y_max,viscosity,V0)
V_analytic = diffusion_analytic(0.5,0.04,V0,dy,viscosity)
plt.figure(figsize=(7,5))
plt.plot(y,V_analytic-V[0.5/dt],'--k',label='small s')

dy = 0.0010
dt = 0.00254
y,V,r,s = diffusion_FTCS(dt,dy,t_max,y_max,viscosity,V0)
V_analytic = diffusion_analytic(0.5,0.04,V0,dy,viscosity)
V_numeric = V[r/2-1,:]

plt.plot(y,V_analytic-V_numeric,'k',label='large s')
plt.xlabel('distance from wall (m)',fontsize=12)
plt.ylabel('velocity difference (m/s)',fontsize=12)
plt.title('difference between numerical and analytic
\n solutions for different \'s\' values',fontsize=14)
plt.axis([0,y_max,-4,4])
plt.legend();


Laasonen implicit scheme

plt.figure(figsize=(6,3))
plt.plot([0,2],[1,1],'k')
plt.plot([1,1],[0,1],'k')
plt.plot([0,1,2,1],[1,1,1,0],'ko',markersize=10)
plt.text(1.1,0.1,'T[n,j]')
plt.text(0.1,1.1,'T[n+1,j-1]')
plt.text(1.1,1.1,'T[n+1,j]')
plt.text(2.1,1.1,'T[n+1,j+1]')
plt.xlabel('space')
plt.ylabel('time')
plt.axis('equal')
plt.yticks([0.0,1.0],[])
plt.xticks([0.0,1.0],[])
plt.title('Laasonen scheme',fontsize=12)
plt.axis([-0.5,2.5,-0.5,1.5]);


Instead of estimating the velocity at time step n+1 with the curvature calculated at time step n, as it is done in the FTCS explicit scheme, we can also estimate the curvature at time step n+1, using the velocity change from time step n to time step n+1:

$s\big(T(x+\Delta x)-2T(x)+T(x-\Delta x)\big) = T(t+\Delta t)-T(t)$

Written in matrix notation, this is equiavlent to

$s\big(T[n+1,j+1]-2T[n+1,j]+T[n+1,j-1]\big) = T[n+1,j]-T[n,j]$

After some reshuffling we get

$-sT[n+1,j+1] + (1+2s)T[n+1,j] - sT[n+1,j-1] = T[n,j]$

This is the Laasonen fully implicit scheme. Unlike the FTCS scheme, the Laasonen scheme is unconditionally stable. Let’s try to write some Python code that implements this scheme. First it is useful for me to go through the logic of constructing the system of equations that needs to be solved. Let’s consider a grid that only consists of 5 nodes in space and we are going to estimate the values of T at the locations marked by the red dots in the figure below. Black dots mark the locations where we already know the values of T (from the initial and boundary conditions).

plt.figure(figsize=(6,3))
plt.plot([0,4,4,0],[0,0,1,1],'k')
for i in range(0,4):
plt.plot([i,i],[0,1],'k')
plt.plot([0,1,2,3,4,0,4],[0,0,0,0,0,1,1],'ko',markersize=10)
plt.plot([1,2,3],[1,1,1],'ro',markersize=10)
for i in range(0,5):
plt.text(i+0.1,0.1,'T[0,'+str(i)+']')
plt.text(i+0.1,1.1,'T[1,'+str(i)+']')
plt.xlabel('space')
plt.ylabel('time')
plt.axis('equal')
plt.yticks([0.0,1.0],['0','1'])
plt.title('first two time steps on a 1D grid of five points',fontsize=12)
plt.axis([-0.5,4.8,-0.5,1.5]);


First we write the equations using the Laasonen scheme centered on the three points of unknown velocity (or temperature) — these are the red dots in the figure above:

$\begin{array}{rrrrrcl} -sT[1,0]&+(1+2s)T[1,1]&-sT[1,2]&+0T[1,3]&+0T[1,4]&=&T[0,1] \\ 0T[1,0]&-sT[1,1]&+(1+2s)T[1,2]&-sT[1,3]&+0T[1,4]&=&T[0,2] \\ 0T[1,0]&+0T[1,1]&-sT[1,2]&+(1+2s)T[1,3]&-sT[1,4]&=&T[0,3] \end{array}$

It may seem like we have five unknowns and only three equations but T[1,0] and T[1,4] are on the boundaries and they are known. Let’s rearrange the equation system so that the left hand side has ony the unknowns:

$\begin{array}{rrrrcrr} (1+2s)T[1,1]&-sT[1,2]&+0T[1,3]&=&T[0,1]&+sT[1,0] \\ -sT[1,1]&+(1+2s)T[1,2]&-sT[1,3]&=&T[0,2]& \\ 0T[1,1]&-sT[1,2]&+(1+2s)T[1,3]&=&T[0,3]&+sT[1,4] \end{array}$

In matrix form this is equivalent to

$\begin{bmatrix} 1+2s & -s & 0 \\ -s & 1+2s & -s \\ 0 & -s & 1+2s \end{bmatrix} \times \left[ \begin{array}{c} T[1,1] \\ T[1,2] \\ T[1,3] \end{array} \right] = \left[ \begin{array}{c} T[0,1]+sT[1,0] \\ T[0,2] \\ T[0,3]+sT[1,4] \end{array} \right]$

This of course can be extended to larger dimensions than shown here.
Now we are ready to write the code for the Laasonen scheme. One important difference relative to what I did in the explicit scheme example is that in this case we only keep the last two versions of the velocity distribution in memory, as opposed to preallocating the full array of nt x ny size as we did before. This difference is not a significant time saver for simple problems like this but once you start dealing with more complicated tasks and code it is not possible and/or practical to keep the results of all time steps in memory.

from scipy.sparse import diags
def diffusion_Laasonen(dt,dy,t_max,y_max,viscosity,V0,V1):
s = viscosity*dt/dy**2  # diffusion number
y = np.arange(0,y_max+dy,dy)
t = np.arange(0,t_max+dt,dt)
nt = len(t) # number of time steps
ny = len(y) # number of dy steps
V = np.zeros((ny,)) # initial condition
V[0] = V0 # boundary condition on left side
V[-1] = V1 # boundary condition on right side
# create coefficient matrix:
A = diags([-s, 1+2*s, -s], [-1, 0, 1], shape=(ny-2, ny-2)).toarray()
for n in range(nt): # time is going from second time step to last
Vn = V #.copy()
B = Vn[1:-1] # create matrix of knowns on the RHS of the equation
B[0] = B[0]+s*V0
B[-1] = B[-1]+s*V1
V[1:-1] = np.linalg.solve(A,B) # solve the equation using numpy
return y,t,V,s


Because this is a stable scheme, it is possible to get reasonable solutions with relatively large time steps (which was not possible with the FTCS scheme):

dt = 0.01 # grid size for time (s)
dy = 0.0005 # grid size for space (m)
viscosity = 2*10**(-4) # kinematic viscosity of oil (m2/s)
y_max = 0.04 # in m
V0 = 10.0 # velocity in m/s
V1 = 0.0 # velocity in m/s

plt.figure(figsize=(7,5))
for time in np.linspace(0,1.0,10):
y,t,V,s = diffusion_Laasonen(dt,dy,time,y_max,viscosity,V0,V1)
plt.plot(y,V,'k')
plt.xlabel('distance from wall (m)',fontsize=12)
plt.ylabel('velocity (m/s)',fontsize=12)
plt.axis([0,y_max,0,V0])
plt.title('Laasonen implicit scheme',fontsize=14);


Just for fun, let’s see what happens if we set in motion the right side of the domain as well; that is, set V1 to a non-zero value:

dt = 0.01 # grid size for time (s)
dy = 0.0005 # grid size for space (m)
viscosity = 2*10**(-4) # kinematic viscosity of oil (m2/s)
y_max = 0.04 # in m
V0 = 10.0 # velocity in m/s
V1 = 5.0 # velocity in m/s

plt.figure(figsize=(7,5))
for time in np.linspace(0,1.0,10):
y,t,V,s = diffusion_Laasonen(dt,dy,time,y_max,viscosity,V0,V1)
plt.plot(y,V,'k')
plt.xlabel('distance from wall (m)',fontsize=12)
plt.ylabel('velocity (m/s)',fontsize=12)
plt.axis([0,y_max,0,V0])
plt.title('Laasonen implicit scheme',fontsize=14);


Crank-Nicolson scheme

The Crank-Nicholson scheme is based on the idea that the forward-in-time approximation of the time derivative is estimating the derivative at the halfway point between times n and n+1, therefore the curvature of space should be estimated there as well. The ‘footprint’ of the scheme looks like this:

plt.figure(figsize=(6,3))
plt.plot([0,2],[0,0],'k')
plt.plot([0,2],[1,1],'k')
plt.plot([1,1],[0,1],'k')
plt.plot([0,1,2,0,1,2],[0,0,0,1,1,1],'ko',markersize=10)
plt.text(0.1,0.1,'T[n,j-1]')
plt.text(1.1,0.1,'T[n,j]')
plt.text(2.1,0.1,'T[n,j+1]')
plt.text(0.1,1.1,'T[n+1,j-1]')
plt.text(1.1,1.1,'T[n+1,j]')
plt.text(2.1,1.1,'T[n+1,j+1]')
plt.xlabel('space')
plt.ylabel('time')
plt.axis('equal')
plt.yticks([0.0,1.0],[])
plt.xticks([0.0,1.0],[])
plt.title('Crank-Nicolson scheme',fontsize=12)
plt.axis([-0.5,2.5,-0.5,1.5]);


The curvature at the halfway point can be estimated through averaging the curvatures that are calculated at n and n+1:

$0.5s\big(T[n+1,j+1]-2T[n+1,j]+T[n+1,j-1]\big) + 0.5s\big(T[n,j+1]-2T[n,j]+T[n,j-1]\big) = T[n+1,j]-T[n,j]$

This can be rearranged so that terms at n+1 are on the left hand side:

$-0.5sT[n+1,j-1]+(1+s)T[n+1,j]-0.5sT[n+1,j+1] = 0.5sT[n,j-1]+(1-s)T[n,j]+0.5sT[n,j+1]$

Just like we did for the Laasonen scheme, we can write the equations for the first two time steps:

$\begin{array}{rrrrrcl} -0.5sT[1,0] & +(1+s)T[1,1] & -0.5sT[1,2] & = & 0.5sT[0,0] & +(1-s)T[0,1] & +0.5sT[0,2] \\ -0.5sT[1,1] & +(1+s)T[1,2] & -0.5sT[1,3] & = & 0.5sT[0,1] & +(1-s)T[0,2] & +0.5sT[0,3] \\ -0.5sT[1,2] & +(1+s)T[1,3] & -0.5sT[1,4] & = & 0.5sT[0,2] & +(1-s)T[0,3] & +0.5sT[0,4] \end{array}$

Writing this in matrix form, with all the unknowns on the LHS:

$\begin{bmatrix} 1+s & -0.5s & 0 \\ -0.5s & 1+s & -0.5s \\ 0 & -0.5s & 1+s \end{bmatrix} \times \left[ \begin{array}{c} T[1,1] \\ T[1,2] \\ T[1,3] \end{array} \right] = \begin{bmatrix} 1-s & 0.5s & 0 \\ 0.5s & 1-s & 0.5s \\ 0 & 0.5s & 1-s \end{bmatrix} \times \left[ \begin{array}{c} T[0,1] \\ T[0,2] \\ T[0,3] \end{array} \right] + \left[ \begin{array}{c} 0.5sT[1,0]+0.5sT[0,0] \\ 0 \\ 0.5sT[1,4]+0.5sT[0,4] \end{array} \right]$

Now we can write the code for the Crank-Nicolson scheme. We will use a new input parameter called ntout that determines how many time steps we want to write out to memory. This way you don’t have to re-run the code if you want to plot multiple time steps.

def diffusion_Crank_Nicolson(dy,ny,dt,nt,D,V,ntout):
Vout = [] # list for storing V arrays at certain time steps
V0 = V[0] # boundary condition on left side
V1 = V[-1] # boundary condition on right side
s = D*dt/dy**2  # diffusion number
# create coefficient matrix:
A = diags([-0.5*s, 1+s, -0.5*s], [-1, 0, 1],
shape=(ny-2, ny-2)).toarray()
B1 = diags([0.5*s, 1-s, 0.5*s],[-1, 0, 1], shape=(ny-2, ny-2)).toarray()
for n in range(1,nt): # time is going from second time step to last
Vn = V
B = np.dot(Vn[1:-1],B1)
B[0] = B[0]+0.5*s*(V0+V0)
B[-1] = B[-1]+0.5*s*(V1+V1)
V[1:-1] = np.linalg.solve(A,B)
if n % int(nt/float(ntout)) == 0 or n==nt-1:
Vout.append(V.copy()) # numpy arrays are mutable,
#so we need to write out a copy of V, not V itself
return Vout,s

dt = 0.001 # grid size for time (s)
dy = 0.001 # grid size for space (m)
viscosity = 2*10**(-4) # kinematic viscosity of oil (m2/s)
y_max = 0.04 # in m
y = np.arange(0,y_max+dy,dy)
ny = len(y)
nt = 1000
plt.figure(figsize=(7,5))
V = np.zeros((ny,)) # initial condition
V[0] = 10
Vout,s = diffusion_Crank_Nicolson(dy,ny,dt,nt,viscosity,V,10)

for V in Vout:
plt.plot(y,V,'k')
plt.xlabel('distance from wall (m)',fontsize=12)
plt.ylabel('velocity (m/s)',fontsize=12)
plt.axis([0,y_max,0,V[0]])
plt.title('Crank-Nicolson scheme',fontsize=14);


Fault scarp diffusion

So far we have been using a somewhat artificial (but simple) example to explore numerical methods that can be used to solve the diffusion equation. Next we look at a geomorphologic application: the evolution of a fault scarp through time. Although the idea that convex hillslopes are the result of diffusive processes go back to G. K. Gilbert, it was Culling (1960, in the paper Analytical Theory of Erosion) who first applied the mathematics of the heat equation – that was already well known to physicists at that time – to geomorphology.

Here I used the Crank-Nicolson scheme to model a fault scarp with a vertical offset of 10 m. To compare the numerical results with the analytical solution (which comes from Culling, 1960), I created a function that was written using a Python package for symbolic math called sympy. One of the advantages of sympy is that you can quickly display equations in $\LaTeX$.

import sympy
from sympy import init_printing
init_printing(use_latex=True)
x, t, Y1, a, K = sympy.symbols('x t Y1 a K')
y = (1/2.0)*Y1*(sympy.erf((a-x)/(2*sympy.sqrt(K*t))) + sympy.erf((a+x)/(2*sympy.sqrt(K*t))))
y


The variables in this equation are x – horizontal coordinates, t – time, a – value of x where fault is located, K – diffusion coefficient, Y1 – height of fault scarp.

from sympy.utilities.lambdify import lambdify
# function for analytic solution:
f = lambdify((x, t, Y1, a, K), y)

dt = 2.5 # time step (years)
dy = 0.1 # grid size for space (m)
D = 50E-4 # diffusion coefficient in m2/yr
# e.g., Fernandes and Dietrich, 1997
h = 10 # height of fault scarp in m
y_max = 20 # length of domain in m
t_max = 500 # total time in years
y = np.arange(0,y_max+dy,dy)
ny = len(y)
nt = int(t_max/dt)
V = np.zeros((ny,)) # initial condition
V[:round(ny/2.0)] = h # initial condition

Vout,s = diffusion_Crank_Nicolson(dy,ny,dt,nt,D,V,20)

plt.figure(figsize=(10,5.2))
for V in Vout:
plt.plot(y,V,'gray')

plt.xlabel('distance (m)',fontsize=12)
plt.ylabel('height (m)',fontsize=12)
plt.axis([0,y_max,0,10])
plt.title('fault scarp diffusion',fontsize=14);
plt.plot(y,np.asarray([f(x0, t_max, h, y_max/2.0, D)
for x0 in y]),'r--',linewidth=2);


The numerical and analytic solutions (dashed red line) are very similar in this case (total time = 500 years). Let’s see what happens if we let the fault scarp evolve for a longer time.

dt = 2.5 # time step (years)
dy = 0.1 # grid size for space (m)
D = 50E-4 # diffusion coefficient in m2/yr
h = 10 # height of fault scarp in m
y_max = 20 # length of domain in m
t_max = 5000 # total time in years
y = np.arange(0,y_max+dy,dy)
ny = len(y)
nt = int(t_max/dt)
V = np.zeros((ny,)) # initial condition
V[:round(ny/2.0)] = h # initial condition

Vout,s = diffusion_Crank_Nicolson(dy,ny,dt,nt,D,V,20)

plt.figure(figsize=(10,5.2))
for V in Vout:
plt.plot(y,V,'gray')

plt.xlabel('distance (m)',fontsize=12)
plt.ylabel('height (m)',fontsize=12)
plt.axis([0,y_max,0,10])
plt.title('fault scarp diffusion',fontsize=14);
plt.plot(y,np.asarray([f(x0, t_max, h, y_max/2.0, D)
for x0 in y]),'r--',linewidth=2);


This doesn’t look very good, does it? The reason for the significant mismatch between the numerical and analytic solutions is the fixed nature of the boundary conditions: we keep the elevation at 10 m on the left side and at 0 m on the right side of the domain. There are two ways of getting a correct numerical solution: we either impose boundary conditions that approximate what the system is supposed to do if the elevations were not fixed; or we extend the space domain so that the boundary conditions can be kept fixed throughout the time of interest. Let’s do the latter; all the other parameters are the same as above.

dt = 2.5 # time step (years)
dy = 0.1 # grid size for space (m)
D = 50E-4 # diffusion coefficient in m2/yr
h = 10 # height of fault scarp in m
y_max = 40 # length of domain in m
t_max = 5000 # total time in years
y = np.arange(0,y_max+dy,dy)
ny = len(y)
nt = int(t_max/dt)
V = np.zeros((ny,)) # initial condition
V[:round(ny/2.0)] = h # initial condition

Vout,s = diffusion_Crank_Nicolson(dy,ny,dt,nt,D,V,20)

plt.figure(figsize=(10,5.2))
for V in Vout:
plt.plot(y,V,'gray')

plt.xlabel('distance (m)',fontsize=12)
plt.ylabel('height (m)',fontsize=12)
plt.axis([0,y_max,0,10])
plt.title('fault scarp diffusion',fontsize=14);
plt.plot(y,np.asarray([f(x0, t_max, h, y_max/2.0, D)
for x0 in y]),'r--',linewidth=2);


Now we have a much better result. The vertical dashed lines show the extent of the domain in the previous experiment. We have also gained some insight into choosing boundary conditions and setting up the model domain. It is not uncommon that setting up the initial and boundary conditions is the most time-consuming and difficult part of running a numerical model.

R. Slingerland and L. Kump (2011) Mathematical Modeling of Earth’s Dynamical Systems

W. E. H. Culling (1960) Analytical Theory of Erosion

L. Barba (2013) 12 steps to Navier-Stokes – an excellent introduction to computational fluid dynamics that uses IPython notebooks

I have blogged before about the geosciency aspects of the diffusion equation over here.

You can view and download the IPython Notebook version of this post from Github.

Questions or suggestions? Contact @zzsylvester

# Rivers through time, as seen in Landsat images

Thanks to the Landsat program and Google Earth Engine, it is possible now to explore how the surface of the Earth has been changing through the last thirty years or so. Besides the obvious issues of interest, like changes in vegetation, the spread of cities, and the melting of glaciers, it is also possible to look at how rivers change their courses through time. You have probably already seen the images of the migrating Ucayali River in Peru, for example here. This river is changing its course with an impressive speed; many – probably most – other rivers don’t show much obvious change during the same 30-year period. What determines the meander migration rate of rivers is an interesting question in fluvial geomorphology.

The data that underlies Google Earth Engine is not accessible to everybody, but the Landsat data is available to anyone who creates a free account with Earth Explorer. It is not that difficult (but fairly time consuming) to download a set of images and create animations like this:

This scene also comes from the Ucayali River (you can view it in Google Earth Engine over here) and it is a nice example of how both neck cutoffs and chute cutoffs form. First a neck cutoff takes place that affects the tight bend in the right side of the image; this is followed by a chute cutoff immediately downstream of the neck cutoff location, as the new course of the river happens to align well with a pre-existing chute channel. The third bend in the upper left corner shows some well-developed counter-point-bar deposits. There is one frame in the movie for each year from 1985 to 2013, with a few years missing (due to low quality of the data).

# Memorable moments and photos from 2013

We are well into 2014 now, but it is not too late to look back at 2013 and pick some of the best moments (which means photos in my case) of the year that just passed.

We started out the year with a trip to Zion and Bryce Canyon National Parks. Although it was fairly cold (especially at Bryce Canyon NP — this park has a much higher elevation overall than Zion NP), we had lots of sunshine and did several day hikes. Visiting these parks in the winter is a great idea — they are a lot less crowded than in the summer, and obviously the landscapes and sights are quite different when they are covered with snow.

Zion National Park is paradise for a sedimentologist: there are endless, top-quality exposures of the Navajo Sandstone, showing all kinds of sedimentary structures characteristic of deposits of wind-blown sand. I have included two examples here; you can find more on my Smugmug site.

Sedimentologically, Bryce Canyon National Park is a bit less exciting than Zion, but this is counterbalanced by the fantastic geomorphology of this place. I haven’t seen Bryce Canyon in the summer, but I wouldn’t be surprised if it was more beautiful when it’s covered with snow.

In February, I went on a ‘business’ trip to Torres del Paine National Park in Southern Chile: I attended a field consortium meeting organized by Steve Hubbard’s group at the University of Calgary. I have been to this area several times before, as it has some of the best outcrops of turbidites (= deep-water sediments) in the world, but I was once again shocked how uniquely beautiful Chilean Patagonia can be.

At the end of the official trip, Zane Jobe (who is blogging at Off the Shelf Edge) and I did a bit of geo-turism: we went to see Glacier Grey and Lago Grey, and then did a day hike in the park to check out the actual Torres del Paine. The rest of the photos are here.

In July, my wife and I took a few days to do some hiking and running in Rocky Mountain National Park. I was struggling with a running injury at that time, but the mountains and the trails acted as efficient tranquilizers. More photos at Smugmug.

In September I attended a research conference on turbidity currents in Italy and Peter Talling showed us some of the classic outcrops of the Marnoso-Arenacea Formation. These rocks are very unique because they were deposited by huge submarine flows that covered the entire basin floor. Always wanted to see them and it was enlightening to get up close to them.

Turbidites of the Marnoso-Arenacea Formation, Italian Apennines. David Piper and Bill Arnott for scale

In October we spent a long weekend in Moab, Utah, to participate in our first trail races, but we also did some hiking. Running the Moab Trail Marathon was an amazing experience (I think I will have to do it again this year); unfortunately I didn’t take a camera with me, as I was trying to focus on running (and surviving the race).

Typical view in Canyonlands National Park

To continue with the theme of ‘national parks in winter’, some friends from California and the two of us wrapped up the year with a Christmas trip to Yellowstone and Grand Teton National Parks. More photos, of course, at Smugmug.

# Exploring grain settling with Python

Grain settling is one of the most important problems in sedimentology (and therefore sedimentary geology), as neither sediment transport nor deposition can be understood and modeled without knowing what is the settling velocity of a particle of a certain grain size. Very small grains, when submerged in water, have a mass small enough that they reach a terminal velocity before any turbulence develops. This is true for clay- and silt-sized particles settling in water, and for these grain size classes Stokes’ Law can be used to calculate the settling velocity:

where R = specific submerged gravity (the density difference between the particle and fluid, normalized by fluid density), g = gravitational acceleration, D is the particle diameter, C1 is a constant with a theoretical value of 18, and the greek letter nu is the kinematic viscosity.

For grain sizes coarser than silt, a category that clearly includes a lot of sediment and rock types of great interest to geologists, things get more complicated. The reason for this is the development of a separation wake behind the falling grain; the appearance of this wake results in turbulence and large pressure differences between the front and back of the particle. For large grains – pebbles, cobbles – this effect is so strong that viscous forces become small compared to pressure forces and turbulent drag dominates; the settling velocity can be estimated using the empirical equation

The important point is that, for larger grains, the settling velocity increases more slowly, with the square root of the grain size, as opposed to the square of particle diameter, as in Stokes’ Law.

Sand grains are small enough that viscous forces still play an important role in their subaqueous settling behavior, but large enough that the departure from Stokes’ Law is significant and wake turbulence cannot be ignored. There are several empirical – and fairly complicated – equations that try to bridge this gap; here I focus on the simplest one, published in 2004 in the Journal of Sedimentary Research (Ferguson and Church, 2004):

At small values of D, the left term in the denominator is much larger than the one containing the third power of D, and the equation is equivalent of Stokes’ Law. At large values of D, the second term dominates and the settling velocity converges to the solution of the turbulent drag equation.

But the point of this blog post is not to give a summary of the Ferguson and Church paper; what I am interested in is to write some simple code and plot settling velocity against grain size to better understand these relationships through exploring them graphically. So what follows is a series of Python code snippets, directly followed by the plots that you can generate if you run the code yourself. I have done this using the IPyhton notebook, a very nice tool that allows and promotes note taking, coding, and plotting within one document. I am not going to get into details of Python programming and the usage of IPyhton notebook, but you can check them out here.

First we have to implement the three equations as Python functions:

import numpy as np
import matplotlib.pyplot as plt
rop = 2650.0 # density of particle in kg/m3
rof = 1000.0 # density of water in kg/m3
visc = 1.002*1E-3 # dynamic viscosity in Pa*s at 20 C
C1 = 18 # constant in Ferguson-Church equation
C2 = 1 # constant in Ferguson-Church equation
def v_stokes(rop,rof,d,visc,C1):
R = (rop-rof)/rof # submerged specific gravity
w = R*9.81*(d**2)/(C1*visc/rof)
return w
def v_turbulent(rop,rof,d,visc,C2):
R = (rop-rof)/rof
w = (4*R*9.81*d/(3*C2))**0.5
return w
def v_ferg(rop,rof,d,visc,C1,C2):
R = (rop-rof)/rof
w = ((R*9.81*d**2)/(C1*visc/rof+
(0.75*C2*R*9.81*d**3)**0.5))
return w


Let’s plot these equations for a range of particle diameters:

d = np.arange(0,0.0005,0.000001)
ws = v_stokes(rop,rof,d,visc,C1)
wt = v_turbulent(rop,rof,d,visc,C2)
wf = v_ferg(rop,rof,d,visc,C1,C2)
figure(figsize=(10,8))
plot(d*1000,ws,label='Stokes',linewidth=3)
plot(d*1000,wt,'g',label='Turbulent',linewidth=3)
plot(d*1000,wf,'r',label='Ferguson-Church',linewidth=3)
plot([0.25, 0.25],[0, 0.15],'k--')
plot([0.25/2, 0.25/2],[0, 0.15],'k--')
plot([0.25/4, 0.25/4],[0, 0.15],'k--')
text(0.36, 0.11, 'medium sand', fontsize=13)
text(0.16, 0.11, 'fine sand', fontsize=13)
text(0.075, 0.11, 'v. fine', fontsize=13)
text(0.08, 0.105, 'sand', fontsize=13)
text(0.01, 0.11, 'silt and', fontsize=13)
text(0.019, 0.105, 'clay', fontsize=13)
legend(loc=2)
xlabel('grain diameter (mm)',fontsize=15)
ylabel('settling velocity (m/s)',fontsize=15)
axis([0,0.5,0,0.15]);
D = [0.068, 0.081, 0.096, 0.115, 0.136, 0.273,
0.386, 0.55, 0.77, 1.09, 2.18, 4.36]
w = [0.00425, 0.0060, 0.0075, 0.0110, 0.0139, 0.0388,
0.0551, 0.0729, 0.0930, 0.141, 0.209, 0.307]
scatter(D,w,50,color='k')
show()


The black dots are data points from settling experiments performed with natural river sands (Table 2 in Ferguson and Church, 2004). It is obvious that the departure from Stokes’ Law is already significant for very fine sand and Stokes settling is completely inadequate for describing the settling of medium sand.

This plot only captures particle sizes finer than medium sand; let’s see what happens as we move to coarser sediment. A log-log plot is much better for this purpose.

d = np.arange(0,0.01,0.00001)
ws = v_stokes(rop,rof,d,visc,C1)
wt = v_turbulent(rop,rof,d,visc,C2)
wf = v_ferg(rop,rof,d,visc,C1,C2)
figure(figsize=(10,8))
loglog(d*1000,ws,label='Stokes',linewidth=3)
loglog(d*1000,wt,'g',label='Turbulent',linewidth=3)
loglog(d*1000,wf,'r',label='Ferguson-Church',linewidth=3)
plot([1.0/64, 1.0/64],[0.00001, 10],'k--')
text(0.012, 0.0007, 'fine silt', fontsize=13,
rotation='vertical')
plot([1.0/32, 1.0/32],[0.00001, 10],'k--')
text(0.17/8, 0.0007, 'medium silt', fontsize=13,
rotation='vertical')
plot([1.0/16, 1.0/16],[0.00001, 10],'k--')
text(0.17/4, 0.0007, 'coarse silt', fontsize=13,
rotation='vertical')
plot([1.0/8, 1.0/8],[0.00001, 10],'k--')
text(0.17/2, 0.001, 'very fine sand', fontsize=13,
rotation='vertical')
plot([0.25, 0.25],[0.00001, 10],'k--')
text(0.17, 0.001, 'fine sand', fontsize=13,
rotation='vertical')
plot([0.5, 0.5],[0.00001, 10],'k--')
text(0.33, 0.001, 'medium sand', fontsize=13,
rotation='vertical')
plot([1, 1],[0.00001, 10],'k--')
text(0.7, 0.001, 'coarse sand', fontsize=13,
rotation='vertical')
plot([2, 2],[0.00001, 10],'k--')
text(1.3, 0.001, 'very coarse sand', fontsize=13,
rotation='vertical')
plot([4, 4],[0.00001, 10],'k--')
text(2.7, 0.001, 'granules', fontsize=13,
rotation='vertical')
text(6, 0.001, 'pebbles', fontsize=13,
rotation='vertical')
legend(loc=2)
xlabel('grain diameter (mm)', fontsize=15)
ylabel('settling velocity (m/s)', fontsize=15)
axis([0,10,0,10])
scatter(D,w,50,color='k');
show()


This plot shows how neither Stokes’ Law, nor the velocity based on turbulent drag are valid for calculating settling velocities of sand-size grains in water, whereas the Ferguson-Church equation provides a good fit for natural river sand.

Grain settling is a special case of the more general problem of flow past a sphere. The analysis and plots above are all dimensional, that is, you can quickly check by looking at the plots what is the approximate settling velocity of very fine sand. That is great, but you would have to generate a new plot – and potentially do a new experiment – if you wanted to look at the behavior of particles in some other fluid than water. A more general treatment of the problem involves dimensionless variables; in this case these variables are the Reynolds number and the drag coefficient. The classic diagram for flow past a sphere is a plot of the drag coefficient against the Reynolds number. I will try to reproduce this plot, using settling velocities that come from the three equations above.

At terminal settling velocity, the drag force equals the gravitational force acting on the grain:

We also know that the gravitational force is given by the submerged weight of the grain:

The drag coefficient is essentially a dimensionless version of the drag force:

At terminal settling velocity, the particle Reynolds number is

Using these relationships it is possible to generate the plot of drag coefficient vs. Reynolds number:

d = np.arange(0.000001,0.3,0.00001)
C2 = 0.4 # this constant is 0.4 for spheres, 1 for natural grains
ws = v_stokes(rop,rof,d,visc,C1)
wt = v_turbulent(rop,rof,d,visc,C2)
wf = v_ferg(rop,rof,d,visc,C1,C2)
Fd = (rop-rof)*4/3*pi*((d/2)**3)*9.81 # drag force
Cds = Fd/(rof*ws**2*pi*(d**2)/8) # drag coefficient
Cdt = Fd/(rof*wt**2*pi*(d**2)/8)
Cdf = Fd/(rof*wf**2*pi*(d**2)/8)
Res = rof*ws*d/visc # particle Reynolds number
Ret = rof*wt*d/visc
Ref = rof*wf*d/visc
figure(figsize=(10,8))
loglog(Res,Cds,linewidth=3, label='Stokes')
loglog(Ret,Cdt,linewidth=3, label='Turbulent')
loglog(Ref,Cdf,linewidth=3, label='Ferguson-Church')
# data digitized from Southard textbook, figure 2-2:
Re_exp = [0.04857,0.10055,0.12383,0.15332,0.25681,0.3343,0.62599,0.77049,0.94788,1.05956,
1.62605,2.13654,2.55138,3.18268,4.46959,4.92143,8.02479,12.28672,14.97393,21.33792,
28.3517,34.55246,57.57204,78.3929,96.88149,159.92596,227.64082,287.31738,375.98547,
516.14355,607.03827,695.8316,861.51953,1147.26099,1194.43213,1513.70166,1939.70557,
2511.91235,2461.13232,3106.32397,3845.99561,4974.59424,6471.96875,8135.45166,8910.81543,
11949.91309,17118.62109,21620.08203,28407.60352,36064.10156,46949.58594,62746.32422,
80926.54688,97655.00781,122041.875,157301.8125,206817.7188,266273,346423.5938,302216.5938,
335862.5313,346202,391121.5938,460256.375,575194.4375,729407.625]
Cd_exp = [479.30811,247.18175,199.24072,170.60068,112.62481,80.21341,45.37168,39.89885,34.56996,
28.01445,18.88166,13.80322,12.9089,11.41266,8.35254,7.08445,5.59686,3.92277,3.53845,
2.75253,2.48307,1.99905,1.49187,1.27743,1.1592,0.89056,0.7368,0.75983,0.64756,0.56107,
0.61246,0.5939,0.49308,0.39722,0.48327,0.46639,0.42725,0.37951,0.43157,0.43157,0.40364,
0.3854,0.40577,0.41649,0.46173,0.41013,0.42295,0.43854,0.44086,0.4714,0.45225,0.47362,
0.45682,0.49104,0.46639,0.42725,0.42725,0.40171,0.31214,0.32189,0.20053,0.16249,0.10658,
0.09175,0.09417,0.10601]
loglog(Re_exp, Cd_exp, 'o', markerfacecolor = [0.6, 0.6, 0.6], markersize=8)

# Reynolds number for golf ball:
rof_air = 1.2041 # density of air at 20 degrees C
u = 50 # velocity of golf ball (m/s)
d = 0.043 # diameter of golf ball (m)
visc_air = 1.983e-5 # dynamic viscosity of air at 20 degrees C
Re = rof_air*u*d/visc_air
loglog([Re, Re], [0.4, 2], 'k--')
text(3e4,2.5,'$Re$ for golf ball',fontsize=13)
legend(loc=1)
axis([1e-2,1e6,1e-2,1e4])
xlabel('particle Reynolds number ($Re$)', fontsize=15)
ylabel('drag coefficient ($C_d$)', fontsize=15);


The grey dots are experimental data points digitized from the excellent textbook by John Southard, available through MIT Open Courseware. As turbulence becomes dominant at larger Reynolds numbers, the drag coefficient converges to a constant value (which is equal to C2 in the equations above). Note however the departure of the experimental data from this ideal horizontal line: at high Reynolds numbers there is a sudden drop in drag coefficient as the laminar boundary layer becomes turbulent and the flow separation around the particle is delayed, that is, pushed toward the back; the separation wake becomes smaller and the turbulent drag decreases. Golf balls are not big enough to reach this point without some additional ‘help'; this help comes from the dimples on the surface of the ball that make the boundary layer turbulent and reduce the wake.

You can view and download the IPython notebook version of this post from the IPython notebook viewer site.

References
Ferguson, R. and Church, M. (2004) A simple universal equation for grain settling velocity. Journal of Sedimentary Research 74, 933–937.

# My job description, in simple words

I couldn’t refrain from playing with the ‘upgoerfive’ tool, an online text editor that only allows you to use the one thousand most common words of the English language. Here is my attempt to describe what I do.

I study how small pieces of rock get together, move along, and settle down to form beds. Water usually moves smaller pieces further than large ones but lots of small and big pieces together can move really fast and build thick beds on the water bottom. I also look at the form of these beds and think about how much space is left in between the small pieces. This space is filled with water and sometimes with other stuff that people like to burn in their cars. It is nice that such beds are often seen in beautiful places around the world that I like to visit.

Here is the link to the original upgoerfive post. You can try it out yourself over here.

# The case for scales, rates, and numbers in geology

A long time ago I used to study geology in a nice city called Cluj, in the middle of that interesting part of Romania known as Transylvania. This was a place and time where and when I learned about quartz, feldspar, species of coral and foraminifera in great detail, heard about sequence stratigraphy and turbidites for the first time, and went on some great geological field trips. Not to mention the half-liter bottles of beer that would be significant components of any decent geological trip or spontaneous philosophical discussion in the evening. The less pleasant part was that many of the classes I took involved brute-force memorization of fossils, minerals, chronostratigraphic names, and formations. Although the geological vocabulary that I picked up was pretty broad and proved useful as a good set of words, terms, and definitions to play with, I forgot many of the details by now. If you asked me what was the difference between granite and granodiorite, I would have to check. I don’t remember at all what fossils are characteristic of the Late Jurassic. And, despite doing some fieldwork myself over there, I cannot remember the stratigraphic nomenclature in Transylvanian Basin; I would have to look it up (probably in this paper).

After college, it took me about one year to realize with convincing clarity that there was a lot left to learn. I went on to grad school on the other side of the planet, at a well-known university. Many of the classes I took over there were – unsurprisingly – quite different; a lot more focus on laws, processes and the connections between geological things than the ‘things’ themselves. It was also there that I started to see the links between geology and physics and math. I picked up quite a bit of math and physics during high school, but then quickly relegated them to the status of “stuff that is rarely used in geology”. At grad school, it dawned on me that numbers and mathematical laws are not only useful in geology, but are in fact necessary for doing good earth science. Maybe I am stating the obvious, but here it goes anyway: geology deals with enormous variations in scale, both in space and time; and it is not enough to say that the river was deep (how deep?), the tectonic deformation was fast (how fast?), the sea-level highstand lasted long (how long?), or the sediment gravity flows were high-energy flows (I am not even sure what that means). One of the most important things I learned was an appreciation for physical and quantitative insight in geology, that is, having at least an idea, a feel for what are the scales and rates involved in the formation of the rocks you are looking at. I cannot say it better than Chris Paola, one of the important and influential advocates of moving sedimentary geology closer to physics and math:

“For the ‘modal’ sedimentary-geology student, it is not sophisticated computational skills or training in advanced calculus that is lacking, but rather the routine application of basic quantitative reasoning. This means things like estimating scales and rates for key processes, knowing the magnitudes of basic physical properties, and being able to estimate the relative importance of various processes in a particular setting. Understanding scales, rates and relative magnitudes is to quantitative science what recognizing quartz and feldspar is to field geology. Neither requires years of sophisticated training, but both require repetition until they become habitual.”

Developing these skills is a lot easier if one is not afraid of tinkering with simple computer programs. Want to really understand what Stokes’ law is about? There is no better way than typing the equation into an Excel spreadsheet or a Matlab m-file and see how the plot of settling velocity against grain size looks like. What about settling in a fluid with different viscosity? Change the variable, and compare the result with the previous curve. High-level programming languages like Matlab or Python* are a lot easier to learn than languages closer to ‘computerese’ and farther from English, and they are great tools for these kinds of exercises and experiments. As somebody interested in stratigraphic architecture, I have become especially fond of creating surfaces that vaguely resemble real-world landscapes and then see how the evolution of these surfaces through time – deposition over here, erosion over there – creates stratigraphy. Complex three dimensional geometry is a lot easier to grasp if you can visualize and dissect it on the computer screen.

Of course, numbers, diagrams and images that come from computer programs are only useful if they demonstrably say something about the real world. Data collection in the field and the laboratory are equally important. But nowadays we often have more data than we wished for, and quantitative skills come handy for visualizing and analyzing large datasets – and comparing them to models.

Not everyone is excited about the growing number of earth scientists who tend to see equations ‘in the rocks’. The logo of the Sedimentology Research Group at the University of Minnesota features the Exner equation carved into a pebble, allegedly as a response to the exclamation “I haven’t seen yet an equation written on the rocks!” There is some concern that many geology graduates nowadays do not get to see, to map and to sample enough real rocks and sediments in the field. Although I think this unease is not entirely unsubstantiated, I wouldn’t want to sound as pessimistic as Emiliano Mutti – one of the founding fathers of deepwater sedimentology – does in the last phrase of a review article:

“This approach raises a problem, and not a small one: in connection with data collection in the field, how many field geologists are being produced in these times of increasingly computerized geology; and how good are they?”

As far as I know, geological field work is still an important part of the curriculum in many departments of geology – as it clearly should be. The number one reason I have become a geologist was that I loved mountains, hiking, and being outdoors in general, way before I started formally studying geology. And I still take every opportunity to go to the field. But I cannot see the growth of “computerized geology” – and of quantitative geology in general – as a bad thing. Does dry quantification take away the beauty and poetry of geology? I don’t think so. Unweaving the rainbow, unfolding a mountain, and reconstructing a turbidity current only add to our appreciation of the scale and grandeur of geology.

* I will let you know later whether this is true about Python…

** I have started writing this post for Accretionary Wedge #38, mostly because I found the call for posts quite inspiring, but haven’t finished it in time. Read all the good stuff at Highly Allochthonous.

# Salt and sediment: A brief history of ideas

Salty weirdness
Salt is a weird kind of rock. At first sight, it behaves like most other rocks: if you pick up a piece, it is hard, it is heavy, and it breaks if hit with a hammer. But put it under stress for thousands of years, and salt will behave like a fluid: relatively small forces can cause it to flow toward less stressful surroundings. This often means it will try to find its way to the surface.

When deposited, sand and mud have lots of pore space filled with water and have relatively low density. However, as they get buried by more sediment, much pore space is lost, both through compaction and cementation. Sediments turn into sedimentary rocks, become harder, and their density increases. In contrast, salt doesn’t have much pore space to begin with; its density will stay the same, regardless of depth of burial. As both salt and sediment are buried to greater depths, an unstable condition develops: lighter salt lying under denser material. In addition, the location of the salt layer in the sediment column is not entirely random: it is in the nature of sedimentary basins to initially place salt at the bottom of the sediment pile. Extensive salt layers usually form early in a basin’s lifetime, when seawaters invade for the first time shallow depressions on a continent that is about to split into two along a rift zone. The Dead Sea is an obvious example that comes to mind.

Layering salt and sediment in this unstable order is a recipe for a spectacular geological show. As salt is trying to find its way to the surface, it forms drop-shaped blobs called diapirs; but also ridges, walls, and salt sheets. Several sheets can connect laterally into a huge salt canopy, a new salt layer that is entirely out-of-place or allochtonous. Salt can also act as a lubricating layer at the base of a thick sequence of sedimentary rocks. But I am rushing ahead a little bit; salt tectonics is such a new – but rapidly growing – science that salt canopies, despite their widespread presence in the subsurface Gulf of Mexico, were not recognized and described until the 1980s.

Tectonics vs. buoyancy, Europe vs. America
Before the beginning of the twentieth century, even with the role that salt played in human history, little was known about how salt domes formed. This was an age of rampant speculation; surface data was scarce because salt does not last very long after exposed as it quickly gets dissolved and washed away by precipitation. Many geologists thought that formation of salt domes didn’t require any significant salt deformation or displacement. But things have changed dramatically in 1901, with the discovery of the Spindletop oil field on top of a salt dome in southeastern Texas. The recognition that oil is often found on top of and around salt domes created a much stronger interest in understanding how exactly salt formations are put in place.

European geologists thought that the main driving force was compression, the force that causes folding and thrusting and builds mountains. In Romania, where the Eastern Carpathians take a sharp turn toward the southwest, salt was found in the cores of oil-bearing anticlines. The contacts with the surrounding rocks were clearly discordant. These are the structures that prompted Ludovic Mrazec, professor of geology at University of Bucharest, to coin the term “diapir” in 1907.

 Mrazec’s explanation of how salt diapirs form. From Barton (1925).

Salt in Germany and Poland also seemed to occur invariably in a compressional setting, in the cores of folds, next to folds that had no salt associated. It seemed obvious that salt was ‘pushed up’ by tectonic forces, and it appeared unlikely that the rise of salt itself was causing the folding.

But the discovery of a multitude of salt diapirs in the Gulf of Mexico made it clear that they can occur far away from any mountains and compressive tectonic forces. The much simpler setting and relative lack of deformation in the Gulf proved informative. “The Roumanian salt-dome geologist possibly may have more to learn from the American salt domes than the American salt-dome geologist has to learn from the Roumanian domes. The occurrence of the American domes in a region of tectonic quiescence suggests that tectonic thrust cannot have the importance postulated by Mrazec” – wrote Donald Barton in 1925.

This was also the time when the density difference between salt and sediment came into discussion. Gravity measurements in the Gulf of Mexico showed anomalies above salt domes that were due to the lower density of salt. It was increasingly recognized that density inversion must play an important role in diapirism, especially where compressive tectonic forces were absent. In addition, by the 1930s geologists have reached a consensus that salt diapirs must somehow punch through the overlying sediment. They seemed to ignore the fact that, as Wade (1931) put it, you cannot drive a putty nail through a wooden board. As mentioned before, salt does behave like a fluid over geological time scales. But how can it penetrate thick layers of hardened sedimentary rock?

A brilliant idea: downbuilding
The solution to this problem came in 1933, from the same Donald Barton who was discussing the differences between European and American salt domes in 1925. He suggested that diapirs can form without much piercement of the sediment above. Instead, once a small dome is initiated, it simply can stay in place, always at or close to the surface, while sediment is deposited around it and the source salt layer subsides: “it is the sediments which move, and not the salt core. The energy requirement (…) is very much less than if there were actual upward movement of the salt.”

 The evolution of salt diapirs through ‘downbuilding’. Salt domes are always close to the surface and diapirism goes hand-in-hand with sedimentation. From Barton (1933).

This was a key insight: it got rid of the “room problem”, the need for moving huge volumes of hard rock out of the way of the rising salt. It also highlighted that salt movement can happen at the same time with sedimentation, a fact that became abundantly obvious later as high-quality seismic data became available. But the concept of ‘downbuilding’ was ignored for the next fifty years.

Animation showing how downbuilding works. Blue represents salt, yellow is sediment. To mimic mass balance for salt (-- what is lost from the source layer must go into the salt dome), the blue area is kept constant through the animation.

The beauty of instabilities

The main reason for conveniently forgetting Barton’s idea was that density inversion between two fluids could be nicely studied in the lab and described with elegant equations. In one of the papers that kicked off this fascination with Rayleigh-Taylor instabilities, Nettleton (1934) used corn syrup and less dense crude oil to visualize diapir-like blobs of fluid in a transparent cylinder and to show that gravity alone, without any help from contractional forces, was enough to generate structures similar to salt domes.

 Less dense crude oil (black) forming diapir-like blobs as rising through higher-density corn syrup (yellow). Redrawn from Nettleton (1934).

One problem with this approach was that oil and syrup can be photographed during deformation, but the transient structures could not be carefully dissected and analyzed later. Materials of higher viscosity were needed for that; however, increasing the viscosity resulted in a density difference too small to get the fluids moving in the first place. The trick was to place the whole experiment in a centrifuge and use the centrifugal force to imitate a larger-than-normal gravitational force. This approach formed the basis of a productive line of research on gravity tectonics in the laboratory of the Norwegian-Swedish geologist Hans Ramberg. The results are probably more relevant to what is happening deeper in the Earth, at higher temperatures and pressures, where most rocks become more similar in behavior to salt.

Modern salt tectonics
By the late 1980s it has become quite obvious that kilometer-thick piles of sedimentary rock cannot be treated as fluids and salt-sediment interaction is more similar to placing and deforming slabs of brittle material on top of a viscous fluid. Seismic from salt-bearing sedimentary basins suggested that the history of salt movement and sedimentation were highly interconnected and Barton’s downbuilding concept was strongly relevant.

Three-dimensional seismic data also showed the variety and complexity of allochtonous salt bodies in salt-rich sedimentary basins. Sandbox experiments with more realistic material properties and ongoing sedimentation during deformation were performed and the results beautifully visualized. The behavior of turbidity currents flowing over complex salt-related submarine topography was investigated. Hundreds of scientific papers were written on salt tectonics, both by industry geoscientists and researchers in the academia.

 N-S cross section in the Gulf of Mexico. Large volumes of the Jurassic Louann salt have been displaced and squeezed into a salt canopy surrounded by much younger sediments. From Pilcher et al., 2011

And there is quite a bit left to explore and understand.

Barton, D. C. (1926) The American Salt-Dome Problems in the Light of the Roumanian and German Salt Domes, AAPG Bulletin, v. 9, p. 1227–1268.

Barton, D. C. (1933) Mechanics of Formation of Salt Domes with Special Reference to Gulf Coast Salt Domes of Texas and Louisiana, AAPG Bulletin, v. 17, 1025–1083.

Hudec, M., & Jackson, M. (2007) Terra infirma: Understanding salt tectonics. Earth Science Reviews, 82(1-2), 1–28.

Jackson, M. (1996) Retrospective salt tectonics, in M.P.A. Jackson, D.G. Roberts, and S. Snelson, eds., Salt tectonics: a global perspective: AAPG Memoir 65, p. 1–28. [great summary of the history of salt tectonics]

Mrazec, L. (1907) Despre cute cu sȋmbure de străpungere [On folds with piercing cores]: Bul. Soc. Stiint., Romania, v. 16, p. 6–8.

Nettleton, L. L. (1934) Fluid Mechanics of Salt Domes, AAPG Bulletin, v. 18, p. 1–30.

Pilcher, R. S., Kilsdonk, B., & Trude, J. (2011) Primary basins and their boundaries in the deep-water northern Gulf of Mexico: Origin, trap types, and petroleum system implications. AAPG Bulletin, v. 95(2), p. 219–240.

Wade, A. (1931) Intrusive salt bodies in coastal Asir, south western Arabia: Institute of Petroleum Technologists Journal, v. 17, p. 321–330, 357–361.

# Where on Google Earth? WoGE #296

Hindered Settling hasn’t hosted a Where-on-Google-Earth in a long time, but WoGE #295 (hosted at Andiwhere’s) had such a range of colors and geological features that I couldn’t refrain from looking for it and, once found it, had to post the solution. So, after a short break in the game (busy week!) here is WoGE #296 — the rules of the game are nicely described over here. I invoke the Schott rule. Posting time is July 8, 2011, 14:00 UTC.

 click on image for larger view

# Stretching the truth: vertical exaggeration of seismic data

If someone showed a photograph of the famous Cuernos massif (Torres del Paine National Park, Chile) like the one below, it would be – probably, hopefully – obvious to everybody that something is wrong with the picture. Our eyes and brains have seen enough mountain scenery that we intuitively know how steep is ‘steep’ in alpine landscapes. The peaks in this photograph just look too extreme, too high if one takes into account their lateral extent.

 The Cuernos in Torres del Paine National Park, Chile, vertically exaggerated by a factor of two.

For comparison, here is the original shot:

 The Cuernos, beautiful, without exaggeration.

Yet this kind of vertical stretching of images is the rule rather than the exception when displaying seismic sections, both on computer screens and in scientific papers. There are two main reasons for this: first, subsurface geometries are often more obvious when vertically stretched and slopes are larger than in real life. Second, more often than not, we have no precise knowledge of the actual vertical scale. Raw seismic reflection data is recorded in time: we are measuring how long it takes before a seismic wave propagates down to a discontinuity and comes back to the surface. This time measure is called ‘two-way traveltime'; in order to convert it to depth, knowledge of the velocity of sound through the rocks is needed:

depth = seismic_velocity * two-way_traveltime / 2
The problem is that the velocity of the wave varies as it goes deeper (usually increases with depth as rocks become ‘harder’); and, unless we are looking at perfect layercake stratigraphy (not that common!), it also changes laterally. So, if we want to look at the actual geological structures, without distortions due to varying velocities, we need to do a depth conversion and we need a ‘velocity model’ that roughly describes the spatial distribution of velocities. Precise velocity measurements often come from wells where depth is well known; less precise estimates can be backed out from the seismic recordings themselves, but the solution is often non-unique and multiple iterations are necessary to build a good velocity- and depth model. As a result, seismic reflection data is often interpreted with two-way traveltime on the vertical axis, without depth conversion; and not knowing the true vertical scale makes it easier to use vertical exaggeration with vengeance.

A recent paper, published in Marine and Petroleum Geology, shows that vertical exaggeration of seismic data is indeed very common. Simon Stewart of Heriot-Watt University has looked through 1437 papers published between 2006-2010 and found that 75% of the papers show seismic displays with vertical exaggeration of a factor larger than 2. Only 12% are shown with roughly equal horizontal and vertical scales.

 Histogram of vertical exaggerations in 1437 papers. From Stewart (2011).

One of the effects of vertical exaggeration is the strong steepening of dips. A 10 degree slope at a vertical exaggeration of 10 becomes an almost vertical drop of 60 degrees; it is hard not to think of these exaggerated slopes as steep slopes, even though they are not that abrupt in reality. Depositional geometries often have very small dips and significant vertical exaggeration is needed to illustrate the overall shapes.

The paper suggests that published seismic sections should be labeled with an estimate of the vertical exaggeration, in addition to the usual horizontal and vertical scales [I am guilty myself of not doing this as it should be done]. Even better, one can go further and create several versions of the figure with different vertical exaggerations. The cross section of a submarine lobe deposit below is a fine example of such a display. Showing only the version that was exaggerated vertically 25 times would suggest that this is a deposit at the base of a steep slope; the 1:1 figure at the top brings us back to reality and clearly shows that this morphology and stratigraphy are both extremely flat.

 Dip section of a submarine lobe deposit, offshore Corsica. From Deptuck et al. (2008).

To see more about scales and vertical exaggeration in geology, check out this recent post at Highly Allochthonous; and the nice illustrations that Matt has put together over at Agile*.

References

Deptuck, M., Piper, D., Savoye, B., & Gervais, A. (2008). Dimensions and architecture of late Pleistocene submarine lobes off the northern margin of East Corsica Sedimentology, 55 (4), 869-898 DOI: 10.1111/j.1365-3091.2007.00926.x

Stewart, S. (2011). Vertical exaggeration of reflection seismic data in geoscience publications 2006–2010 Marine and Petroleum Geology, 28 (5), 959-965 DOI: 10.1016/j.marpetgeo.2010.10.003

# Snorkeling and geology in Kealakekua Bay, Big Island, Hawaii

For a long time, I didn’t think it was worth spending more than an hour on a beach, even the most beautiful ones, unless there were some nice cliffs nearby showing some interesting geology. My views in this regard have changed dramatically about three years ago, when I spent a week on The Big Island of Hawaii, and the hotel where we were staying offered free rental of snorkeling gear. I put on the mask and the fins, trying to remember how this was supposed to work (I did a bit of snorkeling in Baja California many years before that), and put my face into the not-too-interesting-looking waters in the front of the hotel.

I was in for a surprise. The water was far from crystal clear, but I could still see fantastic coral creations lined up along the bay and lots of fish of so many colors and patterns that it felt unreal. Until then I thought that this kind of scenery was hard to see unless you were a filmmaker working for Discovery Channel or a marine biologist specializing in tropical biodiversity. The next day I spotted a couple of green turtles frolicking in the water, clearly not bothered by the nearby snorkelers, and I already knew that I needed to look into the possibility of buying a simple underwater camera.

 Lots of coral, mostly belonging to the genera Lobites (lobe coral) and Pocillopora (cauliflower coral)

Three years later I went back to the Big Island with more excitement about tropical beaches, plus bigger plans and a bit more knowledge about snorkeling. After going through a few well-known snorkeling sites on the west coast, like Kahalu’u Beach in Kona and Two Step near Pu’uhonua o Honaunau park, we got on a nice boat (run by a company called Fair Wind – strongly recommended!) and did some snorkeling in Kealakekua Bay.

 Visibility in Kealakekua Bay is usually very good

 Old wrinkles of pahoehoe lava getting encrusted by algae and corals and chewed up by sea urchins

Kealakekua Bay is difficult to reach; there is no road and no parking lot nearby. You either have to hike in, paddle through the bay in a kayak, or take a boat. I have heard before that this was the best snorkeling spot in Hawai’i, but I think that is an understatement. Unlike all the other spots we tried during the last few years in Hawaii (and that includes several beaches on Kauai and Hanauma Bay on Oahu, the presidential snorkeling site), the water at Kealakekua Bay was calm and very clear, with fantastic visibility.

 Heads of cauliflower coral, with yellow tangs for scale

I will not attempt to describe this whole new world; instead I will let the photographs speak for themselves (as always, more photos at Smugmug). Even better, if you go to the Big Island, make sure that you visit this place with some snorkeling gear.

 Yellow tangs (Zebrasoma flavescens) often congregate in large schools and it is difficult to stop taking pictures of them

When I was at Kealakekua Bay, I didn’t know much about the local geology. The big cliff bordering the bay toward the northwest, called Pali Kapu o Keoua (see image above), shows a number of layered lava flows that belong to the western flank of Mauna Loa; and I suspected that this must have been a large fault scarp, but that was the end of my geological insight. A couple of hours worth of research after I got home revealed that Pali Kapu o Keoua was a fault indeed: it is called the Kealakekua Fault and it has been mapped, along with the associated submarine geomorphological features, in the 1970s and 1980s by U.S. Geological Survey geoscientists.  It turns out that one of the shipboard scientists and key contributors to these studies was Bill Normark (see also a post about Bill at Clastic Detritus). While in California in the late 1990s, I was lucky to get to know Bill and have some truly inspiring discussions with him about turbidites, geology, and wine, so this was a doubly valuable little discovery to me.

So what is the origin of the Kealakekua Fault? The Hawaiian Islands are far away from any tectonic plate boundaries, so there is not a lot of opportunity here for inverse or strike-slip faults to develop. However, the Hawaiian volcanoes are humongous mountains and their underwater slopes are extremely steep by submarine slope standards: gradients of 15-10˚ are common. [This is in contrast by the way with the relatively gentle slopes of 3-8˚ the subaerial flanks of the volcanoes, a difference that – it just occurred to me – has to do something with the different thermal conductivities of water and air. Water is ~24 times more efficient at cooling lavas, or anything for that matter, than air, so once a volcano sticks its head out of the water, basaltic lava flows are pretty efficient at carrying volcanic material far away from the crater, thus building gently sloping shield volcanoes. The same flows are promptly solidified and stopped by the cool ocean waters as soon as they reach the coast.] Slopes that are this steep are also unstable; the underwater parts of these volcanoes tend to fail from time to time and large volumes of rock rapidly move to deeper waters as giant submarine landslides. Seafloor mapping around the islands revealed that the underwater topography is far from smooth; instead, in many places it consists of huge slide and slump blocks.

 Topographic map of the Big Island. Note the location of Kealakekua Fault and the rugged seafloor to the southwest of it, marking the area affected by slides and slumps. This is a map based on higher-resolution bathymetric data collected during a collaborative effort led by JAMSTEC (Japan Marine Science and Technology Center). Source: U.S. Geological Survey Geologic Investigations Series I-2809

Kealakekua Fault is probably part of the head scarp of one such giant landslide, called the Alika landslide. This explains the steep slopes in the bay itself: after a narrow wave-cut platform, a spectacular wall covered with coral – the continuation of the cliff that you can see onshore – dives into the deep blue of the ocean as you float away from the shore. In contrast with submarine landslides that involve well stratified sediments failing along bedding surfaces and forming relatively thin but extensive slide deposits, the Hawaiian failures affect thick stacks of poorly layered volcanic rock and, as a result, both their volumes and morphologic relief are larger (see the paper by Lipman et al, 1988). The entire volume of the Alika slide is estimated to be 1500-2000 cubic kilometers. That is about a hundred times larger than all the sediment carried by the world’s rivers to the ocean in one year! The slides have moved at highway speeds and generated tsunamis. There is evidence on Lanai island for a wave that carried marine debris to 325 meters above sea level; this tsunami was likely put in motion by the Alika landslide*.

You don’t want to be snorkeling in Kealakekua Bay when something like that happens. And it will happen again, it is a matter of (geological) time. Giant underwater landslides are part of the normal life of these mid-ocean, hotspot-related volcanoes.

Reference
Lipman, P., Normark, W., Moore, J., Wilson, J., Gutmacher, C., 1988, The giant submarine Alika debris slide, Mauna Loa, Hawaii. Journal of Geophysical Research, vol. 93, p. 4279-4299.

*tsunamis generated by landslides is a whole new exciting subject that we have no time now to dive or snorkel into.