# Exploring (de)compaction with Python

All clastic sediments are subject to compaction (and reduction of porosity) as the result of increasingly tighter packing of grains under a thickening overburden. Decompaction – the estimation of the decompacted thickness of a rock column – is an important part of subsidence (or geohistory) analysis. The following exercise is loosely based on the excellent basin analysis textbook by Allen & Allen (2013), especially their Appendix 56. You can download the Jupyter notebook version of this post from Github.

## Import stuff

import numpy as np
import matplotlib.pyplot as plt
import functools
from scipy.optimize import bisect
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
plt.rcParams['mathtext.fontset'] = 'cm'

## Posing the problem

Given a sediment column of a certain lithology with its top at $y_1$ and its base at $y_2$, we are trying to find the thickness and average porosity of the same sediment column at a different depth (see figure below). We are going to set the new top $y_1'$ and work towards finding the new base $y_2'$.

plt.figure(figsize=(2,5))
x = [0,1,1,0,0]
y = [1,1,1.5,1.5,1]
plt.text(-0.6,1.02,'$y_1$',fontsize=16)
plt.text(-0.6,1.52,'$y_2$',fontsize=16)
plt.text(-0.6,1.27,'$\phi$',fontsize=16)
plt.fill(x,y,'y')
x = [3,4,4,3,3]
y = [0.5,0.5,1.15,1.15,0.5]
plt.text(2.25,0.52,'$y_1\'$',fontsize=16)
plt.text(2.25,1.17,'$y_2\'$',fontsize=16)
plt.text(2.25,0.9,'$\phi\'$',fontsize=16)
plt.fill(x,y,'y')
plt.plot([1,3],[1,0.5],'k--')
plt.plot([1,3],[1.5,1.15],'k--')
plt.gca().invert_yaxis()
plt.axis('off');

## Porosity decrease with depth

Porosity decreases with depth, initially largely due to mechanical compaction of the sediment. The decrease in porosity is relatively large close to the seafloor, where sediment is loosely packed; the lower the porosity, the less room there is for further compaction. This decrease in porosity with depth is commonly modeled as a negative exponential function (Athy, 1930):

$\phi(y) = \phi_0 e^{-\frac{y}{y_0}}$

where $\phi(y)$ is the porosity at depth $y$ and $y_0$ is the depth where the initial porosity $\phi_0$ was reduced by $1/e$.

This is an empirical equation, as there is no direct physical link between depth and porosity; compaction and porosity reduction are more directly related to the increase in effective stress under a thicker overburden. Here we only address the simplest scenario with no overpressured zones. For normally pressured sediments, Athy’s porosity-depth relationship can be expressed in a slightly different form:

$\phi(y) = \phi_0 e^{-cy}$

where c is a coefficient with the units $km^{-1}$. The idea is that c is a characteristic constant for a certain lithology and it can measured if porosity values are available from different depths. Muds have higher porosities at the seafloor than sands but they compact faster than sands. The plot below show some typical curves that illustrate the exponential decrease in porosity with depth for sand and mud. The continuous lines correspond to the parameters for sand and mud in Appendix 56 of Allen & Allen (2013); the dashed lines are exponential fits to data from the Ocean Drilling Program (Kominz et al., 2011).

c_sand = 0.27 # porosity-depth coefficient for sand (km-1)
c_mud = 0.57 # porosity-depth coefficent for mud (km-1)
phi_sand_0 = 0.49 # surface porosity for sand
phi_mud_0 = 0.63 # surface porosity for mud
y = np.arange(0,3.01,0.01)

phi_sand = phi_sand_0 * np.exp(-c_sand*y)
phi_mud = phi_mud_0 * np.exp(-c_mud*y)

plt.figure(figsize=(4,7))
plt.plot(phi_sand,y,'y',linewidth=2,label='sand')
plt.plot(phi_mud,y,'brown',linewidth=2,label='mud')
plt.xlabel('porosity')
plt.ylabel('depth (km)')
plt.xlim(0,0.65)
plt.gca().invert_yaxis()

c_sand = 1000/18605.0 # Kominz et al. 2011 >90% sand curve
c_mud = 1000/1671.0 # Kominz et al. 2011 >90% mud curve
phi_sand_0 = 0.407 # Kominz et al. 2011 >90% sand curve
phi_mud_0 = 0.614 # Kominz et al. 2011 >90% mud curve
phi_sand = phi_sand_0 * np.exp(-c_sand*y)
phi_mud = phi_mud_0 * np.exp(-c_mud*y)
plt.plot(phi_sand,y,'y--',linewidth=2,label='90% sand')
plt.plot(phi_mud,y,'--',color='brown',linewidth=2,label='90% mud')
plt.legend(loc=0, fontsize=10);


While the compaction trends for mud happen to be fairly similar in the plot above, the ones for sandy lithologies are very different. This highlights that porosity-depth curves vary significantly from one basin to another, and are strongly affected by overpressures and exhumation. Using local data and geological information is critical. As Giles et al. (1998) have put it, “The use of default compaction curves can introduce significant errors into thermal history and pore- fluid pressure calculations, particularly where little well data are available to calibrate the model.” To see how widely – and wildly – variable compaction trends can be, check out the review paper by Giles et al. (1998).

## Deriving the general decompaction equation

Compacting or decompacting a column of sediment means that we have to move it along the curves in the figure above. Let’s consider the volume of water in a small segment of the sediment column (over which porosity does not vary a lot):

$dV_w = \phi dV_t$

As we have seen before, porosity at depth y is

$\phi(y) = \phi_0 e^{-cy}$

The first equation then becomes

$dV_w = \phi_0 e^{-cy} dV_t$

But

$dV_w = A dy_w$

and

$dV_t = A dy_t$

where $y_w$ and $y_t$ are the thicknesses that the water and total volumes occupy respectively, and A is the area of the column we are looking at. So the relationship is equivalent to

$dy_w = \phi_0 e^{-cy} dy_t$

If we integrate this over the interval $y_1$ to $y_2$ we get

$y_w = \int_{y1}^{y2} \phi_0 e^{-cy} dy_t$

Integrating this yields

$y_w = \phi_0 \Bigg(\frac{1}{-c}e^{-cy_2} - \frac{1}{-c}e^{-cy_1}\Bigg) = \frac{\phi_0}{c} \big(e^{-cy_1}-e^{-cy_2}\big)$

As the total thickness equals the sediment thickness plus the water “thickness”, we get

$y_s = y_t - y_w = y_2 - y_1 - y_w = y_2 - y_1 - \frac{\phi_0}{c} \big(e^{-cy_1}-e^{-cy_2}\big)$

The decompacted value of $y_w$ is

$y_w' = \frac{\phi_0}{c} \big(e^{-cy_1'}-e^{-cy_2'}\big)$

Now we can write the general decompaction equation:

$y_2'-y_1' = y_s+y_w'$

That is,

$y_2'-y_1' = y_2 - y_1 - \frac{\phi_0}{c} \big(e^{-cy_1}-e^{-cy_2}\big) + \frac{\phi_0}{c} \big(e^{-cy_1'}-e^{-cy_2'}\big)$

The average porosity at the new depth will be

$\phi = \frac{Ay_w'}{Ay_t'} = \frac{\phi_0}{c}\frac{\big(e^{-cy_1'}-e^{-cy_2'}\big)}{y_2'-y_1'}$

## Write code to compute (de)compacted thickness

The decompaction equation could be solved in the ‘brute force’ way, that is, by gradually changing the value of $y_2'$ until the right hand side (RHS) of the equation is the same as the left hand side (LHS) – see for example the Excel spreadsheet that accompanies Appendix 56 in Allen & Allen (2013). However, we (and scipy) can do better than that; we will use bisection, one the simplest optimization methods to find the root of the function that we set up as RHS-LHS.

# compaction function - the unknown variable is y2a
def comp_func(y2a,y1,y2,y1a,phi,c):
# left hand side of decompaction equation:
LHS = y2a - y1a
# right hand side of decompaction equation:
RHS = y2 - y1 - (phi/c)*(np.exp(-c*y1)-np.exp(-c*y2)) + (phi/c)*(np.exp(-c*y1a)-np.exp(-c*y2a))
return LHS - RHS


Now we can do the calculations; here we set the initial depths of a sandstone column $y_1,y_2$ to 2 and 3 kilometers, and we estimate the new thickness and porosity assuming that the column is brought to the surface ($y_1'=0$).

c_sand = 0.27 # porosity-depth coefficient for sand (km-1)
phi_sand = 0.49 # surface porosity for sand
y1 = 2.0 # top depth in km
y2 = 3.0 # base depth in km
y1a = 0.0 # new top depth in km


One issue we need to address is that ‘comp_func’ six input parameters, but the scipy ‘bisect’ function only takes one parameter. We create a partial function ‘comp_func_1’ in which the only variable is ‘y2a’, the rest are treated as constants:

comp_func_1 = functools.partial(comp_func, y1=y1, y2=y2, y1a=y1a, phi=phi_sand, c=c_sand)
y2a = bisect(comp_func_1,y1a,y1a+3*(y2-y1)) # use bisection to find new base depth
phi = (phi_sand/c_sand)*(np.exp(-c_sand*y1)-np.exp(-c_sand*y2))/(y2-y1) # initial average porosity
phia = (phi_sand/c_sand)*(np.exp(-c_sand*y1a)-np.exp(-c_sand*y2a))/(y2a-y1a) # new average porosity

print('new base depth: '+str(round(y2a,2))+' km')
print('initial thickness: '+str(round(y2-y1,2))+' km')
print('new thickness: '+str(round(y2a-y1a,2))+' km')
print('initial porosity: '+str(round(phi,3)))
print('new porosity: '+str(round(phia,3)))


## Write code to (de)compact a stratigraphic column with multiple layers

Next we write a function that does the depth calculation for more than one layer in a sedimentary column:

def decompact(tops,lith,new_top,phi_sand,phi_mud,c_sand,c_mud):
tops_new = [] # list for decompacted tops
tops_new.append(new_top) # starting value
for i in range(len(tops)-1):
if lith[i] == 0:
phi = phi_mud; c = c_mud
if lith[i] == 1:
phi = phi_sand; c = c_sand
comp_func_1 = functools.partial(comp_func,y1=tops[i],y2=tops[i+1],y1a=tops_new[-1],phi=phi,c=c)
base_new_a = tops_new[-1]+tops[i+1]-tops[i]
base_new = bisect(comp_func_1, base_new_a, 4*base_new_a) # bisection
tops_new.append(base_new)


Let’s use this function to decompact a simple stratigraphic column that consists of 5 alternating layers of sand and mud.

tops = np.array([1.0,1.1,1.15,1.3,1.5,2.0])
lith = np.array([0,1,0,1,0]) # lithology labels: 0 = mud, 1 = sand
phi_sand_0 = 0.49 # surface porosity for sand
phi_mud_0 = 0.63 # surface porosity for mud
c_sand = 0.27 # porosity-depth coefficient for sand (km-1)
c_mud = 0.57 # porosity-depth coefficent for mud (km-1)
tops_new = decompact(tops,lith,0.0,phi_sand_0,phi_mud_0,c_sand,c_mud) # compute new tops


Plot the results:

def plot_decompaction(tops,tops_new):
for i in range(len(tops)-1):
x = [0,1,1,0]
y = [tops[i], tops[i], tops[i+1], tops[i+1]]
if lith[i] == 0:
color = 'xkcd:umber'
if lith[i] == 1:
color = 'xkcd:yellowish'
plt.fill(x,y,color=color)
x = np.array([2,3,3,2])
y = np.array([tops_new[i], tops_new[i], tops_new[i+1], tops_new[i+1]])
if lith[i] == 0:
color = 'xkcd:umber'
if lith[i] == 1:
color = 'xkcd:yellowish'
plt.fill(x,y,color=color)
plt.gca().invert_yaxis()
plt.tick_params(axis='x',which='both',bottom='off',top='off',labelbottom='off')
plt.ylabel('depth (km)');

plot_decompaction(tops,tops_new)


Now let’s see what happens if we use the 90% mud and 90% sand curves from Kominz et al. (2011).

tops = np.array([1.0,1.1,1.15,1.3,1.5,2.0])
lith = np.array([0,1,0,1,0]) # lithology labels: 0 = mud, 1 = sand
c_sand = 1000/18605.0 # Kominz et al. 2011 >90% sand curve
c_mud = 1000/1671.0 # Kominz et al. 2011 >90% mud curve
phi_sand_0 = 0.407 # Kominz et al. 2011 >90% sand curve
phi_mud_0 = 0.614 # Kominz et al. 2011 >90% mud curve
tops_new = decompact(tops,lith,0.0,phi_sand_0,phi_mud_0,c_sand,c_mud) # compute new tops

plot_decompaction(tops,tops_new)


Quite predictably, the main difference is that the sand layers have decompacted less in this case.

That’s it for now. It is not that hard to modify the code above for more than two lithologies. Happy (de)compacting!

## References

Allen, P. A., and Allen, J. R. (2013) Basin Analysis: Principles and Application to Petroleum Play Assessment, Wiley-Blackwell.

Athy, L.F. (1930) Density, porosity and compaction of sedimentary rocks. American Association Petroleum Geologists Bulletin, v. 14, p. 1–24.

Giles, M.R., Indrelid, S.L., and James, D.M.D., 1998, Compaction — the great unknown in basin modelling: Geological Society London Special Publications, v. 141, no. 1, p. 15–43, doi: 10.1144/gsl.sp.1998.141.01.02.

Kominz, M.A., Patterson, K., and Odette, D., 2011, Lithology Dependence of Porosity In Slope and Deep Marine Sediments: Journal of Sedimentary Research, v. 81, no. 10, p. 730–742, doi: 10.2110/jsr.2011.60.

Thanks to Mark Tingay for comments on ‘default’ compaction trends.

# Cutoffs in submarine channels

A paper on sinuous submarine channels that is based on work I have done with Jake Covault (at the Quantitative Clastics Laboratory, Bureau of Economic Geology, UT Austin) has been recently published in the journal Geology (officially it will be published in the October issue, but it has been online for a few weeks now).

The main idea of the paper is simple: many submarine channels become highly sinuous and cutoff bends, similar to those observed in the case of rivers, must be fairly common. The more high-quality data becomes available from these systems, the more evidence there is that this is indeed the case. I think this observation is fairly interesting in itself, as cutoffs will add significant complexity to the three-dimensional structure of the preserved deposits. We have addressed that complexity to some degree in another paper. However, previous work on submarine channels has not dealt with the implications of how the along-channel slope changes as cutoff bends form.

So, inspired by a model of incising subaerial bedrock rivers, we have put together a numerical model that not only mimics the plan-view evolution of meandering channels, by keeping track of the x and y coordinates of the channel centerline, but it also records the elevations (z-values) at each control point along this line. The plan-view evolution is based on the simplest possible model of meandering, described in a fantastic paper by Alan Howard and Thomas Knutson published in 1984 and titled “Sufficient conditions for river meandering: A simulation approach“. The key idea is that the rate of channel migration is not only a function of local curvature (in the sense of larger curvatures resulting in larger rates of outer bank erosion and channel migration), but it also depends on the channel curvature upstream from the point of interest; and the strength of this dependence declines exponentially as you move upstream. Without this ‘trick’ there is no stable meandering. The model starts out with an (almost) straight line that has some added noise; meanders develop through time as a dominant wavelength gets established. The z-coordinates are assigned initially so that they describe a constant along-channel slope; as the sinuosity increases, the overall slope decreases and local variability appears as well, depending of the local sinuosity of the channel. More discrete changes take place at the time and location of a cutoff: this is essentially a shortcut that introduces a much steeper segment into the centerline. The animations below (same as two of the supplementary files that were published with the paper) show how the cutoffs form and how the along-channel slope profiles change through time.

Animation of incising submarine channels with early meander cutoffs.

Evolution of along-channel gradient (top) and elevation profile (bottom) through time.

Meandering lowland rivers (like the Mississippi) have very low slopes to start with, so cutoff-related elevation drops are not as impressive as in the case of steep (with up to 40-50 m/km gradients) submarine channels (see diagram below). As a result, knickpoints that form this way in submarine channels are likely to be the locations of erosion and contribute to the stratigraphic complexity of the channel deposits. Previous work has shown that avulsion- and deformation-driven knickpoints are important in submarine channels; our results suggest that the meandering process itself can generate significant morphologic and stratigraphic along-channel variability.

Knickpoint gradient plotted against the overall slope. Dashed lines show trends for different cutoff distances (beta is the ratio between cutoff distance and channel width).

Finally I should mention that, apart from the seismic interpretation itself, we have done all the modeling, analysis, visualization, and plotting with Jupyter notebooks and Mayavi (for 3D visualization), using the Anaconda Python distribution. I cannot say enough good things about – and enough thanks for – these tools.

# Stratigraphic patterns in slope minibasins

A few months ago some former colleagues and myself have published a paper on the stratigraphic evolution of minibasins that are present on many continental slopes (Sylvester, Z., Cantelli, A., and Pirmez, C., 2015, Stratigraphic evolution of intraslope minibasins: Insights from surface-based model: AAPG Bulletin, v. 99, no. 6, p. 1099–1129, doi: 10.1306/01081514082.) We used a simple model that investigates the interplay between subsidence and sedimentation and helps in the understanding of how stratal termination patterns relate to variations in sediment input and basin subsidence. Conventional sequence stratigraphy focuses on what is happening on the shelf and the shelf edge; and the origin of stratal patterns that characterize the typical ‘slug diagram’ are not trivial to relate to the three main parameters that influence continental-margin stratigraphy: changes in sea level, sediment supply, and subsidence. Our minibasin model is simpler because (1) it only has two parameters: sediment supply and subsidence (the impact of sea level – if any – is included in the sediment supply curve); and (2) sedimentary layers are deposited horizontally during each time step. The basic idea is to investigate the geometries of a system that consists of a depression that deepens through time and sediment is deposited in it with a certain rate. The figures that follow are color versions of the model cross sections that were published in the paper. The originals are actually interactive (you can zoom in and pan if you want to, which is useful if you want to check some details, especially along the basin margins). However, I could not get WordPress to display these in an interactive format; you can check out the interactive version of this post on Github. I have used the excellent, d3.js-based – and recently open-sourced – Plotly to create the interactive plots. A few brief comments follow so that the plots make some sense even if you don’t read the paper.

There are two types of layers in the model: (1) the ones that correspond to the sediment supply that comes into the basin from updip, mostly in the form of turbidity currents – these are shown in yellow color in the plots below; and (2) thinner ones that reflect an overall finer-grained background sedimentation – these are shown in brown color. The latter have a constant thickness across the basin. It is important to note that, while it is true that all the sand is supposed to be in the yellow units, not all of the yellow units are coarse-grained. In fact, slope minibsains have overall a much lower sand content than what a ‘yellow = sand’ assumption would imply in these plots. The blue dots are pinchout points, locations where the yellow layers terminate. In the paper we adopted the idea that ‘onlap’ refers to the case when these termination points move toward the basin edge, and ‘offlap’ corresponds to stratigraphic intervals where the pinchout locations migrate toward the center of the basin.

Sedimentary fills of slope minibasins commonly show strata that either converge toward or onlap onto the basin margins (or show a combination of the two). The first two cross sections are examples of ‘pure’ convergence and well-defined onlap. In the first case (convergence), both sediment input rate and subsidence are constant and have similar orders of magnitude, and, as a result, the locations where sedimentary layers terminate are stationary. The basin never gets very deep, as sedimentation keeps up with subsidence.

Obviously animations are much better at conveying the idea of sedimentation taking place while the basin is subsiding (the background sedimentation is set to zero in the animation; the lower panel is a chronostratigraphic – or Wheeler – diagram):

For the formation of onlap, a relatively deep basin and high sediment input are needed. The model below has the same overall sediment volume as the previous one, but it was deposited in half the time. Subsidence was kept constant, just like in the convergence model.

And here is the animation for the onlap scenario:

The next two plots illustrate what happens if the sediment input varies following a sinusoidal curve with three cycles, while subsidence is constant. In the first model, there is enough accommodation in the basin to keep all the sediment; in the second model, the sediment input exceeds the space available in the basin and each cycle has an upper part during which some of the sediment bypasses the basin. If you zoom in, you can see that both onlap and offlap are present in both cases.

And the animation for the no-bypass scenario:

Next we looked at what happens if sediment input is kept constant, but subsidence varies (following a sinusoidal function). It is tempting to think that the result must be similar to the case of variable sediment input and constant subsidence, but that is *not* the case. In  previous models onlap surfaces overlie condensed sections across the whole basin. In this new model, onlap surfaces are restricted to the basin margin and they correlate to sections with high sedimentation rate at the basin center.

It is a cool idea that stratal patterns might tell you whether variability in sediment supply or deformation is more important. However, if both parameters are changing through time, with similar magnitudes, it is practically impossible to tell what is going on:

This animation gives you an idea how can these two cases transition into each other, as sediment input goes from constant to sinusoidal and out-of-phase with subsidence:

So far all the simulations were based on sinusoidally varying parameters; what happens though if sediment supply is closer to an ‘on-off’ function? You can see the result below: each sedimentary cycle only shows onlap, the offlapping upper part is missing.

The next two plots are attempts to reproduce the stratal patterns seen in two well-studied minibasins in the Gulf of Mexico: Brazos-Trinity Basin 4 and the Auger Basin.

Finally, we looked at the common scenario of multiple linked minibasins, with sediment bypassing the updip basin serving as the sediment input for the basin downdip. The classic idea about this depositional setting is called ‘fill and spill’: the updip basin has to be filled first before any sediment gets into or spills the second basin. A ‘static’ version of this idea, with pre-formed basins and no subsidence during sedimentation, is illustrated in the plot below, with three sediment input cycles:

The animated version looks like this:

However, this static fill and spill is probably the exception rather than the rule: most minibasins keep subsiding while sedimentation is taking place. So, in contrast with the previous model, the large-scale fills of the basins are time-equivalent, while at the scale of individual sediment input cycles the fill-and-spill model still works.

Now let’s animate this:

For more details on these models, see the paper in AAPG Bulletin. Ten related animations are open access and available on the AAPG website. And, again, the plots are a lot more interesting if you can zoom, pan and rescale, so you should really check out the intended format of this post.

# 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

# 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.

# 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.

# The complexity of sinuous channel deposits in three dimensions

The beauty of the shapes and patterns created by meandering rivers has long attracted the attention of many geomorphologists, civil engineers, and sedimentologists. Unless they are fairly steep or have highly stable and unerodible banks, rivers do not like to follow a straight course and tend to develop a sinuous plan-view pattern. The description and mathematical modeling of these curves is a fascinating subject, but that is not what I want to talk about here and now. It is hard enough to understand the plan-view evolution of rivers, especially if one is interested in the long-term results – when cutoffs become important -, but things get really complicated when it comes to the three-dimensional structure of the deposits that meandering rivers leave behind. The same can be said about sinuous channels on the seafloor, created and maintained by dirty mixtures of water and sediment (called turbidity currents). An ever-increasing number of seafloor and seismic images show that highly sinuous submarine channels are almost as common as their subaerial counterparts, but much remains to be learned about the geometries of their deposits that accumulate through geological time.

Using simple modeling of how channel surfaces migrate through time, two recent papers attempt to illustrate the three-dimensional structure of sinuous fluvial and submarine channel deposits. In the Journal of Sedimentary Research, Willis and Tang (2010) show how slightly different patterns of fluvial meander migration result in different deposit geometries and different distribution of grain size, porosity and permeability. [These properties are important because they determine how fluids flow – or don’t flow – through the pores of the sediment.] River meanders can either grow in a direction perpendicular to the overall downslope orientation, or they can keep the same width and migrate downstream through translation. In the latter case – which is often characteristic of rivers incising into older sediments -, deposits forming on the downstream, concave bank of point bars will be preferentially preserved. These deposits tend to be finer grained than the typical convex-bank point bar sediments. In addition to building a range of models and analyzing their geometries, Willis and Tang also ran simulations of how would oil be displaced by water in them. One of their findings is that sinuous rivers that keep adding sediment in the same area over time (in other words, rivers that aggrade) tend to form better connected sand bodies than rivers which keep snaking around roughly in the same horizontal plane, without aggradation.

 Map of deposits forming as river meanders grow (from Willis and Tang,  2010).
 Cross sections through the deposits of two meander bends (locations shown in figure above). Colors represent permeability, red being highly permeable and blue impermeable sediment. From Willis and Tang, 2010.

Check out the paper itself for more images like these, plus discussions of concave-bank deposition, cutoff formation, and filling of abandoned channels.

The second paper (Sylvester, Pirmez, and Cantelli, 2010; and yes, one of the authors is also the author of this blog post, so don’t expect any constructive criticism here) focuses on submarine channels and their overbank deposits, but the starting point and the modeling techniques are similar: take a bunch of sinuous channel centerlines and generate surfaces around them that reflect the topography of the system at every time step. However, we know much less about submarine channels than fluvial ones, because it is much more difficult to collect data at and from the bottom of the ocean than it is from the river in your backyard. The result is that some of the simplifications in our model are controversial; to many sedimentary geologists, submarine channels and their deposits are fundamentally different from rivers and point bars, and there is not much use in even comparing the two. Part of the problem is that not all submarine channels are made equal, and, when looking at an outcrop, it is not easy – or outright impossible – to tell what kind of geomorphology produced the  stratigraphy. In fact, the number of exposures that represent highly sinuous submarine channels, as observed on the seafloor and numerous seismic images, is probably fairly limited. One thing is quite clear, however: many submarine channels show plan-view migration patterns that are very similar to those of rivers, and this large-scale structure imposes some significant constraints on the geometry of the deposits as well.

That being said, nobody denies that there are plenty of significant differences between real and submarine ‘rivers’ [note quotation marks]. A very important one is the amount of overbank – or levee – deposition: turbidity currents often overflow their channel banks as thick muddy clouds and form much thicker deposits than the overbank sediment layers typical of rivers. When these high rates of levee deposition combine with the strong three-dimensionality of channel migration, complex geometries result that are quite tricky to understand just by looking at a single cross section.

 Cross section and chronostratigraphic diagram through a submarine channel system with inner and outer levees (from Sylvester et al., 2010).

One of the consequences of the channel migration is the formation of erosional surfaces that develop through a relatively long time and do not correspond to a geomorphologic surface at all (see the red erosional zones in the Wheeler diagram above). This difference between stratigraphic and geomorphologic surfaces is essential, yet often downplayed or even ignored in stratigraphy. In terms of geomorphology, the combination of channel movement in both horizontal and vertical directions and the extensive levee deposition results in a wide valley with scalloped margins and numerous terraces inside:

 Three-dimensional view of an incising channel-levee system (from Sylvester et al., 2010).

This second paper is part of a nice collection focusing on submarine sedimentary systems that is going to be published as a special issue of Marine and Petroleum Geology, a collection that originated from a great conference held in 2009 in Torres del Paine National Park, Southern Chile.

PS. As I am typing this, I see that Brian over at Clastic Detritus is also thinking about submarine channels and subaerial rivers… Those channels formed by saline density currents on the slope of the Black Sea are fascinating.

Willis, B., & Tang, H. (2010). Three-Dimensional Connectivity of Point-Bar Deposits Journal of Sedimentary Research, 80 (5), 440-454 DOI: 10.2110/jsr.2010.046

Sylvester, Z., Pirmez, C., & Cantelli, A. (2010). A model of submarine channel-levee evolution based on channel trajectories: Implications for stratigraphic architecture Marine and Petroleum Geology DOI: 10.1016/j.marpetgeo.2010.05.012

# Hillslope diffusion

Modeling erosion and deposition of sediment using the diffusion equation is among the important subjects that are usually omitted from sedimentary geology textbooks. Part of the reason for this is that ‘conventional’ sedimentary geology tended to only pay lip service to earth surface processes and was more interested in describing the stratigraphic record than figuring how it relates to geomorphology. Nowadays, a good discussion of stratigraphy and sedimentology cannot ignore anymore what geomorphologists have learned about landscape evolution. (One textbook that clearly recognizes this is this one.)

But let’s get back to the subject of this post. Hillslope evolution can be modeled with the diffusion equation, one of the most common differential equations in science, applied for example to describe how differences in temperature are eliminated through heat conduction. In the case of heat, the heat flux is proportional to the rate of spatial temperature change; on hillslopes, the sediment flux is proportional to the spatial rate of change in elevation. This last quantity of course is the slope itself. In other words,

q = -k*dh/dx,

or

q = -k*slope,

where q is the volumetric sediment flux per unit length, k is a constant called diffusivity, h is the elevation, and x is the horizontal coordinate.

We also know that sediment does not disappear into thin air: considering a small area of the hillslope, the amount of sediment entering and leaving this area will determine how large the change in elevation will be:

dh/dt = -dq/dx,

in other words, deposition or erosion at any location is determined by the change in sediment flux.

Combining this equation with the previous one, we arrive to the diffusion equation:

dh/dt = k*d2h/dx2.

Note that the quantity on the right side is the second derivative (or curvature) of the slope profile. Large negative curvatures result in rapid erosion; places with large positive curvature have high rates of deposition. Through time, the bumps and troughs of the hillslope are smoothed out through erosion and deposition.

The simplest possible case is the diffusion of a fault scarp. The animation below illustrates how a 1 m high fault scarp gets smoothed out through time; the evolution of slope and curvature are also shown. The dashed line indicates the original topography, at time 0. [The plots were generated using Ramon Arrowsmith’s Matlab code].

More complicated slope profiles can be modeled as well; here is an example with two fault scarps:

Note how both erosion and deposition get much slower as the gradients become more uniform.

The simplicity of the diffusion equation makes it an attractive tool in modeling landscape evolution. In addition to hillslopes and fault scarps, it has been successfully applied in modeling – for example – river terraces, deltaic clinoforms, cinder cones, fluvial systems, and foreland basin stratigraphy. However, it is important to know when and where the assumptions behind it become invalid. For example, steep slopes often have a non-linear relationship between sediment flux and slope as mass movements dramatically increase sediment flux above a critical slope value. Also, the models shown here would fail to reproduce the topography of a system where not all sediment is deposited at the toe of the steeper slope, but a significant part is carried away by a river. And that brings us closer to advection; a subject that I might take notes about at another time.

Further reading: 1) The book “Quantitative Modeling of Earth Surface Processes” by Jon Pelletier has a chapter with lots of details about the diffusion equation. 2) Analog and numerical modeling of hillslope diffusion– a nice lab exercise.

# Climbing Ripples I.

Ripples, dunes, cross bedding and cross lamination have always been some of the sexiest subjects in sedimentary geology. They are certainly responsible (in part) for my choice of a certain walk of life that consists of studying dirt. You might say that everything has been already said about ripples and dunes, and you clearly get that feeling if you read some of J.R.L. Allen’s work on the subject (and that can be a lot of reading, by the way) or look at the fantastic multimedia material that David Rubin at the USGS put together. [Of course, there are numerous other authors who have written great papers on the subject, but it is not my purpose here to write a history of bedform sedimentology. Although that would be an interesting subject, if somebody had the time for it.]

However, little of this material gets into the standard sedimentology and stratigraphy textbooks. Maybe rightly so: after all, textbooks are not supposed to include all the details about any particular subject. And maybe there are higher-density issues out there, like whether we should call something a turbidite or a debrite. [Sorry, I could not refrain from typing that].

Take for example climbing ripples. They form when several trains of ripples are superimposed on each other and they seem to ‘climb’, by generating stratigraphic surfaces that are tilted in an upcurrent direction. [Note however that these surfaces are *not* topographic – or time – surfaces; more on that later]. Numerous textbooks and many papers mention climbing ripple cross lamination, but often the explanation is something like “they indicate high rates of deposition”, or “the steepness of the climb and stoss-side preservation are a function of the ratio between suspended-load and bedload”. The question is, what do we *exactly* mean by ‘high rates of deposition’? If we cannot put numbers on it, it is not that informative. Also, by ‘suspended load’, do we mean suspended load concentration? Or deposition from suspended load and bedload, respectively? Those statements are not necessarily wrong, but they do not do justice to the models that have been published many years ago, models that actually have some numbers and equations behind the “conclusion” section.

The key paper that I am talking about is “A quantitative model of climbing ripples and their cross-laminated deposit“, by J.R.L. Allen, published in 1970 in the journal Sedimentology.

The most important relationship that Allen has derived links the angle of climb ζ (see the sketch below) to the rate of deposition M (measured in units of mass over unit time and area), the rate of bedload sediment transport j, and the ripple height H:

tanζ = MH / 2j

This is simply based on decomposing the sediment flux to and through the bed into vertical and horizontal components (plus a relationship between the horizontal sediment transport rate in ripples and the horizontal migration rate of the bedforms). Note that the quantity j refers to the sediment mass that moves through a cross section perpendicular to the general current direction, and does this by being part of the ripples themselves. In other words, there is no direct equivalence between M and suspended load deposition, and j and bedload deposition. Although it is possible that in general suspended load contributes more to M than deposition from bedload, it says nowhere that grains transported within the bedload cannot be deposited on the stoss side of the ripples and thus contribute to the vertical growth of the bed.

Obviously, if the angle of climb is smaller than the dip of the stoss side, there will be no stoss side preservation and the resulting cross lamination will look like in the sketch below (which, by the way, was quite an effort to generate in Matlab; you can easily do this and much-much more with David Rubin’s Matlab code, but I wanted to understand things a little better by coding something simple myself):

This is often called ‘A-type’ (or subcritical) climbing ripple cross lamination, but everybody knows what you are talking about if you “simply” call it climbing ripple cross lamination with no stoss-side preservation.

In contrast, aggradation is much more prominent if the angle of climb is larger than the slope of the stoss side, and in this case deposition takes place on the stoss sides as well, resulting in ‘S-type’ (or supercritical) lamination:

Of course, it says nowhere that the rate of deposition M or the bedload transport rate j must stay constant through time. If the ratio of these quantities changes, the angle of climb will change as well. This sketch shows an example where the rate of deposition M increases through time:

One of the main points of the paper is that there is a fundamental difference between the rate of deposition M and the bedload sediment transport rate j. A rate of deposition larger than zero means that the sediment transport rate within the flow must decrease from an upcurrent position to a downcurrent position; a simple mass balance tells us that this change in the sediment transport rate has to equal the rate of deposition. In other words, the rate of deposition M is a derivative of the sediment transport rate, and as such, does not belong in the same drawer of physical quantities as the bedload transport rate.

Along the same line of thought, Allen emphasizes that climbing ripple lamination says something about flow uniformity and steadiness. A uniform and steady flow can only form a single train of ripples; either non-uniformity or unsteadiness is needed to have climbing-ripple deposition.

That’s it for now; to be continued. It’s time to do my taxes.

Further reading: Brian has a Friday Field Photo and a Geopuzzle on climbing ripples. Here are some pictures and a movie of climbing ripples generated by a turbidity current in a flume.