Characterization of atmospheric pollution by means of stochastic indicator parameters

Pergamon PII:

Atmospheric Environment Vol. 30, No. 22, pp. 3811-3823, 1996 Copyright © 1996 Elsevier Science Ltd Printed in Great Britain. All fights re-~rved 1352-2310/96 $15.00 + 0.00

S1352-2310(96)00083-0

C H A R A C T E R I Z A T I O N O F A T M O S P H E R I C P O L L U T I O N BY MEANS OF STOCHASTIC INDICATOR PARAMETERS GEORGE CHRISTAKOS and DIONISSIOS T. HRISTOPULOS Department of Environmental Sciences and Engineering, University of North Carolina, CB # 7400, Chapel Hill, NC 27599-7400, U.S.A. (First received 23 August 1995 and in final form 9 February 1996)

Abstract--Deriving indicator parameters for determining regions at risk and reducing uncertainties is an important part of atmospheric pollution characterization. In this paper we introduce regional indicator parameters that characterize contamination levels by means of stochastic analysis. The expressions obtained are general, and can be used for different types and distributions of pollutant concentrations. We calculate indicator parameters that represent excess contamination above a threshold level specified by the environmental standards, and we study certain scaling properties of the level-crossing contours. Regional indicator parameters are expressed in terms of the pollutant probability distribution, and are calculated explicitly for sl~cific cases. Certain properties of the indicator parameters related to scale effects are also discussed. An application involving a sulfate deposition data set is discussed, which provides valuable insight regarding stochastic indicator parameters and their importance in environmental policy and decision making. Copyright © 1996 Elsevier Science Ltd Key word index: Atmospheric pollution, random fields, indicator parameters.

1, I N T R O D U C T I O N

Contaminants are introduced into the atmosphere from a variety of human activities, such as exhaust emissions from transportation vehicles, industry, power plants, nuclear weapons tests, etc. (e.g. Abrahamsen,. 1987; Black, 1991). The design of any remedial structure regarding air pollution caused by anthropogenic and biogenic sources (e.g. Clark, 1980) requires a rigorous contaminant level analysis. Contaminant level analysis is concerned with the study of the spatial and/or temporal distribution of atmospheric pollutants. Such an analysis must account for critical pollution thresholds determined on the basis of health and regulatory standards associated with a specific region, thus leading to an accurate classification of regions into substantially contaminated and negligibly contaminated zones and allowing the effective implementation of remedial strategies (e.g. Journel, 1987; Neely, 1994). Measurements and predictions of atmospheric pollutants (e.g. carbon :monoxide, ozone concentrations, or acid depositions) are subject to uncertainty, and, thus, any environmental decisions based upon such measurements are also subject to uncertainty. For example, in order to establish the spatial dispersion of environmental poIlutants the uncertainties and variabilities of deposition surveys should be taken into consideration (e.g. sulfate depositions, see Haas, 1994). Another example may be found in ozone clima-

tology studies, where the failure to account for the uncertainty may adversely influence the planning and implementation of procedures that aim to either minimize the occurrence of, or to completely mitigate, exceedances of health and environmental standards; also, the effects of uncertainty are particularly important in the determination of the diurnal maximum ozone concentration (Vukovich, 1994). A rigorous contaminant level analysis is most efficiently accomplished by means of the stochastic theory of random fields. Identifying and researching stochastic indicator parameters, analytical methods and other tools for determining polluted regions at risk and reducing uncertainties across broad spatial and temporal scales is an important part of contaminant level analysis. Indicator parameters (sometimes also called index functions) have a variety of applications in statistical decision analysis, economics and mining (e.g. Raiffa and Schlaifer, 1961; Kendall and Stuart, 1977; Rivoirard, 1994). In this work contaminant level analysis is performed using random fields. Stochastic indicator parameters of contaminant characterization which provide quantitative assessment of the uncertainties and the risks posed by polluted regions are introduced, and their physical meaning is discussed. These parameters can be expressed either in terms of probability laws or by means of spatial moments of the polluted region; therefore, they can be calculated from the available measurements of the contaminant

3811

3812

G. CHRISTAKOS and D. T. HRISTOPULOS

process. Relations can be established between the indicator parameters in terms of algebraic or ordinary differential equations. Regional parameters related to contaminant spatial correlations and level-crossing contours are also discussed and their scaling properties are examined. The effect of observation scale on contaminant level analysis is also considered. Practical insight is gained by means of an application studying sulfate deposition data in the conterminous U.S.

2. R A N D O M F I E L D C O N C E P T S I N C O N T A M I N A N T

LEVEL

ANALYSIS

In contaminant level analysis, where one is faced with complex atmospheric phenomena, the goal is to obtain simplified models of these phenomena that lead to predictable quantities. A meaningful simplification process should account for variabilities in the spatial distribution of the contaminant and uncertainties in the available data and physical information. In this context the indicator formalism can play an important role, especially as regards the study of contaminant concentration exceeding certain threshold levels. Let X(s), where s = (Sl, s2 . . . . . s,) is the position vector, be a spatial random field modelling a contaminant process within a domain D. The mean value of X(s) is defined by m:,(s) = E[X(s)] = S xf~(s; x) dx

(1)

0

where fx(s;x) is the univariate probability density function (pdf); the non-centered covariance of X(s) is given by

cx(s, s') = E[X(s)X(s')] =

xx'S

X(s) = ~ [1 - tx(S, z)] dz.

(s; x, s';

(2)

where fxx(s; x, s';x') is the bivariate pdf; and the centered covariance of X(s) is 8x(s, s') = c~ (s, s') -- mx (s) m~(s').

(3)

A spatial indicator random field, I~(s, z), is defined in terms of X(s) as (e.g. Journel, 1983) {~

ifX(s)~>z otherwise

(4)

for each fixed real threshold value z (for simplicity, in the following the symbols X and Ix will be used to denote both the random fields and their realizations). Note that while the X(s) is a function of the spatial variable s, the Ix is a function of s and the threshold z. To illustrate the implementation of model (4) in practice consider a contaminated region D and assume that the spatial random field X(s) represents the contaminant concentration (in %,/~g kg-a, #g m-2, etc.) at a location s within D. The threshold concentration value z is selected on the basis of environmental and

(5)

0

Note the distinctly different ways that the contamination processes X(s) and Ix(s, z) are determined in terms of each other. While the (s, z)-parametered excess process Ix(s, z) is defined as a step function of the s-parametered process X(s) and the threshold z (equation (4)), the X(s) is derived by integrating the complement of Ix(s, z) over all possible values of the threshold level z (equation (5)). The contaminant process X(s) is, in general, a continuously valued function, but the excess contaminant process Ix(s, z) is a discretely valued function. The equality Ix(s, z) = Ix(s', z') of the indicator field values is necessary given the equality of the random field values X(s) = X(s'). However, the converse is not true, that is, equality of the indicator field values is not sufficient for equality of the random field values. Furthermore, by virtue of equation (4) for any nonnegative-valued contaminant concentration field X(s), Ix(s, z) <~X(s)/z

o o

Ix(s, z) =

sanitary requirements, and the spatial indicator random field Ix(s, z) models the excess contaminant process, i.e. the part of the contaminant that exceeds the threshold level z. From a stochastic analysis viewpoint these random fields constitute a random field pair {X(s), Ix(s, z)}, in the sense that their stochastic characteristics are mutually dependent (Christakos and Hristopulos, 1996). For example, an immediate consequence of definition (4) is Ix(s, z) = P[-X(s) > zlx]. The meaning of this relation is that the excess contaminant Ix(s, z) is equal to the spatially dependent probability that the original contaminant process X(s) exceeds z at location s given that X(s) = x. It also follows that

(6)

for all s~ D. Hence, the ratio of the contaminant concentration field over the threshold level is an upper bound of the excess process. If the random contaminant field X(s) is stochastically homooeneous, the univariate probability distribution is independent of the location s, i.e. Fx(s, z) = Fx (z), and the bivariate distribution depends only on the space lag h = s - s', i.e. Fxx(s, z; s', z') = Fxx(h; z, z'). It follows that its mean mx is also independent of s and its noncentered covariance depends only on the space lag h = s' - s, i.e. c~(s, s') = cx(h = s' - s); similarly, the centered covariance is given by Yx(h)= cx(h) - m2. Lastly, the indicator random field Ix(s, z) of a homogeneous random field X(s) is also homogeneous. Assuming that the atmospheric random field X(s) is eroodic in the domain D, by definition the statistical moments of X(s) can be calculated in terms of sample means over D. For example, the spatial expression for the mean of X(s) is m~ = S X(s)ds/IOl D

(7)

Characterization of atmospheric pollution by stochastic indicator

3813

where IDI represents the area of the domain/). The ergodic hypothesis is plausible for homogeneous atmospheric fields. In many practical situations the ergodic hypothesis cannot be investigated rigorously. However, it is often a reasonable working hypothesis in view of the lack of additional information (Christakos, 1992). Formulas such as (7) will be used below to express regional contamination indicator parameters in terms of spatial statistics.

level-crossing 3. STOCHASTIC I N D I C A T O R PARAMETERS OF

contours

ATMOSPHERIC C O N T A M I N A N T CHARACTERIZATION

Fig. 1. Schematic of level-crossing contours in a contaminated region D. ® are the regions of contamination exceeding the threshold level z.

3.1. Mathematical definitions and physical

interpretations Using the stochastic tools of the previous section, certain indicator parameters can be established which are useful for a practical characterization of polluted regions. Consider the polluted region D represented in terms of the homogeneous random field pair {X(s), I~(s, z)}, where X(s) denotes the contaminant concentration at location s ~ D and Ix(s, z) denotes the associated indicator landom field. Below we give general expressions for certain contaminant indicator parameters in terms of the random field pair (several other indicator parameters we discussed in Christakos and Hristopulos, 1996). Then we use the ergodic hypothesis to replace stochastic expectations with spatial averages over the polluted domain D. A useful indicator parameter is the relative area of excess contaminatiott (RAEC), Rx(z), defined by

Rx(z) = E[I.(s, z)] = ~ Ix(S, z)fx(x) dx.

(8)

0

threshold z, given by

PO,(z) = E[X(s)Ix(s, ,z)] = ~ xlx(s, z)fx(x) dx. 0

(101 The MEC is a positive, continuous and decreasing function in z. The range of P~(z) is between P~°(0) =mx and P~(oo) --- 0. Physically, in view of the ergodic hypothesis the P~(z) expresses the spatial average of the contaminant concentration that exceeds the threshold level divided by the total area,

P°(z) = ~ X(s) Ix(s, z) ds/IDI.

(11)

D

The mean excess differential contamination (MEDC) is defined as L~(z) = E{[X(s)

-

z]l~(s, z)}

Note that for a homogeneous indicator field the integral in equation (8) is independent of the spatial location. The Rx(z) is a continuous, positive and decreasing function of ze [0, ~ ) --} Rxe [1, 0]. It can be also expressed as Rx(z) == 1 - Fx(z). The inverse function z(R~) can be defined, and it is a continuous and decreasing function :in R~. In practice the range of the threshold level z is usually a finite interval [Zmi., Zm~] with end values that lie in the tails of the pdf, instead of the infinite interval [0, oo]. Physically, the RAEC Rx(z) expresses the fraction of the total area where the contaminant concer~tration X(s) exceeds z, i.e.

The MEDC LO,(z) is a convex and decreasing function in z, whose range lies between L~°(0)= m~ and L~(o0) = 0. At zero threshold z = 0, the derivative of L~(z) is equal to - 1. The integral of L~(z) over the threshold level z~[0, o0] is equal to 0.5(m 2 + a2), where ax is the variance of X(s). Physically, the MEDC expresses the spatial average over the total region D of the difference between the excess contamination and the threshold level, i.e.

R~(z) = ~ t~(s, z) ds/IDI = IOI/IDI

L°(z) = j [X(s) - Z]Ix(s, z)ds/IDI.

(9)

n

where ® is the subregion of D within which Ix(s, z) = 1 (Fig. 1). Thus, the RAEC is important for the accurate classification of geographical regions into substantially contaminated and negligibly contaminated zones and for the effective implementation of remedial strategies. Another contamination indicator parameter is the mean excess contamination (MEC) that exceeds the

= ~ (x - z)lx (s, z)fx(x) dx.

(12)

o

(13)

n

Bounds for the MEDC of a polluted region can be estimated as follows. Consider two contaminant processes X1 (s) and with X2(s) with univariate probability distributions Fx, and Fx2, respectively; means m~, and rex2, respectively; and bivariate distribution Fx~x~. Then,

Fx2(z) <~Fx,(z) if and only if Rx,(z) <~Rx2(z)

(14)

3814

G. CHRISTAKOS and D. T. HRISTOPULOS

and if m~ = mx~ = m then

E[X1 I X 2 ]

if and only if L~,(z) >t L~(z)

= X2

(15) for all z/> 0 (Alfsen, 1971). The meaning of the above result is that provided that the contaminant processes X1 (s) and X2 (s) have the same mean, the constrained mean of X l(S) given the measured value of X2(s) = x2 is equal to the measured value x2, if and only if the MEDC of X~(s) is larger than the MEDC of X2(s) for every threshold level z. Similar inequalities can be established for other indicator parameters and are of interest e.g. in the analysis of the scale effects (Section 3 below). The conditional mean excess contamination (CMEC) is defined as

P°(z) = E[X(s)I x t> z].

(16)

bivariate statistics of the contaminant processes involved. The noncentered indicator covariance (NIC) is a bivariate contaminant indicator parameter defined as el(h; z) = E[Ix(s, z)Ix(s', z)] = ~ S l~(s, z)I~(s', z)Fxx(h; du, du') o

(19)

0

where s ' = s + h; note that due to the stochastic homogeneity of the contaminant process X(s), the el(h; z) depends only on the space lag h. Physically, in view of the ergodic hypothesis the NIC expresses spatial correlations between contaminant concentration values that exceed the threshold level, ci(h; z) = S Ix(s + h, z)Ix(s, z) ds/lD I.

(20)

D

The centered indicator covariance (CIC) ~l(h; z) is defined in terms of the RAEC and NIC as

P°(z) is an increasing function in z. Physically, CMEC

Yi(h; z) = c1(h; z) - R~(z)

represents the spatial average of the contaminant concentration values that exceed the permissible threshold level z over the substantially contaminated subregion O = D; i.e.

where the ~(h; z) represents contaminant fluctuation correlations around the mean. Consider a spatial random field X(s) in R 2 representing e.g. contamination distribution in a horizontal plane or contamination averages in the vertical direction. A level-crossing contour is any contour along which the X(s) crosses from below-the-threshold to above-the-threshold values (see Fig. 1). The specific lenoth (SL) of the level-crossing contours is given by

P°(z) = S g(s)I~(s, z) ds/l®l.

(17)

D

The contaminant indicator dispersion (CID) is defined as

~s(z) = l-l(z)/IOl

Vx = E,[L~(z)]/m~ = S a2(z)dz/m~

(18)

o

where Ez[L~(z)] is the threshold-averaged MEDC and a~(z) is the variance of Ix(s, z). The ~Fx provides a measure of the dispersion of the indicator variance with respect to the threshold level z. As already mentioned, in practice the possible threshold levels z usually lie in a finite interval ['Zmin, Zmax.], and the limits 0 and oo in equation (18) should be replaced by Zm~ and z.... respectively. Then, for a given contaminant pdf, ~ provides a measure of the ratio of the average excess contamination for all possible threshold values z over the mean contamination. The ~F~is a characteristic of the contaminant pdf. As we show in Section 5 below, the CIDs of the Gaussian and the log normal contaminant pdf's are functions of the mean concentration and the standard deviation. The exponential pdfhas a constant CID ~F, = 0.5. The ~F~is zero if and only if the pdf of the contaminant pdf is a delta function. If L~,(z)/> L~,(z), then Wx,/> W,,; in fact, ~ provides a quantitative measure of the change in the standard deviation between the various scale levels (Section 3 below). The RAEC, MEC, MEDC, CMEC and CID are indicator parameters that involve the univariate pdf of the contaminant process. The contaminant indicator parameters to be discussed next involve the

(21)

(22)

where H(z) denotes the expected total length of the level-crossing contours. If L denotes the linear size of the domain D and dn is the fractal dimension of the level-crossing contours (Mandelbrot, 1983), the scaling relations below are true: [DIoc L 2, I-I(g) oc L ~", fs(z) OCL d"-2.

(23)

If dn < 2, the SL fs(z) of the level-crossing contours becomes zero. This constraint restricts the permissible indicator covariances to be used for modelling contaminant processes. The following analytical expression has been obtained for the SL of an isotropic Gaussian random field X(s) (Swerling, 1962):

f~(z) = exp [

(z2~x(0)m~)2lj//2 2x

(24)

where 2x = x/rx(0)/- 5"(0) (the 6~'(0) denotes the second derivative of the covariance of X(s) at zero lag) has length dimensions and is proportional to the correlation length b. Note that equation (24) is valid provided that the first derivative of the covariance is zero at the origin, a condition that is not satisfied in the case of the commonly used spherical covariance (see equation (68) below). Contaminant fields are usually simulated on lattices. The level-crossing contours can be approximated on the lattice with the rectilinear contours formed by

Characterization of atmospheric pollution by stochastic indicator the bonds that join lattice sites with opposite indicator-field values. Below we obtain expressions for the lattice specific length :~(z) on a square lattice with length L = aN, where a is the smallest lattice distance and N 2 is the total number of sites. The sites are labelled by a pair of integers that determine their Cartesian coordinate.,;, so that s~ represents the site with coordinates (ia, ja). We define the bond function t(o), (u)(z) by J'l t(u),

(kl)(Z)

=

EEt,,), (i+,.,~(z) (26)

where the first sum over i and j is over all the lattice sites su, and the second sum over tr and tr' runs over the nearest neighbors of the site s u. For a statistically homogeneous X(s) the average EEt(o), (i+~.j+~,)(z)] is spatially uniform. Hence, the lattice SL can be expressed by the expecl;ation of the four bonds joining a site with its neighbors, i.e. E a, ¢'=

e~(0)(x-~12+xl)----2ex(h)xxx2~ 2 [ ~ (0) - g~(h)] J"

(29) Finally, if the contaminant fluctuations are stochastically isotropic, :b(z)=4

S dx2 ~ dxxf~x(a;xx, x2)/a. (30) z-rex

A summary of the stochastic indicator parameters discussed above is given in Table 1. (A larger list of stochastic indicators can be found in Christakos and Hristopulos, 1996.)

The indicator parameters defined above can be related to each other by means of simple equations. For a given region, these equations provide (a) information regarding the way changes in certain of these parameters affect the others, and (b) the means to calculate each one of the parameters in terms of the others. As is shown in Christakos and Hristopulos (1996), the indicator parameters Rx(z), PD(Z), L~(z) and P~(z) are related by means of the following set of equations:

EEt(o),(/+a,J3(z) + 1

+ t(o~, (~.j+,,)(z)]/2a.

- e~(h)] -~

3.2. Relations between the contamination indicator parameters

a , c r ' = :k 1

+ t(o3, (i.i+a,)(z)]/2L 2

t*b(Z) --

xexp{

- - oo

where s u and su are nearest neighbors on the lattice. The 4(z) is obtained from the expectation of the spatially averaged bond function as follows:

~.,

f~(h; x~, x2) --- [ 2 ~ ( 0 )

(25)

otherwise

i,j

bivariate pdf is given by

z-rnx

if I~(s~,z) + I~(su, z ) = 1

4(z) = a ~'

3815

(27)

The stochastic average in equation (27) can' be evaluated analytically using the E[tt/j). (,+¢d)(z)] =--P[X(sij) >t z, X(s/+,d ) < z]

+ PEX(s 0 < z, X(s,+,,~) >/z]

PV~(z) = R~(z) P°(z)

(31)

L~(z) = P~(z) - zR,(z)

(32)

zR~(z) ~ m'~

(33)

dL~(z)/dz = - Rx(z)

(34)

de°~/dz = z dR~/dz

(35)

dP°,/dPR, = z(Rx) >i 0

(36)

d2p~/dR~ = dz/dRx <~0

(37)

z - mx

=:2

~ dx, z - - rtl x

dL~/dR. = - R~ dz/dR.

S dx2

xf~(h; xl, x2)h=(,,,,, o)

(38)

and

- o~

(28)

wheref~(h; x~, x2) is the bivariate pdf; an analogous expression holds for E[t(o2, (~d+,,)(z)]. The Gaussian

gx dP~/dR. + P~ = z.

(39)

Equations (31)-(39) establish restrictions on the relative changes of the univariate indicator parameters. For example, according to equation (33), the

Table 1. A summary of the stochastic indicator parameters RAEC MEC MEDC

R.(z) e~(z) L~(z)

CMEC

ey(z)

CID NIC CIC SL SL

c,(h; z) ~i(h; z) :,(z) :b(Z)

Wldx

fraction of the total area where contamination exceeds the threshold level z average contamination exceeding the threshold level z over the total region D average difference between the excess contamination and the threshold level z over the total region D average contamination exceeding the threshold level z over the substantially contaminated subregion ® ,- D dispersion of the indicator variance with respect to the threshold level z spatial correlations between contamination values exceeding the threshold level z contaminant fluctuation correlations around the mean expected total length of the level-crossing contours in the domain IDI per unit area expected total length of the transition bonds per unit lattice area

3816

G. CHRISTAKOS and D. T. HRISTOPULOS

average contaminant concentration mx provides an upper bound for zR,`(z). While the derivative of L ° ( z ) (equation (34)) is negative and increasing, the L°~(z) itself is a decreasing function of the threshold z. Equations (36) and (37) imply that P ° ( R x ) is an increasing and concave function in R,`. Equation (38) shows that L ° ( R x ) is an increasing function in R~. The NIC q(h; z) is expressed in terms of the RAEC R,`(z) in the limiting cases of zero and infinite lags; i.e. cl (0; z) = R~(z)

(40)

lira Cl(h; z) = R~(z)

(41)

and h~ao

where h = [hi. At the limit of zero lag, h = 0 equation (21) gives the CIC in terms of RAEC, viz. ?,(0; z) = R,`(z) - R~(z)

(42)

where the indicator variance is given by a ~ ( z ) =

~(0; z). Equation (5) leads to the following relationships between the mean and the covariance of X(s) and the contaminant parameters RAEC and NIC; m,` = S R x ( z ) d z

(43)

and 0

(44)

0

It can be shown (Christakos and Hristopulos, 1995, 1996) that the SL :,(z) is related to the two-dimensional NIC cl (h; z) by lira OAi(h; z)/Oh = - f,(z)/n

(45)

boO

where Ai(h; z) - So cl (h; z) d4~/21r and ~b is the polar angle of the lag vector h. These results are generally valid for homogeneous atmospheric random fields. In the special case of an isotropic NIC, Ai(h; z) = c~(h; z), where h = [hi. Then, the SL is related to the derivative of the NIC at zero lag by :s(z) = - n lim Oq(h; z)/Oh

(46)

h~O

i.e. the limit of E~(z) is related to the behavior of the NIC at zero lag. In Christakos and Hristopulos (1996) it is proven that in the case of a Gaussian SRF X(s) equation (46) is equivalent to equation (24). If ~,(z) and f'~(z') are the SLs of the Gaussian atmospheric fields X(s) and X'(s), respectively, the following scaling relationship is obtained from equation (24): :',(z

a'~/a,` + m"

-

m:, a /a,,) = :,(z)2x/2".

(47)

If X'(s) differs from X(s) only in the correlation length, equation (47) reduces to b'e',(z) = b :,(z).

where ® [.] is the scaling functional that defines the upsealing operation. The corresponding pdf's of the initial and the upscaled random fields are related by d fz(x2) = ~x2 o _:f~,~]f l ( x 0 dxl.

0

Fx(h) = S S F1(h; z , z ' ) d z d z ' .

3.3. Scale effects In practice, contamination characterization depends on the analysis of atmospheric processes at various physical scales, as well as on the establishment of suitable quantitative connections between the results obtained at each one of these scales. This is important, for considerable differences may exist between the physical behavior and the statistics of the same process at different scales. The appropriate choice of a scale should lead to a model that contains sufficient detail to reproduce accurately all the physical effects relevant to a specific problem Without rendering the computations impractical. The change of the contaminant indicator parameters under a change-of-scale operation (e.g. upscaling) can be used to calculate the change of the underlying pdf. Assume that an upscaling operation maps the spatial random field X1 (s) into the field X2 (s) such that X2 (s) = 0 [X l(S)] (49)

(48)

Numerical illustrations of the predicted behavior of the indicator parameters are given in the applications of Section 5 below.

(50)

In case that the ®[.] is differentiable and the ® - a [.] has a continuous derivative, equation (50) reduces to f2(x2) =f~ [O-1(x2)] [ d®-l(x2)/dx 2 l-

(51)

The RAEC difference due to the change in the pdf is given by 6R21(z) -= Rx2(z) - Rx,(z) = S [f~2(v) - f ~ , (v)] dr. z

(52) By taking the derivative of equation (52) we obtain an expression for the pdf change, 6f21(z) --f~2(z) -f~l(z) = - d6R~(z)/dz = dRxl (z)/dz - dRy, 2 (z)/dz.

(53)

Equation (53) permits the calculation of the pdf change due to upscaling in terms of the RAEC, which can be calculated from the available data. Similarly, the pdf change can be expressed in terms of other contaminant indicator parameters as 6f21 (z) = -- d26L°x(z)/dz2 =

--

z- x

d6P~(z)/dz

(54)

where 6L~(z) = L°2 - L~I and 6P~(z) = p O _ p,`~. If the initial pdf and the form of the upscaling functional ® [ . ] are known, it is possible to perform the upsealing process repeatedly. Then, the pdf change at each stage of the upscaling process is calculated by means of equations (52), (53) or (54) above. In Section 5 below, we use a numerical ease study to demonstrate the pdf change due to an upscaling of the sulfate deposition field.

Characterization of atmospheric pollution by stochastic indicator In some cases, scale effects can be quantified in terms of the stochastic indicator parameters (Christakos and Hristopulos, 1966). For example, consider a polluted region D with size IDI that is divided into Nx subregions 011) of the same size I~tl~l. The partitioning is continued so that each of the N1 subregions is divided into N2 smaller sub-subregions u~2) of equal size 1~'2)1 ( < Io(1)1), etc., up to a final partitioning consisting of N~ subregions (u!~)) of equal size Iu(~)l ( < 0(~-1)1 < ... < [u(2)lu(1)l). The MEDCs of the contaminant process X(s) over the subregions of sizes I~t')l (~ = 1, 2 .... ) and the global region of size IDI satisfy the following expressions:

I~(~)1~< IDI lim L~"(z) = L~(z)

L~")(z) decreases in

(55)

I~'"1~ IDI

for all z ~> 0. In other words, the L~,"(z) decreases in Iot~)l such that the larger the subregion size, the smaller the MEDC; as 1:he subregion size increases its MEDC approximates more accurately the actual MEDC L~(z). As a consequence of equation (55), the following is also true in terms of the corresponding CIDs; q~'~' decreases in lim

Io(~)l ~< IDI

V~~' = ~P~

Condition (57) is not sufficient for ergodicity, but it provides a practical test for ergodicity breaking that could be implemented using available measurements. For Gaussian pdf's a theoretical analysis is possible showing that the necessary condition (57) can be expressed as follows: arcsinpx (h)

lim o~dh ID]~

S 0

e-(~-mx)2/(l+sin°)dO/IDl=O" (58)

Provided that the fluctuations of the spatial random field X(s) have a finite correlation length b, equation (58) is true at the limit of an infinite domain size L. This happens because the inner integral in equation (58) is essentially cut off by the correlation length of the field X(s). In practice, equation (58) should be satisfied for ergodic random fields provided that the correlation length is significantly smaller than the domain size. Note that long-range power-law correlations lead to ergodicity breaking. Lastly, in several practical applications where no rigorous ergodicity tests can be applied, ergodicity can still be considered as a working hypothesis that permits the calculation of the contaminant indicator parameters by spatial averaging of single-realization data.

(56)

Since CID is a function of the mean and standard deviation for each scale level, equation (56) offers a quantitative measure of the change in standard deviation between the various scale levels.

4. ERGODICITYANALYSISAND THE ESTIMATIONOF THE INDICATORPARAMETERSIN PRACTICE As was mentioned above, if the homogeneous SRF X(s) is ergodic in the domain D the stochastic indicator parameters can be expressed in terms of spatial averages involving the random field pair {X(s), Ix(s, z)} (see equations (9), (1 I), (13), (17); also (20) and (21)). This makes possible the efficient calculation of these parameters in practice. For example, R,,(z) can be approximated in various ways ranging from the simple, equal weighted, cumulative histogram of the available data to the sophisticated kriging estimators (Journel, 1983). Mathematically rigorous theorems for ergodic random fields are given by Adler (1981), where sufficient and necessary conditions for ergodicity valid for Gaussian random fields are also formulated. In Christakos and Hristopulos (1996) it is shown that a necessary condition for the validity of the ergodic hypothesis in terms of the indicator parameters RAEC and NIC has as follows lim o~c,(h; z) dh/lDI - R2(z) = 0. IDI~oo

3817

(57)

5. CHARACTERIZATIONOF THE SULFATEDEPOSITION DATASET 5.1. Sulfate deposition indicator parameters In an earlier work (Christakos and Thesing, 1993) the N A D P / N T N sulfate deposition data set (National Atmospheric Deposition Program/National Trends Network, 1987) in the conterminous U.S. was studied. Sulfate deposition surveys are important for establishing the spatial distribution of pollutants that are being introduced into the environment (e.g. Bilonick, 1985). Data collection makes it possible to determine the amount or rate of changes in the deposition process. The data were weighted average SO4 concentrations (mg # - 1) for the 1986 winter season at 94 monitoring locations which meet data quality objectives. Weekly samples were collected, and the seasonal weighted mean was determined from the concentrations of these weekly samples. The accumulated mass per unit area (mg m-2) was calculated from the average sulfate concentration and the cumulative depth of precipitation through the period. Table 2 displays summary statistics for the 94 data values. A more detailed discussion of the data set may be found in Christakos and Thesing (1993). Given the limited number and the uneven distribution of the data available it is not possible to construct the probability law of sulfate deposition reliably. It is thus necessary to make some assumption regarding the possible form of the sulfate deposition pdf. A detailed analysis of the behavior of the sulfate deposition indicator parameters is made below assuming

3818

G. CHRISTAKOS and D. T. HRISTOPULOS

a g a m m a pdf. The C I D is then calculated for three other pdf models, viz. the exponential, the Gaussian and the log normal. Comparisons are made and conclusions are drawn regarding contaminant characterization for each one of the possible scenarios. We first calculate the indicator parameters for the sulfate deposition data set assuming that the sulfate deposition X(s) is characterized by the g a m m a pdf

The R A E C is calculated from equation (8), i.e.

Rx(z) = r

#,

F(a)

(60)

where F(a, z/fl) is the incomplete g a m m a function (Abramowitz and Stegun, 1970). The M E C (in mg m -2) can be calculated from equation (10) and is given by

(59) The contamination threshold z is given in terms of the R A E C by the inverse of equation (60),

2 2 where c~ = m,,/trx, fl = a2/m~,, mx = 182.53 m g m -2, ax = 156.968 m g m -2 (see Table 2), and F(.) denotes the g a m m a function.

z(Rx) = d P ~ ( R x ) / d R x = [dP~(z)/dz] (dRx/dz) - 1. (62)

Table 2. Summary statistics of the sulfate deposition data set (in mg SOL- m- 2)

The M E D C is calculated (in mg m -z) from equation (32),

Highest value: 865.365 Lowest value: 2.441 Range: 862.924 Mean: 182.530 Median: 163.603 Standard deviation: 156.968

- zr(~,~)/r(~)

'7

/ r3

0.80

'6

'5~v

/"

II[

•

,

//

o.oo

,'4

r.

S. 0.40

•~

,/

0.20

'

*

0

. 0

0 2 10"

~ 4 10"

0 6 10"

8 10"

1 10~

Threshold

Fig. 2a. Plot of the indicator parameters as a function of the sulfate concentration threshold z, obtained for a gamma sulfate deposition distribution. The MEC, MEDC and CMEC are normalized to a mean concentration of 1 mg m- 2. The values of the normalized MEC and MEDC as well as the values of the RAEC are shown on the left vertical axis. The values of the normalized CMEC arc shown on the fight vertical axis.

(63)

Characterization of atmospheric pollution by stochastic indicator with L~(0)= ms. The CMEC function, po, can be obtained directly from equation (16), or from equation (31) using equations 161) and (63), or from the solution of equation (39) with the boundary condition P ° ( R x = 1) = rex; the result is P°(Rs) = 1~F ( = + 1, ~ ) / F ( , , ~ )

(64)

(in mg m - 2). The indicator parameters R~, p D, L D and p o are plotted in Fig. 2a as functions of the contaminant threshold z. In Fig. 2b we plot z, p o, LxOand p o as functions of R,. Given the value of one of these parameters, Fig. 2a and b readily provides the corresponding values of the remaining parameters. Next the CID ~P~ is calculated from equation (18) assuming various pdf models for the sulfate deposition data and z ~ (0, m). It is worth noticing that for the gamma, the Gaussian and the log normal pdf's the CID is a function of the mean concentration and the standard deviation. In particular, the ~'x of the gamma pdf is calculated numerically in terms of the mean and the standard deviation; for the values of Table 2, it was found that u/~ = 0.443. For the CIDs of the Gaussian and the log normal pdf's analytical

3819

expressions are obtained, i.e. Us = a~/x/~ ms

(65)

°l x ~- 2~x(aff~/2 m~) -- 1

(66)

and

respectively. The exponential contaminant pdf, on the other hand, is characterized by a constant CID (= 0.5) that does not depend on the mean and standard deviation. This property is characteristic of the exponential pdf. For the specific sulfate deposition data set, however, the numerical tPs-values obtained for the four different types of pdfs are very close to each other (see Table 3). There is important information conveyed by these q~s-values. For example, the value ~Px = 0.443 (gamma pdf) means that the expected amount by which sulfate deposition will exceed the threshold level z, averaged over all possible z-levels, is equal to 44.3% of the sulfate deposition mean over the conterminous U.S. (i.e. 80.89 m g m -2 = 0.443 x 182.53mgm-2). This is a considerable percentage that should be taken into consideration in environmental policy making. The same conclusion is valid for the other pdf's. In fact, in the case of the Gaussian

1.2 1 0 "

1.0 10 8"

8.0 10 ~" | ,-

6.0 10 ~"

N

4.0 10 ~"

2.0 10 ~"

_

- . . . . -

0.01 0

0.2

0.4

0.6

0.8

1

RAEC Fig. 2b. Plot of the deposition threshold and the indicator parameters MEC, MEDC and CMEC as a function of the RAEC for a gamma sulfate deposition distribution.

3820

G. CHRISTAKOS and D. T. HRISTOPULOS

Table 3. CID values for the sulfate deposition data PDF

CID

of the bond lengths for a single realization, i.e. I'D(Z) = a

Z

Z

[t.~..+o.,~(z)

i , j ~r,a'= ++_l

Gamma Exponential Gaussian Log normal

0.443 0.5 0.485 0.457

+ qua, ,.J+,')(z)]/2L2"

pdfthe z-averaged MEDC is E,[LO~(z)] = a~/x/~. The last expression means that the expected amount by which sulfate deposition will exceed the threshold z, averaged over all possible z-levels, is proportional to the standard deviation of the sulfate deposition data. 5.2. Numerical simulations of the sulfate deposition levels--comparisons with analytic expressions Next, we calculate the specific length (SL) of the sulfate deposition level-crossing contours for a subset of all the data obtained from 55 sampling stations within a geographical window located over the North-Central U.S. (see Haas, 1990). We investigate certain issues related to the estimation of the levelcrossing contours by means of the analytical expressions (24) and (30) as well as numerical simulations. Following the study by Haas (1990), we assume that the sulfate deposition data are normally distributed within the above window with mean mx~ 170mgm -2, and standard deviation as ~ 145mg m -2. In the same study semivariogram data fits were presented that suggested that the correlation length of the sulfate deposition process is smaller than 50 kin, which is the shortest lag distance that can be resolved. We model the deposition field on a square lattice with N = 300 nodes in each direction using the GSLIB computer software (Deutsch and Journel, 1992). The following covariance models have been used for the sulfate correlations: (i) the Gaussian model ex(h) = ~I exp

- ~

(67)

and (ii) the spherical model a~(1

~x(h)

3h h3 ) -- ~'~ + ~

if h ~ [0, b]

=

(68)

0

ifh>.b

where h is the correlation length. We investigate two correlation lengths b/L = 6 . 6 7 x 1 0 -3 and b/L = 16.67 x 10- 3; note that the scaling relationship (48) can be used to obtain the SL of the level-crossing contours for arbitrary values of the sulfate correlation length. The lattice SL #b(Z) is evaluated by means of three different approaches. The first method is a direct calculation of the analytical expression (30) using the numerical integration routine Nlnt from Mathematica (Wolfram, 1991). The second method is based on the assumption that the sulfate deposition field is ergodic, in which case the stochastic average in equation (26) can be approximated by the spatial average

(69)

The sum in equation (69) can be calculated numerically by direct counting of the transition bonds. The calculations with the first two methods were in good agreement, thus providing a posteriorijustification for the ergodic hypothesis. Finally, the third method is based on a numerical estimation that uses a relationship between the SL 4(z) and the behavior of the NIC at the origin. In the case of a lattice random field the slope of the NIC at the origin is replaced by the NIC lattice difference 6c,(0; z) = [q(a; z ) - ct(0; z)]/a = [])i(0; z) -- vi(a; z)]/a

(70)

where Vl(h; z) is the semivariogram of the indicator field, which can be calculated by means of the GSLIB. The lattice SL can be calculated in terms of the NIC lattice difference using probability analysis, and the following relationship is obtained (Christakos and Hristopulos, 1995, 1996): ~b(z) =

-

-

4tSci(0; z).

(71)

Note that NICs of environmental processes often include a nugget term. These processes can be modelled by adding to covariances such as (67) or (68) above a constant term for near-zero lags. However, if the value of the lattice constant is sufficiently small both ci(a; z) and cx(0; z) should include the nugget term which, therefore, vanishes upon calculating the covariance difference in equation (70) above. Equation (71) represents the lattice equivalent of equation (46), which gives the specific length of the levelcrossing contours for a continuum spatial random field with an isotropic indicator covariance in terms of the covariance behavior at the origin. The SL of the continuum level-crossing contours and the lattice SL are plotted as functions of the threshold in Figs 3-5 below. The SL is calculated by three different methods: (i) by direct counting of the transition bonds on a square lattice, (ii) by the semivariogram method using equations (70) and (71), and (iii) from the analytical expression (24), which is valid for continuum fields. Note that only the lattice SL is given for the spherical covariance, Fig. 5, for which the continuum SL is not well-defined by means of equations (24) and (46). On the basis of Figs 3-5 the following remarks can be made: (1) The SL decreases with increasing correlation length in agreement with the scaling relation (48). (2) The lattice SL obtained by the bond-counting and the indicator-semivariogram methods are in good agreement. (3) The lattice SL, however, overestimates the SL of the continuum level-crossing contours. The

Characterization of atmospheric pollution by stochastic indicator 0

--it-=corltinuum (analytical)

0.20 -

I

0.30'

0.20'

Lattice (bond counting)

/% I

'

m

i

--e--Continuum (analytical)

0

0.50 =

0.40

~ l W . Lattice (semivariogram)

I I

Lattice (bond counting)

==4, =Lattice (semivariogram)

3821

,,..,

I; ,,'

\

i

',,,,,

',x.

o.,o 1 0.050-

0.10'

0.0'

oo ....

I ....

0.1

I ' ' ' ' I ' ' ' ' I ' ' ' "

0.2

0.3

0.4

I ' ' ' '

0.5

'',

I

0.6

Sulfate concentration (mg/m s)

0

I,,,

0.1

. I , , , , i , , , , i , , , , i , , ,

gl

0.2

0.6

0.3

0.4

0.5

Sulfate concentration (mg/m2)

Fig. 3. The continuum and the lattice SL are plotted as functions of the sulfate concentration threshold z. The sulfate deposition field is assumed normally distributed and a Gaussian covariance model is used with correlation length b/L = 6.67 x 10-3.

Fig. 4. The continuum and the lattice SL am plotted as functions of the sulfate concentration threshold z. The sulfate deposition field is assumed normally distributed and a Gaussian covariance model is used with correlation length b/L = 16.67 x 10-3.

discrepancy is due to two factors: (a) the continuum contours are approximated on the lattice by linear segments, and (b) the lattice SL also includes contributions from bonds that belong to open-ended line segments. In such cases, therefore, numerical simulation can lead to inadequate characterization of contamination contour properties. An improved estimate of the continuum SL can be obtained from lattice simulations by using, instead of the lattice estimator 4(z), the expression Es(z) ~ n6ci(0; z) as an approximation of the continuum equation (46). (4) The lattice SL for fields with a spherical covariance is higher than the corresponding SL for the Gaussian covariance with the same standard deviation and correlatior~ length. The continuum SL cannot be evaluated by :means of equation (24), which is not valid for the spherical covariance, as we mentioned above. An analysis of the continuum SL shows that the zero-lag limit in equation (46) does not exist due to the behavior of the first derivative at the origin (see Christakos and Hristopulos, 1996).

tices of size 102 denoted by v~ (i = 1 . . . . . 9 x 102), to obtain the coarse-grained lattice D2. The thus upscaled deposition process X2(s) is defined on D2 as the sublattice average of the original sulfate process X1 (s), i.e.

5.3. Upscaling analysis of sulfate deposition Let Xt(s) = X(s) be the sulfate deposition process described in Section 5.2 above. The original scale of the Xt(s) is defined by the square lattice Dt with 9 x 104 sites. Next we partition Dx into 9 x 102 sublat-

X2(s21) = E Xl($J )/100 $j~'oi

(72)

where i = 1 . . . . . 9 x 102 and s21 denotes the centerpoint of the sublattice ui. The RAEC difference, 6gx(Z ) = Rx2(z) - Rxl(z), and the resulting change in the pdf, 6 f 2 a ( z ) - f x 2 ( z ) - f ~ ( z ) , for the processes Xl(s) and Xz(s) were calculated using equations (52) and (53). The corresponding plots are shown in Fig. 6 as functions of the sulfate concentration threshold z. Note the considerable change in 6Rx(z) (up to 28%) and in 6f21(z) (up to 30%). The spatial averaging of a contaminant process may, therefore, change significantly the probabilistic structure of the original process and alter considerably crucial parameters of environmental decision making and risk assessment. Analytical expressions for the change in 6f2 ~(z) can also be obtained by assuming that under a small change of support the main factor contributing to 6f21 (Z) is the change in the variance of the upscaled process. In view of this assumption we obtain the following expression for the pdf change of the

3822

G. CHRISTAKOS and D. T. HRISTOPULOS Lattice (bond counting)

0

I -,,o - RAECdifference I

-Lattice (semivariogram)

I . . . . . PDF difference I

0.30"

1.0"

"2

0.20"

'0

0.80'

.1= ==

~

I

0.60'

e

0.0'

0.40'

/

1

0.20'

\ • • ,1

,,

0.1

• , i,

0.2

,,

=1 =-,

-0.30

, i = = = , 1 , ~

0.3

0.4

0.5

I

0.6

Sulfate concentration (mg/m =)

Fig. 5. The lattice SL is plotted as functions of the sulfate concentration threshold z. The sulfate deposition field is assumed normally distributed and the spherical covariance model is used with correlation length b/L = 6.67 x 10-3. Gaussian sulfate deposition process:

where r = axl/ax 2 is the standard deviation ratio and can be determined by fitting to equation (73) the 6f21(z) obtained from the deposition data using equaton (53). In Fig. 6 we show a least-squares fit of equation (73) to the 6f21 (z) obtained from the simulation of the sulfate process; the standard deviation ratio thus obtained is given by r = 3.67.

6. CONCLUSIONS Many investigations of atmospheric pollution are closely connected with the study of contaminant concentration exceeding certain threshold levels. In this work stochastic indicator parameters that characterize excess contamination levels were formulated using the theory of random fields, and certain comments were made regarding their importance in environmental policy making. These parameters are related to each other by means of simple algebraic and ordinary differential equations, and can be calculated in terms of the statistical moments of the random fields.

i ,=

-0.10' -6

-0.20'

~ 0.0"

-2 . = .

-8

,l||l,.|,l,l,Wl,,,,l==,,i,,,,i,,,|

0.1

0.2

0.3

0.4

Sulfate concentration

0.5

0.6

0.7

(mg/m ~)

Fig. 6. The RAEC difference tSR21(g) and the pdf difference 6f21(z) that result from the upscaling of the normal sulfate deposition field are plotted as a function of the sulfate concentration threshold z. Both the 6R2t(z) (dashed line with markers) and the 6f21(z) (broken line) are evaluated using the simulation data. The 6f21(z) obtained by means of a least-squares fit of the data to the theoretical approximation of equation (73) is also shown (continuous line). The ergodicity assumption is important in practical applications where the contaminant indicator parameters need to be calculated by spatial averaging of single-realization data. Expressions allowing the calculation of the indicator parameters in terms of the data available are discussed and an ergodicity test for contaminant processes based on indicator moments was suggested. We also discussed the important effect of the physical scale in the calculation of the indicator parameters. A sulfate deposition data set was studied in detail by means of the stochastic indicator parameters, assuming four commonly used probability distributions. The specific length of the level-crossing contamination contours in the continuum and on the lattice have been examined, and it has been shown that the slope of the indicator covariance at zero lag provides a good approximation of the lattice specific length. We have also shown how to use lattice results to obtain estimates of the continuum specific length, and we have pointed out that in the case of the spherical model the continuum specific length is not well-defined due to the irregular behavior of this covariance at zero lag. While in many cases numerical simulation provided good indicator approximations, it led to inadequate estimates of the specific length of

Characterization of atmospheric pollution by stochastic indicator the level-crossing contamination contours. Significant changes in the indicator parameters and the probability distributions have been found for the sulfate deposition data as a result of a small change in the observation scale. Tiffs may lead to considerable changes in crucial parameters of environmental decision making and risk assessment. Acknowledgements--This work has been supported by grants from the National Institute of Environmental Health Sciences (Grant No. P42 ES05948-02) and the Army Research Office (Grant No. DAAL03-92-G-0111), Research Triangle Park, North Carolina.

REFERENCES Abrahamsen G. (1987) Air pollution and soil acidification. In Effects of Atmospheric Pollutants on Forests, Wetlands, and Agricultural Ecosystems (edited by Hutchinson T. C. and Meema K. M.), pp. 321-331. Springer, Berlin. Abramowitz M. and Sl:egun I. A. (1970) Handbook of Mathemotical Functions. ])over, New York. Adler R. J. (1981) The Geometry of Random Fields. Wiley, New York. Alfsen G. M. (1971) Compact Convex Sets and Boundary Integrals. Springer, Berlin. Aoki M. (1967) Optimi2ation of Stochastic Systems. Academic Press, New York. Bilonick R. A. (1985) The space-time distribution of sulfate deposition in the Northeastern United States. Atmospheric Environment 19, 1829-1845. Black F. M. (1991) Control of motor vehicle emissionsthe US experieno~.. Crit. Rev. Envir. Control 21, 373-410. Christakos G. (1992) Random Field Models in Earth Sciences. Academic Press, San Diego, California. Christakos G. and Thesing G. A. (1993) The intrinsic random field model and its application in the study of sulfate deposition data. Atmospheric Environment 27A, 1521-1540.

3823

Christakos, G. and Hristopulos D. T. (1995) Contaminant Characterization in Terms of Stochastic Indicator Parameters. Res. Rep. No. SM/12.95, Department of Environmental Science and Engineering, University of North Carolina, Chapel Hill, North Carolina. Christakos G. and Hristopulos D. T. (1996) Stochastic indicators for waste site characterization. Water Resour. Res. (in press). Clark T. L. (1980) Annual anthropogenic pollutant emissions in the U.S. and the southern Canada east of the Rocky Mountains. Atmospheric Environment 14, 961-970. Deutsch C. V. and Journel A. G. (1992) Geostatistical Software Library and User's Guide. Oxford University Press, Oxford. Haas T. C. (1990) Kriging and automated variogram modeling within a moving window. Atmospheric Environment 24, 1759-1769. Haas T. C. (1994) Prediction of a first and second order nonstationary spatio-temporal process. J. Am. Statist. Ass. preprint. Journel A. G. (1983) Non-parametric estimation of spatial distributions. Math. Geol. 15, 445-468. Journel A. G. (1987) Geostatistics for the environmental sciences. EPA Project No. CR811893, Tech. Report, U.S. EPA, EMS Lab, Las Vegas, Nevada. Kendall M. and Stuart A. (1977) The Advanced Theory of Statistics: Vols 1 and 2. Griffin, London, UK. Mandelbrot B. B. (1983) The Fractal Geometry of Nature. Freeman, San Francisco, California. Neely W. B. (1994) Introduction to Chemical Exposure and Risk Assessment. Lewis Publishers, Boca Raton, Florida. Raiffa H. and Schlaifer R. (1961) Applied Statistical Decision Theory. Harvard Business School, Boston, MA. Rivoirard J. (1994) Introduction to Disjunctive Krigin# and Non-linear Geostatistics. Clarendon Press, Oxford, UK. Swerling P. (1962) Statistical properties of the contours of random surfaces. IRE Trans. Infer. Th. IT-S, 315-321. Vukovich F. M. (1994) Boundary layer ozone variations in the eastern U.S. and their association with meteorological variations: long-term variations. J ffeophys. Res. 99, 16,839-16,850. Wolfram S. (1991) Mathematica: A System for Doinff Mathematics by Computer, 2nd Edn. Addison-Wesley, Reading, Massachusetts.