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.

sylvester_covault_submarine_channel

Animation of incising submarine channels with early meander cutoffs.

 

sylvester_covault_submarine_channel_profile

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.

screen-shot-2016-09-16-at-1-13-17-pm

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.

Forest cover change in the Carpathians since 1985

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

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

bend_area

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:

csik

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

hargita

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.

Reference

P.V. Potapov, S.A. Turubanova, A. Tyukavina, A.M. Krylov, J.L. McCarty, V.C. Radeloff, M.C. Hansen, Eastern Europe’s forest cover dynamics from 1985 to 2012 quantified from the full Landsat archive, Remote Sensing of Environment, Volume 159, 15 March 2015, Pages 28-43, http://dx.doi.org/10.1016/j.rse.2014.11.027.

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]);

Screen Shot 2015-02-02 at 9.03.19 PM

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');

diffusion with FTCS scheme

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();

difference between analytic and numerical solutions

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]);

Laasonen scheme

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 two time steps

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

Diffusion with Laasonen scheme

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

Case when velocities are nonzero on both sides

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]);

Crank-Nicolson scheme

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

Diffusion with Crank-Nicolson scheme

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

Screen Shot 2015-02-06 at 5.13.34 PM

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

Fault scarp diffusion

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

Fault scarp diffusion over long time

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

Fault scarp diffusion, extended domain

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.

Further reading

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:

Screen Shot 2013-08-09 at 5.54.33 PM 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

Screen Shot 2013-08-09 at 5.54.40 PMThe 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):

Screen Shot 2013-08-09 at 5.54.47 PM

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

settling1

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

settling3

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:

Screen Shot 2013-08-09 at 5.49.31 PM

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

Screen Shot 2013-08-09 at 5.34.42 PM

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

Screen Shot 2013-08-09 at 5.49.08 PM

At terminal settling velocity, the particle Reynolds number is

Screen Shot 2013-08-09 at 5.59.14 PM

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

cdvsre

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.

A new way to enjoy photographs

Anybody who takes more than ten photos per year (and everybody has at least a point-and-shoot camera these days) needs a good online photo sharing service. I have been a diehard fan of Smugmug for several years now. I love the elegant, somewhat Apple-like interface, the slick animations, the ability to easily organize and tag photos, the fact that pictures can be displayed at seven different sizes, that it is easy and fast to order high-quality prints, not to mention the significant integration with Google Maps and Google Earth. While I have realized that Flickr seems better equipped for more ‘Web 2.0’ interactivity (maybe largely due to the sheer number of users and photographs), and that there are far more photographs of turbidites on Flickr than Smugmug , I find the Flickr user interface confusing and its design inferior to that of Smugmug, with a lack of style that does not do justice to the zillions of great photos that are out there on the servers.

Having said that, I have recently started to use and appreciate Flickr a lot more. The reason: the updated Apple TV can stream photos directly from Flickr. Television sets with high-definition screens might be a bit ahead of the time due to the limited number of easily (and cheaply) available HD TV programming and movies, but they are perfect for displaying even relatively low-resolution photographs in brilliant colors and surprising clarity. After all, the best HDTVs have a pixel count of 1080 x 1920, and you get more than two megapixels with most digital cameras. Sitting down with a glass of wine and discovering good photographs on a big screen while listening to music is my favorite new pastime and I think it is a lot more enjoyable than browsing photos on much smaller computer screens that usually have a lot of clutter in addition to the photograph.

Now, if Apple was smart and kind enough to put Smugmug on Apple TV as well…

iPod-kritika a la Cluj

Először is, elnézést a kötekedésért. De. A kolozsvári Szabadság cikkezik az új iPod(ok)-ról, és az egy mondatra eső hibák száma kiemelkedő.

A cím: “Újra tarol az iPodE”. Helyesen: iPod.

Első sor: “A Sony Walkmen találmánya volt az, ami alapvetően forradalmasította az utcai zene hallgatást.” Helyesen: Sony Walkman.

“A leghíresebb mp3 lejátszót az Apple készítette el. Az ipod elnevezést azóta generikusan mp3 lejátszókra is használják, akkora sikernek örvendett a Mac zenedoboza.” Az iPod elnevezést kimondottan csak az Apple által gyártott lejátszókra használják, legalábbis angol nyelvterületen.

“A siker receptje két tényező volt: a jól sikerült design és a kiváló minőségű hang visszaadási képesség.” A hangminőség valóban nem rossz az iPod-on, de nem is jobb, mint sok más lejátszón, és gyengébb, mintha CD-t hallgatnál. A siker receptje nem a hangminőségben keresendő.

“Az alma logos lejátszó kezelhetősége azonban messze elmarad a flash lejátszók mögött: mert ezekhez még az Apple által erőltetett iTunes program szükséges.” Igen, az iTunes szükséges az iPod-hoz, de ez nem jelenti azt, hogy az iPod kezelhetősége elmarad a flash lejátszókéhoz képest. Az egyik lényeg az iPodnál az, hogy simán és intuitíven működik az iTunes-zal. Az Apple-nél mindig is fontos volt a szoftver-hardver integráció.

“…a probléma csak az, hogy az Apple szoftvereire jellemzően az iTunes is nehézkes, sőt agresszív szoft.” Bárcsak minden szoftver annyira lenne nehézkes és agresszív, mint az iTunes.

“A második shuffle modell is a külalakjával hódította meg a piacot (pedig a 100 eurós lejátszónak még kijelzője sincs).” Amerikában az iPod Shuffle-t 60 dollárért árulják. Nincs kijelzője, de egyszerűen és nagyszerűen működik. És akkora, mint egy ütés tapló.

“Az iPhone megint csak olyan ketyere, amivel az Apple a mások által kifejlesztett megoldásokat fölözi le.” Az Apple évek óta dolgozik az iPhone érintésérzékeny képernyőjén és annak szoftverén. Senki eddig nem tudta ezt ilyen szinten megcsinálni.

“A funkciók közt szerepel a netkompatibilitás, mp3 illetve videolejátszás, telefon, GPS, érintős képernyő – mindez egy maximálisan 20 gigabyteos merevlemezzel kombinálva.” Az iPhone-on nincs (még) GPS (sajnos). És a merevlemez csak 8 GB-os (sajnos).

Lehetne egyebekkel is kötekedni, de minek.

Diszkléjmer: alulírott megrögzött Apple-barát és boldog iTulajdonos.

Why is photographic metadata not treated properly?

I am looking forward to the day when there will be a good industry standard for adding metadata to digital photographs. Right now, it is a real mess, especially if you want to add new information to the photos that have already been downloaded to the computer.

More to the point: as I mentioned before, I have a small GPS receiver that I like to use whenever I am doing geological field work or I am just hiking (is there a difference? 🙂 ). In addition, I cannot do any of these things (hiking or fieldwork) without having a camera with me and taking pictures all the time. And, when I get home, I try to add the approximate geographic location to at least a number of photographs. This can be done using a variety of software packages that write the geographic coordinates into the EXIF part of the image file, the part that contains all kinds of information about the picture, such as date taken, camera manufacturer, camera model, exposure time, etc.

To combine the geographic information from the GPS with the photos, I am using GPSPhotoLinker, which is by far the best application among those that I have tried (and I have tried many). The nice thing about it is that you don’t have to process the photos one by one; the program has a batch mode that does everything in one go, using the time stamps in the photographs and the time data from the GPS unit. [There are other programs that claim to be able to do this, but GPSPhotoLinker actually works.]

In my case, the problems start when I try to edit the pictures with iPhoto. I am a big fan of iPhoto, I think it is a fantastic photo management program, but apparently it is unable to properly deal with the modified EXIF data. If I embed the geocoordinates before taking the pics into iPhoto, they show up correctly in iPhoto, but then they (at least some of them) get corrupted when I export the images from iPhoto (usually to put them on Smugmug). So the only option is to geotag the photos after they have been edited in and exported from iPhoto. Fair enough, I can live with that. The problem is that with the new version of iPhoto (iPhoto ’08), GPSPhotoLinker cannot write the coordinates to the exported images.

Scheisse mare, as my friend Radu would say (“mare” means big in Romanian, if you want to know). The only workaround I have found is to open the exported images in Adobe Photoshop, and save them again as JPEG files; after this operation – that probably rearranges the EXIF data – GPSPhotoLinker works fine.

Part of the issues are probably rooted in the fact that the EXIF part of the image can be located anywhere in the file. There is a good reason why the longest section of the the Wikipedia article on EXIF is entitled “Problems”.

I hope that things will improve soon. In the meantime, if you have a good workflow or workaround for doing automated geotagging on a Mac, please let me know.

A média jövője jelene


Apróság, de muszáj megemlíteni: azt írja a július hetedikei Szabadság, a “Tudományos-fantasztikus: a média jövője” című cikkben, hogy

A technika-guruk jóslatai szerint nemsokára a mobiltelefon lesz az a műszer, amelyről egyszerre lehet majd videót nézni, zenét hallgatni, internetezni, chattelni – vagy éppen saját tartalmat készíteni és azt feltölteni a világhálóra.

Nemsokára? Ugyan-ugyan. Ez nem a média jövője, és nem tudományos-fantasztikus elmélkedés, hanem valóság.

A cikkben emlegetett mobiltelefonra egy lehetséges példa az iPhone. Ráadásul már az iPhone előtt is voltak “okos” telefonok, amelyek mindenre képesek, amik a cikkíró szerint csak a technika-guruk jóslataiban léteznek. Az más kérdés, hogy engem igazán csak az iPhone óta kezdett a mobiltelefon-téma igazából érdekelni. Ha idővel még GPS-t is tesznek rá, a Google Maps mellé, akkor valóban ez lesz az egyik legizgalmasabb és élvezetesebb szerkentyű.

Tíz ok, ami miatt Macváltozott az életem

Az utóbbi időben amolyan multikulturalizmus alakult ki körülöttem ami a számítógépeket illeti. Nap mint nap három különböző operációs rendszerrel van dolgom. Az irodában egy HP laptopot és egy Linux desktopot kínozok; itthon pedig egy iMac-kel szórakozom, amelyiken a Windows XP-t meg a Unix terminált is elindítom néha. De térjünk a lényegre — itt van tíz ok ami miatt az iMac-en határozottan élvezek dolgozni, a másik kettőn pedig nem mindig. A sorrendnek semmi jelentősége.

1. Csendesebb. Az iMac zajszintje olyan alacsony, hogy legtöbbször alig lehet észrevenni, hogy be van kapcsolva. Ez nem mondható el az HP laptopról. Vagy a többi PC-ről, amivel valaha is dolgom volt.

2. Szebb a külseje. Az iMac olyan, mint egy remekül sikerült kiállítási darab a Modern Művészetek Múzeumából. Az HP pofalemeze tűrhető, de nem nyújt több esztétikai élményt, mint egy műanyag szerszámosdoboz.

3. Update-ek csak ritkán vannak, és ha vannak, akkor egyszerű és gyors az installálás. A gépet általában nem kell újraindítani. Ha elindítom a Windows XP-t az iMac-en, szinte mindig talál egy néhány update-et, és a letöltés után állandóan agitál, hogy a gépet újra kell indítani.

4. Szebb az interfész (van erre jobb vagy szebb magyar szó?). A Mac OS X – ben minden minimalistán elegáns. Tudom, hogy egy számítógépnél nem ez a legfontosabb dolog, de miért ne lenne szép, ha mindennap dolgunk van vele.

5. Kikapcsoláskor és indításkor kb. háromszor olyan gyors, mint az HP laptop (tudom, hogy nem egészen jogos ez az összehasonlítás, de akkor is).

6. Simán és minden bökkenő nélkül együttműködik az iPoddal. Ezt nem mondhatom el az iMac-et megelőző PC-ről (OK, ez lehet, hogy részben az Apple hibája). A zenehallgatás teljesen más műfaj lett mac-esedésem óta: az iPod-jelenség mellett drót nélkül lehet a zenét az iTunes-ból a nappali szoba hangszóróiba küldeni.

7. Stabilabb. Program-lefagyások az iMac-en is előfordulnak, de a vérnyomásom általában kevésbé veszélyeztetett, ha a Mac-en dolgozom.

8. Az iPhotóval fotós igényeimnek 99%-át el tudom intézni; a Photoshop Elements-et csak ritkán indítom el. Az iPhoto ráadásul gyors, szép, és elegáns. A teljes-képernyős képszerkesztés is sokat megér.

9. A Mac OS X-be épített kereső, a Spotlight, össze sem hasonlítható a Windows XP keresőjével. Tudom, létezik Google Desktop és a Windows Vistának van új keresője. A Google Desktop-ot próbáltam, és lassúnak találtam. A Vista-t nem próbáltam, és valószínűleg nem is fogom egyelőre.

10. Egyszerre tudok rajta Mac, Windows, és Unix programokat futtatni. Erre tulajdonképpen ritkán kerül sor, de maga a tudat nem rossz, hogy gyakorlatilag bármi működik.

Geologic maps in Google Earth

One of the best ways to really start understanding the geology of an area is to look at the geologic map in Google Earth. Of course, unless you are interested in San Francisco or some other top notch place, you will not find the geologic map available in kmz or kml format (although you can get the whole US geological map here, and they are working on covering most of the globe).

Before that happens however, with a little patience it is possible to draw your own maps in Google Earth. You can use an image overlay as a starting point, and draw polygons on it after you managed to position it properly. As a quick test (well, actually it took me almost one day to do it), I created a small map that covers part of the southern East Carpathians in Romania, an area where I did some work for my thesis. It is based on the Geologic Map of Romania, 1:50000 scale, that is, one sheet from the series, edited by Murgeanu et al. and published in 1968. Old stuff, but good stuff. And a lot of work.