Spatial variability of saturated hydraulic conductivity and its links with other soil properties at the regional scale

Saturated hydraulic conductivity (K) is a key property for evaluating soil water movement and quality. Most studies on spatial variability of K have been performed soil at a field or smaller scale. Therefore, the aim of this work was to assess (quantify) the spatial distribution of K at the larger regional scale in south-eastern Poland and its relationship with other soil properties, including intrinsic sand, silt, and clay contents, relatively stable organic carbon, cation exchange capacity (CEC) and temporally variable water content (WC), total porosity (FI), and dry bulk density (BD) in the surface layer (0–20 cm). The spatial relationships were assessed using a semivariogram and a cross-semivariogram. The studied region (140 km2) with predominantly permeable sandy soils with low fertility and productivity is located in the south-eastern part of Poland (Podlasie region). The mean sand and organic carbon contents are 74 and 0.86 and their ranges (in %) are 45–95 and 0.002–3.75, respectively. The number of individual samples varied from 216 to 228 (for K, WC, BD, FI) to 691 for the other soil properties. The best fitting models were adjusted to the empirical semivariogram (exponential) and the cross-semivariogram (exponential, Gaussian, or linear) used to draw maps with kriging. The results showed that, among the soil properties studied, K was most variable (coefficient of variation 77.3%) and significantly (p < 0.05) positively correlated with total porosity (r = 0.300) and negatively correlated with soil bulk density (r = – 0.283). The normal or close to the normal distribution was obtained by natural logarithmic and root square transformations. The mean K was 2.597 m day−1 and ranged from 0.01 up to 11.54 m day−1. The spatial autocorrelation (range) of K in the single (direct) semivariograms was 0.081° (8.1 km), while it favourably increased up to 0.149°–0.81° (14.9–81 km) in the cross-semivariograms using the OC contents, textural fractions, and CEC as auxiliary variables. The generated spatial maps allowed outlining two sub-areas with predominantly high K above 3.0 m day−1 in the northern sandier (sand content > 74%) and less silty (silt content < 22%) part and, with lower K in the southern part of the study region. Generally, the spatial distribution of the K values in the study region depended on the share of individual intrinsic textural fractions. On the other hand, the ranges of the spatial relationship between K and the intrinsic and relatively stable soil properties were much larger (from ~ 15 to 81 km) than between K and the temporally variable soil properties (0.3–0.9 km). This knowledge is supportive for making decisions related to land management aimed at alteration of hydraulic conductivity to improve soil water resources and crop productivity and reduce chemical leaching.

The K value depends largely on the pore size distribution (PSD), especially on the share and continuity of relatively large pores (macropores) 9,[16][17][18][19][20] . In a study conducted by Kim et al. 21 , the area of the largest pores explained almost 80% of variability in soil saturated hydraulic conductivity. As shown by Centeno et al. 22 , macro-porosity can be used as a proxy to estimate the spatial variation of K.
Owing to the high sensitivity to pore size distribution influenced by soil texture and management practices, K displays relatively high spatial variability 9,15,20,22,23 . Therefore, knowledge of the spatial distribution of K is essential in selection of the most appropriate localised management practices and amendments to improve water use efficiency in agriculture and to minimise the use and leaching of chemicals [24][25][26][27] . The spatial distribution of soil properties, including saturated hydraulic conductivity, can be assessed by classical and spatial statistics. The classical statistics can adequately analyse variables that are independent of space 28 . However, when the random variation occurs, geostatistical analysis including direct semivariograms and cross-semivariograms is appropriate 29,30 . Semivariograms define the dependence of the difference between values of a given variable on the distance between sampling locations and, hence, the spatial structure of the variation. They aid in designing a sampling setup with an amount of samples required for satisfactory regionalization of soil properties 31,32 . Once various variables are linked, their combined spatial designs can be evaluated by cross-semivariograms. Semivariogram and cross-semivariogram data and maps obtained using krging and cokriging techniques. Cokriging also allows distinguishing time-consuming and/or expensive variables from those that are more easily measured or available in soil databases. When K shows spatial random variation, the use of both classical statistics and geostatistical models are recommended 30 .
Numerous studies on the spatial variability of soil K have been performed to date at a short scale (< 25 m) (e.g. 33,34 ) or a field scale 35 . However, the variability at a larger regional scale is poorly understood, as suitable spatial characterisation of highly heterogeneous K requires a large number of laborious, time-consuming, and expensive direct measurements 9,36,37 . Therefore, the objective of this study was to determine the spatial distribution of K at the regional scale using a semivariogram and a cross-semivariogram and its relationship with soil properties that affect conductivity and can be more easily measured or obtained from existing soil databases.

Materials and methods
Study area and sampling. Study region with an area of about 140 km 2 is situated in a flat area within Łuków Plain, south-eastern Poland. The height differences in the shallow and often wet river (Krzna) valleys do not exceed 10 m. The study region has mostly low productive Podzol soils 38 derived from sandy and sandy loams of glacial origin. About 80% of the area in the region is used in agriculture, with 62.3% and 18.2% of arable lands and grasslands, respectively. Forests cover only 13.5% of the study area, mostly in the south-eastern and western parts. The climate is largely influenced by the western circulation and polar sea air (about 65% days a year). The average annual air temperature in the region is about 7.3 °C. July and January are the warmest and coldest months with respective mean temperatures of 17.7 and − 3.6 °C. The annual amplitude (the difference between max and min) of air temperature is 21.3 °C and 23.4 °C when calculated from the differences of the average temperature of the hottest (July) and coldest (January) months in individual years. The greatest amount of precipitation is recorded in June and July (more than 70 mm) and the lowest values (less than 30 mm) are noted in January, February, and March. This indicates significant predominance of summer rainfall (212 mm) over winter rainfall (83 mm). The sum of rainfall in the growing season (April-September), i.e. 350.9 mm, constitutes 65.4% of the annual total. During the 50 years under consideration, it ranged from 224 to 530 mm.
Tested soil properties. Soil samples were randomly collected from the 0-20 cm layer into cloth bags and 100 cm 3 (5.02 cm diameter and 5.05 cm high) steel cylinders immediately after harvesting cereals (August). Figure 1 shows the spatial distribution of the soil sampling points. The number of individual samples varied from 216 to 228 (for K, WC, BD, FI) to 691 for the other soil properties. At each randomly selected point we collected from a few to a dozen samples around, which we used in the geostatistical analysis.
Soil dry bulk density (Mg m −3 ) was determined with the gravimetric method from the ratio of the mass of soil dried at 105 °C to the soil volume of 100 cm 339 . The gravimetric water content (WC grav ) was determined from the ratio between water mass and mass of the soil after drying using the same cylinders as for determination of dry bulk density. Soil water content was also measured using a TDR meter (WC TDR ) close to the sampling sites. Saturated hydraulic conductivity was measured with the constant head method in soil samples with a volume of 100 cm 3 in a laboratory permeameter (Eijkelkamp Agrisearch Equipments, The Netherlands) 40 . The method assumes that K does not depend on the hydraulic head difference and the soil is the only resistance to water flow. Particle size distribution was analysed using the sieving and hydrometer method 41 . Soil organic matter was determined based on wet oxidation with K 2 Cr 2 O 7 according to Tiurin's procedure 41 . The soil pH (in H 2 O) was determined potentiometrically using a composite electrode. Particle density (Mg m −3 ) was determined with the pycnometric method 42 . The total porosity (m 3 m −3 ) was calculated as a ratio of 1 − bulk density/particle density 43 . Data analysis. Classical statistics. Basic statistics with the mean, standard deviation, coefficient of variation, minimum, maximum, skewness, and kurtosis were calculated for each soil property. Both kurtosis and skewness values of 0 indicate in general symmetrical distribution with a similar right tail (positive) and left tail (negative) of the distribution curve. When one tail is longer than the other, the distribution is asymmetric. As shown by Dahiya et al. 44 , the variability of the soil properties was categorised as low (0-15%), medium (15-75%), and high (> 75%). Pearson correlation coefficients between the soil variables were determined. The results were analysed using STATISTICA 12 PL (StatSoft 2019) and GS + 10 45 . Data normality was assessed using the Cumulative Frequency Distribution 45 . If the distribution was not normal the natural logarithmic and root square transformations were used to ensure the normal or close to the normal distribution. www.nature.com/scientificreports/ Geostatistical methods. Semivariograms and cross-semivariograms. It is assumed that values of soil properties or physical quantities measured in a given point are dependent (similar, correlated) on other points. In this approach, similarity is described by half of the value expected from the difference between the value of Z(x) of the variable at point x and the value of Z(x + h) at a point distant by vector h. A variable whose values correspond to values Z(x) is "regionalised". This variable has a random aspect, which takes into account local anomalies, including a structural aspect reflecting the multiscale trend of the phenomenon (trend). The analysis of this variable consists of identifying the structure of the variation. Three phases of the analysis can be distinguished: preliminary examination of collected data and evaluation of basic statistics, calculation of the empirical variant of the considered regionalised variable, and adjustment of the mathematical model to the course of the empirical variant. The knowledge of the first two statistical moments of random functions is required: the first (expected value) and the second (variance) moment. It is also required that the examined process is stationary, i.e. it does not change its properties when changing the beginning of the time scale. In the case of fulfilment of the stationary process, the random function Z(x) is defined as the second order stationary and then the expected value exists and does not depend on the position, and the experimental semivariogram γ(h) (for a single variable z 1 ) or the cross-semivariogram γ 12 (h) (for two variables z 1 and z 2 ) for distance h are calculated from the following equations 45 : Three characteristic parameters for the semivariograms and cross-semivariograms are distinguished: nugget effect, sill, and range. When the value of the semivariograms increases not from zero but from a certain value, this value is called the nugget effect. It expresses the variability of the examined physical quantity at a scale smaller than the sampling interval and/or accuracy of measurement. A value at which no further increase in the semivariograms is observed (approximately equal to the sample variance) is called sill, while the distance from zero to the point where the semi-or cross-semivariogram reaches 95% of the sill value is called a range. The latter expresses the greatest distance at which the sampled values are auto-or cross-correlated.
For semi-and cross-semivariograms determined empirically, the following mathematical models were fitted using the least squares method 45 : • The linear isotropic model describes a straight line variogram. There is no sill in this model; the range A 0 is defined arbitrarily to be the distance interval for the last lag class in the variogram. The formula used is: • The exponential isotropic model. The formula used for this model is: • The Gaussian isotropic model. The formula used for this model is: In the case of the linear model, there is no effective range A-it is set initially to the separation distance (h) for the last lag class graphed in the variogram. In the case of the spherical model, the effective range A = A 0 . In the case of the exponential model, the effective range A = 3A 0 , which is the distance at which the sill (C + C 0 ) is within 5% of the asymptote. In the case of the Gaussian model, the effective range A = 3 0.5 A 0 , which is the distance at which the sill (C + C 0 ) is within 5% of the asymptote.
The fractal dimension D was determined based on the log-log semivariogram plots using the formula 46 : where: H is the slope of the semivariogram line plotted in the logarithmic system of coordinates.
Kriging. The estimation of values in unmeasured places was conducted using the kriging estimation method. This method gives the best unbiased estimate of the point or block values of the variable under study with minimal variance during the estimation process. The values of the kriging variance depend on the position of the samples in relation to the estimated location, the weights assigned to the samples, and the parameters of the semivariogram model and is described by a linear equation expressed by the formula 45 : where N is the number of measurements, z(x i ) is the measured value at the point x i , z*(x o ) is the value estimated at the estimation point x o , and λ i are the weights. The weights are determined from the system of equations taking into account the condition of non-loadability and efficiency of the estimator, i.e. when the expected value of the difference between the measured and estimated values is zero and the variance of the differences is minimal 45 : Solving the above system of equations, we determined the weights of kriging-λ i . These weights allow also determination of the estimated value z* and its variance from the formula: Ordinary kriging (OK) was used for the estimation, as it gave a good match between the measured and the estimated value. The inverse distance weighting interpolation (IDW) was a worse interpolator, while the ordinary www.nature.com/scientificreports/ cokriging (OCK) did not appreciably improve the estimation compared to OK. Therefore, kriging maps were created using Gamma Design Software GS + 10 45 .

Results
Classical statistics. The statistical parameters of the examined soil characteristics in the studied region are summarised in Table 1 44 (1984), the variability was low for soil bulk density, total porosity, sand content, and pH in H 2 O (CV 8.8-13.8%), medium for CEC, silt content, OC, WC TDR , and WC grav (32.1-67.4%), and high for saturated hydraulic conductivity (77.3%). Skewness, which characterises the degree of asymmetry of the distribution around the mean, was moderate (< 1) for most variables and slightly more positive (< 2) for soil moisture, clay content, and saturated hydraulic conductivity. Silt content and bulk density showed a slight negative asymmetry of (< -1). Kurtosis, which characterises the relative slenderness or flatness of the distribution compared to the normal distribution (zero), was close to zero for most variables. We noted relatively little flattening for sand, silt, and pH (in H 2 O) (< 0 or from − 0.130 to − 0.158), slight slenderness for bulk density and porosity (< 0.053 or from 0.030 to − 0.030), and somewhat higher value for CEC (< 1 or 0.956). Soil moisture, the OC and clay contents, and saturated hydraulic conductivity showed much greater slenderness of distribution (2.218-4.826). The differences between the mean values and the medians for individual variables as well as the values of asymmetry and kurtosis indicate that the studied variables can be described with a normal distribution with fairly good accuracy. Those with greater asymmetry were square-root or natural-logarithm transformed, thus their data distributions were close to the normal distribution (Table 1).
Correlation analysis. The linear correlation coefficients (r) between the considered soil properties are summarised in Table 2 (the values marked in bold are statistically significant at p < 0.05). The saturated water conductivity of the soil was significantly positively correlated with the porosity (0.300) and negatively with the soil density (-0.283). Other significant correlation coefficients were found between sand and silt contents (-0.996) and sand and clay (-0.182). There was no significant correlation between the contents of silt and clay. Soil pH  and positively with the silt content (0.178, 0.168, respectively). CEC was negatively and significantly correlated (p < 0.05) with the sand content (-0.519) and positively correlated with silt, clay, pH, and OC (0.160-0.607). The significant correlation between the gravimetric vs. TDR soil moistures (0.876) indicates suitability of the TDR measurement system, which is widely used as a benchmark for validation of satellite soil moisture products (e.g. 47 ). Soil moisture did not significantly correlate with other soil properties. Soil porosity correlated negatively (p < 0.05) with the bulk density and sand and silt contents (-0.737, − 0.142) and positively with the sand content (0.145).
Geostatistical analysis. The fitted semivariogram models for K and cross-semivariograms for pairs of cross-correlated K and other soil properties are presented in Table 3. In general, there was a good agreement between the theoretical exponential models for all soil properties and the empirical semivariograms, as indicated by the high values of the determination coefficients (R 2 from 0.592 to 0.923) and the sum of squared residuals (RSS) from < 10 −6 to 81.4 depending on soil properties. This agreement for the cross-semivariograms was fairly good in six cases (R 2 > 0.284), and poor in two (R 2 ~ 0.02). The RSS values were small for most models (5.94 × 10 −3 -7.62 × 10 −6 ). In the cross-semivariograms analysis, five soil properties had exponential dependency, four-Gaussian, and one-linear. The presence of nugget effects indicates that the variability of the examined features is smaller than the adopted minimum distance between the measurement points. The sill values of the semivariance are comparable with the values of variance obtained in the classical way (Tables 1, 3), which may indicate lack of clear trends in the data. The sill values of the semivariograms were a derivative of the content of individual textural fractions. The highest sill values were recorded for the sand and silt contents. However, they were lower for saturated hydraulic conductivity and substantially lower for the contents of clay, organic matter, moisture, and pH, cation exchange capacity, bulk density, and total porosity. The range of spatial dependence displayed by the semivariograms was the smallest for pH (0.012°), intermediate for OC, CEC, clay, sand, silt, BD, WC grav , and WC TDR porosity (0.018-0.057°), and the largest for saturated hydraulic conductivity (0.081°). In the case of the cross-semivariograms for pairs of cross-correlated K with intrinsic and relatively stable properties (sand, silt, clay, OC, CEC, pH), the spatial ranges were much larger (from 0.095° to 0.81°) than with dynamic ones, such as gravimetric and TDR soil moistures, bulk density, and total porosity (0.003°-0.009°). According to the classification of Cambardella et al. 48 , the spatial dependences (nugget/sill) for all semivariograms were moderate (0.25-0.75) and those for cross-semivariograms were in general strong (< 0.25). The distribution of the most widely studied soil properties showed anisotropy with orientation mostly from west to east. Only the clay content and CEC showed anisotropy from north to south.
The estimation of the spatial distribution of the studied properties using the fractal theory showed that all soil properties were characterised by high values of the fractal dimensions D > 1.9, indicating a random distribution in the spatial organisation.
Kriging maps. Based on the obtained models, semivariogram parameters, and measured data, regionalscale maps of the soil properties were generated and the estimation errors were calculated (Fig. 2) using ordinary kriging. The estimation errors were 1-2% in the vicinity of the measurement points and up to approx. 10% at the edges of the estimated areas (measurement grids).
In the northern part of the study region, the large island (area) between approximately 52.01° and 52.04° with high saturated hydraulic conductivities from 3.1 to 4.6 m day −1 (Fig. 2) corresponds with the highest sand contents (> 74%) and the lowest contents of silt (< 22%). In the southern part of the region below 52°, the lower  www.nature.com/scientificreports/ saturated conductivities correspond with the lower sand content (< 74%), greater silt content > 24%, and similar clay content (< 2.8%). In general, latitudinal distribution can be observed for the sand and silt contents. Clay, CEC, and BD are distributed in a small island and OC, pH, and gravimetric and TDR soil moistures-in larger islands. Cation exchange capacity (CEC) showed a relatively uniform distribution throughout the study region in small islands of higher or lower values. Organic carbon content (OC) had mostly an island distribution system with a slightly marked meridional distribution. It can be observed that both gravimetric and TDR soil moistures were more similar to that of the organic carbon content than other soil properties. The gravimetric vs. TDR soil moisture distribution was more variable, which implies greater sensitivity. The highest BD values correspond with the highest sand content and the lowest contents of clay and OC.

Discussion
Geostatistical analysis. K exhibited the strongest spatial heterogeneity (coefficient of variation 77.3%) of all the soil studied properties. This can be largely influenced by various management practices, including crop rotation, tillage, liming, or history of land use and other treatments used at the regional scale. These management effects may be associated with changes in total porosity, which was significantly positively correlated with K (p < 0.05) ( Table 2). Particularly important are large and connected (elongated) biological and inter-aggregate pores 9,16,17 . The spatially heterogeneous effect of the various management practices applied at the regional scale on K can mask the impact of intrinsic soil texture, as indicated by the poor and non-significant overall correlations between K and the contents of all textural fractions (Table 2). Probably, there was no significant effect of terrain attributes associated with the landscape level and pedogenetic processes on K, as the flat study area within the Łuków Plain is covered mainly by Podzol soils 49,50 .
The model parameters describing the spatial relationships of K with intrinsic and dynamic soil properties were better for the cross-semivariograms with auxiliary soil properties than for the direct semivariograms, as shown by the appreciably smaller nugget values (C 0 ), the greater ranges (A), and the stronger degree of spatial dependencies (nugget/sill) in the former. These differences were most pronounced when intrinsic (stable) and relatively stable soil properties, including the contents of textural fractions, OC, and CEC, were used as auxiliary variables (C 0 − 0.165 to 0.170 vs. 2.23), A 0.156° (15.6 km) to 0.810° (81 km) vs. 0.081° (8.1 km), nugget/sill 0.001 Table 3. Fitted semivariogram models (SV) for properties data used in the ordinary kriging interpolation method and cross-semivariogram models (CSV) between saturated hydraulic conductivity and other soil properties; 1° corresponds to approx. 100 km. Exp.-exponential model, Lin.-linear model, Gau.-Gaussian model, C 0 -nugget variance, C 0 + C-sill, A-effective range (°) (1° = approx. 100 km), R 2 -coefficient of determination, RSS-root sum square, A z -anisotropy (°-slope angle), D0-fractal dimension.

Properties
Model C 0 C 0 + C C 0 /(C 0 + C) A (°) R 2 RSS A z (°) D0 www.nature.com/scientificreports/ to 0.223 vs. 0.466). This indicates that the intrinsic and relatively stable soil properties were spatially correlated with K although there was no significant linear correlation of each intrinsic or relatively stable soil property vs. K ( Table 2). The smaller nugget values (C 0 ) in the cross-semivariograms compared to the direct semivariograms imply smoother spatial continuity and stronger dependency between neighbouring sampling points 31,51,52 . It is worth noting that the suitability of soil texture data used as auxiliary variables for improvement of the prediction of the spatial K distribution can be enhanced by their worldwide availability in soil geographic databases (e.g. 53 ). It should be underlined that the range values of the cross-semivariograms (in the case of all pairs) exceeded the length and width of the study region (~ 13 × 16 km).
Kriging maps. The kriging maps generated in this study allowed outlining two sub-areas with predominantly saturated hydraulic conductivity (K) > 3.0 m day −1 in the northern part (latitude 52.01-52.06) and < 3.0 m day −1 in the southern part (52.00-51.93) of the study region. As reported by Stryjewski 54 , the K values in the northern part can be classified as high and very high and those in the southern part as fairly high and low. The comparison of the maps in Fig. 2 shows positional similarity between the sub-area with the higher K value and those with the large sand content (> 74%) and the low silt content (< 22%). This similarity can be attributed to the effect of the sand fraction on the abundance of relatively large and connected pores that mostly contribute to high K (e.g. 55 ). This effect can be illustrated by results from a study conducted by Lim et al. 56 , where K of 5.98 m day −1 of coarse sand decreased by 57, 88, and 96% with the successively decreasing sand content in fine sand, loam, and clay textured soils. Our previous studies in the same region along with visual observations showed that limited crop growth and yields were spatially related to higher sand content 50,57 . This crop response in sandier and more permeable zones can be caused by excessive drainage of rainwater resulting in insufficient plant-available water for unsaturated conditions. Furthermore, the drainage contributes to chemical leaching, thereby limiting the availability of nutrients for plants. This explanation can be supported by the significant negative correlation between the sand content and cation exchange capacity ( Table 2). This implies that high K can be an indicator of a low-yielding zone in the studied area with predominance of coarse-textured soils. This is in contrast to fine-textured soils where low K values are indicative of low-yielding zones. For example, in a study conducted by Keller et al. 3 on loam and clay soils with K varying from 0.6 to 25.2 m day −1 , lower saturated hydraulic conductivity was recorded in low-yielding zones than in high-and medium-yielding zones due to the more blocky soil structure in the former. The low yields in fine-textured soils with low saturated hydraulic conductivity often results from water ponding and limited oxygen concentrations for root and shoot growth, especially in wet years 13 . This indicates that the effect of spatial distribution of K on the spatial distribution of soil productivity and other soil functions depends on soil texture and weather conditions during the growing season. Therefore, different threshold K values with respect to productivity and other soil functions should be considered in the case of coarse-and fine-textured soils. The kriging map of K can be useful for the local authorities and agronomy advisers for spatial planning of management practices aimed at reduction of particularly high soil K. In coarse-textured soils, it can be reduced by addition of exogenous organic materials including biochar (as a soil conditioner) 27,56,58 and locally available recycled composted chicken manure or spent mushroom substrate after mushroom harvesting 59,60 . Such organic materials are currently applied in the studied area in research and by some farmers on agriculturally used sandy soils.
Another important agricultural management practice influencing K is crop rotation, including green manure cover crops or intercropping systems 61,62 . This practice is subsidised in several countries, including Poland, to increase soil organic carbon content 63 , improve soil structure 64 , protect the soil surface from raindrop impact 65 , enhance fixation of atmospheric nitrogen in the case of legumes 66 , and improve agricultural productivity 62 . Over a longer time span, conversion of arable land into grassland that can serve as carbon and water storage may be an efficient option (e.g. 67 ).
Also re-compaction of loose soil by traffic leads to reduced K due to a decrease and increase in the large and small pore volumes, respectively 68 . However, this practice needs to be applied with caution to avoid excessive soil compaction and its harmful effect on root growth and crop yield 13,69 . Saturated hydraulic conductivity values ≤ 0.1 m day −1 are used as an indicator of poor soil structure 70 and more recently as threshold values of excessive soil compaction induced by vehicular traffic 71 . The K values in the study region were in general above the thresholds, which may be in part related to the presence of predominantly small farms where relatively light agricultural vehicles and implements are used.

Summary and conclusions
The saturated hydraulic conductivity (K) of the soils in the studied region (140 km 2 ) varied from 0.01 to 11.54 m day −1 and exhibited high spatial variability (CV 77.3%). This variability was higher than that of the contents of textural fractions and organic carbon, cation exchange capacity, soil water content, bulk density, and total porosity (CV 8.9-67.4%) (67.4% for grav. soil moisture and 53.7% for TDR soil moisture). K was significantly (p < 0.05) positively correlated with the total porosity (r = 0.300) and negatively correlated with the soil bulk density (r = -0.283). The spatial prediction (autocorrelation) of the soil properties using single (direct) semivariograms varied from 0.012° (1.2 km) for pH to 0.081° (8.1 km) for K. Areas with larger K values are mainly determined by the proportion of sand and silt and those with smaller K-by the proportion of clay. The range of spatial dependence in the cross-semivariograms between K and sand and silt contents used as secondary variables was smaller (about 15 km) than that of cation exchange capacity, clay content, and organic carbon content used as secondary variables (about 22 to 81 km). However, the range of the spatial dependence prediction in the cross-semivariograms decreased when the dynamic soil properties (soil moisture, bulk density or total porosity) were used as secondary variables (0.3-0.9 km). The suitability of soil texture and organic carbon data www.nature.com/scientificreports/ to be used as secondary variables in cross-semivariograms for predicting spatial cross-correlations of K can be enhanced by their worldwide availability in soil geographic databases. The kriging maps allowed outlining two sub-areas with predominance of K > 3.0 m day −1 in the northern part (latitude 52.01-52.06) and < 3.0 m day −1 in the southern part of the study region. The comparison of the spatial maps indicates that there is positional similarity (agreement) between the sub-areas with the highest K values and the largest sand contents (> 74%). The spatial maps generated in this study can be helpful for the local authorities and agronomy advisers for spatial planning of management practices aimed at reduction of K of permeable and low-productive soils.