Skip to content

Grain Settling Equations and Plots With R

August 18, 2013

Zoltan Sylvester wrote a nice post a couple weeks ago exploring grain settling equations. He plotted up these equations using some Python code. Here, I simply reproduce what Zoltan did but with R code. So, go read Zoltan’s post first. This was more of an exercise in continuing to learn R rather than an exploration of grain settling physics. And, in this case, doing the ‘translation’ is a nice introduction to how Python works.

My R code works (yay), but seems bulky and cumbersome (boo). I’m sure there are ways to streamline it, feel free to comment with suggestions.

This first box of code defines constants, the three equations to be plotted, and then applies those equations to a range of particle sizes (diameters from 0 to 0.5 mm).

# Some constants for the equations
rop = 2650.0 # density of particle in kg/m^3
rof = 1000.0 # density of water in kg/m^3
visc = 1.002*1E-3 # dynamic viscosity in Pa*s at 20C
C1 = 18 # constant in Ferguson-Church equation
C2 = 1 # constant in Ferguson-Church equation

# Stokes
v_stokes = function(d){
  R = (rop-rof)/rof # submerged specific gravity
  w = R*9.81*(d^2)/(C1*(visc/rof))
  return(w)
}
# Turbulent
v_turb = function(d){
  R = (rop-rof)/rof 
  w = ((4*R*9.81*d)/(3*C2))^0.5
  return(w)
}
# Ferguson-Church equation
v_ferg = function(d){
  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)
}

# Define a range (sequence) of particle diameters
d = seq(0,0.0005,0.000001)

# Invoke the function for that range of particle diameters
ws = v_stokes(d)
wt = v_turb(d)
wf = v_ferg(d)

This next box then creates the plot shown below the code.

# Plot diameter (in mm) vs. settling velocity for diameters 0-0.5 mm
plot(d*1000, ws, ylim=c(0,0.15), type="l", lwd=3, col="blue", 
     xlab="Particle Diameter (mm)", 
     ylab="Settling Velocity (m/s)", yaxt="n")
lines(d*1000, wt, col="darkgreen", lwd=3)
lines(d*1000, wf, col="red", lwd=3)
# Axes labels, ticks, annotation, legend
axis(2, at=seq(0,0.14,0.02), cex.axis=0.85, las=1)
abline(v=c(0.063,0.125,0.250), col="gray40", lwd=0.9, lty=5)
text(0.022,0.105, "clay & silt", cex=0.75)
text(0.094,0.105, "vf sand", cex=0.75)
text(0.190,0.105, "fine sand", cex=0.75)
text(0.410,0.105, "medium sand", cex=0.75)
legend(0,0.15,legend=c("Stokes","Turbulent","Ferguson & Church"),
       cex=0.7, lty=c(1,1,1), lwd=c(3,3,3), col=c("blue","darkgreen","red"))
# Plot points from Table 2 of Ferguson & Church
D = c(0.068,0.081,0.096,0.115,0.136,0.273,0.386,0.550,0.770,1.090,2.180,4.360)
W = c(0.00425,0.0060,0.0075,0.0110,0.0139,0.0388,
     0.0551,0.0729,0.0930,0.141,0.209,0.307)
points(D,W, pch=20)

plot1

The next part is the same equations but for an expanded range of particle diameters (up to pebbles) and plotted in log-log space.

# Define a range (sequence) of particle diameters
d = seq(0,0.01,0.00001)

# Invoke the function for that range of particle diameters
ws = v_stokes(d)
wt = v_turb(d)
wf = v_ferg(d)

# Plot in log-log space
plot(d*1000, ws, ylim=c(0.00001,10), xlim=c(0.01,10), type="l", log="xy", 
     lwd=3, col="blue", xlab="Particle Diameter (mm)", 
     ylab="Settling Velocity (m/s)", xaxt="n", yaxt="n")
lines(d*1000, wt, col="darkgreen", lwd=3)
lines(d*1000, wf, col="red", lwd=3)
# Axes labels, ticks, annotation, legend
axis(2, at=c(0.00001,0.0001,0.001,0.01,0.1,1,10), cex.axis=0.85, las=1)
axis(1, at=c(0.01,0.02,0.05,0.1,0.2,0.5,1,2,5,10), cex.axis=0.85, las=1)
abline(v=c(0.0156,0.031,0.063,0.125,0.250,0.5,1,2,4), col="gray40", lwd=0.9, lty=5)
text(0.011,0.0008, "fine silt", cex=0.75, srt=90)
text(0.022,0.0021, "medium silt", cex=0.75, srt=90)
text(0.043,0.0004, "coarse silt", cex=0.75, srt=90)
text(0.09,0.0004, "very fine sand", cex=0.75, srt=90)
text(0.17,0.0004, "fine sand", cex=0.75, srt=90)
text(0.33,0.0004, "medium sand", cex=0.75, srt=90)
text(0.7,0.0004, "coarse sand", cex=0.75, srt=90)
text(1.4,0.0004, "very coarse sand", cex=0.75, srt=90)
text(2.7,0.0004, "granules", cex=0.75, srt=90)
text(7,0.0004, "pebbles", cex=0.75, srt=90)
legend(0.01,6,legend=c("Stokes","Turbulent","Ferguson & Church"),
       cex=0.7, lty=c(1,1,1), lwd=c(3,3,3), col=c("blue","darkgreen","red"))
# Plot points from Table 2 of Ferguson & Church
D = c(0.068,0.081,0.096,0.115,0.136,0.273,0.386,0.550,0.770,1.090,2.180,4.360)
W = c(0.00425,0.0060,0.0075,0.0110,0.0139,0.0388,
      0.0551,0.0729,0.0930,0.141,0.209,0.307)
points(D,W, pch=20)

plot2

Comparing Zoltan’s Python code to my R code is a great way to see how each language works.

I’ve been working with and thinking a lot about the silt range of grain sizes lately … look for a post in the future focused on that.

See all my posts about R here.

-

(You probably see an advertisement right below the post. Sorry to clutter up this content with annoying ads, I don’t like it either, but I would otherwise have to pay for hosting. This option is the least intrusive.)

Friday Field Photo #187: Lava Channel on the Big Island

August 9, 2013

DSCF9164

As an Earth scientist who studies surface processes I’m fascinated by channels. Channels are found all over the Earth’s surface, on land and in the deep sea, and we have mapped ancient and active channels on the surface of other planets/moons. As a sedimentary geologist I typically ponder the types of channels that shape landscapes and seascapes, move sediment across the surface, and serve as long-term repositories of deposited sediment.

This week’s photo is indeed a channel but instead of moving sediment it moved molten rock across the surface of the Big Island in Hawaii. It’s now frozen in time, an abandoned channel, perhaps it will be reoccupied by younger lava flows someday. I don’t know the first thing about lava channels, but it seems that it must solidify from the edges inward, resulting in the down ‘stream’ texture (?).

Happy Friday!

-

This photo and more from the Big Island on my Flickr page.

Using R to Import and Perform Simple Calculation on Multiple Files

July 25, 2013

My first post about using R back in January focused on creating custom plots and led to some interaction with other R users in the comment thread and this post from R expert Gavin Simpson in which he kindly offered a solution (and code) to what I was trying to do. I ended up using a slightly modified version of that code to make a bunch of plots for a paper I’m working on. (In fact, I should be working on that paper right now … yep.)

Anyway, for this post I want to briefly summarize another aspect of R (or most programming languages for that matter) that is quite powerful: automating data processing and/or computational tasks that need to be done for multiple files. This funny graphic* below was going around some months ago and very nicely illustrates what I’m talking about.

geeks-repetitive-tasks

Just a bit of background before I get on to the R code. I’m starting some research that will take me and my current/future students at least a few years, perhaps more, to accomplish simply because of the huge number of samples. I have a literal boatload of deep-sea sediment samples from IODP Expedition 342 (see this post for context) that will be used to better understand the history of deep ocean circulation in response to past climate change. Specifically, we’ll be measuring the variability in grain size of the terrigenous (land-derived, non-biogenic) sediment through time. To the naked eye, all this sediment is mud. But, subtle differences in the mean size of the silt fraction over time can tell us about the relative intensity of long-lived abyssal currents that transported the sediment. (If you want to know more about this approach, including all the assumptions, limitations, and caveats, this is a nice review.)

Read more…

Friday Field Photo #186: Coastal California Turbidites

July 12, 2013

DSCF9895

It’s been two years since we moved from California to Virginia, time flies! We really love where we live now, but there are times where I do miss California. I especially miss exploring the outcrops along the coast. Many of these wave-washed bluffs have gorgeous exposures of bed-scale sedimentary geology, a real delight for those who like sedimentary structures.

The above photo shows some Paleocene turbidites along beach cliffs in the town of Gualala, a few hours north of the Bay Area.

Happy Friday!

Friday Field Photo #185: Las rocas del otro lado del río

June 28, 2013

riozamora

This photo is from the past field season in Patagonia. We wanted to get to these rocks on the other side of the river. My student eventually got there … by going ~30 km downstream to where there was a proper bridge and then driving on some gnarly roads for a couple hours. If only we had jet packs. Drones are all the rage these days … drones shmones, I want jet packs!

Happy Friday!

Summertime research in the lab

June 12, 2013

samples

Summer is in full swing and this summer is all about lab work. An undergraduate researcher and I are in the midst of extracting the terrigenous (land-derived) sediment from marine sediment. We are interested in the grain-size and compositional characteristics of the terrigenous component to better understand the history of a long-lived oceanic current that transported the sediment. Isolating the terrigenous material means we have to get rid of the other components, namely the biogenic (marine microfossils made of carbonate and opal) and authigenic (metalliferous oxyhydroxides that formed in place) material. We are using this method with some tweaks from a helpful collaborator.

I have >1,000 samples from IODP Expedition 342 (see this post for more about that expedition) that will eventually be processed in this way. But, because we have not done this before and are still getting the lab fully equipped, we are currently using ‘practice’ sediment (a chunk from one of the core catchers) to fine-tune the methodology. That is, when we make a mistake — which is inevitable when learning something new — we won’t be sacrificing a ‘real’ sample. This training will pay off in a couple weeks as we ramp up and start processing samples in batches.

While we work on that two of my graduate students are busy crushing, grinding, sawing, drilling, etc. rock samples for their respective Ph.D. projects. We’ve got mineral separation underway for bedrock thermochronology, sample preparation for stable isotope measurement of carbonate rocks, and thin sections being made.

Summertime!

Friday Field Photo #184: Death Valley Dunes

May 31, 2013

death-valley

Here’s a shot of the sand dunes not far from Stovepipe Wells in Death Valley National Park from this past March. Check out this Flickr set for a whole bunch more. Happy Friday!

Follow

Get every new post delivered to your Inbox.

Join 118 other followers