## Abstract

**ABSTRACT** In this study, continuous field models of geochemical landscapes were obtained by interpolating stream sediment geochemical data while discrete field models of geochemical landscapes were obtained by attributing stream sediment geochemical data to their sample catchment basins. This study aimed to: (1) compare and contrast anomaly maps derived from continuous and discrete field models of stream sediment geochemical landscapes; and (2) determine which empirical frequency distributions – those of original point data or those of pixels values in models of stream sediment geochemical landscapes – are more useful in mapping of anomalies in such geochemical landscapes. Anomalies were mapped by using the mean+2*SDEV* (standard deviation), median+2*MAD* (median absolute deviation) and concentration–area (C–A) fractal methods of identifying threshold values in a geochemical data set. The results of the study in the Aroroy gold district (Philippines) highlight the following findings. In mapping of anomalies in either continuous or discrete field models of stream sediment geochemical landscapes, the C–A fractal method performs best, followed by the median+2*MAD* method and then by the mean+2*SDEV* method. Anomalies mapped in discrete field models, compared to anomalies mapped in continuous field models, of stream sediment geochemical landscapes mostly have stronger positive spatial associations with the known epithermal Au deposit occurrences in the study area. Empirical frequency distributions of either the original point data or the pixels values in the models of stream sediment geochemical landscapes are similarly useful in applying the C–A fractal method, but not in applying either the median+2*MAD* or mean+2*SDEV* method, to map anomalies in either continuous or discrete field models of such geochemical landscapes.

- exploratory data analysis
- spatial interpolation
- catchment basins
- concentration-area fractal method
- spatial association analysis
- GIS
- Aroroy (Philippines)

Various geographical phenomena or processes are modelled, based on certain pertinent spatial data, as either continuous or discrete fields. On the one hand, a *continuous field* model is one in which spatial data are portrayed invariably as a *raster* of pixels (or a grid of cells), in which the values in every pixel (or grid cell) change continually from one location to its adjacent locations. Thus, in a continuous field model a rate of change in pixel values per unit distance (i.e. pixel or cell size) can be measured anywhere and in any direction. An example of a continuous field model is a Landsat Thematic Mapper (TM) image, in which pixel values represent spectral reflectance of materials on the Earth's surface. On the other hand, a *discrete field* model is one in which spatial data are portrayed either as a raster map or a *vector* (i.e. polygon) map of mutually exclusive and bounded units. In either a raster or vector map of discrete fields, pixels or polygons representing a discrete field are characterized by the same data value. Thus, in a discrete field model, portrayed as a raster or vector map, a rate of change in values per unit distance cannot be measured anywhere and in any direction because values do not change (at least symbolically) within each unit (i.e. slope is equal to zero), but values change abruptly only across unit boundaries. An example of a discrete field model is a map of lithologic units, each of which symbolizes various geological processes. Now, in a raster map of either a continuous field model or a discrete field model of certain spatial data, every pixel is a model of a *discrete object* (i.e. a sample) of spatial data attributed to the map (*x*, *y*) coordinates of every pixel. However, whereas every pixel in every raster map of either continuous or discrete fields has attribute spatial data, not every pixel in every raster map of discrete objects has attribute spatial data. An example of a discrete object model, but which is neither a continuous nor a discrete field model, is a raster map of geochemical concentration values at sampling points only. This paper deals with mapping of anomalies in raster models of continuous and discrete fields of geochemical landscapes derived from stream sediment geochemical data.

Geochemical landscapes (i.e. spatial distributions of geochemical variables) can be modelled as either continuous or discrete fields based on point geochemical data. On the one hand, a continuous field model of a geochemical landscape is derived invariably by interpolation of geochemical data at sampled locations in order to estimate values at every unsampled location. On the other hand, a discrete field model of geochemical landscape can be derived, for example, by creating a Thiessen or Voronoi polygon (i.e. an area of influence) around every sampled location and then attributing data at individual sampled locations to their associated polygons. Thus, continuous and discrete field modelling of geochemical landscapes involves, respectively, point-to-surface and point-to-area spatial transformations of point geochemical data. The reason for applying either point-to-surface or point-to-area spatial transformation to point geochemical data sets is to facilitate analysis and outlining (i.e. mapping) of anomalous areas because point-symbol representations of geochemical data are inadequate for such purpose. Readers are referred to Bonham-Carter (1994) for a comprehensive introduction to various techniques for spatial transformations of point data.

Modelling of geochemical landscapes based on stream sediment (or water) geochemical data requires more prudence compared to modelling of geochemical landscapes based on, for example, regolith geochemical data. That is because stream sediments (or waters), unlike regolith, are: (a) found only along fluvial environments; and (b) are associated only with materials within their catchment basins. Thus, unlike uni-element concentrations in regolith samples, uni-element concentrations in stream sediment (or water) samples do not strictly represent spatially continuous fields (or variables). Nevertheless, geochemical landscapes can be and have been modelled as continuous fields by interpolating stream sediment (or water) geochemical data if: (a) stream sediment sampling density is considered to be sufficiently high; and (b) uni-element concentrations in stream sediment (or water) samples exhibit positive spatial autocorrelation (i.e. they can be considered regionalized variables). Alternatively, geochemical landscapes can be and have been modelled as discrete fields by attributing or imposing stream sediment (or water) uni-element data (or derivative values representing multi-element associations) to their sample catchment basins (e.g. Bonham-Carter *et al*. 1987; Carranza & Hale 1997; Spadoni *et al*. 2004).

The appropriateness and/or advantage of either continuous or discrete field modelling of geochemical landscapes based on stream sediment (or water) geochemical data would depend largely on the mapping scale of exploration and/or environmental geochemical surveys. On the one hand, continuous field modelling of geochemical landscapes based on stream sediment (or water) data is arguably advantageous in regional-scale geochemical surveys, in which the objective is to delimit broad anomalous areas for further investigations at higher scales. In regional-scale geochemical surveys over large areas (say, western Europe; see Lima *et al*. 2008), discrete field (e.g. catchment basin) modelling of geochemical landscapes based on stream sediment (or water) geochemical data is not totally inappropriate but, because of variations in sizes of sample catchment basins due variations in sampling density and sample distribution, could be both tedious and impractical with respect to the scale(s) of output map(s) (e.g. 1:100 000 or smaller). On the other hand, discrete field (i.e. catchment basin) modelling of geochemical landscapes based on stream sediment (or water) data is arguably advantageous in many cases of district-scale to local-scale environmental and/or exploration geochemical surveys, in which the objective is to establish precisely sources of contamination and/or significant anomalies in individual catchment basins (e.g. Goodyear *et al*. 1996; Ódor *et al*. 1998; Moon 1999). In many cases of district-scale to local-scale geochemical surveys, interpolation of stream sediment (or water) geochemical data may or may not be defensible for reasons stated earlier in the preceding paragraph, whereas catchment basin representation of stream sediment (or water) geochemical data is viable with respect to the scale(s) of the output map(s) (e.g. larger than 1:100 000). The question is: in district- to local-scale geochemical surveys, which model – continuous field or discrete field – of geochemical landscapes based on stream sediment (or water) geochemical data, provided that both models are defensible and viable with respect to output map scales of 1:100 000 or larger, results in better maps of significant anomalies (i.e. maps of geochemical anomalies with stronger positive spatial associations with known occurrences of mineral deposits of the type sought)?

Mapping of anomalies in geochemical landscapes entails dividing a geochemical data set into a number of classes and then judging which of these classes represent, say, presence of mineralization in a study area. On the one hand, the traditional methods applied in mapping of uni-element anomalies include (cf. Levinson 1974; Govett *et al*. 1975; Rose *et al*. 1979; Sinclair 1983; Howarth & Sinding-Larsen 1983): (a) comparison of data from the literature; (b) comparison of data with results of an orientation geochemical survey; (c) graphical discrimination from a data histogram; (d) calculation of threshold as the sum of the mean and some multiples of the standard deviation of data; (e) plotting cumulative frequency distributions of the data and then partitioning the data into background and anomalous populations; (f) constructing normal probability plots of the data and then partitioning the data into background and anomalous populations; and (g) recognition of clusters of anomalous samples when data are plotted on a map. On the other hand, the traditional methods applied in mapping of multi-element anomalies include (cf. Rose *et al*. 1979; Howarth & Sinding-Larsen 1983): (a) stacking uni-element anomaly maps of the same scale on top of each other on a light table and then outlining areas with multiple intersections of uni-element anomalies; (b) determining inter-element correlations, usually by calculating Pearson correlation coefficient; (c) recognizing and quantifying multi-element associations based on their correlations; and (d) mapping of quantified multi-element association scores. Recently, Reimann & Filzmoser (2000) and Reimann *et al*. (2005) compared various statistical methods for analysis of threshold values in order map anomalies in geochemical landscapes, while Reimann (2005) reviewed various techniques for selection of data classes in order to display and interpret the various processes (e.g. significant anomalies) in geochemical landscapes. These workers explain that: (a) mapping of anomalies in geochemical landscapes generally involves applications of various statistical methods for analysing empirical frequency distributions as well as spatial distributions of geochemical data; and (b) neglecting the fact that empirical frequency distributions of many geochemical data do not follow a normal (or log-normal) distribution model can lead to faulty results. Thus, some workers propose to map anomalies by analysing not only empirical frequency distributions of geochemical data but also fractal (e.g. concentration-area (C–A)) properties of geochemical landscapes (Cheng *et al*. 1994) and/or scales of spatial structures of geochemical data (Cheng 1999).

In practise, mapping of anomalies in models of geochemical landscapes generated from measured point geochemical data are based on analysis of empirical frequency distributions of the model values rather than based on empirical frequency distributions of the measured values. This practise requires scrutiny in view of the following arguments. Firstly, mineral explorationists strive to obtain point geochemical data that are representative of the geochemical landscapes they survey. Secondly, accuracy and precision of measured point geochemical data are better controlled and defined than those of model values derived from the point data. Therefore, it can be hypothesized that mapping of anomalies in continuous or discrete field models of geochemical landscapes can be based on empirical frequency distributions of measured point data from which models of geochemical landscapes were derived. If that is so, then the following question about mapping of anomalies in geochemical landscapes based on stream sediment geochemical data arises. Which empirical data frequency distributions – those of the measured values or those of the model values – are more useful in mapping of significant anomalies in continuous or discrete field models of geochemical landscapes?

It is important to provide an answer to each of the two questions stated above because, although continuous or discrete field modelling of geochemical landscapes based on stream sediment geochemical data is justifiable, a prime objective in, say, mineral prospectivity mapping, is to create and integrate layers of spatial evidence (e.g. geochemical anomalies) having good positive spatial associations with known occurrences of mineral deposits of the type sought. This paper describes and discusses the results of a case study, in the Aroroy gold district (Philippines), addressing the two questions stated above by applying the mean+2*SDEV*, the median+2*MAD* and the concentration-area (C–A) fractal methods for separating background and anomalies in a geochemical data set (see details below). The case study area and spatial data sets are the same as those examined by Carranza (2004*a*, 2008) in similar but different case studies of GIS-based modelling of stream sediment geochemical anomalies in the same district.

## CASE STUDY AREA AND GEOCHEMICAL DATA

The Aroroy gold district is located in the northwestern portion of Masbate Island in the Philippines (Fig. 1). In this district, the oldest rocks belong to the Eocene–Oligocene Mandaon Formation (Fig. 1), which consists mainly of andesitic-dacitic agglomerates (Baybayan & Matos 1986; JICA-MMAJ 1986). The Mandaon Formation is intruded by the Miocene Aroroy Diorite, which varies in composition from quartz diorite to hornblende diorite. Unconformably overlying the Mandaon Formation and the Aroroy Diorite are feldspathic wackes belonging to the Early Miocene Sambulawan Formation. Andesitic lithic tuffs (Late Miocene to Early Pliocene Lanang Formation) disconformably overlie the Mandaon Formation and the Sambulawan Formation. The Pliocene Nabongsoran Andesite consists of porphyritic stocks, plugs and dikes that intrude into the series of dacitic-andesitic volcano-sedimentary rocks (i.e. the Mandaon, Sambulawan and Lanang Formations) and the Aroroy Diorite. The Nabongsoran Andesite porphyry intrusions, many of which are not mappable at the map scale of Figure 1, are probably responsible for hydrothermal alteration and epithermal Au mineralization in the intruded rocks (Mitchell & Balce 1990; Mitchell & Leach 1991). Gold, in at least 13 mineral deposit occurrences in the district (Fig. 1), is associated with sulphide (mainly pyritic and/or arsenopyritic) minerals in wide-sheeted and manganese-bearing quartz or silicified veins in generally northwest-trending faults/fractures that cut the volcano-sedimentary rocks.

Stream sediment geochemical data for 135 samples in the district are used. These stream sediment geochemical data, representing a total drainage basin area of about 101 km^{2} (i.e. an average sampling density of one sample per 1.3 km^{2}), are a subset of multi-element geochemical data for more than 2200 stream sediment samples representing drainage basins measuring a total of at least 3 000 km^{2} in Masbate Island (JICA-MMAJ 1986). The quality of the individual uni-element data sets (based on results of replicate analyses; JICA-MMAJ 1986) has been re-examined by Carranza (2004*a*) using robust analysis of variance (Ramsey *et al*. 1992) and uni-element data sets with procedural variances exceeding or approaching 20% of the total (i.e. geochemical + procedural) variance were not used in this case study. Thus, out of 10 elements determined in the stream sediment samples (JICA-MMAJ 1986), only six elements (Cu, Zn, Ni, Co, Mn, As) are studied here.

## METHODS AND RESULTS

The stream sediment geochemical data, the lithological map and locations of known epithermal Au deposits (Fig. 1) were compiled in a geographic information system (GIS), in which the analyses described below were performed.

### Exploratory data analysis

In the initial stages of mapping of geochemical anomalies it is useful to explore the point data sets in order to obtain impressions about data structure, behaviour, local densities, gaps and, of course, outliers (i.e. anomalies). Descriptive statistics and statistical graphics are the useful tools in exploratory analysis of empirical frequency distributions of geochemical data before proceeding with detailed analysis (e.g. to define threshold). Because many uni-element data sets do not follow a normal (or symmetric) frequency distribution model (e.g. Vistelius 1960; Reimann & Filzmoser 2000), Reimann *et al*. (2005) pointed out that ‘data should approach a symmetrical distribution before any threshold estimation methods are applied’. Certain transformations are usually applied in order to ‘normalize’ the values in a uni-element data set (Miesch 1977; Joseph & Bhaumik 1997), but even then most, if not all, transformed uni-element data sets only approximate a normal distribution (e.g. McGrath & Loveland 1992). In this study, asymmetry in the empirical frequency distributions of the individual uni-element data sets (see Carranza 2004*a*) was addressed by transforming the data into natural logarithms (e.g. Miesch 1977; Garrett *et al*. 1980), although this is not a judgment that the individual data sets follow a log-normal distribution model.

Reimann *et al*. (2005) demonstrated that geochemical anomalies are modelled better by application of statistics for exploratory data analysis or EDA (Tukey 1977) instead of by application of classical statistics for confirmatory or probabilistic data analysis. They have shown, for example, that geochemical anomalies based on a threshold value defined from EDA statistics as median+2*MAD* are more robust compared to geochemical anomalies based on a threshold value defined from classical statistics as mean+2*SDEV*. The *MAD* stands for median absolute deviation while the *SDEV* stands for standard deviation. The *MAD* is estimated as the median of absolute deviations of all values from the data median (Tukey 1977). The *MAD* in EDA statistics is analogous to the *SDEV* in classical statistics because both of them are measures of dispersion of data values; and, because the median and mean are also analogous as they are both measures of central tendencies of data values, a threshold value defined as median+2*MAD* is, therefore, also analogous to a threshold value defined as mean+2*SDEV*. Taking the mean+2*SDEV* was, in the early days of exploration geochemistry, one of the three methods proposed by Hawkes & Webb (1962) for defining a threshold value in order to separate geochemical anomalies from background. Recently, Reimann *et al*. (2005) proposed to discontinue the use of the mean+2*SDEV* method because ‘…it was originally suggested to provide a ‘filter’ that would identify approximately 2½% of the data for further inspection…’. In this study, geochemical anomalies based on a threshold value defined as mean+2*SDEV* and geochemical anomalies based on a threshold value defined as median+2*MAD* are compared further in order to prove or disprove the findings of Reimann *et al*. (2005) that the former are more robust than the latter.

Table 1 shows that, in the study area, uni-element threshold values defined as mean+2*SDEV* are greater than uni-element threshold values defined as median+2*MAD*. This implies that a uni-element threshold value defined as median+2*MAD*, compared to a threshold value defined as mean+2*SDEV* for the same element, results in classification of higher number of anomalous samples for that element (Fig. 2). This does not necessarily mean, however, that every uni-element threshold value defined as median+2*MAD* is better than the corresponding uni-element threshold value defined as mean+2*SDEV* in mapping of significant anomalies. Nevertheless, anomalies of As based on a threshold value estimated as median+2*MAD* (Fig. 2C), compared to anomalies of As based on a threshold value estimated as mean+2*SDEV* (Fig. 2A), exhibit better spatial association with the known epithermal Au deposit occurrences in the area. This implies that, at least in the study area, uni-element threshold values defined as median+2*MAD*, compared to uni-element threshold values defined as mean+2*SDEV*, are more useful in mapping of significant anomalies.

The concept of defining uni-element threshold values using either classical statistics or EDA statistics was also applied to define a threshold value in multi-element scores derived via application of principal components (PC) analysis of the uni-element data sets. From the results of PC analysis of the log_{e}-transformed uni-element concentration data sets in the study area ( Table 2), PC3 and PC4 (together accounting for 18% of total variance) were both interpreted as reflecting the presence of epithermal Au mineralization, because both PCs represent As, which is a pathfinder for epithermal Au deposits, and because PC3 and PC4 reflect Cu–As enrichment and As–Ni enrichment, respectively, unrelated to enrichment of Mn (i.e. scavenging effect) in the surficial environments. PC1 (accounting for 62% of total variance) represents overall correlations among uni-element concentrations, which is plausibly due to mixing of weathering products in the surficial environments. PC2 (accounting for about 13% of total variance) represents a Ni–Co–Cu association reflecting lithologic controls, that is antipathetic with a Mn–Zn–As association reflecting metal-scavenging controls of Mn-oxides in the surficial environments. PC5 and PC6 (together accounting for at most 7% of total variance) represent, respectively, Mn–Co and Ni–Mn associations, which, based on the small percentage of explained variance, mostly likely reflect scavenging of Co and Ni by Mn-oxides in some small parts of the study area. Thus, the PC3 and PC4 scores were integrated in order to (a) derive (a) a single variable that represents an As–Cu–Ni association, reflecting presence of epithermal Au mineralization, and (b) a set of multi-element scores for mapping of geochemical anomalies. The PC3 and PC4 scores were integrated into ‘As–Cu–Ni’ scores by first re-scaling both sets of PC scores linearly to the range [0,1], in order to avoid creation of false anomalies from PC3 and PC4 scores that are both negative, and thereafter performing multiplication on the re-scaled variables. Then, a mean+2*SDEV* threshold value and a median+2*MAD* threshold value were defined using the As–Cu–Ni scores. Sample locations with anomalous As–Cu–Ni scores based on threshold score estimated as median+2*MAD* (Fig. 2D), compared to sample locations with anomalous As–Cu–Ni scores based on threshold score estimated as mean+2*SDEV* (Fig. 2B), exhibit better spatial associations with the known epithermal Au deposit occurrences in the case study area.

Thus, the results of the exploratory data analyses of the point stream sediment geochemical data show and support the findings of Reimann et al. (2005) that, in mapping of significant anomalies, threshold values estimated as median+2*MAD* apparently outperform threshold values estimated as mean+2*SDEV*. However, point geochemical anomaly maps (e.g. Fig. 2) only vaguely suggest the spatial extents of anomalous areas or the source catchments of anomalous stream sediment samples. Therefore, in mapping of anomalous areas, it is preferable to perform point-to-surface transformation (i.e. continuous field modelling) and/or point-to-area transformation (i.e. discrete field modelling) of stream sediment geochemical data. In this study, continuous and discrete field modelling of the stream sediment geochemical data were performed using raster GIS operations (i.e. the geochemical landscape models consisted of square pixels or grid cells). A grid cell or pixel size of 10 m was used because: (a) in practise stream sediment samples per sampling location usually consist of several grabs of sediments along a more-or-less 10-m stretch of the stream; and (b) this pixel size translates to map scales larger than 1:100 000 (see Hengl 2006), which are suitable for district-scale mapping.

### Continuous field modelling of geochemical landscapes

The literature in exploration/environmental geochemistry is replete with numerous case studies in which point data of stream sediment uni-element concentrations have been transformed, usually via conventional ‘weighted moving average’ (WMA) interpolation techniques, into continuous fields. Of the various conventional WMA methods, ordinary inverse distance weighting (IDW) and kriging are the most commonly used methods. Ordinary IDW requires some knowledge of spatial correlation and variability of the point uni-element concentration data in order to determine the appropriate size of the ‘moving average’ window (i.e. search kernel) and distance-decay parameters (i.e. limiting distance and weight exponent). Kriging assumes that spatial variability is too complex to be modelled mathematically such that it must be treated as a stochastic process and the interpolation parameters (form of variation, magnitude and spatial association) are analysed via variography. Both ordinary IDW and kriging smooth, to different degrees, local spatial variability of geochemical data. In contrast, multifractal interpolation (Cheng 1999; Agterberg 2001) preserves local variability of spatial data by taking into account both spatial association and scale invariant properties of geographical processes, such as those present in geochemical landscapes. Recently, Lima *et al*. (2003, 2008) demonstrated that applications of multifractal IDW, compared to applications of ordinary IDW and/or kriging, results in more realistic models of geochemical landscapes based on stream sediment/water geochemical data.

In this study, continuous field models of geochemical landscapes based on stream sediment uni-element concentrations were generated via ordinary IDW instead of via kriging or multifractal IDW because: (a) it is simple but can ensure that the output values are more-or-less equal to the input values at original data points (Shepard 1968; Watson & Philip 1985; Farwig 1986; Zuppa 2004); and (b) among these three interpolation methods it is, based on the findings of Lima *et al*. (2008), a moderately satisfactory method for interpolation of stream sediment geochemical data. Interpolation of the individual point geochemical data sets in this study is justified because the average sampling density is quite high and the Moran's and Geary's indices derived per data set by using a lag distance of 1000 m are all positive and are all less than [1] for lag distances of at least 2000 m. In interpolating the individual point geochemical data sets via ordinary IDW, a search radius of 1700 m was used so as to be somewhat consistent with the average stream sediment sampling density. In addition, a weight exponent of [2] was used so as to give more weights to the original point data and thus avoid creating very smooth geochemical landscapes. The results of the ordinary IDW interpolation are satisfactory because, for all the geochemical variables examined, the estimated values are nearly identical to the measured values at all but one or two sample locations (Fig. 3).

In respect of the zone of influence of each stream sediment sample, interpolated uni-element concentrations outside any of the stream sediment sample catchment basins (see below) were masked out in the succeeding analyses in this study. For the purpose of illustration and because of limitations in displaying spatial data distributions in grey-scale, the continuous field models of As and As–Cu–Ni landscape shown in Figures 4A and 4B, respectively, are classified using 25-percentile intervals of values in the As and As–Cu–Ni point data sets. Figures 4A and 4B show that parts of the study area characterized by the upper two quartiles of the As and As–Cu–Ni point data have good spatial associations with the known epithermal Au deposit occurrences.

### Discrete field modelling of geochemical landscapes

If element concentrations in stream sediments are a result of mixing of weathering products of various materials within their catchment areas, then the measured element concentrations in a stream sediment sample represent ‘average’ element contents of surficial sediments within that sample's catchment basin. Based on this supposition, catchment basins of individual stream sediment samples were first created automatically in a GIS by using a digital elevation model and the sample locations as input data. (Readers who are interested in computer-aided algorithms for delineating stream sediment sample catchment basins are referred to Jones 2002). Then, stream sediment uni-element concentrations and the derived integrated As–Cu–Ni scores were attributed to their sample catchment basins in order to model geochemical landscapes as discrete fields. For the purpose of illustration and because of limitations in displaying spatial data distributions in grey-scale, the discrete field models of As and As–Cu–Ni landscapes shown, respectively, in Figures 4C and 4D are classified using 25-percentile intervals of values in the As and As–Cu–Ni point data sets. Figures 4C and 4D show that catchment basins of stream sediment samples associated with the uppermost quartile of the As and As–Cu–Ni point data have good spatial associations with the known epithermal Au deposit occurrences in the district.

### Analysis and mapping of geochemical anomalies in geochemical landscapes

The threshold values estimated as either mean+2*SDEV* or median+2*MAD* of the As and As–Cu–Ni point values (Table 1) are larger than the values corresponding to the 75th percentile of the As and As–Cu–Ni point data (Fig. 4). This implies that, in the As and As–Cu–Ni landscapes, anomalous areas based on threshold values defined either as mean+2*SDEV* (Fig. 5) or median+2*MAD* (Fig. 6) of the As and As–Cu–Ni point data, compared to areas with values greater than 75th percentile of each of the point data sets (Fig. 4), have weaker spatial associations with the known epithermal Au deposit occurrences and, thus, are poor spatial evidence of epithermal Au prospectivity. Similar results were found in mapping of anomalies in the remaining uni-element landscapes based on threshold values estimated as either mean+2*SDEV* or median+2*MAD* of the point uni-element data sets. Nevertheless, the ‘areal’ anomaly maps in Figures 5 and 6, compared to the point anomaly maps in Figure 2, provide better visualization of spatial extents of anomalous areas or the plausible source catchments of anomalous stream sediment samples. However, anomalous areas based on threshold values defined as median+2*MAD* of the As and As–Cu–Ni point data (Fig. 6), compared to anomalous areas based on threshold values defined as mean+2*SDEV* of the As and As–Cu–Ni point data (Fig. 5), apparently show stronger spatial association with the known epithermal Au deposits in the study area. These visual interpretations are quantified via spatial association analysis (see further below).

Experiments were performed to determine if mapping of geochemical anomalies can be improved by calculating threshold values as either mean+2*SDEV* or median+2*MAD* of the pixel values in both continuous and discrete field models of geochemical landscapes. Table 3 shows that the threshold values calculated as either mean+2*SDEV* or median+2*MAD* of the pixel values in either continuous or discrete field models of geochemical landscapes are larger than the corresponding threshold values derived from the individual point geochemical data sets (Table 1). This means, for example, that anomalous areas based on threshold values derived as either mean+2*SDEV* or median+2*MAD* of the pixels values in both continuous and discrete field models of As and As–Cu–Ni geochemical landscapes are smaller than the corresponding anomalous areas shown in either Figure 5 or 6, which are based on threshold values derived as either mean+2*SDEV* or median+2*MAD* of the original As and As–Cu–Ni point data. It is shown below that the anomalous areas shown in either Figure 5 or 6 mostly lack or have weak positive spatial associations with the known epithermal Au deposits in the study area. This means that anomalous areas based on threshold values derived as either mean+2*SDEV* or median+2*MAD* of the pixels values in both continuous and discrete field models of As and As–Cu–Ni geochemical landscapes also lack or have much weaker positive spatial associations with the known epithermal Au deposits in the study area. The results illustrate, therefore, that threshold values derived as either mean+2*SDEV* or median+2*MAD* of pixels values in either continuous or discrete field models of geochemical landscapes, compared to threshold values derived as either mean+2*SDEV* or median+2*MAD* of the original As and As–Cu–Ni point data, do not improve mapping of anomalies.

The larger threshold values derived from the pixel values in either continuous or discrete field models of geochemical landscapes, compared to the threshold values derived from individual point geochemical data sets, are attributed to changes in empirical frequency distributions of the values as a result of the point-to-surface or point-to-area transformations. Figure 7 shows that, compared to the relative frequencies of point values in the stream sediment geochemical data, there are lower relative frequencies of minimum to median (or 50th percentile) pixel values and there are higher relative frequencies of median to maximum pixel values in either the continuous or discrete field models of geochemical landscapes. The empirical frequency distributions of pixels values in the continuous field models, compared to empirical frequency distributions of pixel values in the discrete field models, of geochemical landscapes deviate more from the empirical frequency distributions of the point stream sediment geochemical data. This is a manifestation of the smoothing (or ‘blending’) effect, which is expected, of spatial interpolation. The results, thus, suggest caveats in mapping of anomalies in either continuous or discrete field models of stream sediment geochemical landscapes by defining threshold values based on the empirical frequency distributions of derived field model values rather than on the empirical frequency distributions of measured point values. Moreover, the empirical frequency distributions of either the measured point values or the derived continuous/discrete field model values show varying degrees of skewness (Fig. 7), implying the presence of multiple and intertwined geochemical populations in either the measured or derived data. This is a plausible explanation why the uni-element threshold values calculated as mean+2*SDEV* or median+2*MAD* were inadequate for mapping of anomalies not only in the continuous or discrete field models of stream sediment geochemical landscapes (Figs 5 and 6; see results of spatial association analysis further below) but also in the point stream sediment geochemical data sets (Fig. 2).

Mapping of anomalies in geochemical data sets containing multiple populations can be facilitated by application of, for example, probability plots (Sinclair 1974, 1991). However, the method of threshold estimation using probability graphs was not applied here because it assumes that the data under examination exhibit normal or lognormal distribution, but the original point data and derived model data sets in this study show neither normal nor strictly lognormal distributions (Fig. 7). Many workers have demonstrated that geochemical data invariably show neither normal nor lognormal distributions (e.g. Vistelius 1960; McGrath & Loveland 1992; Reimann & Filzmoser 2000), while several other workers have postulated that geochemical landscapes based on stream sediment geochemical data are multifractals (Bölviken *et al*. 1992; Cheng *et al*. 1996; Cheng 1999; Agterberg 2001; Rantitsch 2001; Li *et al*. 2002; Shen & Cohen 2005). Among these workers, Cheng *et al*. (1996) and Cheng (1999) modelled geochemical landscapes as continuous fields by contouring stream sediment geochemical data and then mapped geochemical anomalies by application of the C–A fractal method (Cheng *et al*. 1994). In contrast, Shen & Cohen (2005) modelled geochemical landscapes as discrete fields by attributing stream sediment geochemical data to their sample catchment basins and then mapped geochemical anomalies by application of a summation fractal method, which they proposed and which is a variant of the C–A fractal method. In another study, Li *et al*. (2003) demonstrated a concentration–distance (C–D) fractal model for direct analysis (i.e. without spatial interpolation) of point geochemical data, although they finally showed an interpolated map of geochemical data with anomalous classes/contours above the threshold derived via the method they proposed. Perhaps they, too, realized that point geochemical anomaly maps only vaguely suggest the spatial extents of anomalous areas so that they opted to portray their anomaly map as a continuous field model of geochemical data. In this study, the C–A fractal method (Cheng *et al*. 1994), which is far more tested than the C–D method (Li *et al*. 2003) and summation method (Shen & Cohen 2005), was adopted in mapping of anomalies in the continuous and discrete field models of stream sediment geochemical landscapes. It can be argued that the overall approach followed by Li *et al*. (2003) to first define a threshold directly from the point geochemical data and then apply this threshold to define anomalies in an interpolated geochemical map is analogous (i.e. ‘reverse’) to one of the approaches discussed here to map anomalies in continuous field models of stream sediment geochemical landscapes via the C–A method (Cheng *et al*. 1994) using classes (or contours) defined from analysis of empirical frequency distributions of point stream sediment geochemical data (see further below).

The concept of C–A fractal method (Cheng *et al*. 1994) can be summarized as follows. For a series of contours of concentration values in a geochemical landscape, concentration contours (** v**) and cumulative areas (

**) enclosed by each concentration contour (i.e.**

*A***(≥**

*A***)) are plotted on log-log graph, with**

*v**v*plotted along the x-axis and

**plotted along the y-axis. A C–A plot characterizes not only the empirical frequency distributions of concentration values but also the spatial distributions and geometrical properties of the features defined by concentration contours in a geochemical landscape. Concentration-area plots can be described invariably by power-law functions, which are depicted as straight lines (or line segments) on a log-log graph. If, on the one hand, a C–A plot can be depicted by one straight line, then it represents a geochemical landscape that can be described by a group of similarly-shaped concentration contours. In this case, the single straight line probably represents a fractal distribution of geochemical background. If, on the other hand, a C–A plot can be depicted by at least two straight-line segments, then it represents a geochemical landscape that can be described by at least two individual groups of similarly-shaped concentration contours. In this case, the rightmost straight-line segment (i.e. highest concentration values) probably represents a fractal distribution of geochemical anomalies, whereas the straight-line segment(s) to the left probably represent(s) a multifractal distribution (or intertwined fractal distributions) of geochemical background. Agterberg (2007) has shown that geochemical anomalies with Pareto distribution plot as a straight-line segment on the right of a C–A plot. Accordingly, the breaks in slopes of straight-line segments fitted through a log-log plot of the C–A relationship represent threshold values of different ranges or populations of concentration values in a geochemical landscape. These different populations would represent different background and anomalous geochemical processes, one of which could be mineralization.**

*A*In the C–A fractal analysis of anomalies in the continuous and discrete field models of stream sediment geochemical landscapes, two sets of experiments were performed by using concentration contours based on: (1) 5-percentile intervals of pixel values in the individual geochemical landscape models; and (2) 5-percentile intervals of values in the individual point geochemical data sets. The two sets of experiments were performed in order to obtain an answer for the second question raised in the introduction. The C–A plots obtained in the two sets of experiments are (due to lack of space) illustrated only for the As and As–Cu–Ni geochemical landscape models (Figs 8 and 9) but the threshold values defined for all the geochemical variables examined are listed in Tables 4 and 5.

The C–A plots obtained by using concentration contours based on 5-percentile intervals of pixel values in the continuous field models of the As and As–Cu–Ni geochemical landscapes (Fig. 8A, B) are strongly similar to the corresponding C–A plots obtained by using concentration contours based on 5-percentile intervals of pixel values in the discrete field models of the As and As–Cu–Ni geochemical landscapes (Fig. 8C, D). The same observations can be made for the results of the C–A fractal analyses of the models of stream sediment geochemical landscapes of the other elements. Accordingly, the sets of threshold values obtained from the continuous field models of stream sediment geochemical landscapes are strongly similar to the corresponding sets of threshold values obtained from the discrete field models of stream sediment geochemical landscapes (Table 4). It appears, however, that most of the threshold values obtained from the continuous field models of stream sediment geochemical landscapes (except for Co) are slightly higher than the threshold values obtained from the corresponding discrete field models of stream sediment geochemical landscapes. This is possibly due to the fact that pixel values in the continuous field models of As and As–Cu–Ni geochemical landscapes have somewhat stronger variances than pixel values in the corresponding discrete field models As and As–Cu–Ni geochemical landscapes (Fig. 7). Nevertheless, the results shown in Figure 8 and Table 4 imply that mapping of anomalies via C–A fractal analysis of both continuous and discrete field models of stream sediment geochemical landscapes are likely to yield strongly similar results.

When concentration contours in both continuous and discrete field models of stream sediment geochemical landscapes are defined according to 5-percentile intervals of values in the individual point stream sediment geochemical data sets, the C–A plots derived from the continuous field models of As and As–Cu–Ni geochemical landscapes (Fig. 9A, B) are also strongly similar to the corresponding C–A plots derived from the discrete field models of As and As–Cu–Ni geochemical landscapes (Fig. 9C, D). The same observations can be made for the results of the C–A fractal analyses of the models of stream sediment geochemical landscapes of the other elements. Accordingly, the sets of threshold values derived from the continuous field models of stream sediment geochemical landscapes are either strongly similar or identical to the corresponding sets of threshold values derived from the discrete field models of stream sediment geochemical landscapes (Table 5). It appears, however, that some of the threshold values derived from some of the continuous field models of stream sediment geochemical landscapes are slightly higher than the threshold values derived from the corresponding discrete field models of stream sediment geochemical landscapes. A plausible reason for this is that pixel values in the continuous field models of As and As–Cu–Ni geochemical landscapes have somewhat stronger variances than pixel values in the corresponding discrete field models As and As–Cu–Ni geochemical landscapes (Fig. 7). Nevertheless, the results shown in Figure 9 and Table 5 imply that mapping of anomalies via C–A fractal analysis of both continuous and discrete field models of stream sediment geochemical landscapes are likely to yield strongly similar results.

Thus, the results of the two sets of experiments (summarized in Tables 4 and 5) imply that, with regard to defining concentration contours in stream sediment geochemical landscapes for the application of the C–A fractal method to identify threshold values, concentration contours defined from empirical frequency distributions of the point stream sediment geochemical data are as useful as concentration contours defined from empirical frequency distributions of derived pixel values in corresponding continuous or discrete field models of stream sediment geochemical landscapes. Consequently, as shown in Figures 10 and 11, maps with strongly similar spatial distributions of background and anomalous populations of As and As–Cu–Ni values can be obtained by using concentration contours defined from empirical frequency distributions of either the measured values in the original point data sets or the derived pixel values in continuous/discrete field models of stream sediment geochemical landscapes. Because of this, further spatial analyses are required in order to answer both the first and second questions stated in the introduction.

### Analysis of spatial association of geochemical anomalies with mineral deposits

Several methods exist for quantifying spatial associations between patterns in two maps (e.g. geochemical anomalies and mineral deposit occurrences). Some of these methods (e.g. calculation of Jaccard's coefficient and Yule's coefficient) involve converting maps into binary patterns, while other methods (e.g. weights-of-evidence (WofE) modelling) are applicable to maps with either binary or multiple patterns (Bonham-Carter 1994). Here, the WofE method (Bonham-Carter *et al*. 1989; Agterberg *et al*. 1990) was applied to quantify spatial associations of the anomalous and background areas in the geochemical landscape models with the locations of known epithermal Au deposit occurrences in the study area. The binary anomaly maps resulting from the application of the mean+2*SDEV* method and the median+2*MAD* shown in Figures 5 and 6, respectively, were evaluated via the binary approach to WofE modelling. The multi-class anomaly maps resulting from the application of C–A method shown in Figures 10 and 11 were evaluated via the cumulative descending approach to WofE modelling (i.e. using maximum to minimum values (or anomalous to background populations) of data) because locations (or occurrences) of mineral deposits are themselves geochemical anomalies and thus stream sediment geochemical anomalies, if significant (i.e. related to mineral deposits), are expected to exhibit positive, rather than negative, spatial associations with mineral deposits. In WofE modelling using either binary or multi-class maps, the type and relative magnitude of spatial association between two map patterns (e.g. geochemical anomalies versus mineral deposits) can be characterized by calculation of a spatial statistic called contrast (*C*) (for details of estimation see Bonham-Carter *et al*. 1989, Agterberg *et al*. 1990). The value of *C* is greater than zero in case of positive spatial association and less than zero in case of negative spatial association. The statistical significance of spatial association can be determined by calculating Studentized *C*, which is the ratio of *C* to its standard deviation. Bonham-Carter (1994, p.323) stated that it is desirable to obtain statistically significant spatial evidence of mineral prospectivity having Studentized *C* values of ≥2 with respect to mineral deposit occurrences.

In WofE modelling with small number of mineral deposit occurrences (e.g. Carranza 2004*b*), quantified spatial associations are likely to be statistically non-significant. This is so in the present case study because there are only 12 epithermal Au deposit occurrences within the sample catchment basins (see, for example, Fig. 11) and particularly because the pixel size used in the analysis (i.e. 10 m) is very small so that the number of deposit pixels is very small compared to the total number of pixels within the sample catchment basins (i.e. 1007704). This problem was overcome by considering the four nearest neighbour pixels of a deposit pixel (i.e. to its north, east, south and west) to be deposit pixels also (cf. Stensgaard *et al*. 2006).

The results of the WofE analysis of the binary anomaly maps obtained from the application of the mean+2*SDEV* and median+2*MAD* methods are shown in Table 6. On the one hand, based on threshold values defined as mean+2*SDEV* of the As and As–Cu–Ni point data sets, anomalies mapped in either continuous or discrete field models of stream sediment geochemical landscapes do not coincide spatially with any of the known epithermal Au deposits in the district. On the other hand, based on threshold values defined as mean+2*MAD* of the As and As–Cu–Ni point data sets, anomalies mapped in either continuous or discrete field models of stream sediment geochemical landscapes coincide spatially with some of the known epithermal Au deposits in the district. These results imply, therefore, that the median+2*MAD* method, compared to the mean+2*SDEV* method, is more useful for mapping of anomalies based on stream sediment geochemical data. In addition, Table 6 shows that anomalies mapped in discrete field models, compared to anomalies mapped in continuous field models, of stream sediment geochemical landscapes, have stronger positive spatial associations with the known epithermal Au deposits in the district. These results imply, therefore, that discrete field models, compared to continuous field models, of stream sediment geochemical landscapes are more useful for mapping of anomalies.

The results of the WofE analysis of the multi-class anomaly maps resulting from the application of C–A method are given in Tables 7 and 8. On the one hand, Table 7 shows that if concentration contours in models of geochemical landscapes are based on 5-percentile intervals of their pixel values, high and/or low anomalies mapped in the discrete field models of As, Cu, and As–Cu–Ni geochemical landscapes have statistically significant positive spatial associations with the epithermal Au deposits, whereas high and low anomalies mapped in only the continuous field model of As geochemical landscape have statistically significant positive spatial associations with the epithermal Au deposits. In addition, Table 7 shows that if concentration contours in models of geochemical landscapes are based on 5-percentile intervals of their pixel values, high background mapped in the discrete field model of Ni geochemical landscape has statistically significant positive spatial association with the epithermal Au deposits, whereas high background mapped in the continuous field model of Ni geochemical landscape has statistically non-significant negative spatial association with the epithermal Au deposits. On the other hand, Table 8 shows that if concentration contours in models of geochemical landscapes are based on 5-percentile intervals of values in the original point geochemical data sets, high and/or low anomalies mapped in the discrete field models of As, Cu, and As–Cu–Ni geochemical landscapes have statistically significant positive spatial associations with the epithermal Au deposits, whereas high and low anomalies mapped in only the continuous field model of As geochemical landscape have statistically significant positive spatial associations with the epithermal Au deposits. In addition, Table 8 shows that if concentration contours in models of geochemical landscapes are based on 5-percentile intervals of values in the original point geochemical data sets, high background mapped in the discrete field model of Ni geochemical landscapes has statistically significant positive spatial association with the epithermal Au deposits, whereas high background mapped in the continuous field model of Ni geochemical landscape has statistically non-significant positive spatial association with the epithermal Au deposits. The results shown in Tables 7 and 8 imply that discrete field models, compared to continuous field models, of stream sediment geochemical landscapes are more useful in mapping of significant anomalies. The results shown in Tables 7 and 8 imply further that, with regard to defining concentration contours in stream sediment geochemical landscapes for the application of the C–A fractal method to identify threshold values, concentration contours defined from empirical frequency distributions of the point stream sediment geochemical data are as useful as concentration contours defined from empirical frequency distributions of derived pixel values in corresponding continuous or discrete field models of stream sediment geochemical landscapes.

## SUMMARY AND CONCLUSIONS

Modelling of stream sediment geochemical landscapes, in order to map significant anomalies, is not a trivial activity. Each method of doing so has its own advantages and disadvantages depending on, for example, the scale of mapping of geochemical anomalies. This study has focused on district-scale mapping of geochemical anomalies in continuous and discrete field models of stream sediment geochemical landscapes because point-symbol representation of stream sediment anomalies only vaguely suggests the spatial extents of anomalous areas or the plausible source catchments of anomalous stream sediment samples. In point anomaly maps, one can surely draw reasonable boundaries between background and anomalous stream sediment samples, but then this would be equivalent to either continuous or discrete field modelling. Continuous field modelling of geochemical landscapes via, for example, ordinary IDW interpolation of stream sediment geochemical data allows mapping of spatial extents of anomalous areas, but, like point-symbol representation of stream sediment geochemical data, it also vaguely suggests the plausible source catchments of anomalous stream sediment samples. In contrast, discrete field modelling of stream sediment geochemical landscapes, by attributing measured or derived geochemical variables to their sample catchment basins, allows mapping of spatial extents of anomalous areas and suggests more clearly the plausible source catchments of anomalous stream sediment samples.

In this case study, the C–A fractal method outperforms the mean+*2SDEV* and median+2*MAD* methods in separating background and anomalies in either continuous or discrete field models of stream sediment geochemical landscapes. That is because, as implied by the C–A plots (Figs 8 and 9), the pixel values in either continuous or discrete field models of stream sediment geochemical landscapes plausibly represent mixtures of at least three populations. It is usually preferable to apply classical statistical techniques, such as the mean+2*SDEV* method, or EDA statistical techniques, such as the median+2*MAD* method, to more-or-less homogeneous subsets of geochemical data in order to properly define threshold values. In contrast, the application of the C–A method does not require subdivision of geochemical data into more-or-less homogeneous subsets in order to properly define threshold values.

The results of the study support the hypothesis in the introduction that mapping of anomalies in continuous/discrete field models of geochemical landscapes can be based on analysis of empirical frequency distributions of measured point geochemical data from which those models were derived. The logic for doing so is that empirical frequency distributions of values in point geochemical data sets would be more trustworthy (i.e. we know more about the precision and accuracy of the measured point data sets) than the empirical frequency distributions of pixel values in the corresponding continuous/discrete field models of geochemical landscapes. However, the results of the study show that empirical frequency distributions of either the measured values in original point geochemical data sets or the pixels values in corresponding continuous/discrete field models of geochemical landscapes can be used as basis for mapping of anomalies via application of the C–A fractal method.

Finally, the results of the study demonstrate that, in district-scale stream sediment geochemical exploration surveys, catchment-based discrete field modelling of stream sediment geochemical landscapes is preferable to continuous field modelling of stream sediment geochemical landscapes. That is because the former respects the zone of influence of each stream sediment sample while the latter does not and because anomalies mapped in catchment-based discrete field models of stream sediment geochemical landscapes, compared to anomalies mapped in continuous field models of stream sediment geochemical landscapes, are likely to exhibit stronger positive spatial associations with mineral deposit occurrences. This finding in the present case study is consistent in both the results of statistical analyses of threshold values based on the original point data sets (Table 6) and the results of C–A fractal analyses of threshold values based on either continuous or discrete field models of geochemical landscapes (Tables 7 and 8).

## Acknowledgments

The author would like to thank Frits Agterberg and an anonymous reviewer for the valuable comments on an early version of the paper.

- © 2010 AAG/Geological Society of London