Multiple discriminant analysisof the DravaRiveralluvial plain sediments
Zoran Peh1 (corresponding author), Robert Šajn2, Josip Halamić1andLidija Galović1
Abstract:
Three discriminant function models are raised and cross-compared in order to distinguish geochemical patterns characteristic for the DravaRiverfloodplain sediments. Based on data representing total element concentrations in samples collected from alluvium (A), terrace (T),and unconsolidated bedrock (B) at the border of a floodplain, four element clusters emergedaccounting for discrimination between the referred groups of sediments. The most prominent is contaminant/carbonate cluster characteristic for alluvium. The other two are: silicate cluster typical for unconsolidated geological substrate (Neogene sedimentary rocks); and naturally dispersed heavy metal cluster separating terrace from the former two groups. Models introducing depth intervals and single profiles as grouping criteria reveal identical sediment-heavy metal matrices. The second important issue of this paper is possibility of reclassification of samples originally assigned to one of the a priori defined groups of sediments, based on established geochemical pattern. The mapped geological units can be reconsidered by the post hoc assignments to a different group if geological border between alluvium and terrace or between terrace and bedrock can not be established geologically with absolute certainty.
Key words: floodplain sediments, heavy metals, multiple discriminant analysis,Drava River, Slovenia, Croatia
1Z. Peh () J. Halamić L. Galović
Croatian Geological Institute,
Sachsova 2, P.O. Box 268, HR-10000 Zagreb, Croatia
E-mail:
Tel.:+385 (1) 6160-753
Fax: +385 (1) 6144-718
2R. Šajn
Geological Survey of Slovenia,
Dimičeva 14, 1000 Ljubljana, Slovenia
E-mail:
1.INTRODUCTION
The modern floodplains are potentialrepositories for a wide range of contaminants among which heavy metals represent by far the greatest threat to human health and environment. The release of metal-rich waste into the fluvial environmentoften combines various anthropogenic sources such as mining operations, industrial and urban wastes, transport, agricultural activity, deforestation, and others. This impact is recorded in the vertical overbank profiles of river alluvium, sometimes with dramatic increase of heavy metal concentration, mainly in the upper part (e.g. Macklin et al. 1994, 2006; Swennen et al. 1994; Swennen et al. 1998; SwennenandVanDer Sluys 2002; Langedall and Ottesen 1998; Ciszewski 2002; Bølviken et al. 2004;Pavlovic et al. 2004; Cappuyns et al. 2006),or,generally,in depth-integrated surficial overbank sediment samples or topsoil developed on floodplains (e.g. Ottesen et al. 1989; McConnel et al. 1993; Ódor et al. 1997; Rognerudet al. 2000; Bidovec et al. 1999; Halamić et al. 2003; Romicand Romic 2003). However, in longer terms, heavy metals may bealso introduced into the soil and hydrosphere by natural processes – weathering, erosion and transport – easily blending anthropogenic signal with a natural one and causing the pattern of metal distribution in alluvial plain sedimentsmore difficult to assess. The problem is magnified by the fact that floodplains join together various distinct geomorphic and hydrologic features which directly influence their geochemical signature. It has long been recognized that dispersion patterns, storage and remobilization of heavy metals in the fluvial system are directly associated with processes of sediment transport (e.g. Foster and Charlesworth 1996; Macklin et al. 2006). Floods regularly rejuvenate the lowest extensive land surfaces by lateral and vertical accretion(overbank deposition) of mixed recent material. Conversely,relatively flat surfaces, or terraces, lying abovethe modern floodplain may escape inundationand thus (largely) preserve the geochemical signature from the time of their origin. Sincetheseare no longer in direct relation to recent hydrological conditions(especially flood events), distinguishing between alluvium and terrace may be focused on difference between “current” and “historic” geochemical patterns which may be useful in mapping young terraces close in altitude to the modern floodplain but showing no visible scarp.This is particularly important since various studies show that, in general, variations in the heavy metal concentration across terraced floodplains are more strongly associated with terrace height than with terrace age (Lambert and Walling19781987: Brewer and Taylor 1997). By the same tokenterraces can be contrasted with bedrock if boundary between floodplain and upland isdifficult to trace due to unconsolidated material and poorly dissected reliefwhich indicate the latter.In this case different geochemical signals can be expected on account of diverse nature of the geological substrate:overbank sediment fillingthealluvial plain versus soil developed on unconsolidated bedrock.It must be also borne in mind that during dry times floodplain sediments, especially on higher terrace levels,are subject to various external influences including soil formation processes which can substantially alter element concentrations(Volden et al. 19961997; Rognerud et al. 2000).
The general aim of this research was to find geochemical contrast between the “recent” and “historic” sedimentary material (alluvium and terrace) deposited on alluvial plain of a great river during its long fluvial history, and both against the older, Upper Miocene/Pliocene weakly consolidated sedimentary rocks (bedrock) underlyingthe floodplain. In this light the presentedstudy can be viewed as a multivariate statisticsvariant of the geomorphological-geochemical approach to provenance analysis,based on the chemical composition of Holocene-Pleistocene fluvial sediments as the main carriers of the heavy metal load. The study draws on the premise that chemical (and physical) properties of the sediment can be used to infer its provenance although the applied method is based on indirect observation rather than direct comparison of floodplain sediments with individual potential sources as in the case of classical fingerprinting technique using multiple discriminant functions (e.g. Molinaroli and Basu 1993: Collins et al. 1997; 1998; Owens et al. 1999; Bottrill et al. 1999). This approach is deemed appropriatedue to different age of floodplain sediments and the lack of precise terrace chronology. Again, it differs from the line using direct comparison of a number of well defined predictors (diagnostic properties) with defined hierarchical relationship (e.g. lithology-mineralogy-geochemistry) such as in canonical correlation analysis where each element, mineral and lithologic unit represents a family of fingerprint properties (Čović 2003; unpublished MSc thesis).A multivariate statistical assessment of primary geochemical pattern distinguishing between a priori defined groups of lithological units (“lithology”) – alluvium, terrace and bedrock – was performedpreferring multiple discriminant analysis (MDA) as the most helpful mathematical tool. Additional discriminant models were constructed in order to examine the influence of depth (“depth-lithology”) and local variations (“profile”) on primary geochemical signature defined by distinctive variable clusters, especially the families of heavy metals.
2.DESCRIPTION OF THE STUDY AREA
The study area lies between 14˚55’ and 16˚55’ E longitude and 46˚129 13 and 46˚38’N latitude. It occupies the large portion of the DravaRiver drainage area spreading from lead-zinc mining and smelting region in northeastern Slovenia (Mežica) to the mostly agricultural lowlands (Drava alluvial plain) of northern Croatia (Fig. 1). The upstream, Slovenian, part is built predominantly of igneous and metamorphic rocks (Palaeozoic) with rare occurrences of sedimentary rocks (mostly Mesozoic) (e.g Mioč and Žnidarčič 1977) while downstream, towards Croatia, the terrain is geologically substantially different, being composed almost solely of sedimentary rocks of the Tertiary and Quaternary age (e.g Mioč and Marković 1998). The youngest, alluvial, sediments along the Drava River are notorious for their increased heavy metal composition which may be either naturally dispersed by Drava tributaries draining the Pb-Zn (Cu, As, Cd) mineral occurrences and deposits hosted in Triassic sedimentary rocks (mostly limestones and dolomites) and Palaeozoic igneous-metamorphic (Palaeozoic) complex, or may draw its origin from anthropogenic input. The latter case is of special concern for the study area because heavy metals are released into the Drava fluvial system for a long time being by erosion of the mine spoils, tailings and contaminated topsoil (Mežica) or in the forms of contaminated aqueous effluents from the metal processing facilities situated in the upstream reaches of the DravaRiver.Diffuse inputs of airborne metal particulates from the furnace stacks (e.g., Prevalje, Ravne, Kidričevo) can also contribute to overall pollution (Šajn 2002; 2003; Vreča et al. 2001; Šajn and Gosar, 2004).
3. FIELD AND ANALYTICAL PROCEDURES
3.1Sampling
During earlier investigations carried out for the purpose of the Basic Geochemical Maps of Slovenia and Croatia, anomalous values for Pb, Zn, As, and Cd were obtained from the chemical analysis of soil samples (topsoil A0-25) collected in the area encompassing Drava River valley and its tributary valleys (Andjelov 1993; Šajn et al. 2000, 2005; Miko et al. 2001; Halamić et al. 2003; 2005). In order to carry the research further for origin and distribution of anomalous elements six profiles lines in Slovenia (Dravograd – DR; Mežica – ME, Radlje – RA; Maribor – MB; Ptuj – PT and Zlatoličje – ZL) and three in Croatia (Sračinec – SR; Belica – BE and Legrad – LE) in a downstream sequence were sampled more or less perpendicularly to the river course before the confluence with the Mura River. The tenth profile line Ormož (OR),placed between ZL and SR was common for both Slovenian and Croatian part and was sampled separately by each research team (Fig. 1). Sampling was carried out according to recommendations of the IGCP and FOREGS task group for the geochemical mapping (Darnley et al. 1995; Salminen et al. 1998).
The total of 47 boreholes(individual profile points) of various depthswasdrilled along the selected profile lines providingas much as 333 samples for chemical analysis. Boreholes were located in a pre-designed scheme set across the river valley to includeeach terrace level on the both sides of the river bed and alluvial sediments in the middle. The peripheral boreholes were located on the bedrock. As the hand-operated drilling device could not go through gravel drilling was halted as soon as the first thick gravel layer was reached. Thesamples were taken at 20 cm intervals as a composite, while the maximum depth varied between 60 cm (PT) and 260cm (DR).
3.2Sample preparation and analysis
Samples were air-dried at temperature <400C for approximately three months to prevent the loss of volatile components and then sieved to the <0.125 mm fraction and finally homogenized in agate mortar. The <0.125 mm fraction was used following some earlier investigations which have shown that high content of most elements, particularly the trace elements, is present in the fine-grained fraction (e.g. Whitney 1975; Gibbs 1977; Förstner and Whitmann 1981;Foster and Charlesworth 1996; Sutherland 1998; Singh et al. 1999; Rhoads and Cahill 1999). Also, this fraction usually contributes to more than 95% of the particles in most samples (Swennen et al. 1998).
All analytical work was carried out at the ACME Analytical Laboratories in Vancouver(Canada). Samples were analyzedby inductively-coupled plasma mass spectrometry (ICP-MS) after total hot 4-acid (HCl-HNO3-HF-HClO4) digestion at temperature of 200C on the set of 41 elements. The digestion was only partial for some Cr and Ba minerals and some Al, Hf, Mn, Sn, Ta and Zr oxides. The volatilization during fuming may have also resulted in some loss of As, Sb and Au.Mercury was determined with cold vapour atomic absorption spectrometry (CV-AAS) after aqua regia digestion.
3.3Accuracy and precision
Accuracy was controlled by certified geological reference materials (DST6/DS6 – ACME Laboratory). For most elements analyzed in reference soil materials accuracy wasfound in the range of ±10% of the certified values, except for Ag (±25%). Precision was determined by repeated analysis of both certified reference samples and randomly selected soil samples (16 duplicated samples).The resulting coefficient of variation was, on average, approximately 5%.
4.DATA PROCESSING
4.1Univariate statistics and data transformation
A suite of 27 elements including 8 major and 19 trace elements were selected as predictor variables incomputing discriminant functions. Summary statistical results forthe whole dataset prior to the multivariate statistical procedure are given inTable 1 (minimum, maximum, median, mean, standard deviation and skewness).Taking into account that most variables are characterized bynon-normal, mostly positively skewed frequency distributions (except Al and Na which are negatively skewed),appropriate transformations were found necessary to get more symmetric distribution for most elements.As a general rule, most geochemical and environmental data do not follow normal (neither lognormal) distribution (e.g. Matschullat et al. 2000; Reimann and Filzmoser 2000; Reimann et al. 2005), which is the outcome of complex non-linear dynamics, feedbacks and thresholds resulting in outliers within investigated system (sample media such as, for example, soils, stream- or overbank sediments) represented by the particular set of variables (e.g. Hugget 1998; Phillips 1999). To stabilize their variance the usual remedy in environmental and geochemical exploration was applied in this work represented by conventional normalization procedure using log-, ln-, and sqrt-transformations,which included all elements except Al and K.Distributional characteristics of the input data were examined by the test of normality (Table 1). According to the applied Shapiro-Wilk’s(S-W) statistical test only two original variables exhibit normal distribution (p>0.05, marked by the asterisk). Even after the different transformation methods were used to approach the normal distribution, normalization procedure failed to improve normality ofa number of elements.In the case where applied transformation produced even greater skew than original (such as with V) the data were rather left untransformed.
4.2 The Strategy of Analysis
Multiple (multi-group) discrimination analysis (MDA) can be particularly usefulmultivariate method when applied to the data scattered within and across variousspatial boundaries such as geological (lithological) units,depth intervals,or profile lines, being composed of the same suite of observed attributes. Its primary purpose is to establish the major sources of between-group differences (Dillon and Goldstein 1984; Rock, 1988)often represented, in environmental studies, by accumulation of heavy metals in various geologic (sampling) media distributed horizontally and vertically from the focal point, or line of dispersion.
Data collected from a number of profile linesset across the investigated part of the Drava valley can be organized in several ways as the grouping criteria can follow different conceptions. The most obvious approach to grouping is based onthe underlying (mostly Quaternary) geology as well as onthe depth intervals from which the samples were collected, after which the two types of divisions (models) can be formed, namely LITHOLOGY and DEPTH. The first of such “a priori” groupings can be accepted as most “genuine”, so that basic differences between the groups can be easily inspected provided that the data variability is great enough to disclose them. It must be noted, however, that sampling procedure skipping the natural boundaries between the soil horizons may decrease the between- against within-group variability rendering the geochemical depth differences more difficult to understand. The “noise” can be partly reduced by combining the two divisions in such a way that newly formed groups more clearly indicate differences between the topand the deeper intervals, for exampleDEPTH-LITHOLOGY. This approach is introduced into analysis to dispel ambiguity between upper intervals sampled in areas further from the most possible source of pollution on the one side and deeper intervals nearer to water course exposed to recurrent flooding on the other.
However, with groups thus formed their mutual relationship can notbe represented meaningfully on relevant field maps since there is not enough mapping points to lay a grid. Also, the profile data are arranged in intervals rendering their representation, unless portrayed in a cross-section, largely abstract. To avoid this problem and, possibly, to focus on a narrower stand such as the single profile, the groups can be reshaped by combining singleprofiles(sampling points, or boreholes) withLITHOLOGY or DEPTH, or both. Thus, each profile pointcan be portrayed as a group of its own, its intervals representing objects (samples). By this rearrangement, the new groups can be created (PROFILE) following the adjusted grouping criterionand represented on a map via the relevant profile locations. However, being arbitrarily created, these groups can introduce a sort of (additional) artificiality into original data. In such a case interpretation must be exercised with great care since differences between larger groups, following the criteria ofLITHOLOGY or DEPTH, can be easily offset by the variations within the smaller ones, localized inindividual profiles.
More than oneMDA can be independently performed following the modes of data grouping discussed above. First and foremost, the data are arranged in such a way that sources of variation between thegeochemistry of underlyinglithological units (LITHOLOGY)– alluvium (A), terrace (T) and bedrock (B) – can be examined (3 groups). The second grouping criterion deals with differences between combination of lithological units and three depth intervals (top, middle and bottom) at which the samples were collected (DEPTH-LITHOLOGY) – including 9 groups altogether which are composed of 20 cm intervals (0-20, 20-40, and 40 to the bottom). The third model is a depth-profile combination which is extended to all three lithological units and contains 47 groups ofprofiles (PROFILE). This last approach was designed to highlight the possible variations within the singleprofiles. The accuracy of the group affiliation can be later easily examined by comparison of the “a priori” (observed) and mathematically computed (post hoc) classification for each object.In all three casesthe same set of objects (N=333) and variables (p=27) is utilized.
5. RESULTS AND DISCUSSION
The results of the MDA are briefly summarizedin the joint table (STATISTICA, Version 6)(Table3) comprising all three of the a priori built exploratory models (cases). Before that the overall significance of their discrimination is tested by theappropriate multivariate tests (Table 2)divulging the vanishingly low associated probabilities (p<0.000), a prerequisite to safely proceed with computing discriminant functions (DF).According to different number of groups (K) in various models the total number of DFsis K-1 or, as with the PROFILE model where the number of groups (K=47) is greater than number variables (p=27), it can not exceed the latter (p-1). Irrespective of statistical significance of the variation between the observed groups (p-level in Table2) which defines dimensionality of the discriminant space, the smaller number of DFs is used to explain the natural variation between groups. The natural variation hidden behind the original data actually served as an effective parsimony criterion reducing the number of discriminant axes to only two or three. Note, however, that in the LITHOLOGY modelonly two DFs exist (K-1) which completelyexplains the difference between the three observed groups.