Evaluation and Spatial Variability of Cryogenic Soil Properties (Yamal-Nenets Autonomous District, Russia)

: Agricultural development in northern polar areas has potential as a result of global warming. Such expansion requires modern soil surveys and large-scale maps. In this study, the abandoned arable experimental ﬁeld founded by I.G. Eichfeld one century ago in Salekhard city (Russian Arctic), located in the polar circle, was investigated. Our aims were to assess the nutritional soil properties and their spatial variability. For spatial assessment and mapping, ordinary kriging (OK) and inverse distance-weighted (IDW) methods were employed. We found that due to long-term agriculture use, the soil cover was represented by a unique Plaggic Podzol (Turbic) that is not typical of the region. The soil was characterized by relatively low soil organic carbon (SOC) content, high acidity and a high content of plant-available forms of phosphorus in the humus-accumulative horizon. The results showed that some properties (pH H 2 O, pH CaCl 2 ) were characterized by large-scale heterogeneity and showed clear spatial dependence. However, some properties (ammonium and nitrate nitrogen, basal respiration) showed a pure-nugget effect, presumably due to experimentation with fertilizer over a long period of time.


Introduction
Agriculture is the most sensitive to climate change among all sectors of the economy. Numerous global climate change reports indicate that the rate of warming in Arctic regions is accelerating [1][2][3]. Despite the many negative consequences for agro-ecosystems (diseases, insect pests, carbon emissions, erosion, etc.), warming in the northern regions can have advantages [4]. These benefits include an expansion of arable land and an increase in crop diversity and yields [5]. For example, farming practices have been carried out for centuries in the circumpolar areas of Russia [6][7][8], as well as in other Nordic countries [9].
Nowadays the majority of agricultural soils are in an abandoned, fallow state due to the winding down of agricultural practices in the 1990s. Only small areas of arable soil were cultivated without interruption in the Arctic zone in the post-Soviet period. Numerous experimental fields were founded according to Vavilov's idea of polar agriculture in the first part of the 20th century [10]. One of them is the experimental field of the Yamal monitoring station (Figure 1). Arable land began to be cultivated in 1932 as an experimental site to study the adaptation of different types of crops to the conditions of the Far North [11]. Here, Soviet scientists were engaged in the development of new frost-resistant varieties of cereals, vegetables and forage crops, and the creation of local agro-technical measures. During the initial development, the adverse effects on yields of factors such as slow soil thawing, high moisture in the spring period, high soil acidity, early frost and the irregular distribution In recent years, Russia has implemented many global projects to foster new development of the Arctic territories [13]. The expansion of farming in the Arctic areas may become more necessary in the future due to the increasing population and depleting soil resources in the southern regions. The successful management and use of arable land require the study of soil properties and their spatial distribution [14]. Thus, up-to-date largescale soil-property maps are also required. The process of creating large-scale soil maps requires field sampling work, while medium-and small-scale maps can be created via cartographic generalization of larger-scale maps.
In the USSR, large-scale soil mapping began to actively develop in the first half of the 20th century. This was necessary to account for and value agricultural land and for the rapid recovery of the economy and agriculture [15]. Almost all agricultural land in the USSR was covered by a network of large-scale soil maps. These cartographic works were completed using the traditional soil-cartography method based on the idea of V.V. Dokuchaev-formed at the turn of the 19th-20th centuries-which stated that under In recent years, Russia has implemented many global projects to foster new development of the Arctic territories [13]. The expansion of farming in the Arctic areas may become more necessary in the future due to the increasing population and depleting soil resources in the southern regions. The successful management and use of arable land require the study of soil properties and their spatial distribution [14]. Thus, up-to-date large-scale soil-property maps are also required. The process of creating large-scale soil maps requires field sampling work, while medium-and small-scale maps can be created via cartographic generalization of larger-scale maps.
In the USSR, large-scale soil mapping began to actively develop in the first half of the 20th century. This was necessary to account for and value agricultural land and for the rapid recovery of the economy and agriculture [15]. Almost all agricultural land in the USSR was covered by a network of large-scale soil maps. These cartographic works were completed using the traditional soil-cartography method based on the idea of V.V. Dokuchaev-formed at the turn of the 19th-20th centuries-which stated that under certain conditions of soil formation factors in a particular area, only a certain soil could form [16]. This approach has formed the basis of soil mapping around the world [17]. After the collapse of the Soviet Union, the economy and agriculture deteriorated, and this also affected soil mapping [18]. Nowadays, digital large-scale mapping of the soil properties of arable land in Russia is predominantly carried out at more southern latitudes [19][20][21][22][23][24].
Under climate change, Siberian territories have become the potential food basket of Eurasia. Thus, the soils of many cryolithozone regions will be involved or reinvolved in agriculture. That is why the investigation of the quality of arable soils of the cryolithozone should be intensified in terms of ecosystem services of agricultural ecosystems. The aims of this study are the following: 1. To evaluate the soil agrochemical properties of an arable plot located in Salekhard city (Yamalo-Nenets Autonomous District, Russia), and 2. to estimate the spatial variability of soil properties using geostatistical methods and create digital maps.

Study Site Description
The study area is an abandoned arable field of the Yamal Agricultural Experimental Station and is located within the city of Salekhard (Yamal-Nenets Autonomous District) ( Figure 1). The district occupies the Arctic, subarctic and temperate climate zones. Salekhard city and the study area are located on the border between subarctic and temperate climates, and are characterized as "cold continental" (Dfc) by the Köppen climate classification [25]. The climate is characterized by short summers and very harsh, long winters with polar nights. The average January temperature is −23.2 • C; the average July temperature is +14.8 • C. The average annual precipitation ranges from 450 to 500 mm, most of which falls in July and August. The number of days with snow cover and persistent frost is up to 200 per year.
The area is characterized by naturally leveled topography with an elevation range of 8-14 m above sea level. No erosion processes have been observed. The parent materials of the region are modern alluvial and lake alluvial Quaternary sediments of loamy and sandy loam composition [26], while that of the study plot is alluvial middle-Holocene sands. Pristine soils in the vicinity of Salekhard city are mainly represented by reference groups of Histic Cryosols, Histic Gleysols and Podzols (without the thick organic horizon) [27,28]. Today, the abandoned field that represents the research object contains motley-grass vegetation cover of Artemisia tilesii Ledeb, Poaceae, Calamagróstis epigéjos, Campanula rotundifolia L., Chamaenérion angustifolium, Crepis nigrescens Pohle, Deschampsia cespitosa, Elytrígia répens, Mathricaria hookeri (Sch. Bip.), Ranunculaceae Juss, Stellāria gramīnea, Taráxacum officinále, Vícia crácca, etc.

Soil Samples and Analyses
The soil sampling and assessment work were carried out under dry conditions in August 2021. One full soil profile pit was excavated. Soil samples were collected in two replicates from each soil genetic horizon, whereas the Ap horizon was divided into 3 parts and samples were taken every 10 cm. For the geostatistical analyses, the soil sampling work was conducted via a stratified simple random-sampling scheme. Most of the studies on digital soil mapping approaches predict soil property from a single topsoil depth (0-30, or various intervals from 0 to 30 cm) [29]. Considering the surface layer as the most vulnerable to change, we estimated the spatial variability of soil properties for a depth of 0-10 cm. A total of 75 soil samples for the analysis were collected from the surface layer ( Figure 1). The exact coordinates of each soil point were identified using a global positioning system (GPS) with an accuracy of ±0.5 m.
All collected samples were air-dried, homogenized and sieved. The chemical analyses were carried out using standard laboratory methods [30][31][32][33]. Values of pH in water and salt suspension were measured using pH −150 m (1:2.5 soil:solution ratio). SOC content was determined via a wet-combustion method (Walkley-Black method). The content of available forms of ammonium (N-NH 4 ) and nitrate nitrogen (N-NO 3 ) were determined using a potassium chloride solution. The content of mobile available phosphorus (P 2 O 5 ) and potassium (K 2 O) was determined using the Kirsanov method. Soil basal respiration is defined as the steady rate of respiration in soil, which originates from the mineralization of organic matter. Basal respiration was evaluated by measuring CO 2 in a sodium hydroxide 0.1 molar solution. The incubation of CO 2 was conducted for 10 days in plastic sealed containers, according to the standard protocol [34,35].

Spatial Assessment and Digital Mapping
The ordinary kriging (OK) method was used to determine the spatial variability of soil properties on the study plot. The basic OK approach is to predict values in the nonsampled points from a weighted average of the observed points [36]. We selected a spherical variogram to fit the models. The formula for OK is defined as: where Z(x 0 ) is the predicted value at the nonsampled location x 0 ; Z(x i ) is the observed value at location x i ; and w i is the weight.
To verify the spatial distribution based on the OK approach, we created maps using the inverse distance-weighting (IDW) method. The IDW is one of the most widely used tools for interpolating various values, including soil properties [37]. The formula for IDW is defined as: where Z(x 0 ) is the predicted value at the nonsampled location x 0 ; x i is the ith data value; h ij is the separation distance between the interpolated values; and β indicates the weighting power.
For the assessment of the spatial correlations between the study properties, we used the nugget/sill ratio (C 0 /(C 0 + C 1 )). A ratio ≤ 25% indicates that the variable has strong spatial correlation; 25-75% indicates a moderate spatial correlation; and a ratio ≥ 75% indicates a weak spatial correlation [38]. Additionally, we applied a cross-validation procedure for evaluating the performance of the OK approach. The root-mean-square error (RMSE) was used to validate the spatial-prediction accuracy. The RMSE is defined as: where O i and P i are the observed and predicted values of the soil properties, and n is the number of samples. For the statistical analyses, OK and IDW interpolations were performed using the "gstat" packages in R 4.0.4 [39] and RStudio (version 1.3.1093, Boston, MA, USA) [40].

Results and Discussion
The long-term agricultural use of this site has resulted in the formation of a unique soil profile (Figure 2), radically different from the typical soils of the region [41]. Over nearly a century of agro-industrial activities has formed a thick (up to 30 cm), brown, light loamy, root-penetrated humus-accumulative (Ap) horizon. The organomineral horizon is underlain by a massive (45 cm) illuvial-iron, sandy loam, gleyic with placic-horizon layers (Bs). Further, the soil texture changes to sandy and unstructured with placic-horizon layers (BCg and Cg@ (Turbic)) with reductimorphic colors and spots and features of cryoturbation. Based on the morphological description, the soil cover of the plot is represented by Plaggic Podzol (Turbic) according to the IUSS Working Group WRB classification [42]. It should be noted that an active permafrost layer on this plot was already detected at a depth of 160 cm Soil Syst. 2022, 6, 65 5 of 14 in previous studies [43]. Nevertheless, in our field studies, a permafrost upper boundary was not detected at a depth of 230 cm, even though studies were conducted in the same periods of the year (August). by Plaggic Podzol (Turbic) according to the IUSS Working Group WRB classification [42]. It should be noted that an active permafrost layer on this plot was already detected at a depth of 160 cm in previous studies [43]. Nevertheless, in our field studies, a permafrost upper boundary was not detected at a depth of 230 cm, even though studies were conducted in the same periods of the year (August). The vertical distribution of Plaggic Podzol (Turbic) soil properties is shown in Figure  3. The pH values over the soil profile on the field of the Yamal agricultural station are distributed heterogeneously. The lowest acidity values are detected in the humus-accumulative horizon (Ap). Below, in the Bg and BCg horizons and parent material, the pH H2O is close to 7. The distribution of the organic-matter content in the soil profile is quite inhomogeneous. The highest content is in the top 30 cm of the soil profile (horizon Ap), where the quantity of soil organic matter reaches 2%. A reduced concentration occurs in the mineral horizons Bs-Cg@, but it should be noted that the processes of removing  The content of key nutrients is one of the main factors reflecting the fertility of soils. The concentration of available phosphorus, exchangeable potassium and mineral forms of nitrogen is especially important to know for the humus-accumulative (arable) horizon, as it is the main plant-feeding and root-containing horizon [44,45]. The soil is characterized by high acidity, and at the same time, in the humus-accumulative horizon Ap, we observe an extremely high content (Figure 3) of plant-available forms of phosphorus (>400 mg/kg) according to the classification [46].
The processes of nutrient migration are also related to the particle size distribution of soils [47]. In our study, the humus-accumulative horizon Ap and the underlying mineral horizons are fundamentally different in this parameter (from sandy loam to sand). This reflects in the vertical distribution of potassium, which gradually increases down the Ap horizon; then, its concentration decreases sharply at the transition to the mineral horizons. The study of denitrification-ammonification processes is especially relevant in polar regions in the context of climate change. The maximum content of N-NO3 (the end-product of the nitrification of nitrogen) is in the upper 10 cm of the soil profile, which, during the growing season, is maximally heated and dried out. The maximum content of NNH4 and N-NO3 is in the organomineral humus-accumulative Ap horizon. According to the content of nitrogen at different depths of the Ap horizon, we can conclude that the nitrification processes proceed differently. As shown in Figure 3, the content of N-NH4 is higher at depths below 10 cm, and the reverse is the case for N-NO3; this may be due to both the depth of penetration of the root systems of plants and the features of aeration. The content of key nutrients is one of the main factors reflecting the fertility of soils. The concentration of available phosphorus, exchangeable potassium and mineral forms of nitrogen is especially important to know for the humus-accumulative (arable) horizon, as it is the main plant-feeding and root-containing horizon [44,45]. The soil is characterized by high acidity, and at the same time, in the humus-accumulative horizon Ap, we observe an extremely high content (Figure 3) of plant-available forms of phosphorus (>400 mg/kg) according to the classification [46].
The processes of nutrient migration are also related to the particle size distribution of soils [47]. In our study, the humus-accumulative horizon Ap and the underlying mineral horizons are fundamentally different in this parameter (from sandy loam to sand). This reflects in the vertical distribution of potassium, which gradually increases down the Ap horizon; then, its concentration decreases sharply at the transition to the mineral horizons. The study of denitrification-ammonification processes is especially relevant in polar regions in the context of climate change. The maximum content of N-NO 3 (the end-product of the nitrification of nitrogen) is in the upper 10 cm of the soil profile, which, during the growing season, is maximally heated and dried out. The maximum content of NNH 4 and N-NO 3 is in the organomineral humus-accumulative Ap horizon. According to the content of nitrogen at different depths of the Ap horizon, we can conclude that the nitrification processes proceed differently. As shown in Figure 3, the content of N-NH 4 is higher at depths below 10 cm, and the reverse is the case for N-NO 3 ; this may be due to both the depth of penetration of the root systems of plants and the features of aeration.
The agriculture development of the Arctic soils leads to a significant transformation of the soil profile. The pristine soil cover of the Yamalo-Nenets Autonomous District Soil Syst. 2022, 6, 65 7 of 14 is characterized by strong acidity (pH H 2 O < 5) caused by the course of podzolization processes [28,48]. Previously, it was also found that abandoned agricultural soils are highly acidic, despite the implementation of melioration measures (the application of ash to the soil) to reduce the soil acidity that occurred in the past [43,48].
According to previous studies, the native soils of adjacent areas are characterized by higher SOC values compared to the studied soils [41,43]. Soils under natural vegetation have a surface layer (O) up to 10-15 cm thick, formed by decomposed organic residues, with SOC content up to 15%; this disappeared as a result of plowing, but in the underlying horizons, the SOC content is comparable.
Native northern soils are not high in nutrients [49]. The concentrations of nutrients (P 2 O 5 , K 2 O and N-NH 4 ) in the study plot are higher than in the background soils. It was previously shown that the content of available phosphorus and exchangeable potassium between abandoned agricultural soils and the background soils of the region can be increased up to 10-fold [50]. However, the degradation of soil fertility in polar ecosystems has not been studied enough, although it has long been proven that the content of available phosphorus closely correlates with soil acidity and the types of fertilizers used [47]. The above-described differences between abandoned and virgin soils can be explained by anthropogenic factors, including fertilizer application.
The general statistics of soil properties (mean, minimum, maximum, standard deviation (SD), coefficient of variation (CV) and median) based on 75 samples at a 0-10 cm depth for geostatistical modeling are shown in Table 1. The studied soil samples have low-to-medium SOC content, and according to pH, the H 2 O values are characterized as moderately-to slightly-acidic reactions. The values of basal respiration range from 14.62 to 30.93 mg CO 2 /100 g/day. The content of available phosphorus (P 2 O 5 ) in the soil ranges from 165 to 1268 mg/kg, which also characterizes the plant availability of phosphorus as high-to-very high. The parameters of the models and their evaluation using the nugget/sill ratio and RMSE values are presented in Table 2. The spherical model was applied for describing the variograms of the studied agrochemical properties. The spatial correlation of the study properties is different. Strong spatial dependence is revealed for pH values and phosphorus, with their nugget/sill ratios varying from 5 to 9%. The study and spatial modeling of SOC content have received particular attention as key indicators of soil fertility. In our study, the nugget/sill ratio of SOC and K 2 O is characterized as moderate (61 and 26%, respectively). Despite the moderate spatial correlation of SOC content, it is close to a weak correlation. The SOC model using the OK approach is also characterized by a relatively high RMSE value. Various studies have reported better results for spatial modeling of SOC and humus content in croplands using the kriging approach [24,[51][52][53], remote-sensing data [54,55] and morphometric parameters of topography calculated using a digital elevation model [20,23].
An assessment of the spatial variation in agrochemical soil properties using the OK method is shown in Figure 4. We do not illustrate CO 2 , N-NO 3 and N-NH 4 variability maps due to the nugget effect, as they suggest that the variation is not spatially structured, and they are not informative. Maps of the spatial distribution of soil acidity and K 2 O are characterized by relatively large-scale heterogeneity and clear spatial dependence. Areas showing a particular property value are few and smoothly continuous. The study and spatial modeling of SOC content have received particular attention key indicators of soil fertility. In our study, the nugget/sill ratio of SOC and K2O is ch acterized as moderate (61 and 26%, respectively). Despite the moderate spatial correlati of SOC content, it is close to a weak correlation. The SOC model using the OK approa is also characterized by a relatively high RMSE value. Various studies have reported bet results for spatial modeling of SOC and humus content in croplands using the krigi approach [24,[51][52][53], remote-sensing data [54,55] and morphometric parameters of topo raphy calculated using a digital elevation model [20,23].
An assessment of the spatial variation in agrochemical soil properties using the O method is shown in Figure 4. We do not illustrate CO2, N-NO3 and N-NH4 variabil maps due to the nugget effect, as they suggest that the variation is not spatially structur and they are not informative. Maps of the spatial distribution of soil acidity and K2O characterized by relatively large-scale heterogeneity and clear spatial dependence. Ar showing a particular property value are few and smoothly continuous. The spatial distribution of properties using the IDW approach is show in Figure  We observe a slightly different spatial mapping of some properties in comparison to O This is a characteristic of the pure-nugget model when using the kriging method. A hough the spatial variability of CO2, N-NO3 and N-NH4 is random, the IDW maps all for a clearer assessment of very short-scale heterogeneity compared to OK. The spatial distribution of properties using the IDW approach is show in Figure 5. We observe a slightly different spatial mapping of some properties in comparison to OK. This is a characteristic of the pure-nugget model when using the kriging method. Although the spatial variability of CO 2 , N-NO 3 and N-NH 4 is random, the IDW maps allow for a clearer assessment of very short-scale heterogeneity compared to OK. The spatial-prediction accuracy according to the RMSE parameter for almost all the properties is slightly higher for OK in comparison to the IDW method, while the prediction results of the pH values are similar ( Table 2). Qiao [56] reported similar results for estimating the spatial distribution of arsenic content in Beijing soils. In this work, the results of cross-validation using the OK and IDW methods were almost the same. The authors concluded that the OK method is more subjective and requires more soil samples. Similarly, there were no significant differences in the accuracy of similar models for the spatial prediction of the physical and chemical properties of soils at a site in Iran [57]. However, the accuracy the phosphorus OK model is markedly superior to the IDW model, which necessitates the verification of these approaches for future research and mapping.
Such different features of the spatial distribution of soil elements, especially at a short distance, are important for the implementation of precision-farming systems. The degree of heterogeneity of the mapped soil properties within the same field can be influenced by various factors. Usually, a strong spatial correlation of soil properties may be due to internal factors, while a weak correlation is due to external (anthropogenic) factors [36,38,58]. The variability of properties in arable areas also depends on the factors of soil formation (microtopography), erosion processes, plowing duration, fertilizer and irrigation systems [52,[59][60][61]. Since the site is not subject to degradation processes and no irrigation measures have been carried out on it, we exclude these factors as leading to the variability. The inhomogeneity of properties can also be explained by the presence of variability on a smaller scale (smaller than the distance between measurements), measurement error and an incorrect choice of point locations [62].
At our experimental site, we primarily relate the absence of spatial correlation to a large number of different fertilization and crop-cultivation experiments in previous years (Table A1). During this period, the field was divided into separate micro-parts, where fertilizers at different doses were applied annually, crops were grown under different precursors, and harrowing and plowing at different depths were carried out. A recent study [63] also links the history of previous agricultural practices (plot division, nutrient-application rates, crop rotations) to poor predictions of soil properties.
Among the studied properties, the highest coefficients of variation were determined for SOC and nutrients (Table 1)  The spatial-prediction accuracy according to the RMSE parameter for almost all the properties is slightly higher for OK in comparison to the IDW method, while the prediction results of the pH values are similar ( Table 2). Qiao [56] reported similar results for estimating the spatial distribution of arsenic content in Beijing soils. In this work, the results of crossvalidation using the OK and IDW methods were almost the same. The authors concluded that the OK method is more subjective and requires more soil samples. Similarly, there were no significant differences in the accuracy of similar models for the spatial prediction of the physical and chemical properties of soils at a site in Iran [57]. However, the accuracy the phosphorus OK model is markedly superior to the IDW model, which necessitates the verification of these approaches for future research and mapping.
Such different features of the spatial distribution of soil elements, especially at a short distance, are important for the implementation of precision-farming systems. The degree of heterogeneity of the mapped soil properties within the same field can be influenced by various factors. Usually, a strong spatial correlation of soil properties may be due to internal factors, while a weak correlation is due to external (anthropogenic) factors [36,38,58]. The variability of properties in arable areas also depends on the factors of soil formation (microtopography), erosion processes, plowing duration, fertilizer and irrigation systems [52,[59][60][61]. Since the site is not subject to degradation processes and no irrigation measures have been carried out on it, we exclude these factors as leading to the variability. The inhomogeneity of properties can also be explained by the presence of variability on a smaller scale (smaller than the distance between measurements), measurement error and an incorrect choice of point locations [62].
At our experimental site, we primarily relate the absence of spatial correlation to a large number of different fertilization and crop-cultivation experiments in previous years (Table A1). During this period, the field was divided into separate micro-parts, where fertilizers at different doses were applied annually, crops were grown under different precursors, and harrowing and plowing at different depths were carried out. A recent study [63] also links the history of previous agricultural practices (plot division, nutrientapplication rates, crop rotations) to poor predictions of soil properties.
Among the studied properties, the highest coefficients of variation were determined for SOC and nutrients (Table 1), where N-NH 4 and N-NO 3 have the highest values (54.16 and 82.62%, respectively). Alternating organic-and mineral-fertilizer applications in recent years of experimentation still determine increased nutrient content (Table A1). Previous studies of abandoned lands in this region have also shown that the nutrient content can remain high, even after 20 years [12]. Thus, we assume that the content and spatial distribution of some nutrients is still chaotic due to the fragmented application of different fertilizers. Such features must be taken into account in the further management of arable land and nutrient application. Since northern soils are less involved in agriculture and the warming climate provides an opportunity to expand arable land in this region, more research is needed regarding agricultural practices.
We also emphasize that improving the accuracy of digital maps is possible by using additional explanatory variables (satellite data, digital elevation models and others). However, for small plots, as in this study, the use of free-access covariates is not always appropriate due to their relatively low spatial resolution. In such cases, spatial data should have high and ultra-high resolutions, which leads to additional costs in purchasing them. Therefore, at the preliminary stage, the use of geostatistical methods to assess the variability and mapping of properties is advisable.

Conclusions
The expansion of northern farming as a result of a milder climate requires the deeper study of local soils and farming practices, including fertilization application. The present study demonstrated an assessment of the agrochemical properties of an abandoned arable plot with a long and unique history in Salekhard city (Russian Arctic). As a result of the involvement of soils in plowing and long-term use, the site has soil cover that differs from other pristine soils in the region. The soil cover is represented by Plaggic Podzol (Turbic) and is characterized by relatively low SOC content and high acidity caused by podzolization processes. We found that nutrient content can remain high, even after 15 years of abandonment.
We investigated the spatial distribution of the studied properties using ordinary kriging and inverse distance-weighted methods in the upper 0-10 cm layer. Our study showed that, in general, the accuracy of the two mapping approaches is comparable. We were not able to estimate a reliable spatial distribution for all of the properties. Some properties were characterized by spatial dependence, while others had a pure-nugget effect. We consider external factors, such as long-term experiments with fragmented fertilizer and crop rotation, as major influences on the random and short-scale distribution of the study properties. Nevertheless, spatial-prediction techniques such as kriging can be used as an initial method for large-scale mapping on a small plot. For effective land management, there is a need for modern and accurate digital maps to assess the spatial distribution of soils and their properties. Further, the accuracy of digital maps could be improved by adding spatial explanatory variables and increasing the number and density of soil samples.  Data Availability Statement: Data may be requested from the corresponding author for an appropriate reason.

Acknowledgments:
The authors are grateful to Gerard Heuvelink from Wageningen University and ISRIC-World Soil Information for their valuable advice on the geostatistical aspects of the study. We also thank the anonymous reviewers for their careful reading of our manuscript and their many insightful comments and suggestions.

1974
Experience for 52 options Experiments on growing 52 varieties of potatoes