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

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 relatively 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 and geologic surfaces that suggest underlying non-linear equations describe their form. Can you think of some examples?

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

In addition, there are many other geoscience phenomena that behave in a non-linear manner. One simple, but important example is accelerating creep, where a hill side starts moving slowly at first, but then faster and faster until it fails catastrophically. Below we explore some geoscience phenomena that are non-linear in character.

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.

One 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, understand the behavior implicit in the equation. Through this 'play' one can come
to appreciate the character of two types of non-linear behavior. You may take the approach of assuming (trusting)
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 sense. 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 a number progression ( **N(0)**, **N(1)**, **N(2)**, .....**N(n) 0**.
If you wanted to do the even positive numbers you start with zero
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**

Where **k** = to a growth constant.

This can also be expressed as follows to calculate the number at a certain increment in the sequence,

**N****(n)=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**.

Some useful basic equations, transformations, equalities:

- for log 10, x = 10^t, then t = log x (you can change x, y and t assignment around as need be).
- log (x*y) = log (x) + log (y)
- log (x/y) = log (x) - log (y)
- log ( x^n) = n * log (x)
- log N = - log(1/N) and ln N = - ln(1/N)
- y = a * x^n transforms into log y = log(a) + n * log(x), which is a straight line form of y = b + m *x. Note that you would use a log-log transformation on both x and y to get a straight line. This is a power function relationship.
- y =b * n^x transforms into log (y) = log(b) + x * log(n), which also can be cast in a straight line form. Note that you would use a semi-log transformation, i.e. on y but not x. The special case is when n = e = 2.718, and this gets into natural logs (ln).
- ln (e^x) = x
- x = e ^ t, then t = ln (x)

The nature of e.

- For a derivation see Vacher (2000) and references therein.
- It is transcendental.
- When differentiated it reproduces itself.
- e^x = 1 + x/1! + x^2/2! + x^3/3! + x^4/4! ...... keep adding terms (thats why its transcendental). ! -> factorial.
- e = limit as N goes to infinity of (1+(1/N)^N.

__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 (**P****o**) 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. In this case, **y **= **N(0)** and **x = t**, and **m =** **(ln (2)/t(1/2))** and **b =** **ln (P****(o))**

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

Then taking the antilog of both sides yields the more general form:

**P = P****o * 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 with precision 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 resulting age 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 .

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.

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.

- Richter scale
**m**- based on corrected amplitude of the**first p-wave arrival**. This is the most common type of body wave magnitude. A size 6 is 10 times as large seismometer response amplitude wise as a 5 and one tenth the amplitude of a 7. - Surface wave magnitude is based on
**Rayleigh wave amplitudes**.**M****s**= log10 A +1.656 * log10(ed) +1.818, where ed = epicentral distance in degrees and A is the amplitude. **m**approximately = .56 ***M****s**+ 2.9- Seismometer amplitude response may be linear equated with accelerations, but not with energy release. Calibrated by explosions the empirical estimate is that energy released = E = 10^5.24 * 10^(1.44*Ms)
- a unit increase in Ms is a 27.5 fold increase in energy release, and about
a 32 fold increase for the Richter scale!
**Interestingly, 80% of the worlds seismic energy is released by the small fraction (<<<1%) that are above 7!!.** - a more accurate measure is the
**moment magnitude**= Mw = [log10(rupture area * average displacement * shear modulus of rock)/1.5]-10.73. This derives from basic physics, and is the theoretical soundest measure of earthquake size. Note that it requires more information for computation. - the Richter magnitude is still generally reported because of broad familiarity with it, simplicity and long history of use.

Isostatic rebound occurs where a load has been taken off the crust of the earth, and the crust, acting as if it is floating, moves upward to come into a new equilibrium position. It does so both slowly and in a non-linear fashion.

The following formula was taken from Davies (1999):

**u = u****o * e^(-t/Ý)** where **Ý = delta p *
g * R / n****m**

and where:

- u = the amount of remaining uplift
- uo = original amount of depression.
- t = time elapsed since unloading.
- delta p = the difference between mantle and continental density .
- R = radius of depression.
- g = gravitational acceleration.
- nm = viscosity of mantle.

Fundamentally it is exponential equation where x is time and y is the amount of remaining uplift, which gets smaller and smaller with time as the pressure disequilibirum at depth that is driving the isostatic response diminishes.

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

**dh/dt = (p****m * 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.__

As the oceanic or any type of floating crust cools its density changes, and thus it will float at a different level. It cools because it starts out as recently crystallized igneous material, and then the heat flows out of it with time. Thus, the profile of the oceanic spreading ridges (with all the bumps and warts taken out) could be thought to reflect a simple history of conductive cooling with time. In words, the basic equation for equilibrium heat flow is that the heat flow is equal to the temperature difference in the direction of heat flow at two points aligned in the direction of heat flow times the thermal conductivity all divided by the distance between the two points.

Combining equations for isostasy and for heat flow for a spreading ridge allows one to model the large scale topogrpahy of the ocean ridge. An equation, taken from Turcotte & Schubert (1982), that describes this model is as follows:

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

This can be shortened to:

**Y depth = k * {(k****m *
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.

Afternote: While not within the scope of this course to pursue in detail, it is interesting to note that this equation does a very good job of describing seafloor topography, there are some consistent departures of reality from the model. For example, the area around the spreading ridge is consistently lower than the equation would predict. Can you think of why? An additional consideration is the phenomena of **non-unique solutions** in geophysics. Could a different physical model produce the same results? While this is a permissible model given constraints, sometimes there are other permissible models that could also explain the reality.

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 * p****m
* 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).

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.

**Here are some other examples:**

**Allometric growth in organisms**: The basic idea is that the "relationship between the size of one part of an organism and the size of another is not linear but follows the general equation y = B * x^a." (Ferguson, 1988). This is apparent for children. This would be a power law function.**Crater/impact flux rate with time**: exponential decay?**Stream profiles and hydraulic phenomena**(Mayer, 1990):**Tertiary or accelerating creep**.**Stream branch development with time**.**Sediment compaction, dewatering and density changes**.

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

Exercise 10

**References:**

- Ferguson, J., 1988, Mathematics in Geology, Allen & Unwin, London, 299 p.
- Mattox, S. R., 1999, An Exercise in Forecasting the Next Mauna Loa Eruption: Journal of Geoscience Education, v. 47, p. 255-260. Describes use of exponential probability equation in the context of predicting volcanic eruptions.
- Mayer, L., 1990, Introduction to Quantitative Geomorphology, Prentice Hall, 380 p.
- Sharma, P. V., 1986, Geophysical Methods in Geology, Elsevier, 442 p.
- Shea, J. H., 2001, Teaching the mathematics of radiometric dating: Journal of Geosciences Education, v. 49, p. 22-24.
- Turcotte, D. & Schubert, G., 1982, Geodynamics, Wiley, 450 p.
- Vacher, H. L., Computational Geology 9 The Exponential Function; Journal of Geoscience Education, v. 48, p. 76.

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.