Multifractal Simulation of Geochemical Map Patterns

NTRODUCTION

The mathematical modeling of geochemical data continues to present an important challenge to earth scientists. Basically, the available data are 2-dimensional realizations of 3-dimensional patterns, which are the results of one or more genetic processes that usually are poorly understood. It may be difficult to project from the observations into the rock mass. In this paper, multifractal modeling is advocated as a methodology to characterize geochemical patterns. If a geochemical pattern is multifractal, its parameters including the multifractal spectrum can be estimated from 2-D map patterns, and projected downward from the surface. This type of extrapolation is similar to extrapolation from 1-D series of observations into 2-D or 3-D space.

Mutifractals arose in physics and chemistry as a generalization of (mono) fractals (Lovejoy and Schertzer, 1990; Feder, 1988). Multifractals can be regarded as spatially intertwined monofractals (Stanley and Meakin, 1988). It is also good to keep in mind that multifractals, as well as monofractals, are self-similar and scale-independent, and that multifractals apply to continuous spatial variability patterns whereas monofractals are for binary patterns.

The relation between monofractals and multifractals was considered by Herzfeld (1993). Herzfeld et al. (1995) showed that the ocean floor could not be modeled as a monofractal. Better results were obtained by means of a multifractal model after incorporation of nonstationary component (Herzfeld and Overbeck, 1999).

In the application to geochemical data, the relation between monofractals and multifractals can be interpreted as follows: suppose that a rock mass is divided into small blocks of equal size and shape and that the average concentration value of 8 chemical element is known for every block. Suppose further that a histogram of the concentration values is created. If the chemical element has multifractal distribution, the blocks corresponding to every class of the histogram approximately form a single fractal with variable fractal dimension (cf. Evertsz and Mandelbrot, 1992). Together, the fractal dimensions of blocks with the same element concentration value form a multifractal spectrum. If the width of the class of the histogram approaches zero, the multifractal spectrum approaches a continuous curve representing its limiting form.

These concepts can be formalized as follows. For convenlence, 1-dimensional space is used but the results can be generalized readily to the map or to 3-dimensional space. Suppose that μ = xε represents total amount of a metal over a line segment of length ε; x is the metal's concentration value. In the multifractal model, it is assumed that (1) μ∞ε^α where ∝ denotes proportionality, and α is the so-called singularity exponent corresponding to concentration value x, and (2) N_α(ε) ∝ ε^-f(α) where N_α (ε) represents total number of line segments of length ε with concentration value x, and f(α) is the fractal dimension of these segments. The multifractal spectrum is a plot of f(α) versusa α.

In computer simulation experiments with discrete segments or cells, the preceding approach can be used without modifications. In practical applications, x and α may be continuous variables. Then the space can be subdivided into small equal-size cells with approximately constant values of x and α. In the histogram method to be discussed later, the frequencies N_α(ε) of cells with the same value of α are used to plot the multifractal spectrum for selected values of ε.

The multifractal model was used in Agterberg (1994) to compare concepts and methods from two different fields: geostatistics and fractal geometry. It was discussed that early geostatistical work by De Wijs (1951) and Matheron (1962) contains basic elements of multifractal modelling. The so-called model of De Wijs is based on the assumption that if a block with average concentration value X is cut in halves with concentration values (1 + d)X and (1-d)X, then d is a fundamental constant independent of block size. Matheron (1962) showed that this model results in the logarithmic De W ijsian semivariogram γ(h)= β·ln (h) where β is a constant that is related to d, and h represents distance between concentration values at equidistant points along a line. Later, Cheng and Agterberg (1996) showed that the logarithmic semivariogram is a good approximation of the multifractal semivariogram (also see later in this paper) for series of element concentration values of adjacent blocks. The concept of self-similarity of geochemical patterns also was illustrated by Krige (1978): moving average maps of gold values from small areas on reefs in Witwatersrand gold fields are similar to moving average maps of gold values from much larger areas.

In Agterberg (1994), results were shown of a 1-dimensional computer simulation experiment for the model of De Wijs. The present paper deals with 2-dimensional multifractal map patterns. Quasiisotropy will be maintained in that the concentration values of series of adjacent blocks in any direction have the same multifractal semivariogram. Working in 2-dimensional space has the advantage with respect to 1-dimensional space in that results for many more blocks can be considered and visualized simultaneously. The purposes of the simulation experiments in this paper are three-fold: (a) convergence testing; (b) effects of artificial distortions of the basic multifractal model, and (c) lognormality and Pareto tails.

(a) Convergence testing. The properties of the model of De Wijs are well-known. They can be readily verified for any realization based on random numbers in a computer simulation experiment. There are two different ways in which a multifractal spectrum can be computed: (1) histogrand method, and (2) method of moments with Legendretransform. It has been: accepted widely that the method of moments produces better results faster than the histogram method. This hypothesis also can be tested by means of computer simulation experiments.

(b) Effects of artificial distortions of the basic multifractal model. Other advantages of computer simulation include that the basic assumptions underlying the multifractal model can be changed. The multifractal model is " stationary" in that its properties are independent of position of points of observation. As in time-series analysis and geostatistical applications, the assumption of stationarity frequently is violated. In geochemical applications, there may be trends (or "drifts") as well as abrupt changes associated with discontinuties (contacts between different rock types). There also is the possibility of superposition of more than a single multifractal model. A few departures from the basic model will be investigated by computer simulation in this paper, in order to see how the multifractal spectrum is affected by these distortions.

(c) Lognormality and Pareto tails. Cheng (1994) had shown that, under very general conditions, existence of a multifractal spectrum implies power-law relationship between cumulative frequencies and concentration values in the tails of element concentration frequency distributions. On geochemical contour maps, the log of area enclosed by an isoline plotted against the log of the concentration value for this isoline generally plots as a straightline both for the highest and the lowest concentration values. Agterberg (1995) adopted this result as a principle to reconciliate lognormality with Pareto tails in some types of mineral resource evaluation problems. For example, although the lognormal model is useful for element concentration values within orebodies as well as size-grade distributions of mineral deposits (cf. Harris, 1984), several authors including Cargill et al. (1981) obtained excellent results by means of power-law models. Likewise, the size of oil pools frequently is approximately lognormal (see e.g. Lee, 1999) but Drew et al. (1982) preferred the use of the Pareto distribution. Computer simulation helps to ilustrate how a lognormal curve can have Pareto high-value tail.

CONSTRUCTION OF MULTFRACTAL MAP PATTERNS

The method described in this section is related to the multiplicative cascade models developed and used by Meneveau and Sreeniva Tasan (1987) and Schertzer et al. (1997). However, the model of De Wijs is taken as the starting point in our applications.

Suppose that a square block is divided into 4 smaller squares as. follows. The block with average concentration value X is cut in half twice, along mutually perpendicular lines through its center, so that the halves have concentration values (1 + d)X and (1-d)X, respectively, where d is the fundamental constant of De Wijs which is independent of block size. During each cut, there is a 50% probability that one of the two rectangular half-blocks obtains the concentration value (1 + d)X. The other half then automatically obtains the value (1-d)X. The result of two successive random cuts along mutually perpendicular lines is that there is a 25% probability that any smaller square obtains the largest concentration value (1 + d)²X. Regardless of which subsquare has the largest value, moving in the clockwise (or anticlockwise) direction, the 4 subsquares have concentration values (1 + d)²x, (1+ d)(1-d)X, (1-d)²x, and (1 + d)(1-d)X, respectively. This process can be repeated indefinitely by subdividing very y smaller block into 4 parts.. Self-similarity and quasi-isotropy are preserved during each iterative step. When the element is a metal, total amount of metal remains the same.

Without loss of generality, X can be set equal to 1. Fig. 1a shows a result for d =0.4 after 2⁷= 128 pairs of cuts, or n = 14 iterations in total. There are 2¹⁴ = 16 384 cells with concentration values ranging from 0.6^l4 = 0.00078 to 1.4¹⁴ = 111.12007. The 4 096 values satisfy a so-called logbinomial frequency distribution. The logarithmically transformed values are binomial and approximately normal in the center. The pattern of Fig. 1a exhibits two features typical of geochemical data: (a) abrupt, discontinuous change, and (b) spatial autocorrelation because neighboring squares have values are relatively close. The pattern has the property of self-similarity in that setting the value of any cell equal to unity and subjecting it to 128 pairs of cuts would yield a similar pattern. The only differences between patterns resulting from different computer simulation experiments are caused by changing the seed, which was set at random at the beginning of each experment in order to generate new series of pseudo-random numbers. The sample size of 4 096 is sufficiently large to permit relatively precise estimation of parameters of the multifractal pattern. These parameters are related to the value of d.

Figure  1.  Three-dimensional plots of two simulated map patterns for element concentration values (overall mean value is equal to 1) obtained by means of a stochastic version of the model of De Wijs (see text for detailed explanations); largest values truncated at upper end. (a). array of 128×128 values of original model; (b). array of version of original model with trend superimposed according to power-law function with maximum at origin exceeding minimum in opposite corner by factor 10; (c). part of array shown in (b). (These are back and white versions of color plots originally obtained with Mathematica 4 software).

DownLoad: Full-Size Img PowerPoint

A second example with d =0.4 and n= 14 is shown in Fig. 1b. It deals with regional deviations from stationarity commonly described as trends or "drifts". A 2-dimensional power-law trend was added to the stationary variablility pattern of Fig. 1a by (a) multiplying all values by a factor equal to 10 at the origin with co-ordinates (0, 0), and equal to 1 at the opposite corner with co-ordinates (128, 128), and (b) readjusting all values by a factor so that the average concentration value remains equal to 1. Because values greater than 2 were truncated in Fig.lb (> 4 truncation was used in Fig. 1a), the pattern of concentration values close to the origin is shown separately with > 30 truncation in Fig. 1c. Later in this paper, it will be shown that the multifractal spectra of the patterns shown in Fig. 1a and 1b are only slightly different. Fig. 1b differs from Fig.la in that the condition of stationarity is clearly violated. The fact that a different set of 4 096 pseudo-random numbers was used for the simulation does not affect the nature of the pattern.

The constant d was set equal to 0.4 in order to facilitate comparison with results for this value previously described in Agterberg (1994). Also, this corresponds to the physical constant controlling multifractal turbulence (Meneveau and Sreenivisan, 1987). Later in this paper results for d =0.6 and n = 20, approximately representing the multifractal distribution of gold in the Mitchell-Sulphurets area of northwestern British Columbia (Cheng, 1994), also will be given.

The procedure used to construct Fig. 1a is not basically different from the 1-dimensional experiment performed by Agterberg (1994). There is one other way in in which the 4 values (1 + d)²X, (1 + d)(1-d)X, (1-d)²X, and (1 + d)(1-d)X can be arranged within a block with concentration value X while preserving quasi-isotropy: starting with 25% probabilty for the largest value, the other 3 values can be assigned at random to the remaining 3 blocks without replacement. A few experiments performed using this other method yielded approximately the same results, and the simpler method used for Fig. 1 was preferred.

CONSTRUCTION OF MULTIFRACTAL SPECTRUM

The method of moments (see e. g. Evertsz and Mandelbrot, 1992) consists of 3 steps. First, the mass-partition function χ_q(ε) of order q and block size (= length of side) ε is plotted against ε on log-log paper. From

where μ_i = x_iε² represents total amount of metal in a cell of size ε labeled i; x_i is the concentration value of the ith cell. The overall concentration value can be set equal to 1 as in the computer simulation experiments illustrated in Fig. 1. In total, there are n(ε) cells. It is convenient to use 2 as base of the logarithms. Setting the side of the smallest block equal to 2 (with n = 2¹⁴ blocks in total), that of the largest block considered equal to 2⁵(with n = 2⁶ blocks), and letting q go from-10 to 10 in steps of 0.5 produced 41 series that were connected by straightlines as in Fig. 2. The fact that the lines connecting the 5 points of each series are straight confirms that the map pattern of Fig. 1a is multifractal with

Figure  2.  Log-log plot (base 2) of mass-partition function versus length of cell side for pattern of 128×128 concentration values of Fig. 1a. Smallest cells have side with log₂ε = 1. Series shown are for integer values of q only. The slopes of the straightlines are plotted in Fig. 3a.

DownLoad: Full-Size Img PowerPoint

where τ(q) is the mass exponent of order q, and ∝ denotes proportionality. The slopes of the straightlines illustrated in Fig. 2 provide estimates of τ(q). The relation between τ(q) and q is shown in Fig. 3a.

Figure  3.  Method of moments continued from Fig. 2. (a). relation between mass exponent τ(q) and q; (b). relation between singularity exponent α(q) and q; (c). relation between multifractal spectrum value f(α) and singularity exponent α. The multivariate spectrum is symmetrical and has maximum value f(α)= 2; the endpoints where f(α)= 0 also are special in that their position does not depend on the number of iterations (n) used for simulation.

DownLoad: Full-Size Img PowerPoint

The singularity exponent α_i can be associated with any value of μ_i= x_iε². In multifractal modeling it is more convenient to work with α than with the concentration value x. The histogram method to be discussed in the next section is based on constructing histograms of α for different values of ε and q. The second step of the method of n moments consists of relating τ(q) to the sngularity exponent α. The relation between t(q) and a can be written in the form

Thus the singularity exponent can be derived from the mass exponent by differentiation. For the results reported here (Fig. 3b), the first derivative was determined by estimating τ(q) for pairs of closely spaced values q±0.001. Each difference between two of these values was divided by 0.002 in order to estimate α(q). Finally, the fractal dimension spectrum f(α) of Fig. 3c, which is the Legendre transform of τ(q) of Fig. 3a, follows from

The points plotted in Fig. 3c coincide with a theoretical limiting curve (f(α) for n→∞) for the model of De Wijs to be constructed in the next section (Fig. 4), showing that the method of moments rapidly leads to a satisfactory multifractal spectrum in this kind of application.

Figure  4.  (a). Histogram method applied to pattern of Fig. 1a with d=0.4 and n = 14; (b). ditto with d= 0.4 and n = 30. In both diagrams (a and b), the larger values belong to the limiting form of f(α) that provides an upper bound and would be reached for n =∞. The difference with the limiting form decreases slowly for increasing n. Note that the method of moments provides estimates (see Fig. 3c) of f(α) for n→ 0.

DownLoad: Full-Size Img PowerPoint

HISTOGRAM METHOD

The methods described in this section are based on results for 1-dimensional series summarized in Evertsz and Mandelbrot (1992). A histogram is constructed for the singularities α_i associated with μ_i = x_iε². The binomial frequencies of the concentration values x_i in the 2-D version : of the model of De Wijs depend on number of iterations n and value of d. Setting n= 14 and d=0.4 as for Fig. 1, gives the results depicted in Fig. 4a. For infinitely large value of n, a limiting multifractal spectrum (f(α) for n→∞) can be derived in analytical form. Values of f(α) from this theoretical curve are shown in Fig. 4a for comparison with the f(α) values for n= 14. There are systematic discrepancies between the two spectra of Fig. 4. However, at the extremes, α_min= log₂ {2/(1 + d)}²= 1.03 and α_max= log₂{2/(1-d)}²=3.47 with f(α)= 0, the two spectra coincide. The other point of equality of the spectra occurs at the center where f(α) = 2.

Figure. 4b shows similar results for n = 30. The histogram values are closer to the theoretical limit values in Fig. 4b, but it is obvious that convergence is exceedingly slow. On the other hand, the 3-step method of moments described in the previous section immediately resulted in f(α) values (Fig. 3c) approximately coinciding with the limit values of Fig. 4 (a and b), illustrating that the method of moments is to be preferred for derivation of the limiting form of the multifractal spectrum (f(α) for n→∞).

Because the realization of the model of De Wijs for a specific value of n is discrete, the multifractal spectrum can be readily interpreted. Each f(α) value represents the fractal dimension of a subset of cells with the same concentration value (cf. discussion in the introduction). Thus box-counting of monofractal binary pattern for any value of x created from a realization such a Is the one shown in Fig. 1a would give a fractal dimension approximately equal to f(α) as in Fig. 4a provided that the number of cells used for box-counting is sufficiently large.

The theoretical multivariate spectrum (f(α) for n→∞) can be used to compute theoretical frequencies of the concentration values x_i for specific values of n and d. These frequencies are not independent of n, and the frequency distribution curve continues to change when n is increased as will be discussed in the last section of this paper. Various scaling and rescaling procedures have to be applied in order to derive the results shown in Fig. 4 and the frequencies, subsequently, to be derived from the limiting form of f(α). The required caleulations are given in detail in the FORTRAN program with brief explanation in Appendix 1.

MULTIFRACTAL SEMIVARIOGRAM

The theoretical multifractal semivariogram satisfies (also see Cheng and Agterberg, 1996, Cheng, 1994, equation 22)

where k = 1, 2, ... represents distance between successive cells measured in multiples of ε, ξ₂(ε) is the noncentered second-order moment obtained by dividing the mass-partition function for q = 2 by number of cells, and τ(2) is the second-order lass ss exponent. Estimates of ξ₂(ε) and τ(2) can be obtained from the linear mass-partition function obtained by adding values of μ_i for successive cells within rows or columns of the pattern shown in Fig. 1a. It is noted that ξ₂(ε) represents the "sill" of this theoretical semivariogram. Experimental semivariogram values also can be estimated from the rows (or columns) of Fig. 1a.

Figure 5 shows the theoretical semivariogram in comparison with 3 experimental semivariograms. The degree of correspondence is satisfactory, keeping in mind that, for this type of application, the sample is too small to obtain relatively precise semivariogram values for longer distances.

Figure  5.  Theoretical form of multifractal semivariogram in comparison with three experimental semivariograms based on 128 rows in patterns similar to the one shown in Fig. 1a. Deviations between experimental semivariograms and continuous curve are relatively large, but there is probably no significant bias.

DownLoad: Full-Size Img PowerPoint

A RTIFICIAL DISTORTIONS OF BASIC MODEL ASSUMPTIONS

Two examples of distortions of the pattern of Fig. 1a will be discussed. First, a constant value of 0.01 was added to all concentration values x. Obviously, this modification is so small that it would not change the appearance of this realization of the model of De Wijs (d =0.4, n = 14). A new mass-partition function was obtained from the values of μ= xε² with preservation of total mass. The log-log plot of the mass-partition function (not shown here) showed some departures from linearity. However, by using sets of points located on straightline segments, the three steps of the method of moments could be completed resulting in the multifractal spectrum of Fig. 6.

Figure  6.  A small value (0.01) was added to all concentration values (cf. Fig. 1a). The multivariate spectrum obtained by means of the method of moments shows differences on the right side in comparison with Fig. 3c.

DownLoad: Full-Size Img PowerPoint

Figure 6 closely resembles Fig. 3c on the left side only. This is the part of the multifractal spectrum that corresponds to the largest concentration values. The addition of 0.01 to all values had most effect on the smallest concentration values represented by the right side of the spectrum. Some asymmetry was introduced and the smallest value of f(α) now is greater than zero. In reality, a situation of this type would arise, for example, if a metal in the rock mass would occur in two forms: : (a) highly concentrated in sulfide crystals satisfying the model of De Wijs, and (b) in small, approximately constant, concentrations in silicate crystals constituting the bulk of the rock mass.

This first example illustrates that it may be difficult to distinguish between minor distortions on the right side of a multifractal spectrum and small systematic departures from the model affecting all values. The pattern used for the second example was shown in Fig. 1b(c). It deals with non-stationarity associated with a superimposed trend. This type of distortion of the basic model which is stationary primarily affects the left side of the multifractal spectrum representing the largest concentration values.

The mass-partition function of the second example (Fig. 7) shows departures from linearity. The lines dipping to the right which correspond to negative values of q are least affected by the distortion. The first three values, for smaller values of ε, fall on straightline segments even when q is close to 10. Taking the slopes of these straightline segments and performing the remaining steps of the method of moments yielded the multifractal spectrum of Fig. 8. Although the trend (systematic changes in element concentration values) is strong, the multifractal spectrum is only moderately changed on the left side that corresponds to the largest values. Comparison with Fig. 3c shows that the limit on the left has become less and that negative values of f(α) were obtained. As in the previous example that resulted in Fig. 6, the differences between Fig. 8 and Fig. 3c are relatively minor.

Figure  7.  Log-log plot (base 2) of mass-partition function versus length of cell side for pattern of 128×128 concentration values of Fig. 1b. Straightline segments through first three points were used only in method of moments. Explanation of symbols same as in Fig. 2.

DownLoad: Full-Size Img PowerPoint

Figure  8.  Multivariate spectrum obtained by means of the method of moments starting from pattern shown in Fig. 1b. In comparison with Fig. 3c, it shows differences on the left side.

DownLoad: Full-Size Img PowerPoint

The real world situation simulated in the second example consists of relatively large concentration values of an element with random fluctuations near a source, with damping of the magnitudes of the fluctuations away from the source. This example was inspired by multifractal modeling of trace element concentration values in soil by Sim et al. (1999) in the surroundings of a smelter in Manitoba that had produced the material absorbed in the soil. At longer distances from the smelter, the trace element concentration values become relatively weak and reflect natural "background" only. A second type of situation to which this example may apply consists of fluctuations in tracer gas (CO) concentration in a room with a single emitting source surrounded by monitors located at different distances. Switzer and McBride (1999) have developed a statistical model for time series analysis of the variable concentrations detected. Closer to the source the mean concentration value remains relatively high and the magnitudes of the random fluctuations are greater. In a discussant paper, Agterberg (1999) suggested to try multifractal modeling in situations of this type.

Although the computer simulation experiments described in this section were few and of limited scope, they strongly suggest that multifractal modeling may continue to provide useful results in situations where basic assumptions of multifractality are obviously violated.

FREQUENCY DISTRIBUTION MODELING

The transformations listed in the FORTRAN program of Appendix 1 not only provide multifractal spectra, they also give frequency distributions of concentration values satisfying the model of De Wijs in 2-dimensional space. As pointed out before, this model has a limiting multifractal spectrum (f(α) for n→∞) which can be represented by means of a continuous curve. On the other hand, any application of the model of De Wijs (see e. g. Fig. 1a) results in discrete values satisfying the logbinomial model. The histogram method is based on these discrete frequency values which can be transformed into a discrete multifractal spectrum (Fig. 4). Each discrete spectrum can be transformed back into the original frequency distribution. By means of the same backward transformation, the continuous limiting form (f(a) for n→∞) of the spectrum can be represented as a (continuous) frequency distribution. Because the transformation depends on n, there is a separate frequency curve for each value of n. These continuous curves more; closely approximate the discrete frequencies when n is increased.

The logarithmic variance of the concentration values satisfies (see e. g. Agterberg, 1994)

This variance continues to increase for larger values of n. Likewise, the largest concentration value goes to ∞. Consequently, a single limiting form of the frequency distribution of the concentration values. for n→∞(similar to the limiting form of the multifractal spectrum) does not exist.

The example to be used for illustration of multifractal frequency distributions of concentration values is based on 1 030 gold values from rock samples of predominantly Jurassic metavolcanics in the Mitchell-Sulphurets area, northwestern British Columbia which measures about 4 km by 6 km. Cheng (1994, Figs. 4.5.2 and 5.2.1b) constructed the lognormal Q-Q plot of these data as well as multifractal spectra for different Au cut-off values (also see Agterberg, 1995, Figs. 3 and 4). All spectra have approximately the same value of α_min = log₂[2/(1 + d)]² =0.65 corresponding to d =0.597. For this reason, d=0.6 will be used in this example with n = 20. Except for ahorizontal line for small Au concentration values representing the 2×10^-9 detection limit, the Q-Q plot essentially is according to a straightline pattern representing a lognormal distribution with σ² (ln X)=8.96. Substitution of this estimate (together with d =0.6) into the preceding equation for the logarithmic variance gives n = 18.94. This estimate is probably too small, because it is likely that rocks with the smallest gold concentration values are underrepresented if the. data set of 1 030 gold values is used to approximate the true frequency distribution of gold in all rocks of the Mitchell-Sulphurets area.

Figure 9a shows discrete and continuous multifractal spectra for the model of De Wijs with d = 0.6 and n = 20. These results are similar to those shown in Fig. 4 with d =0.4 and n= 14, 30. As explained in the preceding paragraph, these new constants were selected on the basis of statistics derived from the distribution of gold in the Mitchell-Sulphurets area. The two frequency distributions of Fig. 9b, which are equivalent to the two spectra of Fig. 9a, closely resemble one another. The frequencies corresponding to f(α) for n→∞ are somewhat higher on the flanks of the curve but in the ceenter and on the tails the two curves approximately coincide. Both frequency distributions are approximately lognormal in the center and Pareto-type on the tails.

Figure  9.  (a). Histogram method illustrated in Fig.4 applied to concentration values with d = 0.6 and n = 20; (b). frequency distribution curves corresponding to the two multivariate spectra shown in (a); frequences of limiting form slightly exceed logbinomial frequencies but difference is zero in the center and at the endpoints; (c). lognormal Q-Q plot of upper bound frequency distribution shown in (b); near the center, the frequency distribution resulting from the model of De Wijs is lognormal, and in the tails it is Pareto-type (according to a power-function). Note that the Pareto tails are weaker than the tails of the lognormal distribution plotting as the straightline.

DownLoad: Full-Size Img PowerPoint

The property that multifractals with f(α_min)= f(α_max) = 0 have Pareto tails for α→α_min and α→α_max was originally demonstrated in Cheng et al. (1994). It is noted that, by plotting frequencies of the largest and smallest gold values on log-log paper followed by straightline fitting, Agterberg et al. (1993) estimated α_min=0.68 and α_max =3.90 from the slopes of the best-fitting straightlines. These estimates are comparable to the values of the endpoints in Fig. 9a. The nature of the Pareto-type departure from lognormality in the tails is clearly demonstrated in the lognormal Q-Q plot shown in Fig. 9c. The straightline in Fig. 9c represents the lognormal distribution in the center. It corresponds to the approximately normal frequency distribution in the center of Fig. 9b. Outward projection of this lognormal curve to the tails illustrates that the Pareto tails of the multifractal frequency distribution are weaker (more rapidly approaching 0) than the tails of the lognormal distribution.

As pointed out in the introduction, the multifractal frequency distribution model with lognormality in the center and Pareto tails provides an explanation for the fact that both the lognormal as well as the Pareto distribution model have been used successfully in quantitative mineral resource estimation, although these two distributions are of different types (also see Agterberg, 1995). The potential implications of this feature are important: the multifractal frequency distribution has weaker tails than the lognormal, but the relative increase in frequency from the high-value tail towards the center is more rapid when the Pareto distribution model is used. This may result in larger estimates of the volumes of lower grade material or larger numbers of relatively small oil pools as originally pointed out by Mandelbrot (1983).

CONCLUDING REMARKS

It would be naive to assume that the multifractal approach advocated in this paper constitutes an accurate representation of natural geochemical patterns. However, this is not the purpose of useful, statistical tools. For example, ordinary frequency distribution analysis frequently has to deal with the fact that the observed data are mixtures from more than a single population. Likewise, the geostatistical approach based on "regionalized" random variables often has to account for trends or "drifts". All approaches have to cope with (existence of strong discontinuities such as contacts between rock formations with different lithologies.

The multifractal model of De Wijs used in this paper to simulate geochemical map patterns is simple but offers an important tool for practical applications. The basic model has the properties of self-si S1milarity (or scale-independence) and stationarity (mean concentration value is independent of location). The multivariate semivariogram has a "sill" and resembles semivariograms frequently encountered in geostatistical practice. It has a symmetrical multifractal spectrum with limiting form that can be derived readily from 1-, 2-, (or 3-dimensional data sets by means of the 3-step method of moments. It has a frequency distribution that is lognormal in the center and has Pareto tails. Some computer simulation experiments described in this paper indicate that distortions of the basic model (non-stationarity trends or mixing of models) have a relatively small effect (on the shape of the multifractal spectrum. An advantage of the multifractal approach with respect to the usual geostatistical approach is that it permits easier frequency distribution analysis. Its advantage with respect to ordinary frequency distribution analysis is that the approach is spatial and allows for the modeling of spatial a autocorre elatio on.

Appendix Ⅰ. FORTRAN computer program for histogram method applied to model of De Wijs.

This program generates multifractal spectra and frequency distributions represented in Fig. 9. The constants d=0.6 and n = 20 are defined near the beginning of the code. The output is written in the file factoria out. The following statistics are computed and printed out in a table which was imported into Excel for creating graphs.

1. xi: row number ranging from 0 to n

2. value(i): concentration value

3. value2(i): log of concentration value (base 2)

4. alpha(i): singularity exponent α ranging from amin to αmax

5. freq(i): frequency of concentration value

6. cumfreq(i): cumulative frequency of concentration value

7. cumfreql(i): log of cumulative frequency of concentration value (base 2)

8. falpha(i): discrete multifractal spectrum value f(α)

9. falphe(i): value of limiting form of multifractal spectrum(f(α) for n-∞)

10. fre(i): frequency of concentration value derived from falphe(i)

program factoria

  real fac(50), value(50), value2(50), alpha(50), freq(50), cumfreq(50)

  real cumfreql(50), falpha(50), falphe(50), fre(50)

  open(6, file= 'factoria.out', status='unknown')

  n= 20

  xd=.6

  eta= (1.- + xd)/(1.-xd)

  amin= alog10(((eta+ 1.)/eta) * * 2)/ alog10(2.)

  amax= alog10((eta+ 1.)* *2)/alog10(2.)

  fac(1)=1

  do 1 i=2, n

  fac(i)=i*fac(i-1)

1 continue

  freq(1) = n

  cumfreq(1) = freq(1) + 1.

  do 2 i= 2, n/2

  freq(i)= fac(n)/(fac(i) * fac(n-i))

  cumfreq(i) = cumfreq(i-1)+ freq(i)

2 continue

  do 3 i= 1, n/2-1

  freq(n/2+ i)= freq(n/2-i)

  cumfreq(n/2 + i)= cumfreq(n/2+i-1)+ freq(n/2+ i)

3 continue

  freq0= 1.

  xi =0.

  xn=n

  value1 =(1.+xd)* * n

  value21 = alog10(value1)/alog10(2.)

  alphal =-value21 * 2./xn+ 2.

  write(6, 1000) xi, valuel, value21, alphal, freq0, freq0, xi, xi, xi, freq0

  do 4 i= 1, n-1

  xi=i

  value(i)= (1. + xd)* *(n-i)*(1.-xd)* *i

  value2(i) = alog10(value(i))/alog10(2.)

  alpha(i)=-value2(i) *2./xn +2.

  constant = alog10(4.)/ alog10(freq(n/2))

  falpha(i) = constant * alog10(freq(i))/ alog10(2.)

  cumfreql(i) = alog10(cumfreq(i))/ alog10(2.)

  al = (alpha(i)-amax)/(amax-amin)

  a2= (amin-alpha(i))/(amax-amin)

  falphe(i)=al * alog10(al * * 2)/alog10(2.) + a2 * alog10(a2* * 2)/alog10(2.)

  fre(i) = alog10(freq(n/2)) * alphe(i)/ alog10(4.)

  fre(i)=2. * * fre(i)

  write(6, 1000) xi, value(i), value2(i), alpha(i), freq(i), cumfreq(i), cumfreql(i), falpha(i), falphe(i), fre(i)

4 continue

  valuen= (1.-xd)* * n

  value2n = alog10(valuen)/ alog10(2.)

  alphan=-value2n* 2./xn + 2.

  cumtreqn = cumtreq (n-1)+1.

  cumfreqln = alog10(cumfreqn)/ alog10(2.)

  zero= 0.

  write(6, 1000) xn, valuen, value2n, alphan, freq0, cumfreqn, cumfreqln, zer+0, zero, freq0

1000 format (f3.0, f12.5, f10.5, f9.5, f8.0, f9.0, f10.5, 2f9.5, f10.2)

  stop

  end