Posing the problem

Given a sediment column of a certain lithology with its top at and its base at , 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 and work towards finding the new base .

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):

where is the porosity at depth and is the depth where the initial porosity was reduced by .

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:

where c is a coefficient with the units . 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).

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):

As we have seen before, porosity at depth y is

The first equation then becomes

But

and

where and 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

If we integrate this over the interval to we get

Integrating this yields

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

The decompacted value of is

Now we can write the general decompaction equation:

That is,

The average porosity at the new depth will be

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

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

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:

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.

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.

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.

There has been a lot of talk in recent months about how bad the deforestation problem is in Romania, especially related to the expansion of the Austrian company Schweighofer. Although there is a lot of chit-chat in the media and on the internet on the subject, data and facts are not that easy to find. Not that easy — unless you make a bit of an effort: it turns out that a massive dataset on forest cover change in Eastern Europe is available for download, thanks to scientists at the University of Maryland who have published a paper on the dataset in the journal Remote Sensing of Environment.

The dataset is based on Landsat imagery collected between 1985 and 2012 and has a pixel size of 30×30 meters. It can be downloaded from this website. I am especially interested in what’s going on in parts of the Carpathians, because that’s where I grew up and went on lots of hikes and through lots of adventures in the 1980s and early 1990s. What follows are a few screenshots that give an idea about what is possible to see with this data.

This first image shows the whole Carpathian – intra-Carpathian area. There are five colors that correspond to different histories of forest cover: black pixels are places without forest during the time of study; green is stable forest; blue is forest gain; red is forest loss; and purple is forest loss followed by forest gain [there are two additional categories in the data, but they are not very common in this area].

The good news is that there is a lot of green in this map, which means that about 28.5% of the Carpathian area was continuously covered by forest since 1985. An additional 3.4% was without tree cover in 1985, but has gained forest cover since then. The forest is gone from 1.5% of the area; and 1.7% has lost and then regained tree cover.

If we zoom in to the ‘bend’ area, which is essentially the southeastern corner of Transylvania, Romania, it becomes more obvious that, although the big picture doesn’t look very bad, some places have been affected fairly significantly over the last three decades:

The black, forest-free patches in the middle are young intramontane basins with no forest. Deforestation looks to be more of a problem in the the Ciuc/Csík and Gheorgheni/Gyergyó basins. Let’s have a closer look at the Csík Basin:

A further zoom-in shows one of the areas with the most forest loss, the western side of the Harghita Mountains:

I am going to stop here; but this dataset has a lot more to offer than I showed in this post. The conclusion from this certainly shouldn’t be that everything is fine; often the issue is not so much the quantity of the forest being cut, but *where* it is being cut. Protected areas and national parks should clearly be green and stay green on these maps; and the Romanian Carpathians could use a few more protected lands, as many of these forests have never been cut (unlike a lot of forests in Western Europe).

I have used IPython Notebook with the GDAL package to create these images. The notebook can be viewed and downloaded over here.

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:

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

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

Solving the first Taylor expansion above for and dropping all higher-order terms yields the forward difference operator:

Similarly, the second equation yields the backward difference operator:

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

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:

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

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

Now the diffusion equation can be written as

This is equivalent to:

The expression is called the diffusion number, denoted here with s:

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

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

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

where 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:

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):

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:

Written in matrix notation, this is equiavlent to

After some reshuffling we get

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:

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:

In matrix form this is equivalent to

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.

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:

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

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

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

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

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.

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 .

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.

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.

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 (click for higher resolution):

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

Update (04/14/2016): If you want to use the animation, feel free to do so, as long as you (1) give credit to NASA/USGS Landsat, (2) give credit to me (= Zoltan Sylvester, geologist), and (3) link to this page. Note that you can see/download the high-resolution version if you click on the image.

Update (04/15/2016): Matteo Niccoli (who has a great blog called MyCarta) has created a slower version of the animation:

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.

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.

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.

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:

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.

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 are three photos that demonstrate how well some animals ‘know’ the composition and grain size of the sediment of their homeland. In all three cases, the salt-and-pepper sand of the background is the result of dark-colored grains of volcanic origin mixing with light-colored fragments of corals and seashells.

Ghost crab on Big Island, Hawaii

Ghost crabs are common on many beaches around the world; below is another one from Costa Rica. Note how the different grain size and sediment color are perfectly captured by the crabs.

Ghost crab near Playa Conchal, Guanacaste, Costa Rica

The third example takes us back to Hawaii, under water, where I managed to get this shot of a flounder while snorkeling in Kona. The only reason I saw this guy was that I got really close and part of the seafloor, which later turned out to be the flatfish, unexpectedly took off.