7.7 SPIRAL TROUGHS ON THE MARTIAN POLAR ICE CAPS 185

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

186 INSTABILITIES

(c)

(b)

(a)

t = 100

t = 10

(d) (e)

200 km

termination

bifurcation

gullwing

gullwing

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

EXERCISE 187

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

8.2 TIME SERIES ANALYSIS AND FRACTIONAL GAUSSIAN NOISES 189

annual peak

1000

Q

10

(m3/s)

S(f )

100

S(f ) ∝ f

1

10

1 0.1

t 8000 12000

4000

0 100

10

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

190 STOCHASTIC PROCESSES

this observed underlying structure. Figure 8.4

10’1

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

10’2

be generated with the Fourier-¬ltering technique.

This technique proceeds in several steps. First, a

10’1

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)

N

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

8.3 LANGEVIN EQUATIONS 191

4.0

b = 1/2

b = 0 (white noise) 4.0

2.0

2.0

0.0

0.0

’2.0 ’2.0

(a) (b)

1000 2000

0

1000 2000

0

b=2

4.0 b = 1

2.0

2.0

1.0

0.0 0.0

’1.0

’2.0

(d)

(c)

1000 2000

0

1000 2000

0

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-