t = 4.8 yr concentration vs. depth

concentration vs. depth 15

10

(a) (b)

diffusivity vs. surface age

Qa7

0.15

D

(cm2/yr)

0.10

Qa6

0.05 Qa3

Qa5 Qa4

0.00

1 10 100

age (kyr)

(c)

Fig 2.13 (a) Diffusion-model results for normalized where C 0 is the initial concentration within a

concentration vs. depth from t = 10 to 100 yr with thin surface layer of thickness dw , z is depth

D = 0.1 cm2 /yr. (b) Cumulative 137 Cs concentration as a

in the soil pro¬le, D is the diffusivity (called

function of depth in four different soil types at Harwell,

D in this section to emphasize that it is not a

England (data from Gale (1964)), along with best-¬t results

hillslope diffusivity), and t is time following de-

for the diffusion model. (c) Plot of diffusivity values vs. surface

position. Equation (2.49) also approximates the

age for the 137 Cs pro¬les on Fortymile Wash alluvial fan.

permeable soil layer as semi-in¬nite. This is an

Horizontal error bars represent the range of age estimates

accurate approximation for man-made fallout

for each surface based on Whitney et al. (2004). Vertical bars

represent the standard deviation of diffusivity values for all of pro¬les, even in a soil with a calcic horizon

the samples collected on each surface. at depth, because radionuclides do not pene-

trate more than several cm into desert soils over

decadal time scales. It should also be noted that

is given by the solution to the diffusion equation

radioactive decay need not be considered explic-

in a semi-in¬nite medium with no-¬‚ux bound-

itly in this analysis because decay does not af-

ary condition at the surface and a fallout mass

fect the relative radionuclide concentration at

C 0 dw input at z = 0 at t = 0 (Carslaw and Jaeger,

different depth intervals, only their absolute val-

1959):

ues. Figure 2.13a gives an example of the diffu-

1 sion model (Eq. (2.49) using D = 0.1 cm2 /yr and

2 /4D t

e’z

C (z, t) = C 0 dw √ (2.49)

πDt t = 10--100 yr).

48 THE DIFFUSION EQUATION

137 Cs

Table 2.1 data and inferred diffusivity values for geomorphic surfaces on the Fortymile Wash alluvial

fan

total 137 Cs (pCi/g) D (cm2 /yr)

Sample ID fraction at 3 cm Geomorphic Unit

Cs-071802-A 0.313 0.8274 0.047 Qa6

Cs-071802-B 0.271 0.5387 0.165 Qa7

Cs-071802-C 0.258 0.8100 0.052 Qa3

Cs-071802-E 0.208 0.7644 0.063 Qa3

Cs-071802-G 0.118 0.9593 0.021 Qa3

Cs-071802-I 0.388 0.9614 0.021 Qa4

Cs-071802-J 0.155 0.6387 0.108 Qa6

Cs-071802-K 0.340 0.9558 0.022 Qa5

Cs-071802-N 0.218 0.9082 0.031 Qa5

Cs-071802-P 0.245 0.9428 0.024 Qa4

Cs-071802-Q 0.205 0.9951 0.011 Qa5

Cs-071802-R 0.237 0.9578 0.021 Qa3

Cs-071802-S 0.285 0.8807 0.037 Qa3

Cs-071802-BB 0.321 0.7943 0.056 Qa5

Radionuclide concentrations measured in the Table 2.1 summarizes the data and results

¬eld are bulk measurements rather than point from Fortymile Wash. To calculate D , the fraction

measurements because the concentration within of total activity at 3 cm was ¬rst computed by di-

different depth intervals is measured. For the viding the activity from 0--3 cm by the total activ-

purposes of extracting model parameters from ity from 0--6 cm. Equation (2.50) was then used to

bulk measurements, it is most accurate to rep- infer the value of the error function argument,

√

resent measured data cumulatively as the frac- equal to z/ 4D t, corresponding to the fraction of

tion of total concentration to a given depth. To activity at 3 cm after 50 yr of diffusion following

compare the diffusion model predictions to this nuclear testing. A table of calculated error func-

normalized cumulative curve, it is necessary to tion values was used for this purpose. The value

integrate Eq. (2.49) to give of the error function argument was then used to

solve for D (column 3 in Table 2.1). Figure 2.13c il-

z

C (·, t) z

d· = erf √ lustrates the relationship between diffusivity val-

(2.50)

C 0 dw 2 Dt

0

ues and geomorphic-surface age on the Fortymile

Wash fan. The vertical bars are the standard de-

Figure 2.13b compares the diffusion model pre-

viation of D values obtained on different pro¬les

diction to 137 Cs pro¬les measured during a 4.7-yr

of the same soil-geomorphic unit. The horizontal

¬eld experiment at Harwell Air¬eld in Oxford-

bars correspond to the age range for each unit

shire, England (Gale et al., 1964). Although these

based on the available age control (Whitney et al.,

data are from a very different climatic regime

2004). The plot exhibits an inverse relationship

than Fortymile Wash, the results illustrate the

between diffusivity and age, with values increas-

precision of the diffusion model in a way that

ing slightly on Qa3 (the oldest surface) relative to

pro¬les in desert soils cannot show because shal-

those on Qa4.

low penetration limits the spatial resolution of

Radionuclides may be redistributed within

desert-soil pro¬les. To obtain the results in Fig-

the soil pro¬le through entrainment and re-

ure 2.13b, the measured 137 Cs pro¬les were ¬t to

deposition by in¬ltrating runoff and/or by

Eq. (2.50). Diffusivity values are lowest for clayey

bulk-mixing processes including freeze/thaw,

and calcareous soils at Harwell, increasing in

wetting/drying, and bioturbation. The rate of

value in sandy soils and peat.

2.2 ANALYTIC METHODS AND APPLICATIONS 49

transport by in¬ltration can be expected to corre- The diffusion equation with radial symmetry

late positively with hydrologic conductivity, and is given by

hence negatively with surface age based on the

‚h 1‚ ‚h

= κr

results of Young et al. (2004). As such, we can con- (2.52)

‚t r ‚r ‚r

clude that transport of radionuclides during in-

If κ is a constant, Eq. (2.52) becomes

¬ltration is a signi¬cant component of transport.

For long time scales, the ¬nite depth of soil

‚h ‚ 2h 1 ‚h

=κ + (2.53)

pro¬les must be considered. In such cases, the

‚t ‚r 2 r ‚r

concentration of radionuclides in the soil pro¬le

The solution to Eq. (2.53) with an initial radi-

will be given by the series solution

ally symmetric topography f (r ) and a constant-

∞

elevation boundary condition h(a, t) = 0 at r = a

dw 2 1 nπdw

C (z, t) = C 0 + sin

π

L n L is given by

n=1

∞

2 J 0 (±n r )

nπ z ’n2 π 2 D t/L 2 2

e’κ±n t

h(r, t) = 2

— cos e (2.51) 2

a J 1 (±n a)

L n=1

a

— r f (r )J 0 (±n r )dr (2.54)

where L is the depth of the soil pro¬le. In some 0

cases, L may be the depth to an impermeable soil

where ±n , n = 1, 2, . . . are the positive roots of

layer in the subsurface. In other cases, it may be

J 0 (±a) = 0 and J 0 (r ) and J 1 (r ) are Bessel functions

more appropriate to take L to be the depth of

of the ¬rst and second kind (Culling, 1963). This

the wetting front (i.e. the maximum distance to

solution utilizes two boundary conditions. First,

which radionuclides are transported by in¬ltra-

at r = 0 Eq. (2.54) implicitly assumes

tive transport).

‚h

=0 (2.55)

‚r r =0

2.2.10 Evolution of cinder cones

Volcanic cinder cones often have a high degree because Eq. (2.54) is a series comprised of even

of radial symmetry and are composed largely of functions of radius (i.e. J 0 (r )). The Bessel func-

tion J 0 (r ) has a derivative of zero at r = 0, so

unconsolidated fallout tephra. As such, solutions

to the diffusion equation in polar coordinates a series comprised of a sum of these functions

can be used to model their evolution. Cinder must also obey that condition. The second bound-

cones are particularly important landforms for ary condition is a constant elevation of zero at

r = a. For simplicity we assumed the elevation

this type of analysis because their initial con-

ditions can often be precisely inferred by mea- of the surrounding alluvial ¬‚at equal to zero,

suring the dip of subsurface bedding planes on but any constant value could be used. In hill-

opposite sides of the crater rim. The initial mor- slope geomorphology, a ¬xed-elevation boundary

phology of pluvial shoreline and fault scarps, in condition applies to a channel that is capable of

contrast, cannot be inferred with comparable pre- transporting all of the sediment delivered by the

cision because the initial mass-wasting phase of hillslope, resulting in neither deposition nor ero-

the scarp typically erases any memory of the ini- sion at the channel location. In the cinder cone

case, a ¬xed elevation at r = a could correspond

tial shape. In addition, 40 Ar/39 Ar dating enables

the eruption age of cinder cones to be precisely to a channel that wraps around the base of the

measured, and hundreds of cones in the west- cone, but this is not a common occurrence. In

ern United States have been dated in this way. most cases, volcanic cones are surrounded by al-

The framework for analyzing volcanic cones with luvial ¬‚ats that do not readily transport mate-

analytic solutions to the diffusion equation was rial from the base of the cone. Equation (2.54)

described by Culling (1963) and Hanks (2000). To can still be used for these cases, but a must be

date, however, no analytic solution has been com- chosen to be much larger than the radius of the

pared to speci¬c natural cones. cone. This way, debris will be removed from the

50 THE DIFFUSION EQUATION

The three integrals that appear in Eq. (2.57)

1

kt/a2 = 0

kt/a2 = 0.0001 can also be written as in¬nite series. In practice,

kt/a2 = 0.0002 however, it is more accurate to evaluate those in-

0.8 kt/a2 = 0.0004

tegrals numerically because the series converge

kt/a2 = 0.0128 very slowly. Therefore, although closed-form ana-

0.6

h lytic solutions for volcanic cones can be written

h0 down, calculating and plotting the solutions re-

0.4

quires some numerical work.

First, we must solve for the roots of J 0 (±a) = 0.

0.2 This is done numerically using a root-¬nding

technique (Chapter 9 of Press et al. (1992) pro-

0 vides many techniques for doing this). The Bessel

0 0.33 0.66

function has an in¬nite number of roots, so how

r/a

many do we need to compute? For our purposes,

Fig 2.14 Plots of analytic solutions to the 3D diffusion the ¬rst one hundred roots are adequate, but

equation for a volcanic cinder cone for a range of times more roots may be required for high-resolution

following eruption.

pro¬les or very young cones. Second, we must

evaluate the integrals in Eq. (2.57) numerically.

To do this, the algorithms in Chapter 4 of Press et

cone and deposited on the surrounding ¬‚at. In

al. (1992) can be used. Appendix 1 provides a pro-

these cases, the boundary condition at r = a will

gram for computing the integrals and perform-

simply serve to maintain the alluvial ¬‚at or pied-

ing the sums in Eq. (2.57).

mont at a constant elevation very far from the

Figure 2.14 illustrates radial pro¬les of Eq.

base of the cone.

(2.57) for the initial cone and at eight sub-

For a volcanic cone of radius rc , crater-rim ra-

sequent times from κt/a 2 = 0.0001 to κt/a 2 =

dius rr , colluvial ¬ll radius of r f , and maximum

0.0128. Both axes are normalized. The y-axis is

height of h 0 (Figure 2.14), f (r ) is given by

normalized to the initial cone height h 0 . As in

§ « all diffusion problems, the solution can be scaled

r f ’rr

⎪ h0 1 + ⎪

if r < r f

⎪ ⎪ up or down in height with no change in the rela-

⎪ ⎪

rc ’rr

⎪ ⎪

⎪ ⎪

⎨h 1+ if r f ≥ r < rr ¬ tive cone shape. The r -axis is scaled to the model

r ’rr

0

f (r ) = rc ’rr (2.56) domain length a.

⎪h 1+ if rr ≥ r < rc ⎪

⎪0 ⎪

rr ’r

⎪ ⎪

⎪ ⎪

rc ’rr

⎪ ⎪ In the cone™s early-stage evolution, the great-

© ⎭

if r ≥ rc

0 est change occurs at the crater rim. This is not

surprising since this is where the pro¬le curva-

ture is greatest. The crater rim is locally trian-

As a technical point, it should be noted that

gular and for early times the crater evolution

the colluvial ¬ll radius r f must be ¬nite in or-

is broadly similar to moraine evolution consid-

der for Eq. (2.54) to apply because otherwise the

ered in Section 2.2.6. For intermediate times,

boundary condition given by Eq. (2.55) is violated.

the position of the crater rim migrates inward.

Physically, r f corresponds to the width of the

This migration is associated with the additional

colluvium that ¬lls the crater shortly following

advective term in Eq. (2.52) that is not present

eruption. Substituting Eq. (2.56) into Eq. (2.54)

gives

∞

r2

2h 0 J 0 (±n r ) 2rr rr

2

e’κ±n t

f

h(r, t) = J 1 (±n r f ) ’ J 1 (±n rr ) + 1 + rc J 1 (±n rc )

±n J 1 (±n a) r c ’ rr r c ’ rr r c ’ rr

2

a2 n=1

±n rr rf rc

+ r 2 J 0 (±n rr )dr + r 2 J 0 (±n r f )dr + r 2 J 0 (±n rc )dr

2

r c ’ rr 0 0 0

(2.57)

2.2 ANALYTIC METHODS AND APPLICATIONS 51

tine et al., 2005). The model solutions compare

(a) 117° W 116° W

well to the observed pro¬les, and best-¬t morpho-

Yucca

Mtn.

logic ages are approximately 400 and 4000 m2 ,

Nevada

Black Cone

respectively. The ratios of the morphologic age

N Ca

ev if Amargosa

ad or

to the radiometric age provide estimates for the

Valley

Funeral a nia

36.5° N Mtns.

time-averaged κ values for these cones: 5.2 and

3.6 m2 /kyr. These values are at the upper range of

Death

Valley

κ values inferred from pluvial shoreline and fault

Black

Mtns.

Lathrop Wells

N 10 km

36.0° N

scarps in the southwestern US (Hanks, 2000).

This approach enables us to estimate how

much height the cone has lost since its erup-

tion. Figure 2.15 suggests that the crater rim of

(b) 1 Lathrop Wells has been eroded by approximately

Lathrop Wells 8% of its original height, or 9--10 m. Black Cone,

0.8 in contrast, has been degraded by approximately

50% according to this approach.

A ¬nal word about Lathrop Wells should be

0.6 kt = 400 m2

h noted. The dating of this cone has been the

h0 Black Cone

subject of considerable debate because it is the

kt = 400 m2

0.4

youngest volcano in the Yucca Mountain region,

the site of the proposed high-level nuclear waste

0.2

repository. Initial attempts to date the cone ge-

omorphically concluded that the cone was no

0

older than approximately 10 ka based on its com-

0.3 0.4

0.2

0 0.1

r (km) plete lack of drainage channels. However, cinder

cones are comprised of coarse (pebble-to-gravel-

Fig 2.15 (b) Observed (thick line) and best-¬t model

sized) tephra. As such, it is likely that little or no

pro¬les (thin line) for a relatively young (Lathrop Wells Cone)

surface runoff can occur until suf¬cient time has

and a relatively old (Black Cone) cone in the Crater Flat

passed for eolian dust to clog surface pore spaces.

volcanic ¬eld, Nevada (location map in (a)).

Therefore, runoff generation on these cones may

be impossible for tens of thousands of years fol-

in 2D problems. Eventually, the crater is ¬lled lowing eruption. In the absence of any radio-

by diffusion with debris and the late-stage cone metric dates for this cone, diffusion-morphologic

morphology is described by J 0 (ar ) with decreas- dating would suggest that the cone is 20 ka to

80 ka assuming a conservative range of κ values

ing amplitude over time. At this point, the expo-

from 0.5 to 2.0 m2 /kyr. This range includes the

nential time-dependent term in Eq. (2.54) has ¬l-

tered all of the high-frequency components in the now-accepted radiometric age but does not in-

topography and only the lowest-frequency term is clude the original estimate of 10 ka. Of course,

signi¬cant in the series. The inset photo in Fig- all morphologic-age estimates suffer from poor

constraints on κ values. However, since κ val-

ure 2.14 shows the pro¬le of an early Quaternary

volcanic cone in the Cima volcanic ¬eld, Califor- ues can be calibrated on a cone-by-cone basis as

nia. The pro¬le shape is similar to the late-stage shown here, they offer great potential for bet-

ter regional calibration of κ values and improved

pro¬les plotted in Figure 2.14.

Figure 2.15 compares the observed topogra- understanding of their underlying climatic and

phy of two cones in the Crater Flat volcanic ¬eld, textural controls.

Amargosa Valley, Nevada, to best-¬t model solu-

2.2.11 Delta progradation

tions. Lathrop Wells Cone has been radiometri-

cally dated to be 77 ka (Heitzler et al., 1999), The application of the diffusion equation to

while Black Cone has an age of 1.1 Ma (Valen- large-scale depositional landforms (i.e. deltas and

52 THE DIFFUSION EQUATION

alluvial fans) is problematic for two reasons. First, U0

sediment

input