Analysis of Spatial Variability of Selected Soil Properties in the Hard Coal Post-Mining Area Sławomir

The result of geomechanical and hydrological degradation caused by the underground extraction of hard coal are the transformations of the environmentally managed terrains and high variability of soil properties occurring in these areas. The analysis of selected soil properties of spatial variability in the post-mining area conducted by means of the kriging technique was presented in the paper. The determined points of empirical semivariogram were described by means of mathematical functions and theoretical semivariograms were plotted. The conducted analysis allowed to plot the maps of variability isolines, which my provide a basis to delineate the boundaries of areas most susceptible to the hydrological degradation of soils. The obtained research results indicate that the applied kriging technique may prove to be a useful tool for determining spatial variability of soils in the areas of hard coal mining operations and allow to delineate the boundaries of the areas most susceptible to soil degradation.


INTRODUCTION
Underground extraction of hard coal causes transformations of the environmentally managed areas.Soils are one of the biotope elements that are the most susceptible to degradation [Klatka et al. 2015].In effect of the extraction activities of mines, the geomechanical land subsidence occurs, leading to the raise of groundwater table and hydrological soil degradation [Kaszowska 2007].A high level of soil moisture content worsens the conditions of plant vegetation and intensifies the expansion of hydrophylic vegetation.In arable lands, a slightly raised groundwater table may cause an increase in agricultural production, but also soil degradation, lowering of the soil quality classes and soil productivity.A necessary change of land use from arable lands to meadow and pasture use is also frequently observed [Strzyszcz 1995, Rosik-Dulewska et at.1999].The extent of damage arising on arable lands and woodlands is clearly spatially diversified, which apart from theanthropgenic pressure, is also affected by the soil properties themselves.As a result of these degradation phenomena, the properties of soils in the areas of extractive industry operations usually reveal a considerable spatial variability.Assessment of these properties should be conducted taking into consideration the recognized character of their variability.Among the methods providing a better estimation of variability character is geostatic kriging technique [Webster and Oliver 1990].The application of mathematical-statistical procedures, as well as determining and analysis of the functions of empirical and theoretical semivariances allow to estimate the value of variable at each point of the investigated area and to plot isoline maps of spatial variability [Marx and Thompson 1987].On the basis of these maps it is possible, among others, to determine the areas prone to various forms of surface degradation and demarcate the boundaries of the most threatened areas [Klatka and Boroń 2008].The method may be also important for the optimization of lands use and increasing their productivity [Sigua and Hudnall 2008].
The aim of the paper was an analysis of the spatial variability of selected soil properties in the 0-25 cm layer in the area subject to degradation due to the underground mining operations.

MATERIALS AND METHODS
The research was conducted in the northeastern part of hard coal mining area, "Ruch-Borynia" coal mine in Jastrzębie-Zdrój within Świerklańska, Zamkowa and Plebiscytowa streets.The geomechanical and hydrological degradation of soils, which occurs in this area is caused by land subsidence between 1 and 2.5 m, in relation to the original position.The investigations were conducted on the basis of a regular measurement grid with 100m side, determined with GPS method.The soil material for laboratory tests was collected at the stabilized measurement points from the depth of 0-25 cm.A total of 30 soil samples were collected in three replications.The location of measurement points was shown in Figure 1.The laboratory tests comprised determining the following soil properties: the particle size distribution, bulk density, the specific gravity, total porosity, organic matter content and soil saturated conductivity coefficient.The particle size distribution was determined by means of aerometric Casagrande's method in Prószyński's modification.Groups and sub-groups of particle size distribution were determined on the basis of PN-R-04033 standard.The physical soil properties (bulk density and total porosity) were determined with Kopecky's method (100 cm -3 cylinders) [Mocek et al. 1997], specific gravity was established by pictometry in distilled water, whereas the soil saturated conductivity coefficients were determined in the laboratory using pressure drop method in the apparatus based on Darcy's law with regulated water pressure and electronic reading of water volume.Determination of saturated conductivity coefficient on this apparatus is based on a linear dependence of the flow on hydraulic gradient at so-called laminar motion, in compliance with Darcy's law [Baver et al. 1972].The organic matter content was determined by means of Tiurin's method in Oleksynowa's modification [Oleksynowa et al. 1987], which relies on humus oxidation with potassium dichromate (Cr 6+ ).

Fig. 1. Localisation of measurement points
The following statistical measures were computed in order to characterize the analyzed lithological parameters and to initially determine the degree of their variability [Hellwig 1993]: maximum and minimum values, median, average value, standard deviation and variability coefficient.Variability coefficient was determined as a variance dispersion measure and standard deviation in the form of a variability quotient of a given feature -standard deviation and average value of a given feature expressed in percent.The results were presented in Table 1.The spatial variability of the studied soil properties was determined with kriging method.The method allows for the assessment of the confidence intervals of statistical estimation, determining the average value in any part of the analyzed area and finding localization for the measurement points, which to the highest extent would narrow the confidence intervals of distributions.Precisely speaking, kriging should be defined as the method optimizing the estimation of spatially correlated value Z, both in the case of its stationarity and in non-stationary cases.If we assume that Z 1 = Z (x 1 ) is the value measured at point x 1 (i = 1,2 ……n), then the problem of point estimation involves determining value Z o at point x o , where value Z was not measured.Through constant change of point x o location, it is possible to determine the entire field of variation of the Z variable [Webster and Oliver 1990].Beside the possibility to estimate the value of Z at each point, kriging technique also allows to: assess the confidence intervals for estimation, determine the mean magnitude of the value Z at any part of the researched area, as well as find the locations for new measurement points, which to the highest extent will narrow the confidence intervals for the calculated distributions [Somorowski 1993].The analysis with the use of kriging method was conducted in two stages.At the first stage, the γ function was determined, in geostatistics referred to as a variogram or semivariogram [Smith et al. 1994], which may defined as: From the semivariogram definition it follows that γ(0) = 0 and γ(h) ≥ 0 for h > 0.
In the case when the assumptions of stationarity are met, the semivariogram contains all the information about spatial variability of the Z value.The assumption of internal stationarity allowed to calculate the value of empirical semivariograms on the basis of a single realization of the process (one measurement series) [Oliver and Webster 1986] as: suggest the following models: spherical, exponential, Gaussian, linear and logarithmic; however, the most commonly used functions belong to the class of so called safe models, i.e. linear and spherical functions.The obtained shape of a semivariogram indicates the kind of spatial correlation of the analyzed variable.There are two basic kinds of semivariograms.The first is characterized by the fact that the value of variance is increasing with growing separation distance of h to some maximum value at which it remains constant with increasing distance.After the initial growth, value γ (h) reaches, for certain separation distance α the value equal to a variable variance.
The separating distance α is called the variogram range and determines the limit of spatial correlation.The Z variable which has the variogram of this kind is not only internally stationary, but itself has the second-order stationarity.In the second type of semivariograms, the variance γ (h) is the function increasing without limits.In this case, the variance is infinite and with growing separation distance h, reaches a certain threshold value determined as C. The threshold value is reached precisely or only asymptotically.If the semivariogram does not pass through the middle of the coordinate system origin, we speak about so called nugget effect, determined as c o .The effect may evidence the fact that if a nugget is collected for sampling, the other samples, even very close ones will differ with their concentrations of a given component [Smith et al. 1994].
For graphic representation of the spatial variability, each of the measuring points was ascribed the position vectors.

RESULTS
According to the Polish Soil Classification [PTG 2011] the investigated soils belong to brown soil order, Eutric Cambisols type, and Eutric Gleyic Cambisols sub-type [Klatka et al. 2015].These are soils revealing the properties of Cambisols, but differing with a higher water content in the profile.The analysis of the obtained results of laboratory tests indicates the occurrence of various sand sub-groups and only in some places lime patches.Therefore, these are light soils, with less than 20% content of floatable particles, high hydraulic conductivity and low water retention.The soils with this mechanical composition are regarded as slightly resistant to degradation processes due to hard coal mining operations [Klatka and Boroń 2008].The specific gravity ranged between 2.4 and 2.66 g•cm -3 and was approximate to the values most frequently registered in the soils of Poland [Zawadzki 1999].The values of bulk density were on average 1.30 g•cm - 3 and were also approximate to the average values for soils, which according to Zawadzki [1999] range between 0.75 and 1.90 g•cm -3 .The total porosity value fluctuated from 44.4 to 52.32%, on average 49.22%.Because total porosity affects, among others, the aerial properties and hydraulic conductivity in soils, it should be stated that at a majority of the investigated points it reached the values favourable for soils and may have influenced the dimension of degradation processes.According to Musierowicz [Zawadzki 1999], in cambisols (brown soils), the humus content in the 0-25 cm layer is between 1.5 and 2.5%.In many cases, the content of organic matter noted in the investigated area exceed these values for brown soils; therefore.in the future reclamation activities one should take into consideration a possible removal of the humus horizon.The soil saturated conductivity coefficients determined with a laboratory method fluctuated from 0.00064 to 0.01287 m•d -1 .The highest values were registered for loose sands.The value of saturated conductivity coefficient plays very important role in hydrological degradation of soils in the areas of mining damage [Klatka et al. 2015, Klatka et al. 2016].The obtained results on an average level of 0.0618 m•d -1 together with a light particle size distribution indicate the studied susceptibility of soils to hydrological degradation.
On the basis of computed variability coefficient, it may be stated than in the 0-25cm soil layer in the investigated area, the most variable spatially was saturated conductivity coefficient, whereas silt and clay fractions content and organic matter content were moderately variable.The lowest variability was noted for bulk density, specific gravity and total porosity.The obtained results are approximate to the findings of other authors [Warrick and Nielsen 1980, Brandyk et al. 1996, Klatka et al. 2015, Kruk et al. 2018].The kriging technique was used to determine the points of theoretical semivariance for each of the analyzed soil properties.For further analysis using kriging technique, only 15 points of the empirical semivariants were used, because the subsequent ones were characterized by the separation distance greater than the distance between the coordinate system origin and the last measurement point.
Most of the plotted empirical semivariograms were smoothed using linear and spherical functions (Fig. 2).The exceptions were the values of theoretical semivariances obtained for the saturated conductivity coefficient, bulk density and specific gravity, which proved impossible to smooth using the above-mentioned models.Generally, for these properties the values of semivariance formed a parallel line to the distance axis.The fact limits the existence of spatial correlation, the limit of which may be lower than the separation distance between the measuring points.The empirical sevariogram plotted for silt fraction content and total porosity was smoothed by means of a linear model.In this case, the semivariance grows along with the separation distance and it is impossible to determine the spatial correlation limit as the value for which the semivariance stabilizes itself.On the other hand, the nugget effect was determined.The empirical semivariograms plotted for the contents of sand and clay fraction were smoothed by the spherical model.In these cases, the nugget effect was observed, evidencing a variability of a given value at the separation distance lower than the distance between the measuring points.The occurrence of this effect may result from the precision of applied measurement method.The limit of spatial correlation, 89.60m and 150.65m, respectively, was also determined for the semivariograms of the discussed soil property.The parameters of models fitted to theoretical variograms were compiled in Table 2.
Using the computed empirical semivariance models, the spatial variability maps were plotted for the studied parameters.Isoline cutting was matched to the random error magnitude -nugget variance.A high value of the variance affected a better smoothing of the spatial distribution.
Many authors investigated the problem of spatial variability of mineral soil properties using geostatic kriging method.L'Opez-Granados et al. [2002] investigated spatial variability of selected chemical properties of soils in two regions of southern Spain in view of selecting a proper kind and doses of fertilizer.Park and Vlek [2002], who studied the spatial variability of soil properties in 64 soil profiles on a slope in Bicknoller Combe, Somerset, UK revealed that recognizing the soil spatial variability gains in importance in ecological modeling.In their investigations, they also characterized the complexity of three-dimensional variances of individual soil properties and researched the possibility of predicting the soil properties distribution using three different regression methods.Wei et al. [2006] used kriging method for the assessment of spatial variability of organic matter content in the soils of northern and eastern parts of China.Stach [1998] determined the spatial variability of the soil arable layer on lithologically heterogenous morain slope.In the areas of reclaimed settling ponds holding soda waste from KZS "Solvay" in Krakow, the problem of variability of the selected sol properties in the insulation layer was investigated by Klatka et al. [2017].The results of the research conducted by the authors mentioned above indicate that the analysis of spatial variability of the studied soil properties is a useful tool allowing to determine the average value of any property in any part of the studied area and to find the location for new measuring points.In their research conducted in the area of the post-mining area at Szczygłowice, Klatka and Boroń [2008] demonstrated that the kriging method and plotted variability isolines may provide a basis for demarcating the boundaries of areas the most susceptible to soil hydrological degradation.
Presently, the agricultural lands in the investigated area, which are the private property of farmers, are characterized by a low level of  applied agrotechnological measures, which undoubtedly affects both their quality and agricultural suitability.In future, one should reckon with an increasing wasteland area accompanied by intensifying processes of soil hydrological degradation.According to Siuta [2007], the most efficient way to counteract ground surface waterlogging and increase rainwater retention in soils, as well as optimizing utilisation of groundwater in soil, is conservation of the existing irrigation networks or construction of a newone.The lack of this type of reclamation actions in a given area may also prove disastrous for some building areas.The threat involves flooding the house foundations and cellars.

CONCLUSIONS
1.The investigated area undergoes intensive effects of the underground hard coal extraction.The geomechanical transformations of the surface caused by the exploitation conducted using so called "roof fall" method led to soil degradation and raised the groundwater table, as well as caused waterlogging of grounds and soils.Apart from the mining activities, the dimension and intensity of this type of degradation were also influenced by the soil properties.Spatial variability of the studied soil properties reveals exceptionally high dynamics, which will intensify in time along with intensification of degradation processes.
2. The conducted geostatical analysis of the examined soil deposits did not determine any theoretical semivariogram model for the saturated conductivity coefficient, bulk density, specific gravity or spatial correlation limit.
It is connected with a considerable random variability of these properties.For silt fraction content and total porosity semivariance was increasing linearly along with the separation distance, and like for the former soil properties, it was impossible to determine their spatial correlation limit.For the content of sand fraction and clay fraction, the limits of spatial correlation were determined.The occurrence of the nugget effect was also observed, which may evidence a variability of a given parameter along a separation distance smaller than the distance between the measurement points.
3. Computations of the statistics for the discussed area should also comprise the character of spatial variability of the investigated soil properties.The variability isolines plotted by means of kriging indicate a relationship of the investigated properties with the parent rock lithology.The conducted research and geostatic computations allow for the statement that the kriging method may prove to be be a tool most useful for determining spatial variability of soils in the areas of hard coal mining activities.However, the reliability of this method is conditioned to a great extent by the form of empirical semivariogram and theoretical variogram fitted on this basis.

Table 1 .
Statistical measures of the soil materials determined in the 0-25 cm layer 2 (ℎ) =1 where: γ (h) -semivariogram value for a separation distance of h, m (h) -number of observations made at a separation distance of h, Z (x i ) -value observed at point x i, Z (x i +h) -value observed at point located at a separation distance of h

Table 2 .
Parameters of the models fitted to empirical semivariances