Modeling with exponential and power law functions in Geosciences

Lecture index: Review of some basics. / Transformations and e. / Radioactive decay and geochronology. / Population growth/dynamics. / Earthquake energy release and magnitude scales. / Isostatic rebound. / Heat flow and spreading ridge topography. / Thermal subsidence in a basin. / Rate laws for chemical reactions. / Other phenomena with nonlinear behavior. / Exercise 10


Geomorphic surfaces.

Nature is full of form, of geometry. These phenomena lend themselves to being described mathematically, and that description can in turn inform understanding of why that form or gemeotry exists. A simple example is the slope of a talus fan and the angle of repose. These slopes are generally straight and a common angle is 30 degrees, the well known angle of repose. Not surprisingly the mathematical description for the strength of granular material is also commonly a linear one (but not for all geologic materials). However, there are many curved geomorphic surfaces that suggest underlying non-linear equations guide their form. Can you think of some examples?

Curved exfoliation joints (and resultant topographic surfaces) in Yosemite National Park.


Reading:

Vacher, H. L., Computational Geology 9 ­ The Exponential Function; Journal of Geoscience Education, v. 48, p. 70-76. You can skim over the part with the derivation of e, but read the rest carefully please.

Section 4-22 Ocean Floor Topography, in Turcotte, D. & Schubert, G., 1982, Geodynamics, Wiley, p. 181-187. At the undergraduate level this can be challenging reading, but you can still pull a lot of information and understanding out of it.


Review of some basics.

The purpose of this exercise is to learn how to play with two forms of non-linear numerical models for some geoscience phenomena. One approach that can be initially taken is that of plug and chug. By plug and chug I mean that you, as a user, may not understand the mathematical derivation of the models, but that you do understand what the input variables are and can manipulate and use the equations to get answers. Through this 'play' one can come to appreciate the character of two types of non-linear behavior. You may take the approach of assuming they are properly formulated mathematically, and then explore how well they describe the real world (or not), and what their behavioral implications are. Always ask your self when you get an answer whether it makes geologic senses. If the equations does not describe the real world particularly well the problem is most likely in the assumptions (or possibly in the input parameters). If the problem is with the assumptions you may need a better or new equation.

Warning - html doesn't have a lot in the way of formating options, and so you should read these equations carefully (and I apologize for the very nonstandard appearance of many of these equations). However, there is some utility in training yourself to recognize the same equation form expressed in different ways.

What's the difference between power law and exponential?

A power law function follows the form that:

y = m * x^b

where x and y are variables and m and b are constants. The "^" symbol means raised to the power of.

An exponential equation follows the form that:

y = m * a ^ (b * x).

Note the very different positions of x. Basically in a power lay function the exponent is constant, and in an exponential equation x the exponent is the variable x. Both forms of equation share the characteristic of being non-linear in form, although log and ln transformations can make them linear.

This Excel graph shows simple examples of linear, power law and exponential equations all with the same constant - 2. Note how the exponential equation "takes off" in comparison to the others - this is a visual demonstration of the power of exponential growth.

Why are exponential equations so useful?

We are taught early on to count in an additive or arithmetic manner. The progression is one, two, three, ... and so on. The next number in the sequence is obtained by adding one to the present number.

N(n+1)=N(n) + 1

Here the subscripts refer to position in the number progression. If you wanted to do the even positive numbers you start with zero and and add 2 instead.

N(n+1)= N(n) +2, where N(0) = 0.

We also learn about counting in different bases, and thus learn that there is an infinite number of counting systems. The Mayans used one that was base 21. In the arithmetic mode note that amount of increase or change is independent of where you are in the sequence. However, in nature the amount of increase or decrease, is often very much a function of where you are in the number sequence, or more simply put, the amount of change in x can be a function of x at that point in time. The population increase is a function of the existing population times a birth rate. So when thinking of a number sequence that describes a natural phenomena a multiplicative incrementation is more useful. The classic example is population growth, but as you will see below many other examples exist. This fundamental equation can be expressed as:

N(n+1) = N(n) * k

This can also be expressed as follows,

N(n+1)=N(0) * k^n

where n is the exponent (equal to x). So this is an exponential equation. The important thing is that the amount of change from number to number in the progression is dependent on the existing N at any point in the sequence. That is simply a much more natural state of affairs. What exists influences how much change occurs.


Transformations and e.

Some useful basic equations, transformations, equalities:

The nature of e.


Radioactive decay and geochronology.

Exponential radioactive decay by half life: We are familiar with the following type of progression. 1000 -> 500 -> 250 -> 125 -> .... where each number in the progression represents the atoms surviving the decay process for a history of integer units of half life. What is the more general expression?

P = P(o) * (1/2)^(t/t(1/2))

The above equation describes how many atoms are left (P) at time t given an initial number of atoms (Po) and a half-life of t(1/2) . Note that (t/t(1/2)) is simply the number of half lives that have passed.

P = P(o) * (1/3)^(t/t(1/3))

The above equation describes how many atoms are left (P) at time t given an initial number of atoms (P(o)) and a third-life of t(1/3), were a third life is the amount of time after which a third of the original population remains. Note that you could write a similar formula for a tenth-life if you so choose. Note also, that if you rearrange you get the fraction (P/Po) of the original material left with time. The equation works for non-integer half-lives.

The equations for half lives can be rearranged to:

ln P = ln (P(o)) - (ln (2)/t(1/2)) * t.

This has a linear form of y = b - m*x, which can be very useful, and so if your data is transformed and plotted in the right way a straight line can result.

We can then let (ln (2)/t(1/2)) = lambda, which is the decay constant.

Then taking the antilog of both sides yields:

P = Po * e ^(-lambda * t).

Some more general information.

A specific example is the Rb-Sr method. It is in some ways the most straightforward of the geochronologic systems. Rb-87 decays to Sr-87, and the basic descriptive equation follows below.

N(Sr-87 at t) = N(Sr-87 at t = 0) + N( Rb-87 at t) * (e^(decay constant * t) -1)

Here N refers to the number of atoms of a specific isotope species.

Note the problem of initial daughter isotope created by the first term on the right side of the equation. If this was equal to zero or we knew it initial value then a simple measure of the ratio of parent to daughter produce would provide an age. Rarely is such an assumption or such information available. What is done to get around this problem is to normalize the equation by the concentration of a stable isotope that should not have changed with time. In this case it is Sr-86. The equation then looks like this:

N(Sr-87 at t)/ N(Sr-86) = N(Sr-87 at t = 0)/ N(Sr-86) + (N( Rb-87 at t)/ N(Sr-86)) * (e^(decay constant * t) -1)

This will be the basis of one of your three problems for this week. This form of the equation has two advantages. First, it is actually easier to experimentally measure the ratio of the isotopes than the absolute amount. More importantly, the line has the form of a straight line equation if y = N(Sr-87 at t)/ N(Sr-86) and x = (N( Rb-87 at t)/ N(Sr-86)). This line is known as the isochron. The line slope is then as follows:

slope = m = (e^(decay constant * t) -1)

and by rearranging this equation it is possible to solve for t, the time the system was closed (some of the transformatios for e given above may be quite useful). What we need then is portions of the rock which started out with different amounts of initial Sr-86 and thus lie along different parts and help define the line, which is known as an isochron. Different mineral separates from the same rock body are usually used to establish points on what is known as a whole rock isochron. Note we get some additional useful information from this analysis, the initial Rb-Sr ratio, which has important implications as to the nature of the source terrane of melting.

A similar thing is done with K-Ar except it is Ar-40 the radiogenic product is normalized against Ar-36. In U-Pb two different decay series are used. Each technique has its own little twist which one needs to learn. The most important thing to learn, however, is what the number means. A critical question is - what is the geologic event that is being dated? Another phrasing of this question is - at what point is the rock, mineral or other system being dated closed?

USGS website - portal to geochronologic facilities/ capabilities - http://geology.cr.usgs.gov/capabilities/gronemtrac/geochron/geochron.html .


Population growth/dynamics

We dealt with this a bit before. The significance from an environmental and complex perspective should be clear, and is explored in a great variety of forums. What systems may be affected by organism population dynamics? Algal blooms and associated anoxic events in bodies of water would be one example. When trying to understand time aspects of how components of the biosphere can alter surface physical conditions (by changing greenhouse gas concentrations, etc.), the exponential response is also important to understand.


Earthquake energy release and magnitude scales

The down and dirty on the Richter Scale and Seismic Moment Scale (mainly taken from Lillie, 1999).

Earthquake magnitude is an attempt to measure the amount of energy released. Typically it is based on the amplitude of the seismic waves. Results have to be standardized and normalized to a specific type of seismometer at a set distance from the hypocenter and for a hard rock substrate. There are different magnitude scales, all with a log character. A size 6 is 10 times as large seismometer response amplitude wise as a 5 and one tenth the amplitude of a 7.


Isostatic rebound

The following formula was taken from Davies (1999):

u = uo * e^(-t/Ý) where Ý = delta p * g * R / nm

and where:

Another formula (from Sharma, 1986) is that the rate of uplift (dh/dt):

dh/dt = (pm * g * h)/ (2 * nm * k) ,

where k is a parameter determined by the wavelength of distortion of the land surface. What influences (or is hidden in) the k factor?

Another way to get a handle on the rate and amount of isostatic uplift is to look at the gravity anomaly caused by lack of isostatic equilibrium.

A note about equilibrium: You can think of two factors that will indicate whether you expect a natural system to attain equilibrium or not. The first is the response time - how long it takes to achieve equilibrium. Often the system only approaches equilibrium but you can think of effective equilibrium. How long does it take for the system to be in 95% or 99% equilibrium? This response time can be compared against the perturbation time scale. If perturbation time is = to or less than response time then the system rarely or never is in equilibrium. On what time scale do glacial pertubuations occur? For glaciation history and considering Milankovich cycles, major changes in loading can occur on a time scale of thousands to tens of thousands of years. Isostatic response time is in tens of thousands of years. This may suggest that in areas of glaciation the crust spends a lot of time out of isostatic equilibrium.


Heat flow and spreading ridge topography.

Shaded relief image of South Atlantic spreading ridge with colors a function of depth.

As the oceanic or any type of floating crust cools its density changes, and thus it will float at a different level. Thus the profile of the oceanic spreading ridges (with all the bumps and warts taken out) could be thought to reflect a simple hidtory of conductive cooling with time. An equation, taken from Turcotte & Schubert (1982), that describes this process is as follows:

Y depth = [( 2 * pm * am * (Tm-To))/ (pm-pw)]* {(km * x )/(pi * uo)}^.5

This can be shortened to:

Y depth = k * {(km * x )/(pi * uo)}^.5 where k = [( 2 * pm * am * (Tm-To))/ (pm-pw)]

Where:

Y depth = depth below the ridge crest elevation.
pm = mantle density, average value of 3250 kg m^(-3)
pw = water density, average value of 1000 kg m^(-3)
am = coefficient of thermal expansion for mantle rocks. Typical value = 3*10^(-5) degrees K^(-1)
Tm = temperature at base of lithosphere in Kelvin (constant T with time) = 1250 degrees Celsius.
To = temperature at top free surface heat is radiating out at in Kelvin. This is pretty close to 0 degrees Celsius for the bottom ocean waters.
km = thermal diffusivity of mantle rocks. Typically 1 mm^2 s^(-1)
x = distance away from ridge center
uo = spreading rate.

Note that x /uo gives the age of the crust. Actual and modeled topographic profiles for specific ridges show a decent fit (better than 85%), with some consistent differences between the model and the observed. The ridge crest areas are not quite as high as the model would predict, and the extreme flanks are higher than the model would predict. Note, also that when thinking about units you will need to have common units for distance and for time in particular for this equation to work. Perhaps change all the distance units to meters, and time units to seconds when using this equation in your exercise.


Thermal subsidence in a basin

The following is mainly taken from Turcotte & Schubert (1982) which provides derivations. It is a modification of the model for cooling continental lithosphere and basin formation. If lithospheric thickening and cooling (a half-space cooling model) is responsible for the subsidence (and a bunch of other assumptions are acceptable) then this equation may work.

Y subsidence = [( 2 * pm * am * (Tm-To))/ (pm-ps)]* (km * t / pi)^.5

The values are as above in the case for the cooling spreading ridge. The one difference being the ps = density of sedimentary rock (typically 2500 kg m^(-3)).

Check units when using this equation. Note that basically this boils down to Y subsidence is a square root function of time. Note that are a whole slew of assumptions behind this that could be violated in any particular case (see Turcotte & Schubert for particulars). For what type or phase of basin might this work? Perhaps the post-rift faulting phase of an extensional basin.

How can you test how well this model describes basin subsidence? Plot the cumulative basin thickness versus time and see if it follows the appropriate line. There is a nice example in Turcotte & Schubert (1982).


Rate laws for chemical reactions

The rate can be described by the equations below:

C = Co * e ^ -(k*t)

ln C = ln Co ­ (k*t)

where Co concentration at the start of a reaction, C concentration at time t during the reaction, k is the rate constant.

If you want to find k experimentally you can measure concentrations with time and then graphically solve for k. Rates of metamorphic or hydrothermal reactions could be investigated in such a framework. Again, think about equilibrium issues.


Other phenomena with nonlinear behavior

Here are some other examples:

This non-linear behavior shows up quite a bit in the real world!!


Exercise 10


References:


Copyright by Harmon D. Maher Jr.. This material may be used for non-profit educational purposes if proper attribution is given. Otherwise please contact Harmon D. Maher Jr.