1. INTRODUCTION

Concentrations of risk elements, particularly Pb, in the soils of central European mountain areas are frequently elevated, which is usually attributed to the long-range transport of contaminants (LRTC) from metallurgy and coal combustion (Steinnes, 2009; Bińczycki et al., 2020; Ciarkowska & Miechówka, 2022). However, the spatial distribution of risk element concentrations in mountain soils is surprisingly variable even over relatively short spatial scales (Szopka et al., 2013; Ciarkowska & Miechówka, 2022), depending on the position on slopes (Jaguś & Skrzypiec, 2019; Ciarkowska & Miechówka, 2022; Demková et al., 2023; Rahmonov et al., 2024) and slope orientation (Ciarkowska & Miechówka, 2022; Miśkowiec, 2022). This spatial variability seems to contrast with the expectation of evenly distributed diffuse contamination attributed to LRTC (Reimann et al., 2008, 2019; Ciarkowska & Miechówka, 2022) and makes attributing the observations exclusively to LRTC questionable. Distinguishing between natural and anthropogenic causes of soil chemistry variability is, however, a challenging task (Reimann et al., 2008, 2018, 2019; Steinnes & Lierhagen, 2018; Bourgault et al., 2022; Ciarkowska & Miechówka, 2022; Adamec et al., 2024), which has not always been reliably performed in previous studies.

The variability of soil risk element concentrations can result from the complex bedrock geology typical for mountain ranges (Ciarkowska & Miechówka, 2022; Hošek et al., 2024), the translocation of elements in soil profiles by pedogenic processes (Waroszewski et al., 2016; Reimann et al., 2019; Bińczycki et al., 2020), and the lateral movement of soil elements/components by gravity or wind (Bourgault et al., 2022; Kellaway et al., 2022; Adamec et al., 2024). Disentangling these processes is complicated by the interrelation between bedrock properties, altitude, slope steepness, wind speed, and vegetation in mountain areas (Korzeniowska & Krąż, 2020; Ciarkowska & Miechówka, 2022); it is not always clear which variables are truly independent and can be considered primary reasons for the risk element distribution. Therefore, the topographic control of emissions' impacts on mountain soils deserves further study. While most preceding papers have focused on the impacts of LRTC, far fewer studies have addressed local contamination from point sources, even though their impact might be easier to understand. Our work specifically targets this gap. 

Many preceding studies of emission impacts on soils and biota have suffered from improper data mining. Particularly incorrect is the evaluation of “raw” concentrations of target elements – results of analyses not corrected for their natural variability. This approach is especially risky in areas with geochemically anomalous bedrock (Hošek et al., 2024), which is typical for many central European mountain ranges. Mathematicians working in geosciences have summarized sufficient arguments for using element ratios instead of “raw” concentrations of risk elements to avoid fundamental statistical errors (McKinley et al., 2016; Greenacre, 2019; Hron et al., 2021). Previous studies (Matys Grygar & Popelka, 2016; Greenacre, 2019; Hron et al., 2021) have shown that only empirically verified and interpretable element ratios should be adopted for reliable data interpretation. In sediments and soils derived from sedimentary rocks with variable quartz and/or calcium carbonate contents, risk element concentrations typically depend on the concentrations of certain lithogenic elements, such as Fe (Hernandez et al., 2003; Sterckeman et al., 2006; Matys Grygar & Popelka, 2016; Aytop et al., 2023; Lučić et al., 2023). The methodology for data mining based on interelement relationships, although still rarely employed in soils, clearly needs further development to combine mathematical rigor with geochemical interpretability, as justified by Greenacre (2019) and Hron et al. (2021). Thus, correct data processing is among the strategic aims of this work.

The target area of this study, the Beskid Mountains at the easternmost part of the Czech-Polish state border (Figure 1), has been impacted by local emissions of acid gases and risk elements from iron and steelmaking in Třinec, the most proximal source (Pavlů et al., 2015), and possibly also from Ostrava as a more distant point source (Matysek et al., 2008). For contamination mapping, the foliage of silver birch, Betula pendula (Roth.), which is suitable for monitoring soil damage caused by Zn emissions and acidification (Bílková et al., 2024). Actually, Zn concentration in soils is worth monitoring, because its content relative to Fe is very high in dust from iron- and steel making (Tsai et al., 2007). Birch can additionally accumulate Mn (Bílková et al., 2024), which is mobilized for plant uptake in acidified soils; iron- and steelmaking produces also considerable amount of acid gas emissions (Tsai et al., 2007). One of the strategic goals of this study is to reveal possible slope transfers of contaminants or other ions mobilized by soil acidification, hypothesized for manganese by Bílková et al. (2024). These slope transfers of ions, while not yet sufficiently understood, could partly explain the yet poorly understood spatial variability in risk element concentrations in mountain soils (Ciarkowska & Miechówka, 2022; Demková et al., 2023). The aims of this work were to map the impact of the Třinec point source on the mountain areas adjacent to the emission center, extending approximately 30 km in the local wind directions, to visualize emission distribution in geochemical maps, and to understand some of the unknown mechanism behind the spatial variability of emission contamination in mountain terrains.

2. STUDY AREA

2.1. Local geology and topography

The flat lowland around Třinec is covered by the Quaternary sedimentary rocks. The mountainous part of the study area features Carpathian flysch bedrock, i.e. sediments deposited in shallow marine environments during the late Mesozoic to early Cenozoic (Samec et al., 2010). Flysch rocks are characterized by alternating finer and coarser silicic strata, with lithologies ranging from mudstones to silty sandstones and conglomerates containing granules or pebbles, primarily composed of quartz. In the target area of this study, claystones are more abundant to the east of Třinec; here, limestones with occasional siderite lenses also occur. During the 18th and 19th centuries, these siderite-rich lenses were locally mined for iron production. Claystones and mudstones, which are finer strata in the flysch sequences, have weathered more easily, while the quartz-richer and coarser flysch strata have remained as residual rock fragments in the soils. In the 19th century, local iron ores and charcoal from the Beskids were replaced by hard coal from the Ostrava area and iron ore imported from abroad, which allowed for a substantial increase in iron and steel production in the 20th century and considerable contamination of that area (Vácha et al., 2015; Matys Grygar et al., 2023; Pavlů et al., 2015).

The characteristics of the flysch dictate the topography of the study area. Hills with mudstone outcrops are lower and have gentler slopes, while higher local peaks typically feature quartzite or sandstone outcrops, with conglomerates found on the steepest peaks. The bedrock thus contains variable proportions of quartz and occasionally carbonates, which primarily control its chemical composition, resulting in stratigraphic and spatial variability (Samec et al., 2010).

The lowland around Třinec is north-west of the Jablunkov Pass (Figure 1). Sampling was conducted along the north-eastern boundary of the Pass on slopes oriented toward the Třinec emission centre as well as in the mountain sites separated from Třinec by local peaks and ridges (Figure 1, Table 1). The subareas of sampling listed in Table 1 have variable distances from and orientations relative to the Třinec contamination source.

2.2. Local contamination of soil environment

The study is focused on an area east and southeast of the iron- and steelworks in Třinec (Figure 1). The entire northeastern part of the Czech Republic has been contaminated by heavy industry, particularly by iron and steel production (Matysek et al., 2008; Pavlů et al., 2015; Bílková et al., 2024) and coal combustion (Vácha et al., 2015). These local emissions peaked in the 1980s when more environmentally friendly technologies were introduced. Agricultural (Vácha et al., 2015; Matys Grygar et al., 2023) and forest soils (Pavlů et al., 2015) in that area have been contaminated by moderately volatile risk elements such as Cd, Pb, and Zn known to be associated with dust of iron and steelmaking (Tsai et al., 2007). Up to now, contamination has mostly been studied in organic horizons of forest soils, while mineral horizons in the Beskids has not yet been mapped, although actually deeper soil horizons are crucial for plants and environment in general.

Prevailing winds in the lowland around Třinec (Figure 1) have transported emissions toward the slopes of the Beskid Mountains in the Czech portion of the study area. Inverse meteorological situations are common in the lowland around Třinec, which contributed to forest dieback on the Beskid slopes in the 1980s, when industrial emissions from Třinec peaked.

3. METHODS

3.1. Sampling and analyses

Soil samples were collected from the top 10 cm of the A horizon after the removal of organic layers (O horizons). The O horizon was not sampled because Mn and Zn ions are mobile in organic matter, and historical loads of those elements has been transferred to the surface of the mineral horizons (Hernandez et al., 2003; Pavlů et al., 2015). Soil samples were sieved through a 2 mm mesh in the field. Residual rock fragments were collected near the soil pits, either from stones scattered on the local surface or directly from the soil pits, similar to the method used by Rambeau et al. (2010). The fragments were crushed in the field to avoid sampling surficial strata of the stones (at least 1 cm of the surface was omitted) that have been exposed to weathering.

Soil and rock fragments were dried, pulverized in a micromill, and analyzed using X-ray fluorescence (XRF) spectroscopy, as described elsewhere (Adamec et al., 2024; Hošek et al., 2024). An Epsilon 3 XRF spectrometer with an energy-dispersive detector (PANalytical, the Netherlands) was used, and single-element calibrations of the XRF analytical signal were performed with a series of certified reference materials, as detailed in previous papers (Adamec et al., 2024; Hošek et al., 2024).

Birch leaves were sampled from heights of 2 to 4 m using telescopic tree pruners (Fiskars, Finland), preferentially from the side of the tree crown exposed to light. The leaves were collected without petioles. Element analyses were performed after drying and pulverizing the samples using a micromill, with the same spectrometer utilized for the analysis of soils and rock fragments. Single-element calibration curves were obtained with certified reference materials and in-house standards, as described in a previous paper (Bílková et al., 2024).

Total number of samples included in statistical analyses (rock fragments, soils, and leaves) was 150. 

3.2. Data processing

Variability in the data series was characterized using quantile-based tools, specifically the median (instead of the mean) and the median of absolute deviations from the median (MAD) instead of the standard deviation. This approach has been recommended for geochemical data due to their non-Gaussian distribution (Reimann et al., 2005).

Interelement relationships were examined as a tool for correcting risk element concentrations for natural variability (Hernandez et al., 2003; Sterckeman et al., 2006).

Spearman rank correlation analysis was used, because it is more robust to outliers and it is not limited to linearity in interelement relationships. The results were presented as the Spearman rank correlation coefficient (ρ) and p-values for the null hypothesis. For predictions, robust regression analysis (RR) was employed, with the median absolute deviation (MAD) used as a relative measure of fit quality. RR was chosen because it is resistant to a certain percentage of outliers, which is a common feature of geochemical datasets (Reimann et al., 2005). RR was effective in finding relationships between target elements (Mn and Zn) and Fe as a reference element. The performance of Fe as a reference element for correcting element contents in rocks and soils for “quartz dilution effect” has already been demonstrated (Sterckeman et al., 2006; Matys Grygar & Popelka, 2016; Lučić et al., 2023).

Visual exploration of the spatial distribution of the geochemical data was conducted using geochemical maps with quartile-classified layers. This functionality is available in ArcGIS Online (Esri, USA). The classified layer is useful for expert-based visual evaluation of spatial patterns, which is advantageous for data with spatial noise, meaning a large fraction of outliers (Matys Grygar et al., 2023; Hošek et al., 2024; Adamec et al., 2024). A general problem with soil geochemical maps is that even in human-impacted areas, the natural variability of bedrocks must be taken into account (Hošek et al., 2024; Adamec et al., 2024). While emission contamination is most probable in topsoils, bioturbation or other disturbances to the soil profile can weaken or erase the contamination signal in the topsoil. The spatial noise can limit the identification of spatial patterns when using conventional geostatistical tools.

The geospatial analyses were performed using the Global Moran’s Index (Global Moran’s I) and the Getis-Ord hotspot/coldspot analysis (Gi*). Both functions are available in ArcGIS Pro (Esri, USA). Global Moran’s I ranges from -1 (indicating a spatially dispersed "chessboard-like" pattern) to +1 (indicating spatially associated values of the examined variable), while a value of 0 indicates no autocorrelation (random pattern). The strength of the autocorrelation in the Global Moran analysis is assessed using the Z-score (higher values indicate stronger spatial patterns) and the P-value, which assesses the probability of the null hypothesis of no spatial pattern. Gi* can be applied to datasets with high spatial autocorrelation to visualize where high and low values of the examined variable are associated to hotspots or coldspots, respectively. Another tool in ArcGIS Pro, geographically weighted regression (GWR), was used to examine spatial patterns in interelement relationships. This method shows the linear regression fit of two variables within a sliding window of local neighborhoods. It is suitable for determining whether relationships with weak correlation coefficients in the entire dataset are stronger in certain parts of the study area, or conversely, whether deviations from these correlations are more frequent in other areas.

4. RESULTS AND DISCUSSION

4.1. Element concentrations in rock fragments and soils

Soils in the target area are depleted in Si relative to the remaining rock fragments, while the concentrations of most other elements, such as Ca, Fe, Mn, and Zn, are comparable to or higher in soils than in the rocks. This observation indicates element-specific selective rock weathering during early pedogenesis. The major elements other than Si exhibit high Spearman rank correlations in both rock fragments and soils (Table 2), which document quartz dilution of all elements except Si as the primary driving force behind that locally fundamental source soil chemical variability.

The soil concentrations of Ca, Fe, Mn, and Zn are correlated with those in the rock fragments (Table 2) – all correlations listed in that table have probabilities higher than 99 %. The relationships between concentrations in soils and rock fragments are graphically represented as histograms of their differences in Figure 2; their shape is clearly non-Gaussian. The histograms reveal that the majority of soil samples are weakly but consistently enriched in Ca, Fe, Mn, and Zn relative to rock fragments, that is consistent with selective rock weathering. Selective weathering is also documented in histograms (Figure 2). The pairs of Ca and Fe concentrations (soil minus rock fragments) show an almost symmetrical distribution, indicating a more or less uniform soil enrichment relative to rock fragments. In contrast, the main populations in the Mn and Zn histograms are more skewed, with elevated heavy tails, indicating that certain portions of the soils are significantly enriched in Mn and Zn compared to the major population (Figure 2). In the histograms for Ca, Mn, and Zn, some minor populations of even more enriched soils are visible in the heavy tails of their distributions, indicating the presence of high outliers in the datasets (Figure 2). Part of the observed soil enrichment, particularly the skewed distributions of the Mn and Zn concentration pairs, can be attributed to emission-related enrichment in some soils.

4.2. Geochemical correction of Mn and Zn soil concentrations

Spearman correlations of element concentrations in soils and rock fragments (Table 2) and the shapes of the histograms (Figure 2) indicate that a substantial portion of Mn and Zn in soils is inherited from the bedrock. To correct soil concentrations of Mn and Zn for the bedrock variability, correlations of Mn and Zn concentrations with Fe concentrations were utilized (Tables 2-3). This correction cannot cover the flysch rock variability, but it is effective to capture geochemical associations relevant for the three major target elements in this study. There are two primary reasons for these correlations: (i) all three mentioned elements (Fe, Mn, and Zn) can occur in the crystal structures of the same bedrock minerals (carbonates, Fe oxides, clay minerals like biotite, or other possible dark minerals), and (ii) in soils, Mn and Zn ions are commonly sorbed on pedogenic Fe oxides and generally all those elements are most abundant in the fine size fractions of soils. Linear correlations were indeed found for both Mn and Zn with Fe in rock fragments and soils using robust regression (Table 3), which documents their inherent geochemical associations of those elements. Soil Fe concentrations could thus be used to predict soil Mn and Zn concentrations, provided that neither emission enrichment nor slope translocation of Mn and Zn had occurred, in a manner similar to what was done by Hernandez et al. (2003) and Sterckeman et al. (2006). Emissions and slope translocations would increase Zn and Mn relative to Fe for several reasons: 1) Zn is more volatile than Fe and Mn in high-temperature processes, leading to its enrichment in emissions relative to Fe (Hernandez et al., 2003; Tsai et al., 2007; Steinnes, 2009); 2) Fe and Mn emissions can often be neglected in the presence of a considerable pool of natural Fe and Mn in soils; 3) Mn is more mobile than Fe in acid soils within catenas (Watmough et al., 2007; Bourgault et al., 2022); 4) Mn is enriched in topsoils relative to Fe due to the accumulation of mobile Mn in tree leaves (Herndon et al., 2015); and 5) Fe(III) ions generally exhibit lower mobility in soils at pH > 3 compared to Mn and Zn (McDaniel et al., 1992; Pavlů et al., 2015).

One advantage of working with element ratios is the possibility to compare local rocks and soils with global backgrounds, regardless of variations in content of organic matter (McKinley et al., 2016), quartz (Matys Grygar & Popelka, 2016), or carbonates (Lučić et al., 2023). In the study area, Zn is correlated with Fe in rock fragments, with a proportionality constant of 28·10-4 in robust regression (Table 3), 31·10-4 in ordinary linear regression, and 28·10-4 for the median Zn/Fe ratio in the dataset. These values are close to the median Zn/Fe ratio of 31·10-4 calculated from analyses of local rocks published by Samec et al. (2010). These ratios are nearly twice the ratio of these elements in global continental rocks (17·10-4, Rudnick & Gao, 2003) and in French forest soils (Hernandez et al., 2003), indicating a certain enrichment of Zn in the local Beskid bedrock relative to offsite references. In soils, robust regression produced a slope for Zn vs. Fe equal to 46·10-4 (Table 3), indicating that local analyzed soils are further enriched in Zn relative to the local rocks, as it is also shown above (Figure 2). This overall enrichment of Zn in soils relative to local rocks can result from emissions; however, a certain contribution from natural processes could also be expected, as Zn is a biogenic element (Reimann et al., 2019) and is weakly accumulated in birch laves (Bílková et al., 2024). Spatial analysis is therefore needed to identify possible anthropogenic influences without interference from natural processes: while emissions should correlate with distance from the Třinec emission center, other processes could affect the entire study area equally. The spatial analysis is presented in the subsequent section.

Weaker correlations were observed for Mn and Fe concentrations in rock fragments and soils (Tables 2-3). Mn/Fe slopes of 103·10-4 and 158·10-4 were found in rock fragments and soils, respectively. These slopes are lower than the Mn/Fe ratio in the Upper Continental Crust (UCC), which is 195·10-4 (Rudnick & Gao, 2003). As with Zn, spatial analysis is crucial for evaluating the potential anthropogenic impacts on Mn/Fe soil enrichment.

4.3. Geochemical maps – spatial patterns of elevated Mn and Zn in soils

The results of the geostatistical analysis of Zn and Mn concentrations in rock fragments and soils are summarized in Table 4. The raw concentrations of Mn and Zn in rock fragments exhibit strong spatial autocorrelation (Table 4), meaning neighboring samples tend to have similar concentrations. Concentrations of Mn and Zn in rock fragments are elevated in a significant portion of the Czech side compared to the Polish sites that are farthest from Třinec (subareas 17 to 19 in Table 1), which could lead to biased evaluations of anthropogenic impacts on soils based solely on raw soil concentrations. These spatial autocorrelations of Mn and Zn in rock fragments are mitigated by normalization with Fe, based on linear regression between Mn (or Zn) and Fe (Table 3). For example, the Zn/Fe concentration ratio in rock fragments shows a weakly disordered spatial arrangement (Table 4, negative Global Moran’s I). The Zn-xFe, Mn/Fe, and Mn-xFe ratios in rock fragments, where x represents the slopes of the linear regressions (Table 3), exhibit a spatially random pattern (Table 4), indicating that the Fe correction effectively eliminated the spatial bias of bedrock Mn and Zn concentrations. A similar Fe correction should thus be applicable also in soils.

Lithologically corrected soil Zn demonstrates high spatial autocorrelation (Table 4). The Gi* analysis (Figure 3B) identified a prominent hotspot in the Ostrý Mountain (subarea 2), along with weaker associations of higher values in subareas 12 to 14 (Figure 3A) and a coldspot in subarea 17 (Figure 3B). However, elevated values of corrected soil Zn are not always identified as hotspots by the Gi* analysis, despite their clear association in the maps with classified layers (Figure 3A). This discrepancy is attributed to significant spatial noise in the geochemical maps, which can arise from the uneven distribution of element concentrations in soils and possibly also in parent rocks. To visualize the excess Zn in soils relative to rock fragments in a spatial context, geographically weighted regression (GWR) was employed.

The GWR in Figure 4 illustrates the relationships between Zn in soils and Zn in rocks in the neighborhood of each point. Green circles in Figure 4 indicate points where soil Zn is proportional to Zn in rock fragments, while purple shades highlight excess Zn in soils compared to rock fragments and blue shades indicate Zn depletion, all relative to neighboring sites. Clustering of purple shades in the GWR (Figure 4) points to Zn contamination based on a different principle than the preceding maps (Figure 3). The high color variability in the GWR map among neighboring sites reflects spatial noise, such as the impact of local translocation of ions or microtopographic variability. The GWR of soil Zn versus Zn in rock fragments (Figure 4), along with the map of classified Zn-xFe ratios (Figure 3A), confirmed a continuous occurrence of soil Zn contamination, particularly prevalent on the Czech side of the study area, specifically in subareas 2, 4, 5, and 12 to 14, all situated on slopes oriented toward Třinec (Table 1). In contrast, relatively low soil Zn concentrations were found in subareas 3 and 6, despite their proximity to the Třinec emission center. Here topography plays role: those sites are at lower altitudes. Visual inspection of the geochemical maps (Figure 3A), Gi* hotspots (Figure 3B), and the associations of magenta shades in the GWR maps (Figure 4) all reveal a similar pattern of soil contamination by Zn.

Lithologically corrected Mn concentrations in soils also show significant spatial autocorrelation (Table 4). Gi* hotspots of corrected soil Mn (Figure 3D) were found in subarea 2, that is the steepest area close to Třinec, and subareas 4, 5, and 7 in slopes of the highest peak in the study area. Elevated soil Mn is found in lower altitudes of ridges or slopes exposed to Třinec, that is indicated by arrows in Figure 3C. All those observations underline topographic control of the Mn distribution in the study area.

4.4. Factors controlling element concentrations in foliage

Macronutrient contents in leaves are regulated by homeostasis; thus, P, K, Ca, and Mg exhibit less variability than Mn and Zn (Table 5). The birch trees in the study area are deficient in K and Mg (Table 5), while the other nutrients are close to the plants' physiological needs. Foliar Mn and Zn concentrations in birch leaves from the study area are an order of magnitude higher than the physiological requirements of birch (Table 5) and, in most subareas, exceed concentrations found in reference areas elsewhere in the Czech Republic and Finland (Figures 5-6). Finland has actually soil nutrition and climatic factors similar to those of the Beskids and thus the data by Hytönen et al. (2014) are suitable as reference for this Beskid study.

Foliar Mn and Zn concentrations in the study area exhibit spatial autocorrelation with probability higher than 99 % (Table 4). Before evaluating the foliar concentrations of Mn and Zn and their spatial ordering, two features were examined: i) the relationships between their concentrations in foliage, rock fragments, and soils, and ii) the interrelations between the difference between foliar and soil element concentrations and soil concentration of Ca (Table 6). The purpose of that evaluation was to assess whether foliar concentrations of elements are simply and directly inherited from soils, or if other factors should also be considered. In our study area, foliar Zn in birch is correlated with soil Zn but is also negatively correlated with soil Ca (Table 6). The impact of Ca deficiency in soils on Zn uptake is demonstrated by the correlations of bioaccumulation factors (BAF, the ratios of foliar to topsoil concentrations) for these elements (Table 7). Therefore, Zn bioaccumulation (the content in leaves relative to that in soils) is amplified by low soil Ca concentrations. This situation is favorable for emission mapping in the surroundings of iron- and steelmaking in Třinec using birch leaves, as both soil Ca deficiency and excess Zn in soils result from the emissions: excess Zn from dust (Tsai et al., 2007) and Ca deficiency in soils caused by nutrient leaching from acidified mountain slopes.

Foliar Mn rarely correlates with total soil Mn concentration (Bílková et al., 2024); instead, it responds more to soil chemistry, including nutrient insufficiency and high soil acidity (Kogelman & Sharpe, 2006; Watmough et al., 2007; Jordan et al., 2019; Turpault et al., 2021; Wills et al., 2023; Bílková et al., 2024). Indeed, foliar Mn in the entire Beskid dataset does not correlate with rock or soil Mn concentrations (Table 6), documenting foliar Mn is not simply inherited from soil Mn. Nevertheless, notably concentrations (Table 6), documenting foliar Mn is not simply inherited from soil Mn. Nevertheless, notably high foliar and soil Mn concentrations are observed in subareas 2, 4, 7, and 12, i.e.in slopes exposed to Třinec emission centers. Nutrient insufficiency plays a certain role in Mn uptake by leaves, as evidenced by the weak negative rank correlation of foliar Mn with soil Ca (Table 6) and foliar K (Table 5), as well as the Spearman correlation of bioaccumulation factors (BAF) for Mn and Ca (Table 7). The elevated uptake of Mn and Zn by birch in the context of insufficient soil Ca could amplify the impacts of soil damage caused by emissions (Ca deficiency due to nutrient leaching). This pattern has also been observed in other tree species elsewhere in world (Kogelmann & Sharpe, 2006; Herndon et al., 2015; Richardson, 2017; Santos & Herndon, 2023). These findings and citations support the proposal of mapping soil acidification damage by foliar birch analyses (Bílková et al., 2024).

4.5. Topographic control of Mn and Zn in soils and leaves and problems of spatial noise in geochemical mapping

Figure 3 shows that the corrected Zn and Mn contents in soils are not homogeneously elevated near the local contamination source. In some sites close to Třinec, the soil enrichment by Mn and Zn is not as significant as might be expected from the spatial proximity, particularly in subareas 1, 3, and 6 (Figures 3A and C). Due to the uneven distribution of contaminants, the plots of soil Mn or Zn versus distance from the source exhibit a “triangular” pattern for all points (Figures 7-8). This probabilistic pattern occurs when the contamination signal is more likely closer to the contamination source, but some sites have not recorded it at all, while most sites have been impacted to varying degrees. Such a pattern cannot be treated by common tools of geochemical data exploration, because the correlation between explained and explanatory variables is not deterministic but rather probabilistic. There are no conventional methods for processing such patterns; specifically, ordinary (least squares) regression is inadequate, as it is suited for unequivocal deterministic patterns. Quantile regression has been tested, because it is based on quantile processing suitable for non-Gaussian data distributions and data with outliers, common for geochemistry data, Quantile regression demonstrated that the triangular pattern exhibits definitive regularity (Tables 8-9, statistical significance of regression is evaluated in Tables 10-11). In particular, quantile regression of non-corrected soil Zn and distance to Třinec showed probability well above 99 % (Table 10). Still the spatial noise makes identification of the distance control of emission load for individual samples not really persuading.

To reduce the impact of spatial noise and account for the probabilistic nature of the contamination record, the dataset was post-stratified to 19 geographic subareas (Table 1), and only the medians for these subareas were plotted as variables (Figures 5-6). There are generally two groups of subareas in Figures 5A-B: the group highlighted by the ellipses in Figures 5A-B shows that soil Zn increases toward the Třinec emission centre for majority of subareas; the subareas in ellipse are primarily located at altitudes close to 600 m a.s.l. Outside the ellipse, there are weakly impacted subareas near Třinec (subareas 1 and 3) and in the Velká Čantoryje Mountain (subareas 7 and 8), consistently with the lack of contamination hotspots in those locations (no Gi* hotspots were identified there, Figure 3C). Similar subarea-specific results were obtained for the geographic controls of foliar Zn concentrations (Figures 5C-D). Foliar Zn maxima were found at distances of 10 to 15 km from Třinec (consistently with findings by Matysek et al., 2008) and altitudes of 600-700 m a.s.l. Conversely, foliar Zn was not elevated in lower altitude sites closer to Třinec (subareas 1, 2, and 13), where high soil Ca concentrations suppress Zn uptake by birch (Tables 6-7, Figure 4). 

The dependence of soil and foliar Mn on distance and altitude was similar to that of Zn. The highest corrected Mn concentrations in soils (Figure 6) were found at altitudes of 600 to 700 m a.s.l., both at distance of 10 to 20 km from Třinec and in the immediate vicinity of Třinec. The relationship between soil Mn concentrations and the distance to Třinec is evident for subareas at these medium altitudes that are appropriately oriented relative to Třinec (all subareas in ellipses in Figures 6A-B). Maximum foliar Mn concentrations were found at similar distances from Třinec and at similar altitudes as the maxima for soil Mn (Figures 6C-D). 

The maximum concentrations of birch foliar Mn at an altitude of approximately 800 m a.s.l. on the Beskid slopes south and southeast of Třinec were reported by Bílková et al. (2024). The novel dataset presented in this work is more representative because of considerably denser sampling than that in Bílková et al. (2024). The novel results in this work thus refined the altitudinal emission load maximum to somewhat lower altitudes of 600 to 700 m a.s.l. Additionally, low soil Mn concentrations were identified in peaks at altitudes between 900 and 1000 m a.s.l.

Emission load maximum at certain altitudes (Figures 6B and D) can be attributed to horizontal (occult) deposition (also called to the “seeder-feeder” effect). Enhanced horizontal emission deposition, e.g. by locally freezing fog or ice accretion have been studied by Doušová et al. (2007) and Voldrichova et al. (2014). Szopka et al. (2013) and Ancora et al. (2021) attributed altitudinal maxima of emission load to direct scavenging (interception) of droplets from fog or the lower parts of clouds. Horizontal and vertical deposition are expected to differ in the distribution of contaminants across the landscape. The vertical deposition (through rain and snow, which are less contaminated) should be spatially more uniform and primarily dependent on total precipitation, but horizontal deposition (such as through ice accretion or smog, which tends to be more concentrated) can be significantly influenced by terrain shape, vegetation, and other factors important for dust or fog retention (Szopka et al., 2013). Understanding to emission load heterogeneity in mountain terrains is poor – and improvement of this state needs more case studies considering local topography.

In the study area in the Beskids, the altitudinal maximum of soil and foliar Mn and Zn, found between 600 and 700 m a.s.l., can result from thermal inversions typical for the valley southeast of Třinec, that is the entry to the Jablunkov Pass (Figure 1). The local inversions would lead to horizontal deposition of the most concentrated contaminants under conditions of poor vertical atmospheric mixing, with weak north-westerly winds transporting emissions from Třinec toward the north-eastern edges of the Jablunkov Pass (wind rose in Figure 1). The digital terrain model (DTM) of the study area southeast of Třinec (Figure 1) shows a major ridge along the Czech-Polish state border, which is important geomorphic boundary in this study. This ridge has a lowest point at 680 m a.s.l., a typical altitude of approximately 800 m a.s.l., and peaks ranging from 900 to 1000 m a.s.l. The ridge would prevent the dilution of Třinec emissions by surrounding air masses during thermal inversions, effectively “focusing” emissions at the altitudinal maximum between 600 and 700 m a.s.l. It is consistent with the altitudinal maxima found in this study (ellipses in Figures 5B and 6B).

4.6. Processing probabilistic impact of emissions on environmental samples by quantile regression

Quantile regression is a robust statistical method first introduced by Koenker & Bassett (1978). It allows for the modelling of the relationship between a response variable and one or more predictor variables at different quantiles of the response variable's distribution. Unlike the standard regression modelling, which focuses on the mean response, quantile regression provides a more comprehensive picture of the relationship between variables by estimating the conditional quantiles of the response variable. In the context of soil geochemistry, quantile regression can be particularly useful. Quantile regression is less sensitive to outliers and can handle heteroscedastic data, and actually presence of outliers and non-normality are encountered in environmental datasets as a rule rather than exception. The application of quantile regression to complex environmental conditions such as site-specific characteristics of soil chemistry or uneven distribution of emissions thus has the potential to improve our understanding of these complex interactions. The subsequent statistical analyses using quantile regression were made using ‘robustbase’ and ‘quantreg’ packages in the R statistical software (R Core Team, 2023). 

Figure 7 shows raw data on soil Zn concentration versus distance to Třinec emission center for the entire dataset. The points are arranged roughly to a triangle, with the most remote samples showing low Zn concentration and growing probability the Zn concentration is increased close to Třinec but still with uncontaminated soils even in sites very close to Třinec. Processing the data by regression does not make sense.

Quantile regression (Table 8) catches this triangular-like data pattern. In Figure 7, the median (red line) shows the central tendency of the dataset, decreasing steadily with growing distance. The 10th and 90th percentiles (blue dashed lines) encompass approximately 80% of the data points, the 25th and 75th percentiles (green dotted lines) show the interquartile range. The lines for different quantiles are not parallel, which suggests that the relationship between Zn and distance varies across the distribution of Zn. The slopes appear steeper for higher quantiles (90th, 75th) compared to lower quantiles (10th, 25th). For lower quantiles (10th), this may indicate negligible or no contamination record and thus mainly natural control of the actual soil Zn concentrations. This natural control is even more obvious in the plot in Figure 8 with Fe corrected Zn concentration as the independent variable and results of its quantile regression (Table 11): here the lines for the lower quantiles have insignificant slopes with small absolute values and intercepts low or even negative and thus showing no control by the distance to Třinec.

It would be possible to examine significance of the quantile regression coefficients using rigorous statistical tests run separately for each quantile. Particularly suitable could be the non-iterative diagonal estimator (NID), the use of which is shown in Tables 10-11. The test using NID show the raw concentration data are dependent on the distance to Třinec, however, the significance of that dependence is weaker in the Fe corrected data. However, still the quantile regression and NID provides proof the distance to Třinec show certain control on soil Zn concentrations. 

The use of quantile regression in this work is to demonstrate there are objective statistical test to confirm probabilistic control of the emission centre in the actual dataset, but perhaps it is not always needed to use rigorous statistical methods, as researchers can also rely more on the sense of the results. In the work in the main body of the article we thus used data post-stratification approach, that is more suitable for visual data evaluation.

4.7. Downslope transport of contaminants in mountain soils and advantage of foliar mapping

Topsoil erosion, inevitable on steep slopes, can also introduce spatial noise into the mapping datasets. This occurs due to the loss of contaminants from erosion-exposed areas of the slopes, often manifested by a thin or missing O horizon that can be apparent even in the field. Consequently, accumulation of airborne and soluble contaminants in footslopes should be expected. Our results indeed demonstrate the downslope transport of Mn across catenas, similar to observations made by others (Bourgault et al., 2022; Bílková et al., 2024): soil Mn concentrations tend to increase toward lower altitudes (red arrows in Figures 3C and 6A-B). In the Velká Čantoryje Mountain peak (subarea 6) and its slopes (subareas 7 and 8), this pattern is also evident for Ca and Mn-xFe in soils (Figure 9). This trend is also noticeable for Zn on the Ostrý and Čantoryje slopes (subareas 2 and 7), where soil and foliar Zn and Mn concentrations are highly variable at short scales, with elevated values being more frequent in footslopes (Figures 3A, 3C, and 4). The resulting spatial noise complicates the use of conventional statistical tools, including spatial autocorrelation analyses by Gi* or Moran’s I. Downslope transport of metal ions has not been considered in most preceding studies, although some have been conducted in very steep terrains (Ciarkowska & Miechówka, 2022; Miśkowiec, 2022; Rahmonov et al., 2024).

Many previous studies on mountain soil contamination have focused on organic soil horizons (Matysek et al., 2008; Steinnes & Lierhagen, 2018) or unspecified topmost soil horizons (Ciawkowska & Miechówka, 2022), where anthropogenic impacts are more pronounced than in mineral soil horizons, making them apparently easier to identify. However, downslope transfer of the litter with contamination, topsoil strata, and mobility of some contaminants complicates already the decision on what soil depth should be sampled to obtain a representative picture of historical contamination. Sampling typical local plants that actively uptake the target risk elements from emissions or have been impacted by acid rains could help: plants can effectively "average" the soil profile through plant roots. The plant uptake is easily characterized by foliar element analyses, as it is demonstrated for Mn and Zn in birch leaves in this work. Such biomonitoring, which is also employed in this paper, should be used more frequently in studies of real environmental impacts of mountain contamination.

5. CONCLUSIONS

Soil contamination studies in mountain areas require the geochemical correction of element concentrations for their natural variability, which results partly from variations in bedrock composition, differential weathering, natural changes in soil-to-bedrock element ratios due to pedogenesis, and the translocation of element ions within soil profiles and across catenas. Our work confirmed the spread of Zn emissions from ferrous metallurgy in Třinec to the forest soils of the Beskids, with maximum contamination observed at approximately 15 km from the source and at altitudes 600 to 700 m a.s.l. Parallel sampling of topsoils and birch leaves proved effective in mapping Zn contamination. The Třinec emissions also increased the concentrations and bioavailability of Mn in topsoils. Birch leaves are suitable for monitoring elevated levels of soil Mn and Zn; their uptake by birch is enhanced by Ca deficiency in soils, which is exacerbated by acid emissions and slope terrains. The probabilistic nature of contamination records and the spatial noise on catenas must be taken into account when mapping the historical impacts of contamination on soils, necessitating further improvements in data processing. Maximum local contamination around the studied point source in Třinec shows an altitudinal peak close to the altitude of the prominent local mountain ridge, attributed to dust scavenging during meteorological conditions of poor pollutant dispersion in the atmosphere (thermal inversions under weak winds). Therefore, the shape of the terrain is an important factor influencing contamination loads in the mountain soils. This work contributes to a better understanding of the uneven distribution of element contaminations in mountain terrains.

Acknowledgment

TMG planned work, participated in fieldwork, and prepared the manuscript, GB provided financial support, performed field GIS, participated in fieldwork, and prepared geochemical maps, JE participated in fieldwork and prepared part of maps, ML performed less conventional statistical analyses. VH participated in fieldwork. The authors thank to Š. Matysová and Č. Matys who helped with field sampling. P. Vorm and M. Maříková (Institute of Inorganic Chemistry CAS) processed all samples in laboratory and performed all XRF analyses. This study was financially supported by the UJEP via student grant for G. Bílková (UJEP-SGS-2022-44-002-3). The manuscript language correction was performed using OpenAI GPT-4o followed by careful reading to correct several mistakes currently typical for that tool. 

REFERENCES

• Adamec, S., Tůmová, Š., Hošek, M., Lučić, M. & Matys Grygar, T., 2024. Pitfalls of distinguishing anthropogenic and geogenic reasons for risk elements in soils around coal-fired power plant: from case study in NW Czech Republic to general recommendations. J Soils Sediments, 24(3), 1274-1288, https://doi.org/10.1007/s11368-024-03726-9
• Ancora, S., Dei, R., Rota, E., Mariotti, G., Bianchi, N. & Bargagli, R., 2021. Altitudinal variation of trace elements deposition in forest ecosystems along the NW side of Mt. Amiata (central Italy): Evidence from topsoil, mosses and epiphytic lichens. Atmospheric Pollution Research, 12(1), 101200, https://doi.org/10.1016/j.apr.2021.101200
• Aytop, H., Koca, Y.K. & Şenol, S., 2023. The importance of using soil series‑based geochemical background values when calculating the enrichment factor in agricultural areas. Environ Geochem Health, 45(8), 6215–6230, https://doi.org/10.1007/s10653-023-01640-6
• Bergmann, W., 1988. Ernährungsstörungen bei Kulturpflanzen (Entstehung, visuelle und analytische Diagnose). VEB Gustav Fischer, Jena, p. 381.
• Bílková, G., Königová, M., Hýlová, V., Elznicová, J., Von Suchodoletz, H., Flem, B. & Matys Grygar, T., 2024. Factors controlling Mn and Zn contents in leaves of silver and downy birch in acidified soils of Central Europe and Norway. Environ Sci Pollut Res, 31(6), 9642–9660, https://doi.org/10.1007/ s11356-024-31837-w
• Bińczycki, T., Weber, J., Mielnik, L. & Asensio, C., 2020. Lead isotope ratios in Podzol profiles as a tracer of pollution source in the subalpine zone of the Karkonosze National Park, Sudety Mts (south-western Poland). Catena, 189, 104476, https://doi.org/10.1016/j.catena.2020.104476
• Bourgault, R.R., Ross, D.S., Bailey, S.W., Perdrial, N. & Bower, J., 2022. Groundwater input drives large variance in soil manganese concentration and reactivity in a forested headwater catchment. Soil Science Soc of Amer J, 86(6), 1553–1570, https://doi.org/10.1002/saj2.20439
• Ciarkowska, K. & Miechówka, A., 2022. Identification of the factors determining the concentration and spatial distribution of Zn, Pb and Cd in the soils of the non-forest Tatra Mountains (southern Poland). Environ Geochem Health, 44, 4323–4341, https://doi.org/10.1007/s10653-022-01201-3
• Doušová, B., Erbanová, L. & Novák, M., 2007. Arsenic in atmospheric deposition at the Czech–Polish border: Two sampling campaigns 20 years apart. Science of The Total Environment, 387(1–3), 185–193, https://doi.org/10.1016/j.scitotenv.2007.06.028
• Demková, L., Bobuľská, L., Árvay, J., Homolová, Z., Michalko, M. & Bálintová, M., 2023. Potentially toxic elements in soil and air along an altitudinal gradient in Tatra National Park. Journal of Geochemical Exploration, 252, 107268, https://doi.org/10.1016/j.gexplo.2023.107268
• Greenacre, M., 2019. Variable Selection in Compositional Data Analysis Using Pairwise Logratios. Math Geosci, 51, 649–682, https://doi.org/10.1007/ s11004-018-9754-x
• Hernandez, L., Probst, A., Probst, J.L. & Ulrich, E., 2003. Heavy metal distribution in some French forest soils: evidence for atmospheric contamination. Science of The Total Environment, 312(1–3), 195–219, https://doi.org/10.1016/S0048-9697(03)00223-7
• Herndon, E.M., Jin, L., Andrews, D.M., Eissenstat, D.M. & Brantley, S.L., 2015. Importance of vegetation for manganese cycling in temperate forested watersheds. Global Biogeochem. Cycles, 29, 160–174, https://doi.org/10.1002/2014GB004858
• Hošek, M., Pavlíková, P., Šoltýs, M., Tůmová, Š. & Matys Grygar, T., 2024. Distinguishing Geogenic Load and Anthropogenic Contribution to Soil Contamination in Mineralised Mountain Landscape of Ore Mountains (Czech Republic) Using Cumulative Distribution Functions. Land, 13(2), 218, https://doi.org/10.3390/land13020218
• Hron, K., Coenders, G., Filzmoser, P., Palarea-Albaladejo, J., Faměra, M. & Matys Grygar, T., 2021. Analysing Pairwise Logratios Revisited. Math Geosci, 53(7), 1643–1666, https://doi.org/ 10.1007/s11004-021-09938-w
• Hytönen, J., Saramäki, J. & Niemisto, P., 2014. Growth, stem quality and nutritional status of Betula pendula and Betula pubescens in pure stands and mixtures. Scandinavian Journal of Forest Research, 29, 1–11, http://dx.doi.org/10.1080/02827581.2013.838300
• Jaguś, A. & Skrzypiec, M., 2019. Toxic Elements in Mountain Soils (Little Beskids, Polish Carpathians). Journal of Ecological Engineering, 20, 197–202, https://doi.org/10.12911/22998993/92694
• Jordan, J., Cernak Sr, R.S. & Richardson, J.B., 2019. Exploring the role of soil geochemistry on Mn and Ca uptake on 75-year-old mine spoils in western Massachusetts, USA. Environ Geochem Health, 41, 2763-2775, https://doi.org/10.1007/s10653-019-00339-x
• Kellaway, E.J., Eimers, M.C. & Watmough, S.A., 2022. Liming legacy effects associated with the world's largest soil liming and regreening program in Sudbury, Ontario, Canada. Science of the Total Environment, 805, 150321, https://doi.org/10.1016/ j.scitotenv.2021.150321
• Koenker, R. & Bassett, G., 1978. Regression quantiles. Econometrica, 46, 33–50, https://doi.org/ 10.2307/1913643
• Kogelmann, W.J. & Sharpe, W.E., 2006. Soil Acidity and Manganese in Declining and Nondeclining Sugar Maple Stands in Pennsylvania. Journal of Environmental Quality, 35(2), 433–441, https://doi.org/10.2134/jeq2004.0347
• Korzeniowska, J. & Krąż, P., 2020. Heavy Metals Content in the Soils of the Tatra National Park Near Lake Morskie Oko and Kasprowy Wierch—A Case Study (Tatra Mts, Central Europe). Minerals, 10(12), 1120, https://doi.org/10.3390/min10121120
• Lučić, M., Mikac, N., Vdović, N., Bačić, N., Dinis, P. & Milačič, R., 2023. Distinguishing between natural and anthropogenic sources of potentially toxic elements in sedimentary materials along the Sava River (Slovenia, Croatia). Applied Geochemistry, 151, 105619, https://doi.org/10.1016/j.apgeochem.2023.105619
• Matys Grygar, T., Elznicová, J., Tůmová, Š., Kylich, T., Skála, J., Hron, K. & Álvarez Vázquez, M.Á., 2023. Moving from geochemical to contamination maps using incomplete chemical information from long-term high-density monitoring of Czech agricultural soils. Environ Earth Sci, 82, 6, https://doi.org/10.1007/s12665-022-10692-3
• Matys Grygar, T. & Popelka, J., 2016. Revisiting geochemical methods of distinguishing natural concentrations and pollution by risk elements in fluvial sediments. J Geochem Explor, 170, 39–57, https://doi.org/10.1016/j.gexplo.2016.08.003
• Matysek, D., Raclavska, H. & Raclavsky, K., 2008. Correlation Between Magnetic Susceptibility and Heavy Metal Concentrations in Forest Soils of the Eastern Czech Republic. JEEG, 13(1), 13–26, https://doi.org/10.2113/JEEG13.1.13
• McDaniel, P.A., Bathke, G.R., Buol, S.W., Cassel, D.K. & Falen, A.L., 1992. Secondary Manganese/Iron Ratios as Pedochemical Indicators of Field‐Scale Throughflow Water Movement. Soil Science Soc of Amer J, 56(4), 1211–1217, https://doi.org/10.2136/ sssaj1992.03615995005600040034x
• McKinley, J.M., Hron, K., Grunsky, E.C., Reimann, C., De Caritat, P., Filzmoser, P., Van Den Boogaart, K.G. & Tolosana-Delgado, R., 2016. The single component geochemical map: Fact or fiction? Journal of Geochemical Exploration, 162, 16–28, https://doi.org/10.1016/j.gexplo.2015.12.005
• Miśkowiec, P., 2022. The impact of the mountain barrier on the spread of heavy metal pollution on the example of Gorce Mountains, Southern Poland. Environ Monit Assess, 194, 663, https://doi.org/10.1007/s10661-022-10316-0
• Pavlů, L., Drábek, O., Borůvka, L., Nikodem, A. & Němeček, K., 2015. Degradation of forest soils in the vicinity of an industrial zone. Soil Water Res, 10(2), 65–73, https://doi.org/10.17221/220/2014-SWR
• R Core Team, 2023. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria, http://www.R-project.org
• Rahmonov, O., Sobala, M., Środek, D., Karkosz, D., Pytel, S. & Rahmonov, M., 2024. The spatial distribution of potentially toxic elements in the mountain forest topsoils (the Silesian Beskids, southern Poland). Scientific Reports, 14, 338, https://doi.org/10.1038/s41598-023-50817-7
• Rambeau, C.M.C., Baize, D., Saby, N., Matera, V., Adatte, T. & Föllmi, K.B., 2010. High cadmium concentrations in Jurassic limestone as the cause for elevated cadmium levels in deriving soils: a case study in Lower Burgundy, France. Environ Earth Sci, 61(8), 1573-1585, https://doi.org/10.1007/s12665-010-0471-0
• Reimann, C., Filzmoser, P. & Garrett, R.G., 2005. Background and threshold: critical comparison of methods of determination. Science of the Total Environment, 346(1–3), 1–16, https://doi.org/ 10.1016/j.scitotenv.2004.11.023
• Reimann, C., Flem, B., Arnoldussen, A., Englmaier, P., Finne, T.E., Koller, F. & Nordgulen, Ø., 2008. Reply to the comments on “The biosphere: a homogeniser of Pb-isotope signals” by Richard Bindler and William Shotyk. Applied Geochemistry, 23(8), 2527–2535, https://doi.org/10.1016/j.apgeochem.2008.05.005
• Reimann, C., Fabian, K., Birke, M., Filzmoser, P., Demetriades, A., Négrel, P., Oorts, K., Matschullat, J., de Caritat, P. & GEMAS Project Team, 2018. GEMAS: Establishing geochemical background and threshold for 53 chemical elements in European agricultural soil. Applied Geochemistry, 88, 302–318, https://doi.org/ 10.1016/j.apgeochem.2017.01.021
• Reimann, C., Fabian, K., Flem, B. & Englmaier, P., 2019. The large-scale distribution of Cu and Zn in sub- and topsoil: Separating topsoil bioaccumulation and natural matrix effects from diffuse and regional contamination. Science of The Total Environment, 655, 730–740, https://doi.org/10.1016/j.scitotenv.2018.11.248
• Richardson, J.B., 2017. Manganese and Mn/Ca ratios in soil and vegetation in forests across the northeastern US: Insights on spatial Mn enrichment. Science of the Total Environment, 581–582, 612–620, http://dx.doi.org/10.1016/j.scitotenv.2016.12.170
• Rudnick, R.L. & Gao, S., 2003. Treatise on geochemistry: composition of the continental crust. In: Rudnick, R.L., Holland, H.D. & Turekian, K.K. (Eds), The Crust, Treatise on Geochemistry, 3. Elsevier-Pergamon, Oxford, pp. 1–64, https://doi.org/ 10.1016/B0-08-043751-6/03016-4
• Samec, P., Kynický, J. & Cheng, X., 2010. Exploration of element contents in the selected rocks of the Slezské Beskydy Mts. (Czech Republic). Acta Mus. Beskid., 2, 2010.
• Santos, F. & Herndon, E., 2023. Plant-soil relationships influence observed trends between manganese and carbon across biomes. Global Biogeochemical Cycles, 37, e2022GB007412, https://doi.org/ 10.1029/2022GB007412
• Steinnes, E., 2009. Comment on ‘‘Geochemical gradients in soil O-horizon samples from southern Norway: Natural or anthropogenic?” by C. Reimann, P. Englmaier, B. Flem, L. Gough, P. Lamothe, Ø. Nordgulen, and D. Smith. Applied Geochemistry, 24, 2019–2022, https://doi.org/10.1016/j.apgeochem.2009.06.009
• Steinnes, E. & Lierhagen, S., 2018. Geographical distribution of trace elements in natural surface soils: Atmospheric influence from natural and anthropogenic sources. Applied Geochemistry, 88, 2–9, https://doi.org/10.1016/j.apgeochem.2017.03.013
• Sterckeman, T., Douay, F., Baize, D., Fourrier, H., Proix, N., Schwartz, C. & Carignan, J., 2006. Trace element distributions in soils developed in loess deposits from northern France. European Journal of Soil Science, 57, 392–410, https://doi.org/10.1111/j.1365-2389.2005.00750.x
• Szopka, K., Karczewska, A., Jezierski, P. & Kabała, C., 2013. Spatial distribution of lead in the surface layers of mountain soils, an example from the Karkonosze National Park, Poland. Geoderma, 192, 259–268, https://doi.org/10.1016/j.geoderma.2012.08.022
• Tsai, J. -H., Lin, K.-H., Chen, C.-Y., Ding, J.-Y., Choa, C.-G. & Chiang, H.-L., 2007. Chemical constituents in particulate emissions from an integrated iron and steel facility. J Haz Mat 147, 111–119, https://doi.org/10.1016/j.jhazmat.2006.12.054
• Turpault, M.-P., Kirchen, G., Calvaruso, C., Redon, P.-O. & Dincher, M., 2021. Exchanges of major elements in a deciduous forest canopy. Biogeochemistry, 152, 51–71, https://doi.org/ 10.1007/s10533-020-00732-0
• Vácha, R., Skála, J., Čechmánková, J., Horváthová, V. & Hladík, J., 2015. Toxic elements and persistent organic pollutants derived from industrial emissions in agricultural soils of the Northern Czech Republic. Journal of Soils and Sediments, 15, 1813–1824, https://doi.org/10.1007/s11368-015-1120-8
• Voldrichova, P., Chrastny, V., Sipkova, A., Farkas, J., Novak, M., Stepanova, M., Krachler, M., Veselovsky, F., Blaha, V., Prechova, E., Komarek, A., Bohdalkova, L., Curik, J., Mikova, J., Erbanova, L. & Pacherova, P., 2014. Zinc isotope systematics in snow and ice accretions in Central European mountains. Chemical Geology, 388, 130–141, https://doi.org/10.1016/j.chemgeo.2014.09.008
• Waroszewski, J., Egli, M., Kabala, C., Kierczak, J. & Brandova, D., 2016. Mass fluxes and clay mineral formation in soils developed on slope deposits of the Kowarski Grzbiet (Karkonosze Mountains, Czech Republic/Poland). Geoderma, 264, 363–378, https://doi.org/10.1016/j.geoderma.2015.08.044
• Watmough, S.A., Eimers, M.C. & Dillon, P.J., 2007. Manganese cycling in central Ontario forests: Response to soil acidification. Applied Geochemistry, 22, 1241–1247, https://doi.org/ 10.1016/j.apgeochem.2007.03.039
• Wills, A., Beier, C.M., Lawrence, G.B. & Drake, J.E., 2023. Foliar chemical composition and respiration rates of sugar maple (Acer saccharum) and American beech (Fagus grandifolia) trees across a gradient of soil acidification. Biogeochemistry, 162, 327–341, https://doi.org/10.1007/s10533-022-01010-x