The study of spatial differentiation of soil properties is undoubtedly one of the most important directions of digital soil cartography and the basis for optimizing agricultural activities, taking into account natural factors. To this end, the scientific and methodological framework of data analysis based on probability theory, known as geostatistics, was employed. This apparatus enables the expression of random phenomena in space and/or time through mathematical and statistical indicators and the classification of the regularities of their distribution.

The extensive range of geostatistical methods, when applied in a dynamic and innovative manner, has enabled the advancement of soil geography, utilizing state-of-the-art technologies and contemporary software solutions. This has facilitated the development of methods for precise, direct, or ”mathematical“ mapping of soil cover and its characteristics [Webster 2007; Krasilnikov & Targulyan 2019]. The development of such methods has enabled the establishment of a new research direction in soil science and soil geography, known as pedometrics [Savin et al. 2019; McBratney et al. 2018].

It is essential to note that direct (mathematical) mapping and geostatistics are not in opposition to the classical approach to land cover mapping. However, as demonstrated by chaos theory, the presence of order and regularity at one hierarchical level does not preclude their absence at another [Esmaeilizad et al. 2024]. The concept of a transition to a ”new soil geography“ was introduced in 2019, representing a synthesis of two divergent approaches [Krasilnikov & Targulyan 2019]. Presently, considerable attention is being directed towards the utilization of machine learning algorithms for the prediction of soil properties. However, methodologies predominantly based on regression matrices overlook the spatial dimension inherent in the data being studied [Heuvelink & Webster 2022; Vaysse & Lagacherie 2017; Wadoux A.M.J.C., et al. 2021]. In certain instances, the integration of machine learning with residual kriging has been shown to yield optimal accuracy [Molchanov et al. 2015]. Concurrently, the analysis of spatial data employing variogram analysis facilitates not only the consideration of data distribution in space but also the identification of patterns and consequences of their differentiation.

Over recent decades, the number of studies in geostatistics and pedometrics aimed at investigating these properties has increased significantly (a summary of these works is given in [Oliver 2010] and [Biswas 2024]). Most of them, however, are limited to the study of the accuracy of the instruments [Mohamed et al. 2023; Soropa et al. 2021; Vaysse & Lagacherie 2017] or to the interpretation of the maps obtained [Suleymanov et al. 2023; Xiao et al. 2023], although several studies dedicated to a comprehensive analysis of the causal relationships underlying the distribution of acidity [Helfenstein et al. 2022] and organic carbon [Szatma´ri & P´asztor 2019] should be noted. The studies face a number of difficulties, the most significant of which are, firstly, the labor-intensive process of obtaining data sets, and secondly, the financial expense associated with soil sample analysis. Concurrently, the selection of soil properties for analysis ought to be informed by the available options for their examination. While moisture and granulometric composition appear to be obvious choices for analysis, their study is very labor-intensive. Conversely, analysis of acidity remains relatively straightforward and does not necessitate substantial financial expenditure. It is also important to note that acidity is one of the most significant soil properties.

The degree of soil acidity exerts a significant influence on the eventual soil fertility, with the process of soil acidification representing a form of chemical degradation [Molchanov et al. 2015]. In the context of leaching water management, soil acidification emerges as a pivotal factor contributing to hazardous agricultural practices in the Republic of Belarus [Klebanovich & Vasilyuk 2003].

The consequences of soil acidification include a decline in the productivity of cultivated plants, a disruption in the growth of natural vegetation, and a deterioration in the viability of aquatic ecosystems. The main causes are deficiencies in Ca, Mg, P, and Mo [Guo et al 2015]. These deficiencies may be exacerbated by a flushing water regime, which depletes ions in the soil.

Although studies that investigate the causes and consequences of acidity’s spatial distribution in soil are important, a significant gap remains in understanding how soil properties vary within the classification units themselves.

The objective of this study is to identify the patterns of spatial differentiation of soil acidity on a detailed scale under different types of land in the Berezina River valley (Belarus), and to assess the degree of transformation of the spatial pattern of acidity distribution under the influence of human agricultural activities.

The following section outlines the general characteristics of the territory. The geographical area under scrutiny is delineated by the geographical station (GS) ”Western Berezina”, situated between the Oshmianskaya and Minsk uplands in the valley of the Berezina River. This territory is situated within the Volozhinsky administrative district of the Republic of Belarus (Fig. 1). The rationale behind the selection of this territory is twofold. Firstly, the high dissection of the relief leads to a high heterogeneity of the soil cover, which in turn contributes to the active course of erosion-denudation processes. Second, the area has undergone significant agricultural development.

The objects of the study are six plots of land, categorized as follows: forest (plots 1 and 2), arable (3 and 4), and meadow (5 and 6). The plots were selected to include all major soil types in the study area (Fig. 2), ensuring thorough representation of hydromorphism and granulometric composition. This approach ensures the representativeness of the selected plots for the study’s objectives.

The sampling process was conducted synchronously at all selected sites throughout May 2023, in accordance with the prescribed methodology (GOST 17.4.3.01-83 ”General requirements for sampling“). The sampling depth was maintained between 0–20 cm. The humus-accumulative horizon capacity in forest plot soils was 10–15 cm thick. Soils on arable land had a 30 cm humus-accumulative horizon. Meadow plots in the flooded part of the floodplain had a 25–30 cm horizon, while those on aeolian [Kurlovich et al. 2024] ridges had a 13–15 cm horizon.

Forest plots of 12.6 and 11.3 ha were selected for analysis, from which 80 and 78 samples were collected. The soil cover of the forest plots consists of sod-podzolic soils (Retisols) and sod-podzolic waterlogged soils (Gleyic Retisols), which more often develop on cohesive and loose loams and light loams. In the northern part of both sites, the vegetation is represented by pine orchard-grass on sod-podzolic automorphic soils. In the southern part, it is represented by mossy spruce, which is confined to sod-podzolic, weakly gleyey soils.

In contrast, the arable plots, measuring 10.6 and 17.6 hectares, exhibit less heterogeneity in soil composition. From these plots, 110 and 102 samples were taken for further analysis. The soils are classified as sodpodzolic, developed on cohesive and loose loams and on light loams. At the time of fieldwork, wheat was sown on the plots. Notably, no crop rotation had been implemented on plot 3 in the previous year. The precursor crops on plot 4 were corn in 2022 and fallow in 2021.

The meadow plots, situated within the maned floodplain, exhibited the most extensive area (30.7 and 44.2 ha, respectively) from which 79 and 85 samples were collected. The soil cover is represented by alluvial soddy-gleyey soils (Gleyic Fluvisols), soddy-gleyey soils (Gleyic Fluvisols) and silty-peaty-peaty-gleyey soils (Hemic/Fibric Histosols).

The soil types present at these locations encompass 42 varieties, indicative of automorphic, semi-hydromorphic, and hydromorphic soils. This distinguishes these sites as unique ”catenas”, illustrating the sequential transition among diverse soil types. This facilitates geostatistical analysis and interpretation of the horizontal distribution of acidity across the predominant soil categories.

Research Methods. This study is based on the use of keys and on the principles of continuity of soil cover properties, probability theory, chaos theory, and the entropy approach [4], which underpin geostatistical interpolators.

The construction of the sampling grid was informed by two criteria: the geomorphological features of the area and the uniformity of plot coverage. This approach was adopted to circumvent the occurrence of both clustering of data and ”empty“ areas devoid of sampling.

The analysis of soil acidity (pHKCl) was conducted in the laboratory using potentiometric titration in accordance with the standards outlined in GOST 26484-85. ”Soils. Method of determination of exchangeable acidity”.

The geostatistical analysis of the obtained results was conducted using ArcGIS ArcMap software, following a complex flowchart that was developed based on the analysis of numerous studies [Lark 2012; Samsonova & Meshalkina 2020; Biswa 2024] and empirical research [Klebanovich et al. 2021; Klebanovich et al. 2018] (Fig. 3).

The distinguishing characteristic of the proposed geostatistical analysis algorithm is its integration of theoretical geostatistics with a decision-making algorithm for modeling the studied phenomenon in solving applied geographical problems. Additionally, it offers the possibility of using it to create separate applications for geostatistical analysis.

The initial block delineates the primary phases of spatial heterogeneity analysis. Spatial heterogeneity (f(r)) is defined as the sum of three components: macrocomponents (m(r)), mesocomponents (s(r)) and microcomponents (e(r)) [Yakushev et al. 2010]. The radius vector (r) defines the position of an arbitrary point in the chosen coordinate system (x; y), i.e., r = r (x, y). The statistical distribution of the data is determined by the center of distribution (mode (Mo) and median (Me)), variant dispersion (variation (V) and dispersion (D, σ2)), and distribution shape (asymmetry (A) and excess (E)). The deterministic trend (macrocomponent) is also identified, in addition to variogram analysis and graphical presentation of results.

The primary processing and analysis of the obtained statistical results were carried out in Microsoft Excel software products, where the data were checked for normality of distribution according to the method of E.A. Pustylnik [Pustylnik 1963] (Eqs. 1 and 2). E. A. Pustylnik’s method was chosen because Microsoft Excel uses its formulas to calculate the skewness and kurtosis coefficients. To maintain uniformity, errors were also calculated using his method.

(1)

(2)

where ОА – asymmetry coefficient error; ОЕ – excess error; n – sample size.

The second block is devoted to the most important component of geostatistical analysis: variography, which is the selection of the mathematical function of the variogram to the variogram constructed from real data (empirical variogram). The calculation of the main parameters of the variogram is then undertaken. The nugget-effect (e(r)), representing the ”noise“ component, is defined as the unaccounted variability. The partial sill (s(r)), which corresponds to the variance included in the model, is another key element. The sill (s(r)+e(r)), representing the maximum level of variance at the variogram rank level, is also a crucial factor. The variogram range (R, m) is defined as the distance at which points cease to correlate with each other. Lag size (L, m) is the average distance between sampling points; residual variance (nugget to sill ratio) (D, %) is the fraction of noise and interpolation accuracy.

It encompasses supplementary operations implemented solely when required: trend removal, data transformation, and declustering. The term also includes anisotropy (An, °), nearest neighbor point, and cross-validation (cross-validation) determinations.

The accuracy of variogram fitting was checked by cross-validation, where each point in turn is removed from the calculation, and the remaining points are used to calculate the predicted value at the location of the removed point (the process is repeated for all input points) and the calculation of interpolation errors (Eq. 3–5):

(3)

where МЕ – average forecast error;

the average difference between the measurement and the predicted value;

(4)

where RMSE – forecast mean square error.

(5)

where RMSS – standardized root mean square error.

The following table (Table 1) presents the indicators of descriptive statistics based on the results of laboratory analyses of soil acidity in the selected plots.

The primary distinction among forest, arable, and meadow (floodplain) plots can be identified using the coefficient of variation. Forest plots exhibit variation levels of 11.9 and 16.5%, while arable plots demonstrate the lowest values of 9.52 and 8.74%. The meadow (floodplain) plots, on the other hand, exhibit a coefficient of variation of 14.8%.

Site 1 is distinctly different, where asymmetry and kurtosis reach very high values of 2.40 and 8.16, which are an order of magnitude higher than their errors (0.41 and 0.76, respectively). A slight right-sided asymmetry is observed in plots on arable land, thus showing a greater proportion of high pH values, which is due to the application of chemical inputs. Regarding kurtosis, except for the first plot, slight leptokurtosis is evident in plot 4, while plot 6 shows slight platykurtosis.

During macrocomponent determination, no trends were observed in plots 1 and 6, while plots 2 and 4 exhibited a second-order trend and plot 5 displayed a first-order trend; these were taken into account when selecting variograms (Fig. 4). The selection of variograms was carried out on the basis of forecast errors and residual variance. The primary variogram types analyzed are outlined in Table 2.

There is a high degree of similarity in variogram error rates; however, a substantial discrepancy is observed in residual variance, which was used as the primary criterion for model selection. An exponential variogram was selected for the first forest area, and a pentaspherical variogram was selected for site 2. The spatial distribution of acidity in plots 4 and 5 is described by exponential and pentaspherical variograms, respectively. The soil acidity of the meadow (floodplain) plots is primarily characterized by the exponential variogram (Fig. 4).

The analysis of soil acidity in the forest plots reveals a decrease in dispersion at distances of 140–180 m (more clearly expressed in the variogram of the second plot), which corresponds to the length of the slope in plot 1 and the distance between ravines with temporary watercourses in plot 2. The increase in variance at distances of 220–280 m and the subsequent decline (more clearly evident in the variogram of the first plot) indicates the presence of general deterministic dependencies emerging at distances greater than 400 m.

Variograms of soil acidity on arable plots exhibit a single distinction, namely the presence of anisotropy at plot 4. Otherwise, variograms demonstrate a smooth increase in dispersion at small distances, oscillations in the central part, and a decrease at the end.

It is also noteworthy that the variograms of sites 2 and 4 are almost identical, a phenomenon attributable to their close proximity to each other and, consequently, to geochemical processes.

It is pointed out at site 5 that anisotropy is observed, along with a sharp jump in variance at distances of 130–170 m, which corresponds to the width of the floodplain mane in the central part of the site; this emphasizes significant differences between acidity in the flooded part of the floodplain and on the elevation.

The experimental variogram at these sites shows similar behavior, evidenced by a slight decrease in variance at medium distances (370–500 m). The decrease is attributable to similar natural conditions at sampling points at these distances (points near the terrace and the mane, or flooded areas separated by the mane).

The main parameters obtained from the variograms are presented in Table 3.

The obtained results allow for the determination of the main geostatistical differences among the sites. First, note that the average distance between points (L) for sites within the same group differs insignificantly, allowing comparison of the results.

The sill values for forest soil acidity (s(r)+e(r)) are observed in the plots, ranging from 0.9 to 1.3, for arable soils from thousandths and hundredths of a fraction (0.006–0.01), and for meadow (floodplain) soils from tenths of a fraction (0.53–0.54).

To confirm these differences, a measure of ”variation per meter“ was derived, since the threshold reflects the variance and the range the distance, it is possible to calculate the variance attributable to one unit of distance. The calculation was performed as follows: 1) The ratio of threshold (maximum variance) to range (distance (m) of autocorrelation) was calculated – ”variance/meter“; 2) From ”variance/meter“, the root (standard deviation dissection) was extracted – ”mean square deviation/meter“; and 3)”mean square deviation/meter“ was divided by the arithmetic mean of plot pH – ”variation per meter“. This indicator allows moving from comparisons of absolute values to relative values, thus removing the influence of site size. The values obtained for this indicator confirm the previously described differences: variation per meter is approximately 2% for forest soils, 0.1–0.2% for arable soils, and about 1% for meadow soils.

The observed low values of residual dispersion indicate the high accuracy of the obtained cartograms (Fig. 5).

Almost all cartograms demonstrate a high level of heterogeneity in acidity. In forest plots, values ranging from 4.0 to 8.5 are evident. The first plot is characterized by a smaller spread of pHKCl values (4.0–7.0) and a standard error of interpolation ranging from 0.20–0.37 (2.79 ha (16.8 % of the total plot area) to 3.36 ha (20.2 %)). The remaining part of the plot (10.5 ha) has another area of 0.37–0.50 (63.0 %). The second forest plot has a variety value range of 0.15–0.95. Only 5.7% (0.64 ha) of this area has an error of less than 0.56, while the remaining areas with values from 0.92 to 0.95 account for 26.7% (3.02 ha).

There is significantly lower differentiation of soil acidity, both within and between fields, in arable land (from 5.0 to 7.5). Due to the smaller scatter of values, prediction errors at site 3 do not exceed 0.43 (1.03 ha / 9.11 %). Most of the area has an index between 0.27 and 0.42 (7.41 ha / 65.5 %), while the remainder belongs to values below 0.27 (0.18 ha / 25.4 %). A similar situation is observed at site 4, where only 2.5 % (0.44 ha) of the entire study area has a prediction error greater than 0.45, with a maximum error of 0.63. The minimum – 0.04 units, with an average range of 0.22 to 0.44 (16.75 % to 95.2 %).

The meadow plots (5 and 6), which are located on the maned floodplain, are characterized by high soil acidity heterogeneity. With an average pHKCl value of 5.0–6.0 across both sites, two areas of approximately equal size can be distinguished: 1) areas with acidic and moderately acidic reactions (pHKCl 4.5–5.5), and 2) areas with slightly acidic and close-to-neutral reactions (pHKCl 5.5–6.5). The RMS error of prediction is similar at both sites. The respective areas are 38.1 % (11.9 ha) and 44.5% (20.4 ha) between 0.47 and 0.55. The minimum (0.20–0.34) and maximum (0.55–0.70) values at both sites account for less than 1.5 ha.

The primary statistical processing indicates a substantial differentiation of acidity at all sites, despite relatively modest values of the coefficient of variation (8.74 to 14.8%), which can be attributed to the logarithmic values of H+ ion content. However, when assessing the values of the coefficient of variation qualitatively, it is important to recognize the difficulty in determining which coefficients to consider as large and which as small [3]. Given that pH values are already logarithmic scale, it can be posited that these values may indeed signify a medium to high degree of heterogeneity in acidity.

To address the anomalously high kurtosis index for acidity values at site 1, a rather uncommon method of data transformation was employed, namely the Box-Cox method of the second degree (used for exponential distribution). Conversely, no such transformation was necessary for plot 2. Plots 3 and 4 show right-skewed distributions, necessitating logarithmic transformation of the data to approximate normality. The floodplain plots (5 and 6) demonstrate approximately normal distributions and thus do not require transformation.

It is quite typical to observe trends in arable plots (3 and 4), which are often a consequence of human impact on the soil, such as fertilizer application and plowing. The key plot 3 polynomials indicate an increase in values in the central part, with an increase followed by a decline in the north-south and west-east directions (more pronounced in the north-south direction), which corresponds to a parabolic dependence. A similar polynomial describes the spatial distribution of acidity at site 4; however, in the west-east direction, there is a decrease in values in the central part. The acidity of meadow plot 5 is characterized by a linear dependence, reflecting an increase in pHKCl values from west to east (from terrace to channel) and from south to north (from peaty soils to alluvial soddy and soddy waterlogged soils).

The selection of the variogram is central to any geostatistical analysis, and it is this that has a significant impact on the final interpolation result. As demonstrated in this study, an exclusive focus on standard criteria can yield skewed and inconsistent results. While differences in errors as small as hundredths or thousandths may appear negligible, disparities in residual variance prove to be of paramount importance to the ultimate outcomes. For instance, at site 1, the RMS normalized error between variograms differs by a mere 0.02 (all other indices being equal), while the differences in residual variance reach 30–40%. If the Gaussian variogram is chosen, the cartogram is found to be 50.9% different from the actual one.

Therefore, utilizing a carefully curated selection of models, variogram analysis enables the clear determination of the distances at which analogous natural processes and idiosyncrasies in geochemical conditions manifest. These conditions pertain to the redistribution of chemical elements within the landscape. Forest plots demonstrate a decline in dispersion at distances of 140–180 m, which corresponds to the length of the slope (plot 1) and the distance between ravines with temporary watercourses at plot 2. The arable plots, on the other hand, demonstrate an increase in variance at short distances (70–80 m) and a decrease at longer distances (250 m at plot 3 and 370 m at plot 4). This indicates strong differences between adjacent points and similarity between distant ones, reflecting the disturbances in spatial structure caused by anthropogenic influence on acidity indices. Consequently, these experimental variograms delineate the anthropogenic influence on acidity heterogeneity, with the ”natural component“ of heterogeneity incorporated only at the conclusion. The experimental variograms for meadow (floodplain) sites (5 and 6) demonstrate a pronounced increase in variance at distances of 130–170 m, which corresponds to the width of the floodplain mane in the central region of the site. This variation highlights the substantial disparities between acidity levels directly in the flooded portion of the floodplain and those in areas of higher topography. The slight decrease in variance at medium distances (370–500 m) is attributable to the similarity of natural conditions at sampling locations at these distances (points near the terrace and the mane or the flooded areas separated by the mane).

The distribution of рНКСІ values in forest plots exhibits striking variations. Plot 1 is distinguished by the prevalence of medium acid soils (рНКСІ 4.5–5.0) with isolated areas of strongly acid soils (рНКСІ < 4.5) in the northwestern part of the study area, restricted to sod-podzolic weakly gleyey and gleyey soils. Conversely, the southeastern and eastern regions of the territory are predominantly characterized by acidic soils (рНКСІ 5.0–5.5), though there are areas exhibiting soils with reactions approaching neutral (рНКСІ 6.0–6.5) and slightly alkaline (рНКСІ 7.0–7.5), confined to a hill with soddy-carbonate soils (fig. 1.A.).

Site 2 is characterized by a wide variety of acidity groups, from strongly acidic (рНКСІ 4.00–4.50) soils to medium alkaline (рНКСІ 7.50–8.00), which is explained by the cumulative heterogeneity of vegetation. In the northern part of the site, vegetation consists of an eagle-grass pine forest on sod-podzolic automorphic soils, while in the southern part, there is a mossy spruce forest. The presence of all acidity groups is observed at this site, indicative of cumulative heterogeneity in vegetation. The northern region is characterized by the presence of eagle-grass pine forest on sod-podzolic automorphic soils, while the southern region is dominated by mossy spruce forest. The latter is confined to sod-podzolic, weakly gleyey soils.

The presence of ravines with temporary watercourses contributes to the transfer of particles from adjacent arable land. The bedrock is hypothesized to be composed of carbonate moraine, which is susceptible to erosion, resulting in the removal of the upper soil horizons. The outcrop of carbonate rocks has led to the formation of moderately alkaline soils (pH > 7.50, up to 8.17) at Plot 2, which is atypical for sod-podzolic soils. Furthermore, the leaching of fertilizers from the adjacent arable land (Plot 4) likely contributed to localized alkalinization, further enhancing the acidity heterogeneity at this site. On plots of arable land, soil acidity has demonstrated a lesser degree of differentiation, both within individual fields and between fields. The majority of the area under both plots is characterized by soils with slightly acidic pH levels (рНКСІ 5.5–6.0) and those approaching neutrality (рНКСІ 6.0–6.5). Site 3 also contains an array of acidic soils in the northern part of the site, adjacent to the woodland and aeolian till. A contour of elevated рНКСІ values (7.0–7.5) was identified in the center, as noted in the trend analysis. In general, the arable plots exhibit acidity levels consistently above 5.0, with the presence of micro-contours, alkaline soils (0.98 ha, constituting 9.17% of the plot area) and a contour of near-neutral soils (29.3 ha, accounting for 27.6%) that do not necessitate lime application. However, these micro-contours and near-neutral soils remain undetected during the determination of lime requirement, potentially leading to additional material costs and a decline in crop yield on microplots with рНКСІ >6.

The topography of meadow lands is characterized by the presence of plots situated on the floodplain, in conjunction with forest areas. These regions are distinguished by a pronounced degree of differentiation in soil acidity. The topography of these regions is predominantly level, with the presence of areas characterized by acidic and medium-acidic soils, which are confined to floodplain ridges. In these areas, the processes of acidification are attributable to the predominance of moisture movement in a downward direction, resulting in weakened hydrological links with the river. Additionally, there are areas of poor aeolian sands. Other notable features include areas of weakly alkaline and neutral soils with close groundwater tables, as well as areas remote from the channel, where peat formation processes occur and medium- and strongly acidic soils are formed.

The complex spatial distribution of acidity in the forested areas near the “Western Berezina” GS is driven by vegetation heterogeneity (eagle-grass pine forest in the north and mossy spruce forest in the south) and topography. Ravines with temporary watercourses facilitate the transfer of particles from adjacent arable land, and hypothetical outcrops of carbonate moraine occur where the upper soil horizons have been eroded. In the meadow areas, acidic and moderately acidic soils are confined to floodplain ridges and areas distant from the river channel where peat formation occurs. In contrast, weakly alkaline and neutral soils are found in areas with a high groundwater table. The low spatial variation of acidity in arable land results from a homogenizing anthropogenic impact.

The spatial differentiation of soil acidity is more accurately assessed through geostatistical indicators that clarify the classical ideas about the formation of acid-alkaline geochemical barriers and transformation of natural landscapes under human influence. The lowest degree of soil acidity homogeneity is observed on forest lands with insignificant anthropogenic transformation and is characterized by a dispersion range from 0.90 to 1.30 (± 5%) and relative variation per meter of 1.73–2.14%; the increase of anthropogenic impact on meadow soils is characterized by the increase of anthropogenic impact on the soils. he human impact on the environment includes 40-50 years of liming, as well as plowing and applying basic nutrients at these sites and across the country. This leads to the flattening of natural geochemical processes and the prevention of soil acidity’s natural heterogeneity.

In conjunction with geostatistical indicators, the presented cartograms demonstrate that conventional notions concerning the specific acidity group or range of рНКСІ values attributed to particular soil types necessitate re-evaluation and elucidation through the utilization of contemporary geoinformation and geostatistical analysis methodologies.