Progress in the Reconstruction of Terrain Relief Before Extraction of Rock Materials—The Case of Liban Quarry, Poland

Abstract

Open pit mining leads to irreversible changes in topographical relief, which makes a return to the original morphology virtually impossible. This is important for quarries that were part of former mining areas. This research presents an innovative approach to the reconstruction of the relief of anthropogenically transformed land on the example of Liban Quarry in Cracow, where operations began before 1873 to 1986. The basis for the reconstructed area was a Topographic Map of Poland with a scale 1:10,000 from 1997, from which a set of data was obtained to perform spatial analyses. The estimation was conducted using the ordinary kriging method, enabling a reconstruction of the morphology of the studied area and presenting it in the form of a hypsometric map and a digital elevation model. The correctness of the modelling was verified by cross-validation and a kriging standard deviation map (SD_OK). These revealed low values of estimation errors in the places without contour lines on the base map. The comparison of the obtained maps and model with a Tactical Map of Poland with a scale 1:100,000 from 1934 indicated great similarities. The highest interpolation error value was recorded in the part of the pit where the difference between the actual and reconstructed elevation was about 30 m on average. In the exploited part, the SD_OK did not exceed 0.52 m, and in the entire studied area, it reached a maximum of 0.56 m. The proposed approach fulfilled the assumptions of reconstruction, as the analysis revealed elements matching the historic relief in both forms of presentation of the topography of the quarry, on the obtained hypsometric map and on the tactical map. Our study is among the very few in the world concerning the application of geostatistics in the restoration of the relief of land transformed by open pit mining activities.

1. Introduction

Open pit mining contributes to irreversible changes in the relief of land. Hydrographic conditions are altered, and soil and vegetation degradation processes occur [1,2,3,4,5]. For this reason, a return to the original morphology becomes practically impossible. This is highly important in the case of anthropogenic areas, in particular, former mining areas, including quarries established in large numbers until the mid-20th century [6,7,8]. Remote sensing techniques were not developed at that time to the extent that would allow for accurate archiving of data pertaining to the original relief, i.e., before the start of the excavation of the rock material. The only materials available were maps, mainly topographical maps with various scales, usually too small to be able to precisely reproduce the former relief of the land affected by exploitation.

The literature describes many examples of locations that have been transformed by mining activity [1,9,10,11,12,13,14,15,16,17,18,19]. However, the number of studies referring to different methods of reconstructing the relief of areas transformed by mining seems to be insufficient [7,8,20,21,22]. The approaches taken so far in this scope mainly focus on the use of photogrammetric methods based on archival satellite images, or aerial and UAV images [7,8,22], as well as the vectorisation of old topographic maps [20]. Such studies are extremely useful in geomorphological analyses, in the assessment of human impact on the environment, and thus in planning the reclamation of anthropogenic areas, as well as in assessing the value of future investment projects [8,20]. However, none of the proposed solutions have so far applied geostatistical tools to reconstruct a relief transformed by anthropogenic pressure.

In Poland, a large group of terrain formations created as a result of the extraction of rock materials occurs in Cracow. These include the inactive Jurassic limestone quarries. The changes in the land relief have not been analysed so far. This is especially true of Liban Quarry, which is part of the Krzemionki mining area, occupying a fragment of Krzemionki Hills. Only a few cartographic materials refer to the morphology of Liban Quarry and its surroundings. These include the Geomorphologic Map of Cracow with a scale of 1:50,000 [23] and the Tactical Map of Poland with a scale of 1:100,000 (1934) [24]. However, the changes in the relief of the section of Krzemionki Hills, resulting from the creation of an extensive mining working—Liban Quarry—have not been studied so far. There are also no studies related to the digital visualisation of the original relief of the quarry in question, i.e., before the commencement of rock extraction.

Therefore, the aim of the study was to reconstruct the relief of the area currently occupied by the former Liban Quarry from the period before the start of the mining operations (before 1873), using geostatistical tools.

2. Data and Methods

2.1. Study Area and Object of the Study

The former Liban Quarry (50°02′10.37″N; 19°57′24.26″E) with the surface area of 13.68 ha, is located in the district of Krakow, Podgórze, near the Krakus Mound from the 7th century (269 m a.s.l.), which is the culmination of the Krzemionki Hills [25] (Figure 1). Together with the Krzemionki mining site, it covers a surface area of 21.11 ha. The size of the entire analysed area, including the surroundings of Krzemionki, amounts to 48.38 ha in total.

The beginning of limestone mining in Liban Quarry probably dates back to the 14th century (northern part) and was carried out on a larger scale starting around 1873 (southern part) [26]. The exploited deposit consisted of thick-bedded Jurassic rocky limestone, with both vertical and horizontal cracking, without clear bedding. The deposit was also recorded to contain diversiform flint nodules and fossils. Mining operations were discontinued on 30 June 1986 [27].

The orientation of the main axis of the working runs from N to S. Liban Quarry is surrounded by vertical walls that, with rock shelves, reach a maximum height of 40 m (near the Krakus Mound). Some debris cones have formed at their feet. The floor of the working is a rock bed with shallow cavities where water is retained [28].

2.2. Cartographic Data

The analysis of the source material concerning Liban Quarry included archival documentation and archival and contemporary cartographic materials of the studied area, with various scales [29,30,31,32,33,34,35,36]. A review of historic maps revealed that the relief of the analysed area was presented on them using the old map drawing technique (hachuring)—hand drawn lines in the direction of the steepest slope that represent a relief [37]. Therefore, these maps did not contain data on the absolute elevations and could not be used to determine the denivelation within the studied area and thus visualise its relief. For this reason, the Tactical Map of Poland with a scale of 1:100,000 (1934) [24] was used, as its topographical basis with a scale of 1:25,000 dates back to before 1873, i.e., to the period before the start of mining operations in the area under analysis. The absolute elevations charted on this map indicate unequivocally that the relief of the examined section of Krzemionki Hills on both sides of the Krakus Mound (271.3 m a.s.l.) (southern and western) remained unchanged (Figure 2a). The insufficient scale of the map prevented obtaining the precise spatial data required for a geostatistical modelling of the historic relief occupied by the now defunct Liban Quarry. It was, however, used as a reference standard to verify the obtained map and model using the ordinary kriging method (OK). Eventually, a contemporary Topographic Map of Poland with a scale of 1:10,000, sheet Kraków–Wola Duchacka (1997) [25] was used as a basis, and an attempt was made to reconstruct and digitally visualise the terrain relief from the second half of the 19th century that is occupied by the now defunct Liban Quarry.

The border of the mining area was delineated as based on the map of the Krzemionki mining area with a scale of 1:1000 [27]. This included two former borders of the area dated 15 November 1960 and 7 January 1987, within which the analysed working is located. At the same time, the border of Liban Quarry was delineated for the purposes of the study by using available sources and cartographic materials. These included Zarychta [28], an orthophotomap (Figure 1), a topographic map (Figure 2b) and a digital elevation model (DEM) (Figure 3a,b). Their analysis revealed small discrepancies in the course of the quarry border that result from the different scales of the materials used.

In order to obtain information on the present relief of the studied mining area in addition to field observations, a digital elevation model (DEM) analysis was implemented. This was developed on the basis of detailed airborne laser scanning that are recorded by the device, named Light Detection and Ranging (LiDAR) spatial data with a resolution of 1 m obtained from airborne laser scanning which was then used to assess the relief on a small scale. Furthermore, a comparative analysis of the current and reconstructed relief was undertaken to evaluate the usefulness of the modelling carried out in order to obtain a visualisation presenting the morphology of the studied area as transformed under the influence of mining works (Section 3).

2.3. Exploratory Spatial Data and Semivariogram Analysis

Geostatistical modelling was preceded by the analysis of descriptive statistics and by generating a histogram with a QQ plot to verify their distribution. A Box–Cox transformation was then carried out to provide a more constant variance over the entire analysed area. For the transformed dataset, descriptive statistics were recalculated in order to compare them with the original dataset.

The presence of spatial autocorrelation in the obtained dataset was examined on the basis of the generated lagged scatterplots for all pairs of $z\left(u_{\alpha }\right)$ and $z\left(u_{\alpha }+h\right)$, the measured $z$ attribute in the location $u_{\alpha }$, separated by the $h$ vector and grouped by distance for each lag interval [38]. The spatial continuity of the studied phenomenon (i.e., elevation) across the entire analysed area was described as a semivariogram function that belongs to the group of measures of spatial variability. It was calculated as half of the variance, the square root difference between the two locations $z\left(u_{\alpha }\right)$ and $z\left(u_{\alpha }+h\right)$, distant by vector $h$ [39], as in Equation (1):

$$\gamma \left(h\right)=\frac{1}{2N\left(h\right)}{\displaystyle \sum }_{\alpha =1}^{N\left(h\right)}{\left[z\left(u_{\alpha }\right)-z\left(u_{\alpha }+h\right)\right]}^2$$

(1)

where $N\left(h\right)$ is the number of the pairs of points.

On this basis, a description of spatial autocorrelation was made for locations where no elevation data were obtained. As a result, a theoretical model was obtained that was approximated using a theoretical semivariogram function—a Gaussian model with a practical range $a$ [38], Equation (2):

$$\gamma \left(h\right)=1-exp\left(-\frac{3h^2}{a^2}\right)$$

(2)

The correctness of the approximation was verified by cross-validation (CV), calculating the following errors: mean error (ME) (Equation (3)), root mean square error (RMSE) (Equation (4)) and root mean square standardized error (RMSSE) (Equation (5)):

$$\mathrm{ME}=\frac{{{\displaystyle \sum }}_{\alpha =1}^n\left(Z^{*}\left(u\right)-z\left(u_{\alpha }\right)\right)}{n},$$

(3)

where $Z^{*}\left(u\right)$ is the predicted value from the CV, $z\left(u_{\alpha }\right)$ is the observed value and n is the number of sampling points.

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{\alpha =1}^n{\left(Z^{*}\left(u\right)-z\left(u_{\alpha }\right)\right)}^2}{n}}$$

(4)

$$\mathrm{RMSSE}=\sqrt{\frac{{{\displaystyle \sum }}_{\alpha =1}^n{\left[\left(Z^{*}\left(u\right)-z\left(u_{\alpha }\right)\right)/{\sigma }^{*}\left(u_{\alpha }\right)\right]}^2}{n}}$$

(5)

where ${\sigma }^{*}\left(u_{\alpha }\right)$ is a standard prediction error for the location $u_{\alpha }$.

Values close to zero in the case of ME and RMSE, respectively, and one in the case of RMSSE, provide information about the correctness of the process of modelling the spatial variability of the studied phenomenon. Descriptive statistics and the Box–Cox transformation were carried out in Past3 software [40], while the initial geostatistical analysis was undertaken utilising the gstat package [41,42] in the R software [43].

2.4. Ordinary Kriging

For the estimation of the values in places where no data were recorded, the ordinary kriging method (OK) was applied. This is considered the best linear unbiased estimator (B.L.U.E.). The linearity of the estimator results from the fact that its estimation takes into account linear, weighted combinations of the available data. The unbiased nature is due to an average error of zero. At the same time, kriging is a method that minimises the error variance because it gives the most accurate results [39].

The OK estimation enabled the estimation of the values in places where no data were obtained, i.e., also within the central part of the examined pit. The values estimated in the OK were calculated from the formula below (Equation (6)):

$$Z_{OK}^{*}\left(u\right)={\sum }_{\alpha =1}^{n\left(u\right)}{\lambda }_{\alpha }^{OK}\left(u\right)Z\left(u_{\alpha }\right) \mathrm{with} {\sum }_{\alpha =1}^{n\left(u\right)}{\lambda }_{\alpha }^{OK}\left(u\right)=1$$

(6)

where ${\lambda }_{\alpha }^{OK}$ is the weight assigned to the $n\left(u\right)$ random variables $Z\left(u_{\alpha }\right)$.

As a result, the OK estimation enabled generating standard deviation (SD_OK) maps of the performed interpolation.

Correlation and determination coefficients showing the combination of the estimated and measured values were calculated using the Statistica 13.1 (StatSoft Inc., USA) software. The geostatistical modelling was performed via applying commercial software: ArcGIS 10.7.1 (ESRI, Readlands, CA, USA) and Surfer 16 (Golden Software, Golden, CO. USA).

3. Results

3.1. Present Relief

The current visualisation of the relief of the studied area, generated from the LiDAR data, was verified by field observations carried out in 2018 and 2019. On comparing the image of the obtained hypsometry of the studied area (Figure 3a,b) with the description of the pit contained in the literature (Section 2) and field observations (Figure 3d), it was found that the lowest floor part of the quarry is located at an elevation of approximately 212–214 m a.s.l. In turn, the Podgórski Cemetery bordering the quarry from the east is located at an elevation of 235 m a.s.l. and the anthropogenic embankments adjacent to its western wall are located at an elevation of 225 m a.s.l. The highest point in the studied area is the Krakus Mound, reaching 269 m above sea level and located to the north-east of the quarry pit. The areas south and south-east of Liban Quarry reach an elevation of 250 and 255 m a.s.l., respectively. The maximum denivelation within the entire area including the Krakus Mound is 58 m (Figure 1 and Figure 3).

On comparing the results of the field observations with the provided analysis of the hypsometric map and DEM obtained from the LiDAR, it was shown that the present relief of Liban Quarry and its adjacent areas does not differ from that visualised on the basis of the data obtained from the airborne laser scanning (LiDAR) (Figure 3d).

3.2. Reconstruction of the Relief before the Commencement of Mining Operations

As a result of the digitisation of all the contour lines from the topographic map with a scale of 1:10,000 [25] within the examined fragment of the Krzemionki Hills, a set of spatial data containing 1378 points was collected (Figure 4), excluding the elevation point describing the Krakus Mound, which is an anthropogenic object. The obtained point cloud used in the statistical analyses and geostatistical modelling enabled the development of a hypsometric map with the errors of the estimation carried out with the application of the ordinary kriging (OK) and digital elevation model (DEM).

From the analysis of the descriptive statistics (

Table 1. Selected descriptive statistics for the Krzemionki mining area and its surroundings.

Descriptive Statistics	Measurement Data	Box–Cox Transformation
N	1378	1378
Minimum	220	33.49
Maximum	256.25	36.59
Mean	240.16	35.23
Standard error	0.26	0.02
Variance	91.23	0.66
Standard deviation	9.55	0.81
Median	241.25	35.33
Skewness	−0.24	−0.26
Kurtosis	−6.06	−6.08

), histogram and the QQ plot (Figure 5a,b) obtained on the basis of the collected spatial dataset, it has been shown that the data meet the criteria for the application of ordinary kriging (OK). At the same time, they were transformed (Box–Cox) in order to ensure a constant distribution of variances over the entire analysed area (Table 1).

In order to check the presence of spatial autocorrelation in the studied dataset over a distance of 240 m, lagged scatterplots for pairs of points located in 12 lag intervals every 20 m were used (Figure 6). Due to their decreasing correlation and increasing semivariance with the distance, a semivariogram function was used to describe the spatial variability.

In the next stage, an experimental semivariogram was calculated for which the maximum lag distance was 84 m with the number of lags, i.e., 12. The theoretical semivariogram—a Gaussian model (Figure 7a)—was then developed. The criterion for selecting an appropriate theoretical semivariogram model was based on the course of the experimental semivariogram and its parabolic behaviour at the origin [44]. The theoretical model selected on this basis described best the topographic elevation of gently undulating hills [38].

The analysis of the experimental semivariogram showed that due to the presence of anisotropy, the spatial variability was not identical in all directions. The lowest variability with a 60-meter radius was observed in the direction 87.54° (W–E) (Figure 7b). Perpendicularly to it, but in the N–S direction, the maximum variability of the examined phenomenon was defined, with a radius of 30.47 m. At the same time, at a distance of 24 to 60 m, a hole effect was also recorded, the presence of which is related to a limited number of measuring points in the given direction (Figure 7c). The analysis of the semivariograms showed that the spatial autocorrelations fade away over a distance of 60 m for the maximum variability in the W–E direction (

Table 2. Parameters of theoretical semivariogram models.

Semivariogram Parameters
Semivariogram
	Model	Gaussian
	Nugget effect	0.0003
	Sill	0.0016
Anisotropy
	Major range (m)	60
	Minor range (m)	30.47
	Direction	87.54
	Anisotropy factor	1.97
Search ellipse
	Maximum neighbors	20
	Minimum neighbors	4
	Sector type	4 sectors with 45° offset
	Angle	87.54
	Major semiaxis	60
	Minor semiaxis	30.47

).

The correctness of the semivariogram modelling was verified with cross-validation (CV) (

Table 3. Cross-validation (CV) error values.

Type of Errors	Code	Value
Mean Error	ME	0.004
Root Mean Square Error	RMSE	0.28
Root Mean Square Standardized Error	RMSSE	0.98

). The values of the ME and RMSE errors obtained ranged from 0.004 to 0.28 and were close to zero. This indicates a lack of bias and, at the same time, small differences between the predicted value and the measured value. In addition, the value of the calculated RMSSE reached 0.98, and thus close to one, indicating a slight overestimation at the sampling points. Overall, the analysis of errors indicated a correct approximation of the theoretical model to the generated experimental model, thus no discernible difference between the predicted value and the measured value was observed.

For the collected dataset, graphs describing the distribution of the measured values in relation to the predicted values and calculated errors (residuum and standardised error) were compiled and compared (Figure 8). On this basis, it was assessed that the predicted values were characterised by low scattering in relation to the measured values (Figure 8a). Similarly, in the case of the graphs presenting the residuum error and standardised error compared with the measured values, points were concentrated around a defined regression line (Figure 8b,d). However, from the QQ plot analysis (Figure 8c) presenting the values of the standardised errors in relation to the corresponding quantiles as derived from the standard normal distribution of the data, the values of errors have a distribution close to normal, their predictions are unbiased and the predictions values are close to those obtained from the sampling points.

The OK estimation was performed using the parameters of the search ellipse, defined for four sectors with an offset of 45°, with a minimum (4) and a maximum (20) number of neighbours. The remaining parameters were read from the modelled course of the theoretical semivariogram described for the defined anisotropy (Table 2).

As a result of the geostatistical analyses, a hypsometric map and a digital elevation model (DEM) were obtained for the studied Krzemionki mining area, together with its surroundings for the period when the site was not subject to an anthropogenic transformation (Figure 9a,b).

The analysis of the reconstructed hypsometric map and DEM showing the original relief of the studied area (Figure 9a,b) showed that before the commencement of limestone mining, the location of today’s Krzemionki mining area, together with Liban Quarry and its adjacent areas, had a corrugated surface, with denivelations reaching up to 38 m. The highest part of this area was the Lasota Hill (256 m), located in its eastern part, near the then non-existent Krakus Mound. Similar high values of elevation, up to 255 m above sea level, were displayed by the area of the present working and the south-eastern part of the studied area. The areas with lowest elevations, from 218 to 219 m a.s.l., were located in the north-eastern part. Before 1873, the reconstructed area was separated from the culminations of the Lasota Hill, which included the Krakus Mound with a narrow anticline necking located at 248 m above sea level (Figure 9a,b).

For the final verification of the resulting models (Figure 9a,b), the elevation values in the locations contained in the original dataset were extracted. This allowed an assessment of the correctness of the estimation. As a result, two sets of data—measured and predicted values—which were the basis for calculating the correlation (r) and determination (r²) coefficients, were compiled (Figure 10).

With the use of ordinary kriging (OK), an interpolation error map was generated, i.e., the kriging standard deviation map (SD_OK), which was used to evaluate the uncertainty of the performed estimation, especially in the places where the reconstruction was carried out. This confirms the normal distribution of the obtained errors. The largest estimation error, with a value of ±0.51 m, was applied to the surface inside the present working and decreased as it moved away from the pit (Figure 11a). The rest of the area showed a standard deviation of 0.1–0.4 m for the kriging.

3.3. Present and Reconstructed Relief Design

Two terrain profiles, A–B and C–D with lengths of 750 m and 820 m (Figure 12), respectively, were developed for a detailed analysis of the relief of the quarry and its surrounding areas under consideration, and were generated on the basis of the LiDAR data and reconstructed on the basis of the ordinary kriging (OK) estimation (Section 3). They illustrate the courses of absolute elevations along the points with the greatest interpolation error, with the A–B profile passing at about 99% through the anthropogenically transformed area in the N–S direction (Figure 12), while the C–D profile passes at slightly more than 50% through both the present and reconstructed areas, in the SW–NE direction (Figure 12). In order to verify the difference in absolute elevation between the present and the reconstructed model, the lowest point in the central part of the floor of the working was determined. It was located at the intersection of both delineated terrain profile lines (A–B and C–D) and its elevation was 212.3 m a.s.l. After transferring the location of the point to the reconstructed model, its elevation value of 249.36 ± 0.51 m a.s.l. was read, including the SD_OK value obtained from the calculated standard deviation map (Figure 11) corresponding to this location. The difference between the elevation obtained from the LiDAR data and the estimated elevation reached the value of 37.6 m ± 0.51 m, including the SD_OK.

In the A–B profile that underwent the reconstruction (wherein a lack of elevation data was also noted), it was found that the reconstructed area rose from 221 to 254 m a.s.l. at a distance of 0 to 320 m, then fell to 240 m a.s.l. at a distance of 320 to 450 m, where for another 150 m it remained at a constant elevation of about 240 m a.s.l., and in the final section of its course, it rose at a slight angle to 250 m a.s.l. (Figure 12c).

In turn, the elevation section in the C–D profile was characterised by a different type of relief course compared to the A–B profile (Figure 12c,d). The reconstructed area in the C–D profile stretched over at a distance of 180 to 540 m, with the remaining area conforming to a large extent to the present relief (Figure 12d). At the distance of 0 to 176 m, there were no significant differences in elevation between the two analysed models of the C–D profile—the present and the estimated (reconstructed). In the initial section of this profile (up to 100 m), a decrease in elevation from 247 to 238 m a.s.l. was recorded. In the next section, from 100 to 220 m, a gradual increase in elevation to 255 m a.s.l. was observed. In the next section, up to 380 m (covering the reconstructed area), a slight decrease in elevation to 246 m a.s.l. was recorded, while in the further section, from 380 to 540 m, a gradual increase in elevation to 255 m a.s.l. was recorded. Further on, the profile in its course within the section from 540 to 690 m covered a fragment of the Lasota Hill with an elevation of 255–253 m a.s.l., where the anthropogenic object Krakus Mound is currently located (the formation was not included in the geostatistical modelling due to its anthropogenic nature). In the final section of the C–D profile, at a distance of 690–820 m, a fall in elevation to 224 m a.s.l. was recorded, which is connected with the presence of the northern slope of the Lasota Hill (Figure 12d).

Changes in the relief of the studied area were also analysed by means of histograms presenting the quantitative share of the studied variable (elevation) in a raster image (Figure 3c, Figure 9c and Figure 11b). Histograms were generated on the basis of prepared maps (Figure 3, Figure 9 and Figure 11). Their analysis revealed that, currently, Liban Quarry shows terrain elevations ranging from 212 to 268 m a.s.l. The low frequency of elevation values above 255 m a.s.l. indicate the presence of the Krakus Mound, which is not included in the geostatistical modelling. A distinctly high frequency of elevations in the range from 212 to 215 m a.s.l. comes about due to the occurrence of the mining pit. The remaining elevation values, in the range 216–254 m a.s.l., refer to the points where the relief has not been transformed over the years (Figure 3c). The described changes in the relief of the land are reflected in the histogram developed for the reconstructed model. There are no outliers found on the model referring to the floor of the mine and the Krakus Mound. The frequency of the elevation values on the histogram for the reconstructed model is in a similar range from 218 to 255 m a.s.l (Figure 9c) as compared to the histogram describing the present relief (Figure 3c). The histogram that refers to the ordinary kriging standard deviation map (SD_OK) shows the frequency of the standard deviation errors that result from the OK estimation. In this case, the most common is the error group in the range 0.48–0.51 m a.s.l. and refer to the floor of the pit undergoing the reconstruction process. However, the most numerous group are the minimum values of the standard deviation of the kriging that are from 0.11 to 0.19 m a.s.l.

4. Discussion

This work presents a way to reconstruct the relief of Liban Quarry in Cracow on the basis of a topographic map with a scale of 1:10,000 [25] by means of applying geostatistical tools. The solution used to recreate the anthropogenically transformed relief is different from that presented in the global literature. Until now, the reconstruction of relief has been carried out mainly with photogrammetric methods based on, e.g., archival oblique aerial photographs [8,45]. Moreover, historical models have been created primarily through the vectorisation of the contour lines on maps [20,41] with simple interpolators, such as bilinear interpolation [46] and triangulated irregular networks (TIN) [7].

Our proposed approach on the anthropogenic geomorphology related to the reconstruction of quarry relief takes advantage of the methods of geostatistics. Furthermore, it is based on the vectorisation of locations in an area that has remained unchanged throughout the entire period of mining activity, while taking into account locations with missing data. Thanks to this undoubtedly novel approach, it was possible to define a model of spatial variability in which the presence of anisotropy in the W–E direction (87.54°) and a small hole effect in the N–S direction of maximum variability were identified. The obtained model was described by means of a theoretical Gaussian semivariogram function. The same Gaussian model was applied by Tripathi et al. [47] to describe the electrical conductivity of a saturated soil extract (ECe). However, they did not find the presence of anisotropy because the dataset of 132 soil samples proved to be too small. The presence of directionality of the phenomenon can be determined when the dataset contains over 300 samples [48]. In our case, 1378 points were used for the reconstruction, so the anisotropy was detected correctly. In addition, Tripathi et al. [47] did not find the hole effect. This effect, related to the locations with missing data, was identified in our spatial variability model. The structural analysis procedure was verified by cross-validation (CV), calculating the mean error (ME), root mean square error (RMSE) and root mean square standardised error (RMSSE). For the Gaussian semivariogram used, their values were, respectively, 0.004 (ME), 0.28 (RMSE) and 0.98 (RMSSE) (Table 3). These were higher for the ME and RMSSE and lower for the RMSE, compared with the same ME (–0.11), RMSE (1.68) and RMSSE (0.70) values for the electrical conductivity of a saturated soil extract (ECe) obtained by Tripathi et al. [47]. At the same time, the high RMSE value (1.68) for ECe at the maximum predicted value for this parameter (>13) caused the misclassification error to be 12.92% and the correctness of the obtained map was at 87.08% [47]. In the case of Liban Quarry reconstruction, the misclassification error had the value of 0.2%, so the correctness of the map and the model reached a high level of 99.80%. The difference in the correctness of the obtained maps based on the same Gaussian model, between the one developed for ECe by Tripathi et al. [47] and the one generated for the purposes of the reconstruction of Liban Quarry amounted to 12.72%. Furthermore, a RMSSE value significantly below one (0.70) may have caused a significant overestimation of the model for the ECe obtained by Tripathi et al. [47]. However, the low values of the calculated ME and RMSE errors (close to zero) indicate that the semivariogram functions used in the Liban Quarry reconstruction were properly selected and fitted.

A correctly performed geostatistical modelling allowed the application of estimation by ordinary kriging (OK). As a result, cross-validation (CV) enabled the verification of the correctness of fit of the selected models of the semivariograms to the elevation prediction both in places where the value of the examined variable was recorded and in places where it was not. This was very important when choosing the type of interpolator. Therefore, OK seems to be the best fit for estimating variables in locations not included in the sampling strategy, which is also confirmed by studies such as Tripathi et al. [47] and Yan et al. [49]. This should be considered as a proof of the correct choice of ordinary kriging (OK) as an estimation method, assuming that the kriging standard errors were dependent on the prediction, which was confirmed by the normal distribution obtained from the residuum (Figure 8c).

The evaluation of the model of the reconstructed Liban Quarry relief was based on the correlation (r) and determination (r²) coefficients calculated between the estimated data and that obtained from the topographic map which, in the end, amounted to 0.99. This emphasises the strong relationship between the obtained models of the present and reconstructed relief. At the same time, the values of both coefficients indicated small differences between the reconstructed model and the point cloud obtained in the process of vectorisation and subjected to the geostatistical modelling. The obtained coefficient of determination (r²) was 31.51% higher than the results obtained by Pandey and Kumar [7], who tested other interpolators, including triangulated irregular networks (TIN). Based on the kriging standard deviation map (SD_OK) for Liban Quarry, it was found that increased SD_OK prediction values occurred in locations not included in the sampling strategy, i.e., within the fragments of the present Liban Quarry pit. On the other hand, the areas directly surrounding the analysed area, which at the same time constituted the basis for the reconstruction of the mined-out part of the quarry, were characterised by lower SD_OK values. This result is confirmed by a variogram structured analysis and its CV verification. The maximum SD_OK error value was estimated at 0.51 m, in the central part of the pit. At the same time, as the distance from it increased, the SD_OK error values decreased due to the large dataset obtained from the areas around the working, as well as by the lack of measurement points inside the present pit. However, the error values obtained were very low in relation to the total reconstructed elevation. Significantly higher values of the standard deviation (SD), although for vertical differences between the two point clouds, were obtained by Midgley and Tonkin [45]. They reconstructed the topographic surface of a glacier (Svalbard, Norway) from archival oblique aerial photographs using the structure-from-motion (SfM) photogrammetric technique. The resulting RMSE values for the elevation variable were 14.9 m for the mountain sides, 8.5 m for all the “stable” areas and 2.9 m for the outwash zone and bedrock, respectively [45]. However, as the authors noted, when using photogrammetric techniques, one should bear in mind the limitations resulting from the quality or resolution of aerial photographs. This is particularly important when archival photographs are analysed, which may generate an inaccurate topography, especially when reconstructing locations with a limited field of vision. Therefore, SfM can be recommended to reconstruct an old topography with a vertical error accuracy of up to 10 m [45]. Such limitations were not noted in the case of the reconstruction of the Liban Quarry carried out using geostatistical methods. At the same time, the vertical errors described here with the SD_OK proved to be even more accurate by 14.39 m and 2.39 m, in relation to the highest and lowest RMSE (Figure 11) when compared with those obtained by Midgley and Tonkin [45].

The reliability of the obtained relief model of Liban Quarry (Figure 9b) is confirmed by the archival Tactical Map of Poland with a scale of 1:100,000 from 1934 [24] (Figure 2a). Its comparison with the model revealed that, to a very large extent, it predicts the measured values of the analysed variable, namely the absolute elevation. Despite the lack of base points within the present pit of the Liban Quarry, the obtained course of the contour lines was similar to the one depicted on the 1934 map. This was possible after defining the spatial variability describing the quarry and its surroundings on the basis of which the reconstruction using the OK method was performed (Figure 9a,b; Table 2 and Table 3).

Data obtained from archival maps are sometimes burdened with errors and great simplification, especially in terms of morphology. One example is the Krakus Mound, which was described on the topographic map with only a symbol. On comparing its relief to that obtained from the LiDAR data, differences can be observed. Despite the low precision of archival maps in relation to new cartographic and remote sensing products (based on LiDAR), they still often remain the only source of information about the historic morphology and relief of anthropogenically transformed sites [46]. Therefore, in order to present the differences between the obtained models, two topographic profiles were generated in the course of the study (Figure 12), together with histograms showing the quantitative share of the frequency of the elevation values (Figure 3c, Figure 9c and Figure 11b). Moreover, profiles and histograms were determined for the terrain surface generated from both the LiDAR data and reconstructed with the histogram for the SD_OK.

The profiles were characterised by a high degree of conformity with one another. The A–B profile has reflected the contrast in the course of the reconstructed and contemporary surface quite well, in particular, as it extends in the N–S direction along the entire length of the anthropogenically transformed area. This profile was treated as experimental because it reflected two types of relief, i.e., historical and present. The C–D profile has formed a section that crosses the A–B profile in the central part of the quarry, thus supplementing the A–B profile in the diagonal direction (SW–NE). At the same time, it illustrates the course of both the transformed areas and areas not subjected to anthropogenic pressure. The C–D profile was considered to be the reference profile as it referred to the areas with known absolute elevations from which the data were collected, as well as those for which they were not extracted from the base map. When comparing the models, current and reconstructed (Figure 3b and Figure 9b), to the A–B profile, it can be seen that the historical relief was once characterised by the presence of a hill, which does not exist today. It was replaced by the Liban Quarry workings. Similarly, when comparing both models to the C–D profile, it was found that the prepared relief provided a complete picture of the reconstructed area in which a small, former anticline formation that probably formed the area in the central part of the pit became visible. This is also reflected in the Tactical Map of Poland from 1934 [24], which may be an element connecting the non-existent hill revealed in the A–B profile with the Lasota Hill. On this basis, it can be concluded that in the area of the present pit, in places without elevation data, i.e., along the entire length of the A–B profile and in the central part of the C–D profile along the section from 180 to 540 m, excluding the Krakus Mound (presented in general on the topographic map), the relief has been correctly estimated, taking into account the greatest interpolation error ±0.51 m (Figure 11). A similar evaluation of the topography using surface profiles was carried out by Midgley and Tonkin [45]. It referred to the profiles of the surfaces of the Midtre and Austre Lovénbreen glaciers, which were reconstructed on the basis of archival oblique aerial photographs and data from LiDAR. However, the authors did not perform a comparative analysis of the reference profile against the profile of the fully reconstructed glacier surfaces.

The conclusions obtained from the analysis of the topographic profiles were confirmed by the histogram analysis (Figure 3c, Figure 8c and Figure 10b), which quantify the height ranges for the current and reconstructed relief and kriging standard deviation (SD_OK). It should be noted, however, that the resulting SD_OK range is small, as with the elevation differences reaching about 30 m, a maximum error of 0.5 m is not noticeable (Figure 11b). The evaluation of the topography using histograms was also carried out by Jaskulski and Nowak [20]. The authors reconstructed the topography of the Bełchatów Coal Mine and compared histograms in terms of the slope of the reconstructed and present surface.

In summary, the lack of ground or aerial elevation measurement data from before mining can make it significantly more difficult to assess the accuracy of the reconstructed relief based on the available topographic map [7]. However, the use of the data obtained from the topographic map (1:10,000) [25] in our study, as a base for the reconstruction of the former relief within the boundaries of the analysed area, proved to be the right solution. Despite the lack of precise information about the former relief of the studied area, i.e., before the quarry was opened, but due to the application of appropriate and advanced calculation methods, it was possible to reconstruct the original relief of the Liban Quarry within the boundaries of the studied Krzemionki mining area. The use of geostatistical methods in our study has revealed their potential and usefulness in studies related to the reconstruction of anthropogenic areas that were subject to open pit mining, which negatively affects the topography of areas covered by its activities. Although the reconstruction of the relief using geostatistics will not achieve the same high level of accuracy as that applying photogrammetric techniques, such as aerial stereoscopic or oblique aerial photographs [7,45,50], it may turn out to be an alternative and accurate method (assuming that the resulting error oscillates at 0.5 m) that is also useful in situations where there is no access to, e.g., aerial photographs.

Contemporary research involving the reconstruction of former relief in anthropogenic areas is important, among others, for understanding the dynamics of changes in the relief and their assessment in the temporal context. In addition, they provide a basis for developing strategies for restoring the landscape to the condition before the commencement of anthropogenic activities. The use of innovative methodological approaches in our study can serve as a model and become useful in the reconstruction of other anthropogenic areas, especially those affected by open pit mining.

5. Conclusions

Mining operations are the cause of large transformations in topography. One form of such activity is open pit mining. This contributes to irreversible changes in the landscape. In such cases, the attempt to reconstruct the relief appears to be an important measure that should be considered more widely, despite the often insufficient, incomplete or even missing data related to the former relief. The results of such research may be helpful in geomorphological analyses, and thus in planning the reclamation of the aforementioned areas or for future investment projects. In addition, they provide an opportunity to develop strategies for restoring the landscape to the condition from before the commencement of anthropogenic activities.

The geomorphologically diverse Liban Quarry within the Krzemionki mining area can be considered as a model object for the research on changes in the relief of areas after open pit mining activity. It is characterised by the presence of vertical rock walls and a nearly horizontal rock floor with cavities where small water reservoirs are formed. On both sides of the quarry, there are also anthropogenic embankments, and in the bottom of the quarry there are debris cones.

The relief of the studied area (which included Liban Quarry, as well as the Krzemionki mining area, together with the adjacent areas before the start of the mining activities in the quarry) was visualised as a hypsometric map and DEM based on a set of point data obtained from a topographic map with a scale of 1:10,000. The estimation was carried out with the use of ordinary kriging (OK), one of the geostatistical interpolators. This has enabled the authors to reproduce the relief of a section of Krzemionki Hills from the period before the start of mining operations. Despite insufficient data on the former morphology of the area, but with the help of appropriate and advanced geostatistical modelling methods, it was possible to fully reconstruct this relief. The result is a reconstructed original morphology of the analysed quarry located within the boundaries of the Krzemionki mining area. This is the first study of its kind relating to the reconstruction of a quarry relief in an anthropogenically transformed area. The innovative methodological solution as applied may serve as a model and be useful in the reconstruction of other anthropogenic areas, in particular those after open pit mining activities.