Geostatistical Tools to Assess Existing Monitoring Network of Forest Soils in a Mountainous National Park

Abstract

Environmental changes in national parks are generally subject to constant observation. A particular case is parks located in mountains, which are more vulnerable to climate change and the binding of pollutants in mountain ranges as orographic barriers. The effectiveness of forest soil monitoring networks based on a systematic grid with a predetermined density has not been analysed so far. This study’s analysis was conducted in the Stolowe Mountains National Park (SMNP), SW Poland, using total Pb concentration data obtained from an initial network of 403 circle plots with centroids arranged in a regular 400 × 400 m square grid. The number and distribution of monitoring plots were analysed using geostatistical tools in terms of the accuracy and correctness of soil parameters obtained from spatial distribution imaging. The analysis also aimed at reducing the number of monitoring plots taking into account the economic and logistic aspects of the monitoring investigations in order to improve sampling efficiency in subsequent studies in the SMNP. The concept of the evaluation and modification of the monitoring network presented in this paper is an original solution and included first the reduction and then the extension of plot numbers. Two variants of reduced monitoring networks, constructed using the proposed procedure, allowed us to develop the correct geostatistical models, which were characterised by a slightly worse mean standardised error (MSE) and root mean squared error (RMSE) compared to errors from the original, regular monitoring network. Based on the new geostatistical models, the prediction of Pb concentration in soils in the reduced grids changed the spatial proportions of areas in different pollution classes to a limited extent compared to the original network.

1. Introduction

Soils play a key role in forest ecosystems where they bridge the abiotic and biotic parts of the ecosystems [1,2]. First of all, soils enable plant growth, but they also participate in the cycling and binding of water, carbon and xenobiotics [3,4,5]. The dynamic development of human civilisation in recent centuries has contributed to adverse changes in the soil environment. Increasing industrialisation and urbanisation have contributed to higher pollutant loads entering soils [6,7,8,9]. In addition, the intensification of agricultural and forest management, often handled in an unbalanced manner, leads to the degradation of soils by nutrient leaching, acidification and erosion [10,11].

Wise management of the environment requires a recognition of the dynamic of contemporary threats to ecosystems, including soils. Therefore, monitoring programmes have become standard [12,13,14], and their scale and scope vary in relation to the intensity of natural and anthropogenic transformation, society awareness and ecological policy [2,15,16,17,18]. The results of monitoring studies support the sustainable management of soil resources and programmes for nature conservation, particularly in areas strongly exposed to unfavourable environmental conditions or substantially transformed by humans, as in the case of large-scale monoculture forests [19,20,21,22,23]. Today, Geographic Information System (GIS) technologies, recently combined with remote sensing, support nearly all environmental inventories and monitoring programmes [24,25,26]. GIS is widely applied from the research planning stages, through to the fieldwork, analysis and final presentation of the results and predictions [27]. The intense development of GIS technologies has provided an increasing number of tools to enhance research effectiveness. The final results of GIS application are numerical databases and thematic maps that present the spatial variability and time trends of the analysed phenomena [28,29,30,31].

The unfavourable transformation of environmental quality continues to be a serious problem in Poland. Hence, the National Programme for Environmental Monitoring was implemented in 1992 [32] with forest monitoring as a key module in the monitoring of biological resources [33]. Unfortunately, protected natural areas, which are excluded from normally managed forests, are not involved in the programme to a satisfactory degree. In contrast, it is believed that ongoing observation of environmental changes should be conducted primarily in national parks [3], particularly those located in the mountains that are especially vulnerable to climate change [34] and pollution retained by mountain ranges as orographic barriers [35]. Thus, independent local monitoring networks have been created in the national parks in conjunction with a managed forest monitoring network [36]. Forest monitoring is needed to track environmental changes. Unfortunately, in Poland the effectiveness of existing networks that use regular grids of an arbitrarily decided plot density, has not been analysed with respect to forest soil monitoring. As far as monitoring studies of environmental components are concerned, it is crucial to adopt an appropriate sampling strategy. The overriding aim in selecting an appropriate structure for the monitoring network, in which the sampling will be carried out, is to ensure that the data obtained from it is sufficiently representative [37]. However, this choice is also conditioned by the number of parameters being monitored and their spatial variability, as well as by the efficient use of time and human and financial resources [38]. Besides, a properly designed network structure should include a minimum number of monitoring plots to guarantee representativeness for a given parameter under study while maintaining the assumed confidence level [39]. The correctness of the developed monitoring network is also determined by the adopted sampling design within its structure. There are many such schemes, e.g., judgmental sampling, simple random sampling, systematic sampling, line-transect random sampling, line-intercept random sampling, composite sampling, cluster sampling, stratified sampling, multi-stage stratified random sampling, and conditioned Latin hypercube sampling. They differ in terms of their sampling strategy and efficiency. Among the most commonly used sampling designs, the judgmental sampling method stands out. It is an effective sampling method that uses the researcher’s knowledge, site history, and field observations, but it has the disadvantage of statistical bias. Moreover, this can lead to an uneven distribution of sampling areas, necessitating additional sampling to ensure site coverage. Another method is simple random sampling, which allows statistically unbiased estimates of mean and variability. This technique is most useful in areas with low variability of the analysed feature as well as a low probability of the occurrence of hotspots. The irregular distribution of sampling areas in this method may result in undesirable densification. In turn, one of the more commonly used techniques is systematic sampling. It is carried out in a square, rectangular, or triangular grid system, or transects. It provides statistical impartiality, and compared to simple random sampling, is simpler to describe and less prone to field error, which usually results in better accuracy of the mean. Also, it is used to detect hotspots. Unfortunately, systematic sampling usually gives the same error margin as simple random sampling. A more complex method is stratified sampling, which is mostly used for large areas characterised by a heterogeneous population that can be divided into fairly homogeneous and non-overlapping subgroups. Each layer may have a different sampling design and sample density, in which case they must be analysed separately. The advantage of this sampling method is its high accuracy due to greater precision in estimating the mean and variance, and the ability to calculate reliable estimates for subgroups. Nevertheless, reliable knowledge of the research area is necessary for its correct application. Due to its specificity, this method may require more complex statistical analysis. A particular variation of stratified sampling is conditioned Latin hypercube sampling (cLHS). It is an effective method of sampling from the variability of the feature space of multiple environmental covariates. A sample is drawn from the covariates, so that each variable is maximally stratified. This provides a close representation of the original distribution of the environmental covariates with relatively small sample sizes. Compared to simple random sampling, stratified random sampling and sampling along the principal components, cLHS enables the most accurate reproduction of the original distribution of the environmental co-variates. This method minimises the number of samples while maintaining optimal representativeness [40,41,42].

Monitoring networks that are designed and used in research should be evaluated for possible improvement of their effectiveness [43,44]. For this purpose, classic statistical, as well as geostatistical methods are used. For example, the classical kriging technique supplemented by Moran’s I analysis [45] or extended by Bayesian modelling methods [46,47] is used to optimise the density and distribution of monitoring plots. Kriging is also used to develop new indexes to determine the ratio of sampling efficiency to performance (RSEP) or sampling density–density of soil samples index (DSSI), or to assess prediction accuracy–prediction accuracy index (CEPA) [48]. Bearing in mind the correctness of the evaluation of the monitoring network, the methods used for this purpose should also be subject to review in terms of their effectiveness.

The present work aimed to examine the effectiveness of a regular monitoring network to assess forest soil Pb concentrations to increase sampling efficiency in subsequent research in the SMNP. The analysis applied geostatistical tools (geostatistical analyser and network density analyser) to survey the number and distribution of observation plots in terms of the accuracy and correctness of imaging the soil properties’ spatial distribution. The analysis was conducted using data for the total concentration of lead in soils obtained from the monitoring network of forest soils in the Stolowe Mountains National Park, SW Poland.

2. Materials and Methods

2.1. Research Area

Stolowe Mountains National Park (SMNP) is one of 23 national parks in Poland (and also one of the country’s 7 mountain parks). Covering an area of 63.4 km² in the Sudeten Mountains, SW Poland (Figure 1), SMNP was established to protect the landscape and natural phenomena connected with its unique geological structure [49]. The central part of the area, including the high plateau and hilltops, is constructed of Upper Cretaceous sandstones and mudstones, whereas the south-western part comprises Carboniferous granites and Permian sandstones [50]. The area is situated between 391 and 919 m a.s.l., with an average altitude of 681 m a.s.l. Mean annual air temperature decreases with altitude from 6.5 °C at 500 m a.s.l. to 4 °C at 750–919 m a.s.l. The mean annual precipitation increases with elevation, ranging from 750 to 920 mm. The prevailing winds blow from the west and south-west directions, which also determines the direction of the inflow of polluted air masses [51]. The prevailing soils are Dystric/Eutric Cambisols and Stagnic Luvisols/Alisols, featuring a loamy and silty texture, and Albic Podzols with prevailing sandy texture classes. These dominant soils are locally complemented by Leptosols, Stagnosols, Gleysols, Histosols and Fluvisols [52]. Forest communities cover ca. 89% of the SMNP [53] with still-prevailing spruce (Picea abies) forming the conifer monocultures (typically identified as the Calamagrostio villosae-Piceetum association). The share of beech, fir, sycamore and other species is gradually increasing as a result of the on-going conversion of conifer monocultures to mixed and broadleaved stands [54]. Both the monocultures and their conversion noticeably affect the physicochemical properties of forest soils and humus type as well as the organic carbon stock [10,55].

2.2. Field and Laboratory Research Methodology

The grid of monitoring plots for the SMNP area was developed in 2005–2006 [57], comprising 403 circular plots, 373 of which were located in the forests. Each soil monitoring plot had a diameter of 16 m, and plot centroids were arranged in a regular, 400 × 400 m grid. For further analysis, this network was named the “primary monitoring network” (PMNet).

The first soil inventory in the monitoring network was conducted in the period 2010–2012. Ten replicates (primary samples) were collected in each plot using a stainless-steel frame for the forest litter and a stainless-steel gouge auger for the mineral samples (from the depths of 0–10 and 10–20 cm). Primary samples from the entire plot were mixed, separately for respective layers, and the representative sample for analyses was separated by subdividing. Samples were air-dried, crushed and sieved. The physicochemical analyses were performed in the fine-earth fraction (<2 mm) and involved particle-size distribution, pH in distilled water, soil organic carbon, base cation concentration and the “near-total” (i.e., aqua-regia extractable) concentration of trace elements (Pb, Zn and Cu). A detailed description of the analytical procedures and a general characterisation of the results have been presented in previous papers [10,52].

In the present study, the analysis of the quality of the monitoring network was based on the total concentration of Pb in the 0–10 cm topsoil layer. Pb was selected for analyses because its concentration in the study area was the most variable (among the elements under analysis) and exceeded 100 mg kg⁻¹ in several plots, which signals potential soil contamination according to the legal regulations for nature conservation areas [58]. Pb concentrations were reported in the range of 3.5–219 mg kg⁻¹, with a mean value of 46.6 mg kg⁻¹ (Table 1). The highest Pb concentrations were mainly recorded on the highly elevated surfaces, particularly in the western part of the area, whereas the lowest concentrations were found in the eastern and northern parts, the latter irrespectively of the altitude and parent rocks [52].

2.3. Geostatistical Analyses

The analysis of the quality of the monitoring network was carried out using the geostatistical tools included in the ArcGIS software platform, version 10.5. Regarding the visual presentation of spatial data, a DEM was used. In the geostatistical analyses, the monitoring plots of the primary monitoring network, as well as the modified networks, were represented by their centroids. PMNet was designed to determine the direction and rate of changes occurring in the forest environment [57]. Therefore, bearing in mind that soils are a selected element of this environment, the analyses were performed, attempting to keep its main structure based on a square grid.

In the first step of the analysis, a geostatistical model was developed for the spatial distribution of Pb in the 0–10 cm soil layer in the PMNet grid. The modelling was conducted using a geostatistical analyser and the kriging technique, which provides a linear unbiased estimate for spatial variables [59], according to Equation (1):

(1)$\mathrm{Z}\left({\dot{\mathrm{x}}}_0\right)={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}\mathrm{Z}\left({\mathrm{x}}_{\mathrm{i}}\right)$

where Z(x₀) is the predicted value at location x₀, Z(x_i) is the measured value at location x_i, n is the number of sites within the search neighbourhood used for the estimation, and λ_i is the weighting function.

To ensure that the estimation is unbiased, the weighting values must sum to one, and the estimation errors should be as small as possible, according to Equation (2):

(2)${\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\lambda }}_{\mathrm{i}}=1$

Based on the model obtained, the subsequent steps of PMNet structure analysis were carried out using a network density analyser. The aim was to examine the number of monitoring plots and their distribution in terms of their accuracy and correctness compared to the obtained spatial distribution of Pb content in soils. The analysis steps were as follows:

• 1.. Modification of the PMNet structure by removing monitoring plots that jointly met the following conditions:

  • (a). According to the network density analyser, they showed the lowest probability of Pb concentration > 30 mg kg⁻¹.

  • (b). The geostatistical analyser assigned them with the lowest weights at the estimation stage when generating a prediction map of Pb spatial distribution in the PMNet.

The adopted threshold value of 30 mg kg⁻¹ refers to the geochemical background for Pb concentration in the soils developed from the rocks prevailing in the area under study [60]. The new monitoring network structure obtained in this way was named the “reduced monitoring network” (RMNet).

• 2.. Modification of the RMNet structure by adding new monitoring plots in areas that jointly met the following conditions:

  • (a). In the prediction of spatial distribution for PMNet, the Pb concentrations in these areas exceeded 100 mg kg⁻¹.

  • (b). In the geostatistical model for PMNet, the areas generated the largest errors.

In this case, the adopted threshold value of 100 mg kg⁻¹ refers to the maximum permissible Pb content in uncontaminated soils according to [58]. The monitoring network structure thus created was named the “reduced and extended monitoring network” (REMNet).

• 3.. Development of new geostatistical models for Pb concentration in the RMNet and REMNet grids and performing predictions of Pb spatial distribution for these grids. Due to time constraints, it was not possible to return to the field and perform additional analysis at new sampling points. Therefore, Pb concentrations from the PMNet prediction map were assigned to the new monitoring plots during the development of the geostatistical model for REMNet.

Optimisation of all developed models was carried out in a geostatistical analyser using manual adjustment of their parameters. For each geostatistical model, the curve describing it was adjusted against the experimental semivariogram by selecting optimal values for its sill, range, and nugget. The lag size was determined using the average nearest neighbour module included in the spatial statistics analyser of the ArcGIS platform. The correctness and accuracy of the geostatistical models were verified using the cross-validation method, which evaluates the “consistency” of the estimations performed with the measurement results. Minimisation of the burden of the mean standardised error (MSE) and root mean squared error (RMSE) was adopted as the main criterion for the optimisation of model parameters [45,61].

3. Results

3.1. Primary Monitoring Network (PMNet)

The PMNet grid comprised 403 monitoring plots arranged in a regular square grid. Due to the lack of technical feasibility of soil sampling in certain areas, only 387 plots were included in the analysis (Figure 2). The geostatistical model for PMNet was developed using the simple kriging method, for which the spatial dependence was best described by a stable curve algorithm. The autocorrelation range for Pb concentration exceeded 1000 m. The nugget value was 0.101, which may suggest the presence of small-scale spatial variations. In turn, the semi-variance values for PMNet indicate the highest spatial autocorrelation as compared to the other networks (RMNet, REMNet).

The autocorrelation strength was also indicated by the nugget-to-sill ratio (RNS index), which determines the ratio of the nugget value to the sill value of semi-variance. RNS is computed as (C₀/(C₀ + C)), where C₀ is the nugget value, and C is the sill value of semi-variance.

As reported by [62,63], this ratio shows strong, moderate and weak spatial dependence in the value ranges <0.25, 0.25−0.75 and >0.75, respectively. In the case of the PMNet network, it was 0.36 (Table 2), which indicates a moderate spatial autocorrelation of Pb.

In assessing the quality of the model by cross-validation, the value of the error rate MSE was found to be close to zero and the RMSE was close to 1 (Table 3), which according to [64] attested to the correctness of the model. The prediction of the spatial distribution of Pb concentration in the 0–10 cm soil layer, developed on the basis of the geostatistical model for PMNet, is presented in Figure 3.

According to the prediction, the highest Pb values (>100 mg kg⁻¹) were recorded in the central and western parts of the SMNP. They occurred mainly in the highest elevated areas of the mountains (i.e., in the region of Mt. Skalniak and Mt. Szczeliniec Wielki).

On the other hand, most locations showed Pb concentrations close to the median (Table 1). These areas covered mainly the eastern part (Batorowo area) as well as the northern part, but they also appeared in the central region of the area. In the eastern region, Pb concentrations did not exceed 10 mg kg⁻¹.

3.2. Reduced Monitoring Network (RMNet)

Based on the criteria defined in Section 2.3, the PMNet structure was reduced by 30 monitoring plots. Finally, the total number of plots in RMNet was 357. Most plots were removed from the northern and eastern region (Figure 4).

As in the case of PMNet, the geostatistical model for RMNet was developed using the simple kriging method. Spatial autocorrelation was best represented by a stable curve algorithm. The autocorrelation range was slightly larger than that in the model for PMNet and was more than 1160 m, but the values of nugget and partial sill of semivariance increased. Removal of the selected monitoring plots resulted in an increase in the mean distance between them and a decrease in spatial autocorrelation as against the PMNet. Similar to the PMNet model, the RNS value indicated, according to [62,63], a moderate spatial autocorrelation of Pb concentration (Table 2). Based on the MSE and RMSE errors, it was found that despite the reduction in the number of monitoring plots, the model was correctly executed (Table 3).

The prediction of spatial Pb distribution was similar to that developed for PMNet (Figure 5). There was a decrease in the range of surfaces with a predicted lead concentration of 10–20 mg kg⁻¹ (mainly in the eastern and northern regions) and >125 mg kg⁻¹ (Mt. Szczeliniec Wielki), whereas the areas with a predicted Pb concentration of 25–50 mg kg⁻¹, close to the mean Pb concentration, slightly increased (Figure 5).

3.3. Reduced and Extended Monitoring Network (REMNet)

The modification of the RMNet structure according to the assumptions set out in Section 2.3 resulted in the addition of 20 new monitoring plots. Their number and location were determined by the network density analyser, taking into account an additional criterion that excluded the overlapping of the plots. The distribution of the new plots was characterised by irregularity. The largest number of new plots was located in the highest elevations (i.e., in the vicinity of Mt. Skalnik and Mt. Szczeliniec Wielki). The REMNet grid created in this manner finally comprised 377 plots (Figure 6) (i.e., 10 less than the PMNet). The simple kriging method was also used for modelling, but in this case, the spatial autocorrelation was best described by the Gaussian curve algorithm. The autocorrelation range and the values of the nugget and partial sill were similar to the model made for the RMNet. As with the RMNet, modifying the network structure resulted in an increase in the mean distance between monitoring plots and a decrease in spatial autocorrelation. Spatial autocorrelation was noticeably lower compared to the PMNet and slightly higher than the RMNet.

Similar to the models for PMNet and RMNet grids, the value of the RNS index (Table 2) indicated a moderate spatial autocorrelation of Pb concentration [62,63]. The MSE and RMSE errors (Table 3) indicated correct modelling [64].

The prediction of the spatial Pb distribution for REMNet increased the share of plots with Pb concentrations ranging from 25 to 50 mg kg⁻¹ in the central part of the area and decreased the share of areas with contents ranging from 10 to 25 mg kg⁻¹ in the eastern part (Figure 7). The prediction of the spatial Pb distribution did not bring significant changes compared to the map developed for RMNet (Figure 5); however, it came closer again to the spatial variability reflected within the PMNet (Figure 3).

4. Discussion

The proposed concept of the evaluation and modification of a monitoring network is an original solution. As mentioned in the introduction, the procedures presented by other authors are different and based on individual solutions, however, they have a similar objective, which is the optimisation of monitoring networks used in environmental studies in terms of correct imaging of changes that may occur in broader temporal and spatial perspectives [65,66].

The geostatistical analyses conducted in this study indicated that the initial PMNet structure allowed for the effective identification of Pb distribution in SMNP soils. The spatial autocorrelation range for this element exceeded 1000 m, despite a large number of factors potentially affecting the soil variability in mountainous areas. This finding could be explained by the poor mobility of Pb in soils or its strong affinity to organic matter [4,9,16], leading to its accumulation in the topsoil layers, particularly in the forest soils. The geostatistical model for PMNet, applying the cross-validation evaluation, developed in this study with respect to RMNet and REMNet models, showed mean standardised error and root mean squared error values closest to the optimal ones (Table 3). However, it did not reach them, most likely due to the distribution of monitoring plots in a rigid square grid, where one plot represents 0.16 km². As indicated by other authors [15,65,67], the representativeness of such large plots may be questionable in regions of high relief variability.

The main idea of the new monitoring network structures was an optimisation oriented towards minimising the number of monitoring plots while maintaining the correctness and accuracy of soil parameter imaging. Achieving this goal using geostatistical methods may translate into better economic performance and shorter time of study [68]. Both network modernisation concepts (RMNet and REMNet) meet the adopted requirements. In both cases, the number of monitoring plots was reduced (RMNet by 30, REMNet by 10), and the geostatistical models developed for them were characterised by fine correctness. The MSE and RMSE errors are close to the optimal values reported by other authors, but slightly worse compared to the parameters obtained for the initial PMNet (Table 3). The higher values of the partial sill and nugget for RMNet and REMNet are due to higher values of semivariance in the spatial autocorrelation models, potentially as a result of the change in the structure of these networks caused by reducing the number of monitoring plots under analysis, which led to an increase in the average distance between them [67,69]. The increase in the value of error(s) as far as the REMNet is concerned, may also be due to the fact that Pb contents in the points added at the last stage were not derived from the new sampling and analysis, but rather, were determined based on the original prediction map. It cannot be excluded that samples from the added points would have revealed higher local contamination, which would have affected the prediction.

The spatial variability of Pb concentration imaged in RMNet and REMNet, despite minor differences, reflected spatial dependencies that compared well to the PMNet. In the map developed for the REMNet, surfaces with high Pb concentration (>100 mg kg⁻¹) are enlarged, suggesting poorer representativeness of some monitoring plots in the PMNet structure. Therefore, the use of the REMNet in future monitoring studies in this study area may be more beneficial for the recognition of the spatial variability of Pb concentration, but will slightly reduce workloads. An unavoidable disadvantage of this proposal is the necessity to designate new monitoring plots outside the original locations. This requirement will reduce the possibility of comparing the results in successive monitoring studies, and thus, may weaken the statistical significance of possible time-related trends.

It is worth noting that one of the alternative solutions for optimising the existing network structure with the use of geostatistical tools could be redesigning it from scratch. As far as the modelling process is concerned, the application of the kriging method makes it possible to recognise the spatial variability structure and to use it to produce the sampling design [70]. As indicated by [71], the maximum standard error of a kriged estimate is a reasonable measure of the goodness of the sampling design. In this case, using a regular triangular grid as the sampling scheme usually minimises the maximum standard error. In turn, for each assumed maximum standard error, the optimal density can be determined based on the semivariogram for the variable. However, it must be emphasised that the concept for the evaluation of the existing monitoring network presented was not intended to introduce a completely new network structure. It aimed to optimise the network with as little interference as possible of its original regular grid structure. As mentioned in Section 2.3, the network was also designed to observe other environmental components (tree stands, vegetation, lichens, etc.). Thus, the introduction of a completely new network structure just for soils would hinder the implementation of monitoring as well as the linking and interpretation of data related to all of the components assessed.

Moreover, only one parameter relating to soil contamination was applied in the current study. As far as the optimisation of the monitoring network is concerned, it is necessary to analyse at least a few key parameters relating, for example, to macroelement concentration and organic carbon pools, which are relevant for semi-natural forest reconstruction and climate-change mitigation policies, respectively. It will only be possible to present the economic calculation for the proposed method, and thus, determine the validity of its application, after performing the analyses for a larger number of parameters.

It cannot be ruled out that geostatistical analyses will provide not one, but different variants of the monitoring network structure, the merging of which will not be substantively justified. In such an event, it would be more reasonable to use distinct modifications of the network structure, separately for selected key parameters, developed according to the methodology presented for the creation of the RMNet structure. It is possible that in the case of the more mobile soil parameters, their spatial distribution might strongly vary over time, thus reducing the potential of the presented approaches (RMNet and REMNet). The RMNet approach will require more field and laboratory work but this will yield more reliable analysis results.

5. Conclusions

The use of geostatistical tools to optimise the structure of the original mountain forest soil monitoring network, which is established based on a regular square grid, resulted in a reduction in the total number of monitoring plots. Different variants of reduced monitoring networks, developed using the proposed procedure, make it possible to build correct geostatistical models characterised by slightly worse MSE and RMSE errors versus the original, regular monitoring network. The prediction of Pb concentration in soils based on the new geostatistical models changed the spatial proportions of areas in different pollution classes to a limited extent compared to the original regular network. Although the application of the geostatistical tools gave positive results, the final proposal for the optimisation of the monitoring network should include additional analyses performed using a similar method for other key soil parameters, particularly macronutrients and organic carbon pools, which are crucial for the achievement of monitoring objectives as well as for forest management and protection policies in the national park.

Author Contributions

All authors contributed to the development of the idea and authoring of the paper: Conceptualisation P.J. and C.K.; methodology P.J.; writing—original draft P.J.; writing—review and editing C.K. All authors have read and agreed to the published version of the manuscript.

Funding

This work was co-financed by the Ministry of Higher Education and Science of Poland, project MNiSW N R09 0029 04/2008 and a statute project of the Wrocław University of Environmental and Life Sciences, Institute of Soil Science and Environmental Protection from a subsidy of the Ministry of Education and Science of Poland.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The study did not report any data.

Acknowledgments

The authors wish to thank A. Bogacz, B. Gałka, B. Łabaz and J. Waroszewski, who were included in the development of the soil studies in the Stolowe Mountains National Park.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Figure A1

Experimental semivariogram (symbols) and models (line) for total concentration of Pb in the topsoil layer—PMNet.

[Figure omitted. See PDF]

Figure A2

Experimental semivariogram (symbols) and models (line) for total concentration of Pb in the topsoil layer—RMNet.

[Figure omitted. See PDF]

Figure A3

Experimental semivariogram (symbols) and models (line) for total concentration of Pb in the topsoil layer—REMNet.

[Figure omitted. See PDF]

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Location of the study area on a digital elevation model (DEM) map. DEM data acquired from light detection and ranging (LIDAR) [56].

Figure 2. Location of the 387 plots of the primary monitoring network (PMNet).

Figure 3. Prediction map of the spatial Pb distribution, developed by simple kriging for the primary monitoring network (PMNet).

Figure 4. Location of 357 plots (red circles) in the reduced monitoring network (RMNet). Yellow circles represent removed monitoring plots.

Figure 5. Prediction map of the spatial Pb distribution, developed by simple kriging for the reduced monitoring network (RMNet).

Figure 6. Location of 377 plots (red and blue circles) in the reduced and extended monitoring network (REMNet). Blue circles represent added monitoring plots.

Figure 7. Prediction map of the spatial Pb distribution, developed by simply kriging for the reduced and extended monitoring network (REMNet).

Table 1

Statistical parameters characterizing variability of total Pb concentration in the soils of SMNP (data source [52]).

N	Minimum	Maximum	Mean	Median	Standard Deviation
	mg kg−1
387 *	3.5	219	46.6	41.0	28.6

* Due to the lack of technical feasibility of soil sampling in certain areas, only 387 of the 403 plots were included in the analysis.

Table 2

Semivariogram model parameters for analysed monitoring networks.

Monitoring Network	Model *	Number of Samples	Nugget	Partial Sill	Range, m	RNS
			(C0)	(C0 + C)	(A0)	C0/(C0 + C)
PMNet	stable	387	0.101	0.182	1002	0.36
RMNet	stable	357	0.301	0.621	1162	0.33
REMNet	gaussian	377	0.362	0.510	1130	0.42

* Graphs of semivariograms in Appendix A.

Table 3

Model errors for analysed monitoring networks.

Monitoring Network	MSE	RMSE
	Mean Standardized Error	Root Mean Squared Error
PMNet	−0.014	1.009
RMNet	−0.025	1.000
REMNet	0.023	1.011