Influence of site-classification approach on geochemical background values

2.1 Sampling and information about sites in the study area

Soil sampling sites were located in the inner and outer peri-urban parts of Vilnius in various directions from the centre (Figure 1).

Figure 1

Location of soil sampling sites in peri-urban part of Vilnius and information about the underlying Quaternary deposits, relief elevation and land use. Information on Quaternary deposits was taken from the Lithuanian Geological Survey [17]. The index in the upper row indicates their origin and age. Origin: lg – proglacial glaciolacustrine, lgt – marginal glaciolacustrine, lg eg – englacial glaciolacustrine, g – basal till, gt – marginal till, f – proglacial glaciofluvial, ft – marginal glaciofluvial, a – alluvial, v – aeolian. Age: II md – Middle Pleistocene Medininkai, III gr – Late Pleistocene Grūda, III bl – Late Pleistocene Baltija and IV – Holocene. The lower row contains site ID, lithology of deposits and relief elevation (m). Lithology: cla – clay; tll – till; different types of sand, i.e. ssc – silty-clayey, scl – clayey, sfi – fine, sme – medium, sva – various and sgr – gravelly. The lower row is in boldface when no information about soil type and lithology (texture) is available in the soil database [18].

Although the sites were not in built-up areas, but rather in forests or meadows, diffuse source influence is inevitable as the sites are not far from roads or settlements. The total number of sites was 25, each of which included 5 sub-sites resulting in 125 samples. The minerogenic soil sample in each sub-site was composite, i.e. formed from five sub-samples taken from topsoil to a depth of 25–30 cm excluding the forest litter. The sub-samples were mixed in a plastic bucket in the field and a portion was taken to the laboratory. Selection of the sampling sites was based on the underlying Quaternary deposits, taking into account their influence on soil texture and elemental composition. Later, using median values of LKS-94 coordinates measured by GPS in five sub-sites, information on these deposits was adjusted based on digital maps from the Geological Survey of Lithuania: (i) for outer peri-urban sites, on the Vilnius Quaternary geological map (scale 1:50,000) [17] and (ii) for inner peri-urban sites, on the Quaternary map of Vilnius (scale 1:16,000) compiled by R. Guobytė.

Information on soil taken from the Geoportal database [18] appeared to be incomplete, i.e. it was absent in the seven sites indicated in boldface type in Figure 1. In other sites, soil types and the respective first-level WRB soil reference groups [19] are presented in Table 1.

In these sites, there is also information on texture (lithological) classes for surface soil and for soil parent material (the latter is in Table 2), which is based on two classifications [20]: (i) the old one (used before 2000) by Kachinskii with 12 classes according to fraction <0.01 mm and (ii) the new one by FAO and modified in Lithuania in 1998 with 15 classes according to the percentages of sand (0.05–2 mm), silt (0.002–0.05 mm) and clay (<0.002 mm) fractions.

2.2 Sample preparation, determination of soil organic matter, texture and total element contents

All samples were sieved through a nylon sieve to <2 mm fractions, which the authors used for analysis: (i) estimation of soil organic matter percentage by the loss-on-ignition (LOI) method with combustion at 550°C [21]; (ii) grain-size measurement by the Fritsch Laser Particle Sizer Analysette 22 MicroTec Plus [22] with calculation of percentage of sand (63–2,000 µm) (SAF), silt (2–63 µm) (SIF) and clay-sized (<2 µm) (CLF) fractions [19]; and (iii) determination of total contents of chemical elements by energy-dispersive X-ray fluorescence Spectro Xepos equipment which uses Turboquant for the pressed pellets calibration method [23].

Additional sample preparation before chemical analysis was the following: (i) part of a <2-mm fraction was milled using an MM 400 mixer mill to obtain two 4 g portions of powdery material, each of which was homogenised with 0.9 g of Licowax binder and (ii) two parallel pressed pellets of 32 mm diameter were prepared from these mixtures using a PP15 press. Total contents of Al, Ca, Fe, K, Mg, Na, P, S, Si, and Ti (ten major elements) and Ba, Cr, Cu, Mn, Ni, Pb, and Zn (seven PHEs) were used in this study. The joint uncertainty of their contents due to sample preparation and measurement was variable in the different samples, but the arrangement of the major elements according to the median value of the relative standard deviation in the paired pellets (%) divided by the square root of 2 (Med_RSD) was as follows: Fe(0.9) < K, Si(1.0) < Ca, Al(1.2) < Ti(1.5) < P(1.7) < Mg(2.1) < Na(2.2) < S(3.6). The sequence of PHEs according to Med_RSD was the following: Mn(1.0) < Zn(1.1) < Ba(2.0) < Pb(2.8) < Ni(3.4) < Cr(7.3) < Cu(8.5). The contents of the elements were adjusted using the re-calibration curves compiled by analysing samples of the “International soil-analytical exchange” program organised by Wageningen University.

2.3 Statistical data treatment

MB values were estimated by the Med ± 2MAD method in Microsoft Excel. STATISTICA 9 software was used to calculate the values of Spearman’s rank correlation coefficients (r_S) and their significance, to perform hierarchical cluster analysis by Ward’s method using 1 − r_S distance (for variables) or Euclidean distance (for sites), to test the non-parametric hypotheses using the Mann–Whitney U test and to perform Kruskal–Wallis analysis of variance and multiple comparisons of mean ranks for several groups (their number can be up to ten). Computation of post hoc comparisons of mean ranks of all pairs of groups (explained in STATISTICA 9 help) is based on the publication by Siegel and Castellan [24]. First, the z values for each comparison between groups u and v are computed according to the formula:

where R with subscripts are the average ranks for groups u and v, n_u and n_v are the numbers of observations in them and N is the total number of observations across all groups. Then p values (two-sided significance levels with a Bonferroni adjustment) associated with each comparison are computed as follows: p=pz⋅k⋅k−1, where pz=2⋅1−PrZ<z and k is the total number of groups compared.

Ethical approval: The conducted research is not related to either human or animal use.

3.1 Three site-classification approaches using available qualitative information

We selected three different natural factors (approaches) which due to their presumed indirect relationship with the increased clay-sized or fine fractions in topsoil can be used for classification of sites into LGs: P approach, i.e. according to parent material (underlying Quaternary deposits); S approach, i.e. according to soil type; and L approach, i.e. according to lithology (texture) of the soil parent material.

Four LGs for each site classification seemed reasonable, because to classify using S approach, there are three groups of soil types (Table 1) and one unidentified group of sites (bolded in Figure 1). Not only Arenosols (Table 1) but also unidentified soils are on sandy Quaternary deposits. We presumed that the latter group has the lowest level (LG1) and the code S0 was assigned to it (zero emphasises uncertainty), and attributed Arenosols to the group with the lowered level (LG2, code S2). Since Albeluvisols (Retisols) contain more sand and less clay in the upper horizons (A and E) than Luvisols [25], they were classified to the group with the elevated level (LG3, code S3), while the Luvisols were classified to the group with the highest level (LG4, code S4).

As the surface soil texture in sites appeared to be very uniform (mainly sandy loam), the texture of the soil parent material was used to classify by L approach (Table 2). Similar to the classification according to the S approach, the unidentified group sites were attributed to LG1 (code L0).

Great variability in Quaternary deposits near Vilnius is obvious in the sampling sites (Figure 1): there are 4 ages, 9 origin types and 15 combinations as well as 8 different lithologies. To classify sites into four LGs according to the P approach, we took into account mainly the generalised origin and lithology (Table 3): (i) glaciolacustrine (P4 or P2 depending on lithology), (ii) glacial (P3), and (iii) other, i.e. glaciofluvial, alluvial, aeolian (the latter are represented by sand, so more detailed information about origin and age was used to subdivide them, presuming that alluvial or aeolian deposits, also relatively younger proglacial glaciofluvial, should be attributed to LG1 with code P1, while relatively older marginal glaciofluvial deposits were attributed to LG2 with code P2).

3.2 Two site-classification approaches using quantitative geochemical and textural data

To distinguish LGs by geochemical approach (G approach), three major elements, Al, K and Ti, were used as variable classifiers based on our previous research [15,26] and presuming that (i) the contents of PHEs should be correlated with them and (ii) they should be correlated with the CLF and possibly fine fraction (<63 µm; [FIF]) percentage. If the closure effect of geochemical data is not taken into account, seven significantly (p < 0.01) inter-correlated major elements (ME₇) are observed, i.e. Al, K, Ti, Fe, Mg, Ca and S, but the association of Al, K, and Ti is the closest (Table 4 and Figure 2). The content of each PHE has a significant (p < 0.01) correlation with all ME₇ members, but the association with Al, K, Ti is relatively greater, because (i) the contents of Cr, Ba and Mn are most of all correlated namely with Al, K, Ti; (ii) the contents of Zn, Ni, Cu, Pb have relatively high correlation with at least one of these three major elements (Table 4).

Figure 2

Cluster analysis dendrogram indicating association of Al, K and Ti.

Ignoring the closure effect of geochemical and grain-size data, the contents of all ME₇ members have a significant (p < 0.01) correlation with both CLF and silt fraction (SIF) as well as with FIF and LOI percentages (Table 4). Almost all ME₇ members except S have a higher correlation with CLF or SIF than with LOI; further, the correlation with SIF is higher than that with CLF. Our selection of only three variable classifiers was based on their grouping in the dendrogram (Figure 2).

Hence, the site classification into four LGs by G approach was based on the clustering way of Ward’s method and the Euclidean distances between the relative median contents of Al, K and Ti (division was by overall median; Figure 3). Ward’s method with city block distances, which was ranked higher [27], provided identical classification.

Figure 3

Four LGs of sites distinguished using the geochemical approach (G approach) for site classification and their comparison with those obtained by the other four site-classification approaches: T (textural), L (lithological), S (soil type), P (parent material). Indices after the letters indicate the presumed level of LGs: 1 or 0 – lowest (LG1), 2 – lowered (LG2), 3 – elevated (LG3) and 4 – highest (LG4).

The significant correlation of all PHEs with both CLF and SIF was the basis for finding out by the same clustering way the LGs of the textural approach (T approach), according to the mean values of CLF and SIF (Figure 4).

Figure 4

Four LGs of sites distinguished using T approach for site classification. The number after the first dash indicates the average FIF (%); after the second dash it indicates the median LOI (%) in the site. The codes after the third dash indicate the WRB textural class [19] according to average percentages of SAF, SIF and CLF determined by laser diffraction in sites: LS – loamy sand (number of sites n = 3), SL – sandy loam (n = 9) and SiL – silt loam (n = 13). The code after the fourth dash indicates WRB textural class [19] corresponding to that in the soil database for the soil parent material: ?? – no information (n = 7, L0), S – sand (n = 5, L2), SL – sandy loam (n = 7, L3), L – loam (n = 5, L4) and CL – clay loam (n = 1, L4). The code after the fifth dash indicates the lithology of the Quaternary deposits (Figure 1).

3.3 Particle size fractions and LOI in LGs distinguished by different site-classification approaches

Four LGs distinguished according to each site-classification approach have analogous similarities according to the parameters of particle-size and LOI (Supplementary Table S1). The most important of these is the increase in CLF, SIF, FIF and LOI as the index of LGs grows (Figure 5).

Figure 5

Increase in average percentages of influential grain-size variables (a, b and d) and soil organic matter (c) in four LGs distinguished by five site-classification approaches (see Figure 3). Error bars are based on standard deviation.

All LG2 groups are sufficiently homogeneous, i.e. with relatively low coefficients of variation (CV < 33.3%) of SIF, CLF and LOI, but LGs of other levels have exceptions (Supplementary Table S1). The G and T approaches provide relatively lower CV values in LG1 groups and enable the greatest number of homogeneous LGs to be distinguished: only their LG4 groups are heterogeneous according to CLF. But other approaches also distinguish heterogeneous LG4 groups, because site 20 in them is anomalous. It is the only one with a clay loam texture and clay lithology of underlying deposits (Table 2); and it has much higher CLF in comparison with other sites attributed to LG4 groups. The U test and multiple comparisons test applied to grain-size variables showed different ability of site-classification approaches to separate adjacent LGs (Supplementary Table S2).

3.4 Median background values: non-differentiated and differentiated by each approach

The NMB values of ME₇ and PHEs divide the range of DMB values into two intervals (Tables 5 and 6): those of LG1 and LG2 belong to the lower interval and those of LG3 and LG4 belong to the upper one (with only a few exceptions). In contrast, the DMB values of Si are higher in LG1 and LG2 but lower in LG3 and LG4.

The DMB values of ME₇ and PHEs usually increase the higher the index of LGs. But only G approach has clear trends with the growth of the index of LGs: (i) increase of DMB values of all ME₇ and decrease of Si; (ii) increase of DMB values of all PHEs. Other approaches have exceptions, especially the S approach.

The DMB values of PHEs taken from the study of Zhao et al. [2] or from “Geochemical atlas of Lithuania” [29] (Table 6) represent different soil textures, so one may expect them to be somewhat similar to the DMB values obtained by L approach. Although their direct comparison is not substantiated due to different analytical methods and different numbers of LGs, some interesting details can be seen. All the DMD values of PHEs in sandy topsoil taken from “Geochemical atlas of Lithuania” [29], which are based on total contents like in our research, are higher than the DMB values in the L0 and L2 groups, which are also sandy. This is likely due to a different method, i.e. optical emission spectrophotometry, or different qualitative site classification. Although the DMB values in sandy topsoil taken from Zhao et al. [2] are based on partial aqua regia contents and should be lower than the total contents, the results in Table 6 contradict this. For Cu, Zn and Pb, they appear to be higher not only than in the L0 and L2 groups but also than in “Geochemical atlas” [29]; for Ni, the DMB value exceeds the respective values in the L0 and L2 groups; for Cr, it exceeds the value in L0. Possible reasons are differences in the underlying deposits and higher DSI in various types of land use in England and Wales.

3.5 Geochemical differences between adjacent LGs and homogeneity of LGs

Three findings concerning the differences between adjacent LGs in the contents of ME₇ and PHEs (after elimination of their anomalies) and in LOI (Supplementary Table S3) are as follows: (i) all approaches reveal many differences between LG3 and LG2; (ii) most approaches reveal differences between LG2 and LG1; and (iii) G approach is most effective to reveal differences between LG4 and LG3. Further, for Al, K and Ti, the G approach results in a minimum sum of robust coefficients of variation (RCV) across LGs (Supplementary Table S4).

3.6 Scoring and ranking of results of five site-classification approaches

The scoring and ranking were done by taking into account the following ten differentiating variables (Dvar): SIF, CLF, LOI and the contents of ME₇. Three scores were calculated based on the features of LGs obtained by each classification approach (Supplementary Table S5): (i) trend score reflecting the consecutive increase in the central values of Dvar in LGs with a growing index; (ii) separation score indicating the overall degree of significant differences in the values of Dvar between adjacent LGs; and (iii) heterogeneity score estimating the overall variability of Dvar within LGs. Three analogous scores were calculated in aforementioned way for seven PHEs.

Trend score is the number of Dvar with consecutive increase in average values (for SIF, CLF and LOI) or DMB values (for ME₇) in LGs with the growth of their indices. Its maximum possible value is 10 (number of Dvar). The arrangement of approaches according to decreasing trend scores G&P > T&L > S can be adjusted by adding the respective scores of PHEs; this provides the sequence G > P > T > L > S.

Separation score takes into account not only the number of Dvar (10) but also three pairs of adjacent LGs and three possible values of STR code-variable (statistical test results obtained by comparing adjacent LGs): STR = 2, if the difference is significant (p < 0.05) according to the multiple comparisons test; STR = 1, if it is insignificant (p ≥ 0.05) according to the multiple comparisons test but significant (p < 0.05) according to the U test; STR = 0, if the differences are insignificant (p ≥ 0.05) according to the U test. The maximum possible value for each Dvar is 2 × 3 = 6, so the separation score of the approach is <60. The arrangement of approaches according to decreasing separation scores is G > T > P > L > S and is rather similar to that according to PHEs: T > G > L > S > P.

The heterogeneity score is the sum of the CV values of CLF, SIF and LOI in four LGs (Supplementary Table S1) plus the sum of the RCV values of ME₇ in all LGs after elimination of anomalies (Supplementary Table S4). The arrangement of approaches according to increasing heterogeneity score, i.e. T < G < L < S < P, as well as respective scores according to PHEs, shows that quantitative approaches provide greater overall homogeneity of LGs than qualitative and that P approach based on Quaternary maps gives the most heterogeneous LGs.

A higher trend rank or separation rank was given to the approach with a higher trend score or separation score, while a higher homogeneity rank was attributed to the approach with a lower heterogeneity score (Table 7). The sequence of approaches according to descending average ranks G > T > P > L > S indicates that the S and L approaches can be treated as low-exactness and P approach as medium-exactness, meanwhile the G and T approaches are obviously high-exactness.

3.7 Similarities in site classifications

When different classification approaches are used, soil from the same site is not always attributed to groups of the same level, and the number of sites in them usually differs (Figure 3). But some LGs distinguished using different approaches (S0 and L0, S2 and L2 and G1 and T1) are identical according to the sites; therefore, the respective DMB values are also equal (Tables 5 and 6). Significant (p < 0.01) r_S values between the indices of LGs in sites show sufficient similarity in different classifications. According to descending r_S values, the pairs of approaches can be arranged as follows: G&T (0.94) > T&L (0.86) > G&L, L&S (0.85) > G&P (0.81) > G&S (0.80) > T&S (0.78) > T&P, L&P (0.75) > S&P (0.73).

3.8 Similarity in background ranges and grain-size variables in LGs of the same level

The similarity in background ranges of ME₇ and PHEs in LGs of the same level obtained using five site-classification approaches was tested by the multiple comparisons test (Supplementary Table S6). In almost all cases where significant (p < 0.05) differences were revealed, a small (sometimes even significant) increase was observed in at least one of the influential grain-size variables, i.e. clay, silt or fine fractions. The arrangement of LGs according to the increasing number of significant differences in the contents of each of 14 elements is LG2 < LG1 < LG3 < LG4 (Supplementary Table S7). The analogous arrangement is according to PHEs, but slightly different according to ME₇: LG2 < LG3 < LG1 < LG4. For PHEs, the number of significant differences in all LGs is higher in comparison to that of ME₇.

3.9 Influence of MB values on ASI_M values

When different MB values are used, the number of sub-sites with elevated (>1.3) values of CI or EF calculated using Al as a reference element is rather high and any consistent pattern is hardly visible, most probably due to accidental waste. But the ASI_M values in five sub-sites, i.e. the median CI (CI_M) or the median EF (EF_M), reduce accidental influence and more properly characterise sites (Table 8).

Since the study sites were selected to determine background, those with at least one PHE having ASI_M > 1.3 will be called problematic (CI-problematic or EF-problematic), because the reason for the higher values is unknown.

According to the number of CI-problematic sites when using NMB values, these sites occur more often for Ni, Cr, Zn and Cu than for Mn, Ba or Pb.

In comparison with the NMB values, the DMB values of Ni, Cr, Cu and Zn obtained by all approaches reduce the number of CI-problematic sites (except Cu using L approach). But the G approach gives the lowest number of CI-problematic sites, and T approach reveals a lower number of CI-problematic sites of Ni and Zn than the L, S or P approaches. There are fewer EF-problematic sites than CI-problematic sites, even when normalisation by Al is done using NMB values; this is especially obvious for Zn. Analogous reduction is also when the DMB values estimated by various approaches are used (for Ni using G approach the numbers are equal). Each approach for background differentiation of Ni, Cr, Cu and Zn is useful in comparison with NMB values, because it results in a lower number of EF-problematic sites. However, the advantage of G approach is that for Ni and Cr, the number of EF-problematic sites is the lowest.

As concerns the other three PHEs (Mn, Ba and Pb), not all approaches ensure that using respective DMB values, the number of their CI-problematic sites will be lower or at least equal to those obtained using NMB values. But the G approach seems to be the most suitable for their reduction or at least stabilisation.

4.1 Why were such aims and tasks chosen?

The regional geochemical variability in soil metal contents reflecting both natural and anthropogenic sources was demonstrated in e.g. ref. [16]. Due to the inevitable influence of natural factors, including parent material and complexity of Quaternary deposits near Vilnius described by R. Guobytė [30], we presumed that this area also has high local geochemical variability and therefore background differentiation is necessary. According to Wu et al. [31], it is better to use methods with prediction for on-site contamination assessment, because methods that estimate TB values are unsuitable. Since TB values have drawbacks [15,32], we used MB values but excluded methods based on prediction, because variable predictors in the study by Wu et al. [31] are attributed to different factors; further, this or other algorithms, e.g. ref. [33,34], are complicated. The word “approach” in our research means background differentiation factor but not background determination method, which, according to Matschullat et al. [10], can be either empirical or theoretical. Our research can be attributed to the integrated method that was acknowledged later [35] when samples are taken in relatively pristine (i.e. peri-urban) areas, but analytical results are subjected to statistical calculations (i.e. Med ± 2MAD method). Although even in strongly urbanised areas, the influence of geological background has been found [36], and Johnson and Demetriades [37] emphasise the role of parent material and soil types in urban geochemical mapping, the determination of B values in such areas is complicated there due to the extremely variable soil resulting from many anthropogenic factors such as site age, human impact, land use, building rubble and municipal waste [38]. Urban topsoil in Klaipėda [15] also demonstrated that eastern and western parts of the city differ in underlying Quaternary deposits and topsoil Al, K and Ti contents, but due to anthropogenic influence on the sandy part, it was impossible to predict the VB values of some PHEs or correctly assess their DMB values by the Med ± 2MAD method.

Peri-urban soil in Vilnius underlain by Quaternary glacigenic or post-glacial deposits is suitable for experiments with DMB values obtained using different site-classification approaches. It does not matter that some sites are within the city boundaries, because urban areas can also be understood as built-up [37]. The main thing is that the sites are far from the centre and industrial districts with many pollution sources [39,40]. The uppermost layer of the Quaternary cover differs not only according to origin but also according to age, i.e. it is at the boundary of two glaciations: the last one being Nemunas (with older Grūda stage and younger Baltija stage) and the penultimate one being Medininkai [41]. Hence, representation of all ages and origin types needs many more samples.

We excluded the approach based on soil texture and organic matter used in ref. [2,29], because it is based on different factors, i.e. soil with high organic matter content is not a textural class [19]. We studied only minerogenic soil, since peaty or organic soil needs separate research, as it may be greatly influenced by anthropogenic pollution, as shown in the study of peri-urban ombrotrophic bog records [42].

4.2 Insufficiency of NMB values and explanation of approach ranking

In comparison with the median elemental contents of agricultural topsoil from ten European countries [28], the study areas are characterised by 1.2 times higher NMB value of Si and therefore lower NMB values of almost all (except K) other major elements (Table 5) and all PHEs (Table 6). A comparison with agricultural soil from Lithuania [28] also shows relatively lower NMB values of most elements related to clay or fine fractions. The sequence of elements according to descending ratios of medians S(2.2) > Ca(2.1) > Cu(1.9) > Cr(1.7) > Ni, Fe, Mg (1.6) > Zn(1.4) > P, Al (1.3) > Ti(1.2) > Ba(1.1) indicates the prevalence of light texture in the topsoil. But these NMB values are obviously unable to reveal increased DSI in sites that have low natural background level, as can be seen in the Ni example. Since its NMB value is in the middle of the range of DMB values (Table 6), all CI > 1.3 values of Ni appear to be in sites from LG3 or LG4, while low CI values are found in sites from LG1 or LG2 (Figure 6). This figure clearly demonstrates the advantage of differentiated background values in comparison with non-differentiated.

Figure 6

Values of contamination index (a) and enrichment factor (b) of Ni in sub-sites and sites using non-differentiated median background values and differentiated median background values obtained by G approach for classification of sites (LGs are in parentheses).

Although L and S approaches provide greater overall homogeneity of LGs of differentiating variables than P approach, they have the lowest trend and separation scores, especially S approach. Its lowest average rank is in agreement with the finding of Zhao et al. [2] that the major soil taxonomical group explains the lower percentage of variation in the contents of PHEs in comparison with soil texture.

The rather low average rank of L approach in comparison with T approach is explained not only by incomplete information from the soil database but also by the horizontal and vertical variability of the soil texture. Generalised information on soil texture in maps can mask its horizontal variability. The vertical variability of texture is obvious from the soil database: in topsoil, it is very uniform; while in soil parent material, it is more variable. In only 22% of study sites, the topsoil texture is similar to that of the parent material; while in about 61% of sites, the parent material is heavier. The vertical variability in soil texture in Luvisols from Northern Poland due to eluviation–illuviation (lessivage) processes was demonstrated by Switoniak et al. [43]. These researchers also found a rather uniform texture of the uppermost layer, and in most cases, a heavier texture in deeper horizons. Leaching of clay particles to illuvial horizons of Luvisols is without chemical degradation [44]. The vertical variability of soil texture is also obvious from discrepancies between the two texture classes, i.e. in topsoil determined by laser diffraction and in soil parent material indicated in the soil database (Figure 4 and Supplementary Table S8). Of course, the discrepancies can also be explained by systematic deviation in laser diffraction data in comparison with traditional pipette method data. The underestimation of clay and sand fractions in favour of silt fraction when using the laser diffraction method was demonstrated by Kun et al. [45]. Therefore, using this method, we obtained silt loam instead of clay or loam as indicated in the soil database and in seven sites where the soil texture class was absent from the soil database, but where sand was expected according to Quaternary deposits, it appeared to be loamy sand, sandy loam or silt loam. Despite this, laser diffraction and the pipette methods provided a similar estimation of sand, silt and clay fractions, since it was confirmed by significant inter-correlation of the respective results [45].

The higher average rank of P approach in comparison with low-exactness site-classification approaches is due to its highest trend rank and higher separation rank (Table 7). But it is close to the average rank of L approach, indicating rather good correspondence between the textural classes given for soil parent material and the lithology of the Quaternary deposits (Figure 4). Presumably, grain-size differences between soil parent material and underlying Quaternary deposits are lower than that between topsoil and soil parent material. Effective control of soil texture by underlying rocks has been demonstrated by Camara et al. [12]. The high heterogeneity score of P approach is most probably due to the existing differences between the topsoil and underlying Quaternary deposits as well as difficulties in classifying them into a few LGs.

An increasing trend seems to be an optional requirement for determination of DMB values, especially in countries with metalliferous mineralisation when classification is done using mixed approaches which include parent material as one of the factors [1]. But in Lithuania, this requirement is natural and useful, since a thick cover of glacigenic Quaternary deposits predetermines the decisive role of soil clay and fine fractions and a regular increase in the contents of most chemical elements. But the increasing trend and the separation of DMB values are not always clearly expressed even using textural approach, e.g. in six LGs of minerogenic soil presented by Zhao et al. [2] (Table 6). The following reasons are possible: (i) too many LGs distinguished during field assessment with small textural differences between some of them and (ii) the use of the aqua regia digestion. The increasing trend of DMB values in LGs distinguished according to minerogenic topsoil texture is better revealed for total contents of elements when the number of LGs is lower. For example, this can be seen in “Geochemical atlas of Lithuania” [29] for Ni, Cr, Cu, Zn, Mn, Ba (Table 6), as well as for Ti; but in comparison with the present research, the separation of adjacent LGs estimated by ratios of respective DMB values is lower. A possible reason is that information on soil texture in atlas was qualitative, i.e. taken from soil maps.

The advantage of T approach compared to L approach is the direct estimation of the percentage of topsoil fractions in samples by laser diffraction. Texture classes in topsoil samples also appeared to be rather uniform and unsuitable for classification into LGs, but clustering of continuous data enabled a higher ranked classification. The average rank of T approach is closer to G approach compared to P approach, because both G and T approaches are based on the topsoil measurements. The higher average rank of G approach compared to T approach is in agreement with the opinion of Dung et al. [46] that normalisation (the aim of which is similar to site classification) by geochemical methods is better than granulometric methods, because similar-in-texture samples can be different according to their mineralogical composition. The drawback of G approach is the difficulty in interpretation of LGs without other information. Qualitative information should never be ignored, even if it is inexact or incomplete. But to interpret the differences between G2 and G1, the grain-size and LOI data were necessary as they show a higher percentage of fine fraction and organic matter in G2 (see Figure 3 and labels in Figure 4).

4.3 Why was it useful to study similarities of site classifications?

The similarities of site classifications by using different approaches indicate that, in general, all of them are acceptable for background differentiation. But the high-exactness approaches (G and T) are preferable as they have higher average ranks according to the essential peculiarities of the soil matrix. Besides, these approaches are most similar. This finding enables us to presume that both will provide the most similar DMB values. The greatest similarity of the high-exactness approaches occurs because both are based on measurements of the same topsoil samples, i.e. determination of Al, K, Ti or clay and silt fractions. This similarity is in accordance with the correlation of Al with the CLF shown by Dreher et al. [14]. It is quite natural taking into account the variables used by researchers for normalisation [46], i.e. fractions <2 and <63 µm for granulometric normalisation and Al, K and Ti for geochemical normalisation.

A rather high similarity of low-exactness classification approaches (L and S) occurs because both are based on the same information source and their LG1 and LG2 contain the same sites. The lowest similarity between classifications based on soil type and Quaternary maps occurs because they use totally different sources of information and because there is a lack of soil-type information in the seven sites. The greater similarity of G approach to other approaches (r_S ≥ 0.80) in comparison with T approach is in favour of G approach. L approach based on the lithology (texture) of the soil parent material is most similar to both quantitative approaches, but its lower trend rank and especially its separation rank indicate that this approach can be used with caution.

Greater similarity of classifications ensures a lower number of significant geochemical differences between LGs of the same level, especially according to ME₇ (Supplementary Tables S6 and S7). Most similar high-exactness classification approaches have no significant differences according to ME₇ in all LGs. Low-exactness approaches are also rather similar according to ME₇ in three LGs but differ in LG4 according to Ti. However, if different exactness approaches are compared, the number of significant differences increases. For ME₇, all significant differences are observed when comparing LGs distinguished by different exactness approaches. In these cases, many more significant differences are also observed in the DMB values of PHEs. Consequently, different exactness site-classification approaches will provide many more different values of assessment indices.

4.4 Should DMB values be adjusted by additional normalisation?

The finding that most PHEs have a lower number of EF-problematic than CI-problematic sites in peri-urban areas used for background estimation enables us to presume that normalisation by Al can sometimes be a useful tool for adjustment of DMB values. The background in each LG becomes slightly variable and if EF values better reduce the influence of soil grain size as the main natural factor, a comparison of sites according to those values can be more reasonable than according to CI values. But properly choosing between two assessment indices for PHEs, i.e. either EF or CI becomes extremely important, because it enables more reasonable revealing of the influence of other factors (natural or anthropogenic).

We selected Al as the candidate reference element for normalisation, since it was one of the variable classifiers. When choosing between CI and EF, we took into account two restrictions for calculation of EF: (i) PHE should be correlated with the reference element used for normalisation [47] and (ii) PHE should have higher variability than the reference element [8]. In the urban area of Klaipėda, there were only a few PHEs which could be normalised, because they did not satisfy at least one of the restrictions [15]. Meanwhile in peri-urban Vilnius with lower anthropogenic influence, all PHEs appeared to have significant (p < 0.01) correlation with Al, but the magnitude of r_S values was different. Four PHEs, i.e. Ni, Cr, Cu and Zn (in comparison with Mn or Pb), show a higher correlation with Al (Table 4); further, they have higher RCV (31.0, 50.1, 32.2, 27.0%) than Al (RCV = 19.5%). Although Mn in comparison with the four aforementioned PHEs has relatively lower r_S with Al, its RCV (21.4%) also slightly exceeds that of Al. Therefore, these five PHEs will be represented by EF values as assessment indices. The expediency of Ni normalisation is obvious, since the difference between EF values obtained by G approach and those obtained without background differentiation is lower than the difference in the respective CI values (Figure 6). But as the RCV values of Ba and Pb (13.3%) are lower than those of Al, they will be represented by CI values.

4.5 What do the values of ASI_M show?

The most important is that in peri-urban areas, these values are an additional tool to check the quality of site-classification approaches. After site classification by G approach, neither clay nor silt fraction has any influence on ASI_M values of PHEs or values of the respective additive accumulation index Z_a, as is obvious from insignificant Spearman’s correlation coefficients r_S (Supplementary Table S9). Other site-classification approaches do not ensure total independence of ASI_M values on fine fraction components. The low-exactness classification approaches retain the greatest number of significant r_S (5), the medium-exactness P approach preserves their lower number (2) and the high-exactness classification T approach has only one. Further, for Ni, Cr, Cu and Zn, the number of EF-problematic sites revealed by G approach (Table 8) is not only lower than or equal to EF-problematic sites revealed using NMB values, but it is also either the lowest or rather low in comparison with other approaches. This fact is an additional indication that G approach reduces the influence of fine fraction best of all. But Mn demonstrates the opposite feature. The number of its EF-problematic sites obtained using G approach is higher in comparison not only with the respective numbers of EF-problematic sites revealed by other approaches or without background differentiation but also with the number of its CI-problematic sites.

We expected that the optimal site-classification approach and therefore the more precise DMB values as well as the proper choice between CI or EF would help us to reveal, according to values of assessment indices, the peri-urban sites with elevated DSI. But it is not strange that we could not fulfil this task, since DSI “usually leads to sites that are relatively uniformly contaminated” [4]. In some sub-sites, high values of assessment index seem to show contamination by local accidental waste (despite that sample was a composite from five sub-samples), but a comparison of sites according to values of ASI_M confirms that DSI is rather uniformly distributed, except for some problematic sites (Supplementary Table S10).

Interpretation of the processes responsible for higher ASI_M values in these sites is complicated. The difficulties in distinguishing between natural and anthropogenic sources with the help of EF were mentioned by Reimann and De Caritat [8] who stated that these indices can be elevated due to many reasons. Background differentiation in peri-urban minerogenic soil followed by optional normalisation reduced the influence of fine fraction, but not of other factors. Although the ASI_M values obtained by G approach are optimal, they do not necessarily reveal anthropogenic influence, since they can be related to some natural processes. For example, two EF-problematic sites are obtained by G approach for Ni (27 and 20) and two for Cu (26 and 20). We can speculate that sites 26 and 27 indicate increased DSI from the Vilnius–Kaunas highway (Figure 1), but we cannot state that site 20, which is far from pollution sources, is influenced by anthropogenic activity. Elevated values of assessment index most probably reflect the fact that even G approach provides insufficiently precise DMB values for this site. But the absence of EF-problematic sites of Cr and Zn, as well as of CI-problematic sites of Ba, enables us to presume that DSI is low in peri-urban territories. Some doubt arises only concerning one CI-problematic site of Pb (site 25, which is relatively close to the highway).

Another problem is the interpretation of the greatest number (5) of EF-problematic sites of Mn obtained by G approach (Supplementary Table S10). We can only say that although research in the central part of Vilnius [26] has shown that Mn correlates with both natural and anthropogenic factors, it is hardly believable that its EF-problematic sites in peri-urban areas indicate anthropogenic anomalies. First, the number of urban–industrial pollution sources of Mn is relatively low [48]. Second, concentration coefficients of Mn in dust from vents and stacks of some pollution sources in Vilnius are much lower than that of Ni, Zn, Cu and Cr [39]. Third, although Mn is a component of methylcyclopentadienyl manganese tricarbonyl, research along the highways in Canada [49] as well as near the Vilnius–Kaunas highway [50] does not show its significant accumulation, especially taking into account that the distance of the present study sites from roads is >15 m. The fact that Seleznev et al. [34] even use Mn as a conservative element for prediction of VB values of PHEs in urban puddle sediments is in favour of the natural origin of Mn, but revealing the natural reasons for elevated EF values of Mn needs separate research.

Since the G approach applied to sites with minerogenic topsoil reduces the influence not only of grain size but simultaneously of organic matter (Figure 5), very few other natural factors can be mentioned as possibly influencing ASI_M values, e.g. pedogenic processes, natural land use types or differences in relief elevation. Investigation of pedogenic processes is impossible when composite samples are considered. It is also complicated to reveal the influence of relief elevation, since its range is narrow (111–199 m) and most sites are at a great distance.

Hence, we tested only land use influence presuming that in our study this factor can be treated as natural, because we sampled only from forests and meadows. Publications that discuss land use influence present diagrams or tables which contain either the contents of PHEs [51] or the values of their assessment indices, e.g. PI [7]. Sometimes both the contents and EFs are given [52] and their diagrams differ. Our search, with the help of the U test for significant differences, according to both median contents in sites and ASI_M values also demonstrated that the results are not identical. Significantly (p < 0.01) higher median contents of Cr, Ni, Cu, Zn and Ba in meadow sites in comparison with forest sites are obviously predetermined by a significantly higher percentage of clay, silt and fine fractions. Meanwhile, the median values of CI in sites (CI_M) obtained after background differentiation reduce the influence of grain size; therefore, the differences between CI_M values are mainly insignificant (p ≥ 0.5), especially when DMB values are estimated using G approach for site classification (Supplementary Table S11). Similar results were obtained by Sun and Chen [7] when comparing the PI values of Cu, Ni, Pb, Zn and Cr in several land use types: only a few differences in PI values were significant (p < 0.05), but not those between forest and grass land use. The formulas of PI and CI are identical, but the difference is that the background differentiation used by Sun and Chen [7] was according to provinces and it is not clear whether they sufficiently reflect the natural peculiarities of parent rocks, topsoil texture and geochemical composition. But researchers’ conclusion – that the impact of land use on heavy metal risks is insignificant – is in accordance with our results. Our results concerning CI_M of Ni, Cr, Cu, Zn (slightly higher in meadows) and Mn (slightly higher in forests) are similar to those based on data from “Geochemical atlas of Lithuania” [29]. Of course, when DMB values obtained by G approach are used and additional normalisation by Al is applied to five PHEs, a significantly higher value for EF of Cu in meadows appears, but the ratio of median EF_M value in meadows to respective value in forests is low (1.13).

4.6 Some problematic aspects for future investigations

Although our results seem logical, further experiments are necessary taking into account the problem that grain-size and geochemical data are compositional (closed). To treat such data, the log–ratio approach was introduced by Aitchison [53] who proposed “additive log–ratio” (arl) and “centered log–ratio” (crl) transformations. Also the “isometric log–ratio” (ilr) [54] was demonstrated as useful in robust principal component analysis [55]. Despite statisticians’ huge input to the theoretical development of the compositional data analysis (CoDA) approach, according to Filzmoser et al. [56], this is “only recently being followed by a significant adoption” in geochemistry.

If the CoDA approach is selected for our data treatment, the inter-correlations of chemical elements and their correlation with grain-size fractions will be different. This will lead to new interpretation, changes in the most suitable major elements for site classification, another decision about expediency of normalisation and most suitable reference element and different values of median assessment indices. Possibly, their interpretation, especially as concerns Mn anomalies, will be easier. It may even be possible to distinguish sites with elevated DSI or to better reveal and explain different natural processes, provided that fine fraction influence is more effectively reduced. We presume, however, that when using the CoDA approach for our data treatment, the DMB values obtained after site classification using other major elements demonstrating a strong association with clay-sized or silt fractions will be quite similar to that obtained when classifying using relative median contents Al, K and Ti. This presumption can be done first of all because the latter site classification is similar to that using the relative median contents of all major elements (r_S = 0.870). Our present results also cannot be neglected, considering the findings of Fačevicova et al. [57] when studying the Devonian/Carboniferous boundary and comparing the results obtained by the standard approach with those obtained by the compositional approach and revealing that they have “synergic effects:” neither of them has prevalence and both can be further used in geochemical investigations. If so, the results of our research can be useful for future investigations of geochemical approaches based on inter-correlated major elements revealed using transformations that take into account the closure effect of data.