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

Saturated hydraulic conductivity (SHC) is a key property for evaluating soil water movement and quality. Most studies on spatial variability of SHC have been performed soil at a field or smaller scale. Therefore, the aim of this work was to assess (quantify) the spatial distribution of SHC at the commune scale and its relationship with other soil properties, including intrinsic sand, silt, and clay contents, relatively stable organic carbon, cation exchange capacity (CEC), dynamic 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 commune (140 km 2 ) 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–228 (for SHC, 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, SHC 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 mean SHC was 2.597 m day –1 and ranged from 0.01 up to 11.54 m day –1 . The spatial autocorrelation (range) of SHC in the single (direct) semivariograms was 0.081° hand, the ranges of the spatial relationship between SHC and the intrinsic and relatively stable soil properties were much larger (from ~15 to 81 km) than between SHC and the dynamic soil properties (0.3-0.9 km). This knowledge is supportive for making decisions related to land management aimed at reduction of hydraulic conductivity and chemical leaching and improvement of soil water resources and crop productivity.

hand, the ranges of the spatial relationship between SHC and the intrinsic and relatively stable soil properties were much larger (from ~15 to 81 km) than between SHC and the dynamic soil properties (0.3-0.9 km). This knowledge is supportive for making decisions related to land management aimed at reduction of hydraulic conductivity and chemical leaching and improvement of soil water resources and crop productivity.
Keywords: saturated hydraulic conductivity, intrinsic and dynamic soil properties, communescale variability, geostatistics, kriging maps 1. Introduction Saturated hydraulic conductivity (SHC) determines the maximum capacity of soil to transmit water, pathways of water movement partitioning precipitation and irrigation water into surface runoff and retention in the soil 1,2 , and the soil water dynamics in the soil profile 3 .
High SHC leads to rapid water infiltration and drainage 4,5 and reduced time for attenuation of dissolved agrochemicals before entering ground waters 6 , whereas low SHC increases surface runoff and erosion 7,8 . Thereby, SHC helps farmers to apply an appropriate amount of irrigation water 9 . Furthermore, SHC affects soil aeration capacity 10,11 influencing nutrient transformations and uptake by plants 12,13 . Due to the numerous contributions, SHC is often used as a measure of soil physical quality (e.g. 12 ). Also, it is a key parameter in mathematical models for predicting soil hydraulic behaviour 2, 14 , 15.
The SHC 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 SHC.
Owing to the high sensitivity to pore size distribution influenced by soil texture and management practices, SHC displays relatively high spatial variability 9,15,20,22,23 . Therefore, knowledge of the spatial distribution of SHC 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 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 description of soil properties 31,32 . Once various variables are linked, their combined spatial designs can be evaluated by crosssemivariograms. Cross-semivariogram data and maps obtained using the cokriging technique allow distinguishing time-consuming and/or expensive variables from those that are more easily measured or available in soil databases. When SHC shows spatial random variation, the use of both classical statistics and geostatistical stochastic models is recommended 30 .
Numerous studies on the spatial variability of soil SHC 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 scale is poorly understood, as suitable spatial characterisation of highly heterogeneous SHC requires a large number of laborious, time-consuming, and expensive direct measurements 9,36,37 . To overcome these constraints, we have estimated the spatial distribution of SHC at the commune 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.

Study area and sampling
Trzebieszów commune with an area of about 140 km 2 is situated in a flat area within Łuków Plain, Podlasie region, south-eastern Poland. The height differences in the shallow and often wet river (Krzna) valleys do not exceed 10 m. The commune has mostly low productive Podzol soils 38 derived from sandy and sandy loams of glacial origin. About 80% of the commune area is used in agriculture, with 62.3% and 18.2% of arable lands and grasslands, respectively. Forests cover only 13.5% of the commune 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 of air temperature calculated from the monthly averages of January and July is 21.3°C, and 23.4°C when calculated from the differences of the average temperature of the hottest and coldest month 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 a 0-20 cm layer into cloth bags and 100 cm 3 Kopecky cylinders immediately after harvesting cereals (August). The number of individual samples varied from 216-228 (for SHC, WC, BD, FI) to 691 for the other soil properties. Figure 1 shows the spatial distribution of the soil sampling points. 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 3 39 . The gravimetric water content (WCgrav) 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 (WCTDR) 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 . Particle size distribution was analysed using the sieving and hydrometer method 41 . Soil organic matter was determined based on wet oxidation with K2Cr2O7 according to Tiurin's procedure 41 . The soil pH (in H2O) 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 .

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 .

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 in 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 or spatial 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 z1) or the cross-semivariogram 12 (ℎ) (for two variables z1 and z2) for distance h are calculated from the following equations 45 : 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 A0 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 Ait 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 = A0. In the case of the exponential model, the effective range A = 3A0, which is the distance at which the sill (C + C0) is within 5% of the asymptote. In the case of the Gaussian model, the effective range A = 3 0.5 A0, which is the distance at which the sill (C + C0) 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 where N is the number of measurements, z(xi) is the measured value at the point xi, z*(xo) is the value estimated at the estimation point xo, 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: Regular 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 regular cokriging (OCK) did not appreciably improve the estimation compared to OK.

Classical statistics
The statistical parameters of the examined soil characteristics in the studied commune are summarised in Table 1 Mg m -3 and 0.424, 0.308, 0.524 m 3 m -3 , respectively. As in the study conducted by Dahiya et al. 44 (1984), the variability was low for soil bulk density, total porosity, sand content, and pH in H2O (CV 8.8-13.8%), medium for CEC, silt content, OC, WCTDR, and WCgrav (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 H2O) (<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 squareroot 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

Geostatistical analysis
The fitted semivariogram models for SHC and cross-semivariograms for pairs of crosscorrelated SHC 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 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. 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.

Kriging maps
Based on the obtained models, semivariogram parameters, and measured data, commune scale 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 commune, 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 commune below 52°, the lower 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 moisturesin larger islands. Cation exchange capacity (CEC) showed a relatively uniform distribution throughout the commune 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.  Table   2). The smaller nugget values (C0) 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 SHC 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 commune (~13×16 km).

Kriging maps
The kriged maps generated in this study allowed outlining two sub-areas with predominantly  Figure 2 shows positional similarity between the sub-area with the higher SHC 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 SHC (e.g. 55 ). This effect can be illustrated by results from a study conducted by Lim et al. 56  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 SHC 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 SHC values with respect to productivity and other soil functions should be considered in the case of coarse-and fine-textured soils.
The kriged map of SHC can be useful for the commune authorities and agronomy advisers for spatial planning of management practices aimed at reduction of particularly high soil SHC. In coarse-textured soils, it can be reduced by addition of exogenous organic materials including biochar (as a soil conditioner) 27,56,58Error! Reference source not found. 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 commune area in research and by some farmers on agriculturally used sandy soils.
Another important agricultural management practice influencing SHC 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 SHC 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 SHC values in the commune area 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 (SHC) of the soils in the studied commune (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). SHC 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° 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 to be used as secondary variables in cross-semivariograms for predicting spatial cross-correlations of SHC can be enhanced by their worldwide availability in soil geographic databases. The kriged maps allowed outlining two sub-areas with predominance of SHC >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 commune. The comparison of the spatial maps indicates that there is positional similarity (agreement) between the sub-areas with the highest SHC values and the largest sand contents (>74%). The spatial maps generated in this study can be helpful for the commune authorities and agronomy advisers for spatial planning of management practices aimed at reduction of SHC of permeable and low-productive soils.