<< . .

. 36
( : 51)

. . >>

perature, T 0 , and a positive feedback for larger

occur before an incipient scarp is suf¬ciently cation to Mars, here we focus on the excitable
steepened for subsurface dust to be exposed regime.
and the maximum ice-surface temperature to be The solution of Eq. (7.52) in 1D starting with a
narrow pulse of localized sublimation (T = 0.25
reached. The model behavior is not sensitive to
the precise value of „i as long as it is small com- in the pulse, T = 0 elsewhere) is given in Figure
pared with „ f . The threshold temperature for sub- 7.19b. The pulse initiates a self-sustaining train
limation, relative to the equilibrium and maxi- of waves that propagate with the trough trailing
mum temperatures, is between 0 and 1 and varies the temperature wave. The width of the troughs

as a function of latitude. At 75—¦ latitude, the max- is given by » ≈ κ„ f , or 13 km for the parameters
imum summer temperature of a nearly ¬‚at, high- characterizing the Martian polar ice caps. This
albedo slope is about ’10 —¦ C (Howard, 1978). The value is in good agreement with observed trough
temperature difference between equator-facing widths on Mars. This behavior contrasts with that
troughs and nearby ¬‚ats is about 40 —¦ C (Howard, of previous models, in which troughs widened
1978). Using these values, T 0 , or 0 —¦ C, is equal to inde¬nitely without reaching a steady state. An-
0.2. We assume this value applies to the entire other similarity with the observed topography
ice cap for simplicity. Equation (7.52) can be fur- of spiral troughs is trough asymmetry: both in
ther simpli¬ed by scaling time to the value of „ f . the model and on Mars, equator-facing scarps are
With this rescaling, the model is reduced to three about twice as steep as pole-facing scarps. Code
independent parameters: κ = 175 km2 , „i = 0.05, for implementing the 1D spiral trough model is
and T 0 = 0.2. given in Appendix 6.
Equation (7.52), without the directional de- Spiral formation in Eq. (7.52) is illustrated
pendence in f (T , h), is one example of the in Figures 7.20b--7.20e for a circular region
FitzHugh--Nagumo (FHN) equations, a well- with a permanently-frozen region near the pole
studied class of equations that describe excitable and starting from random initial conditions. A
media. The FHN equations have been applied to shaded-relief image of a polar-stereographic DEM
predator-prey systems (Murray, 1990), electrical of MOLA topography is given in Figure 7.20a for
conduction in nerve ¬bers (Winfree, 1987), chem- comparison. In these calculations we have used
ical oscillators (Fife, 1976), and electrical trans- the Barkley (1991) approximations to Eq. (7.52),
mission lines (Nagumo et al., 1965). The FHN equa- which use a semi-implicit algorithm for the cu-
bic function in f (T , h) in order to gain computa-
tions come in several forms that all reproduce
solitary waves in 2D (i.e. h(x)) and spiral waves tional ef¬ciency and investigate spiral evolution
in 3D (h(x, y)) for a range of parameter values. over long time scales. The early-time evolution
Excitable media are composed of two variables: of the model is characterized by the initiation,
an ˜˜activator™™ (T in Eq. (7.52)) and an ˜˜inhibitor™™ growth and merging of spirals. It should be noted
(h in Eq. (7.52)). Triggering or excitation of the that heat conduction in the neighborhood of an
medium occurs when a positive feedback is ini- equator-facing scarp may facilitate the sublima-
tiated above a threshold value of the activator. tion of nearby regions, even if those regions face
Local triggering may initiate the triggering of ad- toward the pole. The late-time evolution of the
jacent zones through diffusive spreading of the model is characterized by continued spiral merg-
activator. The value of the inhibitor increases dur- ing, alignment of troughs to face the equator,
ing excitation and eventually returns the system and a gradual increase in trough spacing as ˜˜de-
to the unexcited state through the negative cou- fects™™ are removed from the system. Mature spi-
pling between h and T in Eq. (7.52). The FHN equa- rals face the equator at low latitudes but curve
tions exhibit both excitable behavior, character- towards the pole at higher latitudes, just as those
ized by wide swings in h and T when „i „f, on Mars do. In the model, this poleward orienta-
and oscillatory behavior, characterized by smaller tion is required for the maintenance of steady-
¬‚uctuations when „i ≈ „ f (Gong and Christini, state, rigidly rotating spiral forms. By developing
2003). Although both regimes may have appli- poleward orientations, high-latitude scarps face


t = 100
t = 10

(d) (e)
200 km


wave-like undulations
t = 1000 t = 2000

Fig 7.20 Numerical model results and north-polar
troughs, bifurcations, and terminations (Figures
topography. (a) Shaded-relief image of Martian north-polar ice
7.20a, 7.20c, and 7.20d). The steady state of the
cap DEM constructed using MOLA topography. The
model after t = 2000 is a set of seven rigidly-
large-scale closeup indicates examples of gullwing-shaped
troughs, bifurcations, and terminations. Highest elevations are rotating spirals. The sense of spiral orientation
red and lowest elevations are green. (b)“(e) Shaded-relief has an equal probability of being clockwise or
images of the model topography, ’h, for (b) t = 10, (c)
counter-clockwise in the model depending on
t = 100, (d) t = 1000, and (e) t = 2000 starting from
the random number seed used to generate the
random initial conditions. The model parameters are
initial conditions. The spiral orientation of Fig-
L = 250 (number of pixels in each direction), x = 0.4,
ure 7.20e is opposite to those on the north pole
T0 = 0.3, „i = 0.05, „ f = 1. In the Barkley approximations to
of Mars, but 50% of random initial conditions
Eq. (7.52), T0 is combined into two parameters, a = 0.75 and
b = 0.01. The model evolution is characterized by spiral
merging and alignment in the equator-facing direction. A
steady state is eventually reached with uniformly rotating
spirals oriented clockwise or counter-clockwise depending
on the initial conditions. For color version, see plate section.
Modi¬ed from Pelletier (2004b).

perpendicular to the equator, thereby reducing
their ablation rates and decreasing their migra-
tion speeds. The decreased migration speed off-
sets the smaller distance required to complete
a revolution at high latitudes, enabling spirals
to rotate at a uniform angular velocity and pre-
serve their steady-state form. The model also re-
Fig 7.21 Vegetation bands of Ralston Playa, Nevada. Dirt
produces a number of secondary features of po-
road provides scale. Individual bands are 1“3 m wide.
lar troughs on Mars, including gullwing-shaped

lead to spirals that rotate in the same sense of the western US and are comprised of brushy, salt-
tolerant vegetation. These vegetation bands most
as those on Mars (i.e. clockwise away from the
likely form by a positive feedback in which vege-
pole). This model behavior suggests that a uni-
tation growth promotes mounding of the surface
form sense of spiraling need not represent an
beneath the vegetation. Mounds trap occasional
asymmetrical process of trough evolution. In-
runoff into the playa, which promotes more vegeta-
stead, it could be the result of a minor asymme-
tion growth. Bands evolve because a line of shrubs
try in the initial pattern of sublimation, as in the
oriented along a contour pond the most water for
model. the least biomass. Develop a model for the forma-
tion of these bands based on a simpli¬ed model for
the positive feedback between vegetation, mound-
ing, and ponding. Actual vegetation bands follow
Exercise contours closely and are more closely spaced on
more-steeply dipping portions of the playa. Is your
7.1 Figure 7.21 illustrates the curious phenomenon of
model consistent with these observations?
vegetation bands. These bands form on some playas
Chapter 8

Stochastic processes

observed statistical distributions of rainfall
8.1 Introduction events, for example, can be used to generate syn-
thetic storms that match the statistical behav-
ior of actual storms. In cases where there is sig-
Thus far we have considered models in which
ni¬cant uncertainty in model input parameters,
the future behavior of the system can be de-
Monte Carlo methods can be used to determine a
termined using equations and boundary condi-
range of model outputs corresponding to a range
tions known at some initial time. In many Earth
of input parameters. Monte Carlo methods are
surface systems, however, the future behavior of
based on deterministic models; they work by run-
the system cannot be predicted with certainty,
ning many successive examples of a deterministic
either because the system behavior is sensitive
model using different input parameters sampled
to small-scale processes that cannot be fully re-
from prescribed probability distributions.
solved and/or because there is signi¬cant uncer-
tainty in model input parameters. The climate
system is a good example of a system that de-
8.2 Time series analysis and
pends on small-scale processes (i.e. turbulence)
that cannot be fully resolved in any numerical
fractional Gaussian noises
model. As a result, climate and weather models
are inherently limited in their ability to predict
Climatological and hydrological time series pro-
the details of future climate and weather pat-
vide many good examples of phenomena that re-
terns. Soil permeability is a good example of a
quire a stochastic modeling approach. Time se-
model input parameter with signi¬cant uncer-
ries data for rainfall and discharge at a point in
tainty. Soil permeability depends on the detailed
space are the result of complex processes, includ-
structure and composition of the soil at a range
ing convective instabilities in the atmosphere. Al-
of spatial scales, making an exact determination
though some elements of the climate system can
of permeability very dif¬cult over large areas. As
be quanti¬ed deterministically, convective insta-
such, numerical models that require soil perme-
bilities are very dif¬cult to model more than a
ability as an input (e.g. models for runoff, in¬l-
few days into the future. The occurrence of con-
tration, aquifer recharge, etc.) are limited in their
vective instabilities is one of the fundamental
precision no matter how ¬nely they resolve the
limitations on accuracy of weather forecasts.
underlying physics of the problem.
Figure 8.1 plots the daily mean discharge
In some applications where deterministic
in the Gila River near the town of Winkel-
models cannot predict the future behavior of
man, Arizona from 1942 to 1980. Discharge in
a system precisely, stochastic models can be
this case is plotted on a logarithmic scale be-
useful for understanding the range of possible
cause of its highly skewed distribution, with peak
system behaviors. Stochastic models based on

annual peak

S(f )
S(f ) ∝ f


1 0.1
t 8000 12000
0 100
10 10 10
(days, 1942“1980)
f (1/day)
Fig 8.1 Stream¬‚ow of the Gila River near Winkelman,
Fig 8.2 Power spectrum of the Gila River discharge plotted
Arizona, from 1942“1980.
in Figure 8.1.

values many times larger than average values. because they are controlled by a single meteoro-
Many kinds of datasets that arise in climatolog- logical event. At longer time scales, a seasonal
ical and Earth surface studies involve spatial or cycle can be expected and is apparent in Fig-
temporal series data. Figure 8.1 is a time series ure 8.1. At still longer time scales, interannual
because it is a quantity that varies as a function and decadal variations in ocean temperatures (re-
of time. Spatial series data, such as topographic sulting from known periodic processes such as
pro¬les, can be analyzed using many of the same ENSO and also aperiodic variations) lead to clus-
tools used to analyze time series data. Time se- ters of wet years and clusters of dry years. Hurst
ries data can be characterized by their probabil- et al. (1965) quanti¬ed this multi-scale autocor-
ity distributions and autocorrelation functions. relation pattern in hydrological data sets, now
The probability distribution of a dataset quanti- called the Hurst phenomenon. In this section we
¬es its mean value and its variability from the explore a particular type of stochastic model, i.e.
mean. The Gila River near Winkelman is charac- the fractional Gaussian noise, that reproduces
terized by a mean stream¬‚ow value of approxi- the autocorrelation structure common to many
mately 100 m3 /s. The mean value is not very rep- hydrologic and climatic time series.
resentative of the typical stream¬‚ow in this case, The autocorrelation structure of a time series
however, because the river also has extended dry can be quanti¬ed in a number of ways, but the
periods with discharge values less than 10 m3 /s power spectrum is one of the simplest. The power
and occasional large storm events that trigger spectrum is de¬ned as the square of the coef¬-
discharges greater than 1000 m3 /s. This river, like cients in a Fourier series representation of the
most in the southwestern US, has a large variance time series. The power spectrum quanti¬es the
in stream¬‚ow values due to the episodic nature variance of a time series at different frequencies.
of precipitation in semi-arid environments. The Figure 8.2 plots the power spectrum of the Gila
autocorrelation function of a time series quan- River stream¬‚ow dataset on logarithmic scales.
ti¬es how strongly points in the series are cor- This spectrum shows a prominent seasonal cycle
related over different time scales. In the case of (i.e. a large peak at a period of 1 year) superim-
river discharge, signi¬cant correlations can be ex- posed on a broad spectrum that decreases with
pected over many time scales due to the nature larger frequencies. Aside from the seasonal cy-
of rainfall and hydrologic response in a large cle, the background spectrum follows a power-
law relationship with frequency given by S( f ) ∝
river basin. During an individual ¬‚ood event
f ’1/2 . This relationship turns out to be nearly
that lasts several days or more stream¬‚ow values
can be expected to have signi¬cant correlation universal for stream¬‚ow datasets worldwide. In a

this observed underlying structure. Figure 8.4
plots examples of three synthetic time series with
power spectra given by S( f ) ∝ f ’β , with β = 0,
1/2, and 1. Fractional Gaussian noises are simply
S(f ) time series that have Gaussian probability distri-
(yr) ’1/2
S(f ) ∝ f butions and power spectra have a power-law de-
pendence on frequency f . As the magnitude of
the power-spectral exponent increases, the cor-
relations between adjacent points in the series
become stronger. Fractional Gaussian noises can
be generated with the Fourier-¬ltering technique.
This technique proceeds in several steps. First, a
10’3 10’2 100 101
Gaussian white noise sequence is generated. Sec-
f (1/yr)
ond, a discrete Fourier transform of the Gaus-
Fig 8.3 Average of power spectra of monthly average sian white noise sequence is computed. The re-
sulting Fourier spectrum will be ¬‚at, i.e. β = 0
discharge data from 636 gage stations with the annual
periodicity ¬rst removed.
in the expression S( f ) ∝ f ’β . Except for the sta-
tistical scatter, the amplitudes of the Fourier co-
ef¬cients, |Ym | will be equal. Third, the Fourier
study of 636 gage stations, Pelletier and Turcotte coef¬cients are multiplied by a power-law func-
(1997) found that, when the annual peak was tion of frequency:
removed, the average power spectrum followed
m 2
Ym =
a nearly straight line on log--log scales with a Ym (8.1)
slope of ’0.50 (Figure 8.3). The data for this anal-
The exponent β/2 is used because the power
ysis were chosen from all complete records with
durations greater than or equal to 512 months, spectrum is proportional to the square of the
including 636 records in the analysis. Pelletier Fourier coef¬cients. Fourth, an inverse discrete
and Turcotte (1997) took advantage of the large Fourier transform is computed using the new
number of available stations to investigate the coef¬cients Ym . This technique is called Fourier
possible regional variability of the power spec- ¬ltering because Eq. (8.1) has the effect of re-
tra. All of the regions exhibited the same power- moving, or ¬ltering out, some of the variability
at high frequencies (if β > 0) or low frequencies
law power-spectral form with an average expo-
nent of ’0.52 and a standard deviation of 0.03, (if β < 0). Two-dimensional fractional noises can
indicating little variation. Proxy data for hydro- also be constructed (Figure 8.5). These noises have
logic datasets (e.g. tree ring data in semi-arid en- applications to stochastic models of topography,
vironments) also exhibit that same spectrum out soil moisture, and the porosity structure of sed-
to time scales greater than 1000 years. Pelletier imentary basins. The details of this technique,
and Turcotte (1997) also showed similar power- together with codes for its implementation, are
spectral behavior in precipitation time series at given in Appendix 7.
time scales greater than 10 years, indicating that Fractional noises have many applications in
the correlation structure in stream¬‚ow data are geosciences (Pelletier and Turcote, 1999). For ex-
primarily a result of correlations in the climate ample, generating synthetic fractional noises can
system on decadal time scales and greater, while help to assess the likelihood of future hydrolog-
shorter-term correlations are at least partially the ical events that depend on the cumulative dis-
result of surface and subsurface ¬‚ow pathways. charge over a series of years. A drought of a given
Thus far we have shown that hydrological duration, for example, is the result of several con-
data sets have a common underlying autocorre- secutive years of below-average ¬‚ow. The correla-
lation structure. It is useful for many applica- tion structure of hydrological time series plays a
tions to generate synthetic data sets that match very signi¬cant role in the frequency of drought

b = 1/2
b = 0 (white noise) 4.0



’2.0 ’2.0
(a) (b)
1000 2000
1000 2000
4.0 b = 1


0.0 0.0

1000 2000
1000 2000
on a stochastic diffusion model of convective
Fig 8.4 Examples of fractional Gaussian noises with β = 0,
transport. Fluctuations in precipitation, river dis-
1/2, 1, and 2.
charge, and tree ring growth are related in com-
plex but direct ways to variations in temperature
occurrence. If hydrological time series were com- and/or humidity.
pletely uncorrelated, then the probability of 10 The simplest approach to the turbulent trans-
years of below-average stream¬‚ow would be 1/210 port of heat and water vapor vertically in the at-
or 1/1024. However, the fact that hydrologic data mosphere is to assume that the ¬‚ux of heat or
sets have positive correlation over a range of time water vapor concentration is proportional to its
scales makes the probability of droughts, espe- gradient. The mixing of heat energy and humid-
cially long droughts, far more frequent (Pelletier ity will then be governed by the diffusion equa-
and Turcotte, 1997). tion. Experimental and theoretical observations
support the hypothesis that diffusive transfer ac-
curately models vertical transport in the atmo-
8.3 Langevin equations sphere. Kida (1983) has traced the dispersion of
air parcels in a hemispheric GCM and found the
dispersion of those air parcels to obey the diffu-
If a broad range of climatological and hydrolog-

<< . .

. 36
( : 51)

. . >>