Full Terms & Conditions of access and use can be found at
https://www.tandfonline.com/action/journalInformation?journalCode=gags20
Archives of Agronomy and Soil Science
ISSN: 0365-0340 (Print) 1476-3567 (Online) Journal homepage: https://www.tandfonline.com/loi/gags20
Geostatistical approach as a tool for estimation of
field capacity and permanent wilting point in
semi-arid terrestrial ecosystem
Tülay Tunçay, Oğuz Başkan, Ilhami Bayramın, Orhan Dengız & Şeref Kılıç
To cite this article: Tülay Tunçay, Oğuz Başkan, Ilhami Bayramın, Orhan Dengız & Şeref Kılıç (2018) Geostatistical approach as a tool for estimation of field capacity and permanent wilting point in semi-arid terrestrial ecosystem, Archives of Agronomy and Soil Science, 64:9, 1240-1253, DOI: 10.1080/03650340.2017.1422081To link to this article: https://doi.org/10.1080/03650340.2017.1422081
Accepted author version posted online: 29 Dec 2017.
Published online: 20 Jan 2018. Submit your article to this journal
Article views: 180
View related articles
ARTICLE
Geostatistical approach as a tool for estimation of
field capacity
and permanent wilting point in semi-arid terrestrial ecosystem
Tülay Tunçaya, Oğuz Başkana, Ilhami Bayramınb, Orhan Dengızcand Şeref Kılıçd
aSoil Fertilizer and Water Resources Central Research Institute, Ankara, Turkey;bAgricultural Faculty Department
of Soil Science and Plant Nutrition, Ankara University, Ankara, Turkey;cAgricultural Faculty, Department of Soil
Science and Plant Nutrition, Ondokuz Mayis University, Samsun, Turkey;dEngineering Faculty Department of
Environmental Engineering, Ardahan University, Ardahan, Turkey
ABSTRACT
The understanding of spatial and seasonal variability in soil water retention properties is critical for careful management of soil water in agricultural production in semi-arid regions. The main objectives of this study were to spatially predict and prepare distribution maps of the field capacity and permanent wilting point in semi-arid terrestrial eco-system region by using regression kriging and cokriging, and kriging models in order to predict soil water budget with limited data. Capability of the models was compared, including the use of descrip-tive statistics, semivariogram and cross-semivariogram. For this aim, 287 disturbed and 167 undisturbed soil samples (0 – 20 cm) were collected in grid system and used for this research. Overall results showed that regression kriging produced smaller mean absolute error (MAE) and mean square error (MSE) than kriging and cokriging forfield capacity, but cokriging was superior for the interpolation of permanent wilting point in terms of smaller MAE and MSE obtained at the sam-pling point. Consequently, the final maps and calculated results can also be used in the decision processes for land and water manage-ments and soil conservation practices by authorities, as well as by farmers for irrigatedfields in semi-arid areas.
ARTICLE HISTORY
Received 9 August 2017 Accepted 22 December 2017
KEYWORDS
Field capacity; permanent wilting point; regression kriging; cokriging; geostatistics
Introduction
The characterization of spatial variability is crucial to both the determination and monitoring of soil hydraulic properties and processes, especially in large catchments with limited data available to conservation planners and decision makers pursuing sustainable natural resources use. The determi-nation of the plant available water as a difference between the field capacity (FC) and permanent wilting point (PWP) to characterize soil water retention properties is relevant to ecological processes, environmental issues and particularly, hydrological models and management of irrigation practices. The understanding of spatial and seasonal variability in soil water retention properties is critical for precise water inputs into agricultural production in semi-arid regions (Pachepsky and Rawls2004; Wösten et al.2012). However, it is impractical to measure FC and PWP continuously at each sampling point across the watershed or agriculturalfields due to the expensive and time-consuming nature of the process. For this reason, scientists have successfully developed interpolation models for estimation of soil hydraulic characteristics by using easily available soil data such as bulk density, clay, sand and silt content, organic matter etc. Studies have showed that these auxiliary variables used for
CONTACTOrhan Dengız odengiz@omu.edu.tr Agricultural Faculty, Department of Soil Science and Plant Nutrition, Ondokuz Mayis University, Samsun, Turkey
VOL. 64, NO. 9, 1240–1253
https://doi.org/10.1080/03650340.2017.1422081
interpolation were significant for the prediction of soil properties. Furthermore, many scientists have reported that some interpolation techniques, such as cokriging and regression kriging which use auxiliary variables, enhance the interpolation accuracy of the relationship between the variable of interest and the auxiliary variable (Bishop and McBratney2001; Khosla et al.2002; Mzuku et al.2005; Liu et al.2006). Moreover, soil water retention properties reflect interaction among the various soil properties such as bulk density and particle size distribution. That resulted into indirect estimation of soil water retention characteristics by pedotransfer functions, which are also widely developed and used (i.e. Wösten et al.2001; Minasny and Hartemink2011; Miháliková et al.2014). Among researchers in environmental science, geography, ecology and other earth sciences, spatial interpolation is one of the best known techniques for defining the regional variability of soil characteristics (Anderson2001; Meng et al.2010). Numerous indirect interpolation methods are available to estimate soil hydraulic properties such as ordinary kriging (Bishop and McBratney 2001; Igbal et al. 2005; Lloyd 2005), cokriging (Lark2000; Webster and Oliver2001) and multiple regression kriging (hybrid model).
Ordinary kriging is a geostatistical interpolation technique that is described by the acronym BLUE– ‘best linear unbiased estimator’. It is the ‘best’ because it aims at minimizing the variance of errors,‘linear’ because its estimates are weighted linear combinations of the available data, and ‘unbiased’ since it tries to have the mean residual or error equal to zero (Isaaks and Srivastava1989; Mesic2016). Besides, cokriging involves correlation between the primary and auxiliary variables for optimal estimation of the variable of interest (Goovaerts 2000; Adhikary et al. 2017). Regression kriging is a spatial prediction technique that combines the regression of the dependent variable on auxiliary variables with kriging of the regression residuals. It is mathematically defined as universal kriging and kriging with external drift (McBratney et al.2000; Hengl et al.2004; Hengl 2007). It is used to explore the spatial correlation between the response variable and the non-spatial correla-tion between the dependent variable and auxiliary variables (Hengl et al.2004,2007; Sullivan et al.
2005; Miháliková et al.2015,2016), especially in soil science.
Applied effectively in other disciplines over recent decades, the combination of kriging and hybrid interpolation techniques for the evaluation of the spatial variability of soil properties is partly new in the soil science community because regression kriging combined classical statistical technique with geostatistical methods. The accurate estimation of the spatial variability of soil hydraulic or soil water retention properties is essential in the hydrological modelling of watersheds in semi-arid region so that policy makers can develop plans and strategies, particularly at the broad landscape level.
The careful management of soil water is crucial for sustainable agricultural production in arid and semi-arid regions. Determination of soil hydraulic properties in large catchments has been used in some modelling studies (ordinary kriging, hybrid techniques) with limited data; moreover, it is impossible to measure FC and PWP continuously on a large scale due to high costs and the amount of time and resources it would consume. Given that scenario, soil texture is an important parameter in predicting soil hydraulics in arid and semi-arid areas. Changes of the soil texture e.g. variable sand content, affect interpolation techniques.
Against this background, the specific objectives of the current research in a semi- arid terrestrial ecosystem region were to: i) assess the performance of selected interpolation techniques by comparing measured and predicted soil properties for the study area and ii) prepare maps of the spatial distribution of FC and PWP. Thus, it is substantial to detect the best performing model among different estimation methods for prediction of soil water retention properties and soil water budget at large scale in terrestrial ecosystem of the semi-arid areas.
Material and methods Description of the study area
This study was carried out on an area of approximately 64 km2(8 km × 8 km) on the Altınova State Farm (29,608.6 ha, located between eastern longitudes 31° 39ı 20ıı – 32° 49ı 55ıı and northern
latitudes 58° 39ı20ıı- 51° 41ı54ıı(Figure 1) in the middle of Central Anatolia, Turkey. The Altınova State Farm is located approximately 189 km far from Ankara and about 126 km far from Konya in the Great Konya Basin.
According to data from the State Farm’s meteorological station over the 1990–2011 period, the mean annual precipitation was 302.3 mm, the total potential evapotranspiration was 1,296 mm and the mean annual temperature was 12.8°C. Soil temperature and soil moisture regimes at the study site were classified according to Soil Survey Staff (2015) as mesic and aridic, respectively.
The pedological and geological properties of the soils and sediments, and the formation and diagenesis of carbonates under the lacustrine environment in the Great Konya Basin, have been investigated (Driessen1970; Vergouwen1981). The ground base of the study area was composed of magmatic rocks, including gabro, diyabaz, diyolite and serpentine. The deposition of marly
limestone bands, clay marl and gravel occurred across the study area in the Neogene period and that offluvial sand and gravel deposits occurred in the Qauternary period (TİGEM1997).
The study area, Altınova State Farm, is used for agriculture: 100% of the total area is covered by crops, principally wheat, barley, triticale, alfalfa and corn. The most common crop of the study area is wheat, followed by barley and triticale. The fallow rotation system of farming is employed to sustain crop production across the study area.
According to soil classification system of WRB (2014), ten different soil series belonging to the Vertisols, Cambisols and Calcisols great groups were found in the study area and classified into five subgroups (Calcic Vertisol, Vertic Cambisol, Haplic Calcisol, Aric Cambisol and Petric Calcisol). In addition, the Petric Calcisol (42.7%), Aric Cambisol (23.9%) and Haplic Calcisol (21.3%) subgroups account for approximately 88% of the soils. The most prevalent soil series were the Odabaşı (26.6%), Altınova (15.7%), and Kap (12.3%) series. Odabaşı and Kap soil series were classified as Petric Calcisol while Altınova series was classified as Calcic Vertisol in the study area.
Soil sampling and laboratory analysis
Surface soil samples were collected at 287 locations at 500 m spacing at the test site. This part of the study area includes most of the soil of the whole state farm territory. A total of 287 disturbed and 167 undisturbed soil samples (100 cm3) were collected with a core sampler to investigate their physical and chemical characteristics, including pH, electrical conductivity (EC), organic matter (OM), lime content and soil texture, and hydrophysical characteristics such as bulk density, field capacity (FC) and permanent wilting point (PWP), from 0–20 cm soil depth. Both disturbed and undisturbed soil samples were taken from each point of grid system infield (Figure 2). The disturbed soil samples were air-dried and passed through a 2 mm sieve.
The pH and EC (mS m−1) were determined from a prepared mixture of 1:2 soil to distilled water with a pH electrode (Soil Survey Laboratory 2004), the lime content (CaCO3,%) with a Sheibler
calcimeter (Soil Survey Laboratory Methods Manual2004), and the organic matter content (OM, %) with the Walkey-Black method (Richards1954). The sand (S, 0.05–2 mm), silt (Si, 0.002–0.05 mm), and clay (C, <0.002 mm) contents (%) of the soil samples were determined with the hydrometer method (Bouyoucos1951), based on the USDA texture categories. Additionally, the FC (0.033 MPa) and the PWP (1.5 MPa) were determined with the porous plate method (Klute1986). The soil dry bulk density, (BD), Mg m−3, was also determined by dividing the total dry mass of the sample by the volume of the soil core (Blake and Hartge1986) for undisturbed samples.
Descriptive statistics, correlation and regression analysis
Descriptive statistics of soil properties, including minimum, maximum, mean, standard deviation, coefficient of variation, skewness and kurtosis, were calculated. The relationship between FC and PWP and other soil properties was performed with simple correlation analysis; parameters that showed a high correlation (p < 0.05) were selected as potential independent variables to be included in the stepwise regression and geostatistical analysis. Kolmogorov-Smirnov (K-S) tests of normality were applied to the measured soil characteristics, prior to ordinary kriging, cokriging and regression kriging. Moreover, prior to regression analysis, multicollinearity and autocorrelation were checked among the selected variables.
Geostatistical analysis
Geostatistical techniques use the spatial dependence between observed samples to predict values at unsampled points. The semivariogram model that expresses spatial dependence between observed values separated by h distance has been proved to be a useful basic tool for the interpolation of many soil properties. The validity of the interpolation methods and accuracy of
the models used in this study was evaluated by cross validating the observed values Z(Ui) and estimated values Z*(Ui) in terms of mean absolute error (MAE) and mean-squared error (MSE) (Isaaks and Srivastava1989; Eldeiry and Garcia2009).
MAE ¼1nXn i¼1 Z Uð Þi j Z U i ð Þj (1)
WhereZ(Ui) is the observed value at location i, Z*(Ui) is the estimated value at location i, and n is the sample number.
MSE ¼n1Xn
i¼1
Z Uð Þ Zi ð ÞUi
ð Þ2
(2)
Interpolation techniques can be divided into three categories: statistics (e.g. regression tree and multiple regressions), geostatistics (e.g. ordinary kriging and universal kriging) and hybrid techniques (e.g. cokriging and regression kriging) (Robinson and Metternicht2006; Zhu and Lin
2010). This study focused on comparison of ordinary kriging (OK), cokriging (CK) and regression kriging (RK) interpolation methods spatial prediction of FC and PWP for a semiarid region. These models:
Ordinary kriging
In geostatistics, OK is the most frequently used interpolation procedure (Matheron1969; Webster and Oliver2001). ZOK x 0 ð Þ ¼Xn i¼1 wi:Z xð Þi (3)
WhereZOKðX0Þis estimated value at the unsampled location at X0, n is the number of samples
in a search neighborhood,Z xð Þobserved values at points Xi i, wi denotes kriging weights which are chosen in order to minimize the estimation error variance is the sampled location.
Cokriging
Cokriging is a versatile statistical technique for estimating one or more variables of primary data when both primary (estimated data which needs time, cost and labor) and secondary (auxiliary) data are available. When using CK, auxiliary data does not have to be obtained at the same location. Some soil physical parameters such as sand, clay and bulk density were used for the estimation of FC and PWP in this study. This process is an estimation method that minimizes the variance of estimation error by taking into consideration the spatial correlation between primary and secondary data (Goovaerts 1997; Kravchenko and Bullock
1999; Meng et al. 2013). Z CKð Þ ¼X0 Xn i¼1 λi:ZiþX m j¼1 λj:Uj (4)
WhereZCKð Þ is estimated value at the unsampled location at XX0 0,λ and λjare the weightings
assigned to the primary and secondary data, respectively, n and m denote the numbers of primary data Zi (FC and PWP obtained from 167 soil samples) and secondary data Uj (sand % and BD obtained from 287 soil samples).
Regression kriging
Regression kriging, one of the hybrid interpolation techniques, is derived from the sum of regres-sions between interest (primary) variable and secondary variable(s) and kriging of residuals obtained from the regression (Wackernagel2003; Hengl et al.2004).
Z RKð Þ ¼X0 Xp k¼0 βk:qkð Þ þX0 Xn i¼1 Wi:e x 1 ð Þ (5)
WhereZRK ð Þ is the RK estimate at unsampled location at XX0 0,βkande xð Þ are the regression1
coefficients and regression residuals respectively, Wi denotes kriging weights using variogram of residuals,p and n denotes the number of predictors or auxiliary variables and number of samples, are the values of the auxiliary variables, which are the sand percentage and bulk density value of the studied soil in this case. RK was performed with Ordinary Least Square (OLS) and kriging of residuals was performed with OK in order to generate regression coefficients and residuals. Moreover, stepwise regression analysis was used in RK analysis while selecting auxiliary variables (automatically selecting sand values for FC and bulk density values for PWP measured across study area) for inclusion in the OLS equations.
Results and discussion
The relationship between thefield capacity and permanent wilting point and physical and chemical properties
The descriptive statistics of FC and PWP and auxiliary variables such as OM, BD, lime content, pH, EC, sand, clay and silt contents of soil samples are shown inTable 1. The physical and hydrophysical soil properties in the study area were likely affected by the high sand content which ranged from 22.7 % to 84.7 % in all great groups and soils series. In this case, similar results were reported by Pachepsky et al. (2001), they found the coefficient of determination of the relationship between the soil water content and clay content approximately 2.3 times less than that of relationships between the soil water content and sand or silt contents.
The Vertic Cambisol genetic soil group showed the highest variability in all the observed characteristics. Minimum values were determined for the soil series in Başkuyu, whereas the soil series in Odabaşı include maximum values. These soil series were classified as Haplic Calcisol and Petric Calcisol, respectively. Minimum skewness was determined for OM, calculated from 287 sample point values, whereas lime content had maximum skewness. The coefficient of variation (CV) for all soil characteristics except of pH, BD and OM (Table 1) showed values higher than 0.15. These differences are mainly due to soil samples being collected from the different soil series. This high CV is mainly due to heterogeneous distribution of soil properties and it is not possible to explain this variability with classical statistical methods. It is known that classical statistical methods do not take the spatial dependence of soil properties into con-sideration because they assume the observations to be independent in spite of their spatial dependence.
A number of physical and chemical soil properties of the study area, were used to predict the FC and PWP. The Kolmogorov-Smirnov statistical test was applied to test normality (Table 1). Log and square root transformation was applied to several soil properties – except of clay content, FC, EC and OM to complement the normal distribution.
Correlation analysis
Correlations among different soil characteristics were tested by using Pearson’s correlation coeffi-cient (R). FC and PWP values showed a high correlation (p < 0.01) with sand, clay and BD values across the study area, and a non-significant correlation with lime, OM and EC. Although silt and pH values were related to FC and PWP values, they were not selected as independent variables for regression analysis due to their low determination coefficient (R2) (Table 2). The high sand content of the study area might have affected the results of the correlation analyses. High sand content, and thus, higher BD values, have the closest relationship with the FC and PWP, thus they were selected as the main factors affecting values of the FC and PWP. In this study, 75% of the 287 disturbed soil samples include more than 67% of sand.
Regression analysis
The FC and PWP parameters were selected as dependent variables and BD, sand and clay were selected as independent variables to estimate FC and PWP according to the correlation analysis results. Probable multi-collinearity between auxiliary variables was tested before the analysis. The regression model and coefficients of determination estimated for FC and PWP are shown inTable 3. A stepwise regression technique, with a lower number of independent variables, was utilized to estimate values of the FC and PWP. It selected sand content for estimating FC in the regression analysis, with statistically significant dependence (R2= 0.44, p < 0.01). Likewise, using a lower number of variables for PWP provided the best prediction model when using the BD variable with R2 = 0.52 (p <0.01).
Table 1.Descriptive statistics of the soil properties based on genetic soil groups.
Variables N Min Max Mean SD CV Skewness K-S
Calcic Vertisol S (%) 20 23.26 62.75 38.99 7.21 0.18 1.42 1.018 Si (%) 20 18.71 29.24 25.05 2.41 0.09 −0.69 0.519 C (%) 20 18.54 48.68 34.94 5.81 0.16 −0.83 0.842 BD (Mg m−3) 20 1.24 1.46 1.32 0.04 0.03 1.62 1.052 FC (%, 0.033 MPa) 13 25.79 34.78 28.90 2.71 0.09 0.67 0.523 PW (%, 15 MPa) 20 14.56 19.13 16.64 1.53 0.02 0.47 0.828 pH 20 7.02 7.84 7.63 0.16 0.19 −2.97 1.299 EC (dS m−1) 20 0.21 0.44 0.34 0.66 0.61 −0.23 0.846 Lime (%) 20 0.79 17.66 7.62 4.66 0.28 0.68 0.592 OM (%) 20 1.07 3.70 2.39 0.69 0.29 0.26 0.443 Vertic Cambisol S (%) 15 42.42 84.73 56.06 11.27 0.20 1.12 0.725 Si (%) 15 2.84 27.89 17.68 7.85 0.44 −0.75 0.582 C (%) 15 11.93 33.87 26.24 5.99 0.23 −0.89 0.877 BD (Mg m−3) 15 1.35 1.57 1.40 0.05 0.04 1.67 0.767 FC (%, 0.033 MPa) 7 10.38 39.36 22.47 9.48 0.42 0.46 0.692 PW (%, 15 MPa) 7 3.01 17.89 11.72 5.74 0.49 −0.71 0.710 pH 15 7.13 7.81 7.63 0.19 0.02 −1.74 1.269 EC (dS m−1) 15 0.11 0.38 0.26 0.71 0.27 −0.28 0.401 Lime (%) 15 1.18 21.05 6.63 6.32 0.95 1.17 0.859 OM (%) 15 1.47 3.50 2.41 0.61 0.25 0.35 0.470 Haplic Calcisol S (%) 62 22.71 63.06 41.95 8.66 0.20 0.40 0.710 Si (%) 62 16.07 46.71 25.45 4.77 0.18 1.19 0.848 C (%) 62 11.90 55.55 32.58 7.34 0.22 0.68 0.648 BD (Mg m−3) 62 1.22 1.44 1.34 0.04 0.03 0.14 0.811 FC (%, 0.033 MPa) 36 15.57 39.13 27.73 5.78 0.20 −0.04 0.621 PW (%, 15 MPa) 36 10.08 21.63 15.54 2.40 0.15 −0.41 0.876 pH 62 6.99 8.18 7.73 0.16 0.02 −1.27 0.978 EC (dS m−1) 62 0.17 0.48 0.30 0.82 0.26 0.65 0.808 Lime (%) 62 0.32 30.44 9.35 7.35 0.78 1.18 1.211 OM (%) 62 1.11 3.77 2.203 0.61 0.27 0.22 0.476 Aric Cambisol S (%) 72 25.06 55.48 37.87 5.80 0.15 0.35 0.854 Si (%) 72 6.07 32.24 23.48 4.13 0.17 −0.96 0.426 C (%) 72 26.27 58.54 38.63 5.77 0.14 0.78 0.519 BD (Mg m−3) 72 1.23 1.40 1.30 0.03 0.02 0.08 0.285 FC (0.033 MPa) 30 25.02 40.59 30.26 3.89 0.12 0.74 0.761 PW (15 MPa) 30 13.54 23.17 18.20 2.38 0.13 0.21 0.941 pH 72 6.82 7.98 7.60 0.18 0.02 −1.94 0.003 EC (dS m−1) 72 0.14 0.60 0.36 0.81 0.22 0.08 0.623 Lime (%) 72 0.63 22.08 7.41 5.31 0.71 0.68 0.279 OM (%) 72 0.25 3.58 2.20 0.57 0.26 −0.25 0.860 Petric Calcisol S (%) 118 26.59 67.84 43.24 6.69 0.15 0.30 0.490 Si (%) 118 4.64 45.97 25.63 4.62 0.18 −0.21 0.904 C (%) 118 20.00 55.08 31.12 5.29 0.17 0.74 0.701 BD (Mg m−3) 118 1.26 1.45 1.34 0.03 0.02 0.27 1.123 FC (0.033 MPa) 81 19.17 36.90 27.33 3.75 0.13 0.40 1.140 PW (15 MPa) 81 8.10 24.14 15.69 2.36 0.15 0.43 0.860 pH 118 7.44 8.12 7.72 0.11 0.01 0.54 0.971 EC (dS m−1) 118 0.15 0.53 0.29 0.73 0.24 0.88 1.093 Lime (%) 118 10.3 37.54 11.39 7.45 0.65 0.79 1.096 OM (%) 118 1.14 3.85 2.60 0.69 0.26 −0.38 0.842
S sand, Si silt, C Clay, BD Bulk density, FCfield capacity, PWP permanent wilting point, EC electrical conductivity, Lime calcium carbonate content, OM organic matter, SD standard deviation, CV coefficient of variation, K-S Kolmogorov- Smirnov Z (tested at the 0.05 probability level)
Prediction of thefield capacity and permanent wilting point using the three different interpolation techniques
In this study, prior to geostatistical analysis, FC and PWP were also checked for drift, trend and anisotropy, tested with semivariogram models that were used to determine the degree of spatial variability, and then formed in four different directions to check geometric anisotropy. All the model procedures were used with the isotropic semivariograms for the OK, CK and RK techniques, and calculated omni-directionally for all variables.
Both the semivariogram models and coefficients for the OK, CK and RK techniques are provided inTable 4andFigure 3. In the semivariogram and cross semivariogram models for OK, both FC and PWP were fitted with the spherical model; in the spherical models for CK and semivariogram models for RK, both FC and PWP were fitted by using exponential models. Due to high sand content and bulk density which were negatively correlated with both FC and PWP, the cross-semivariogram models produced for cokriging assumed negative values. Jabro et al. (2010) also documented that cross-semivariogram modeling showed a strong, negative spatial interdepen-dence between soil penetration resistance and gravimetric water content in a grassland environ-ment. Structural variance values for semivariogram models ranged from 0.523 to 0.999. In addition, the estimated spatial dependence for values of FC and PWP ranged from 947 m to 5214 m in the semivariogram model. The R2values inTable 4show which modelsfit the experimental
semivar-Table 4.Semivariogram models of FC and PWP– ordinary kriging, regression kriging and cokriging. Variables Semivariogram Model C0 C0+C1 A (m) C1/C0+C1 RSS R2 FC(OK) Spherical 4.74 22.27 1773 0.78 1.65 0.98 FC(RK) Exponential 0.00 0.93 2571 0.99 5.27E-03 0.95 PWP(OK) Spherical 0.44 7.46 947 0.94 0.1 0.96 PWP(RK) Exponential 0.38 0.81 5214 0.52 3.48E-03 0.97 S Exponential 0.01 0.32 1527 0.96 9.24E-07 0.98
BD Spherical 8.0E-04 1.7E-03 1785 0.52 3.65E-09 0.98
CK
Cross-semivariogram models
FCxS Spherical −0.28 −0.63 1703 0.55 9.16E-04 0.97
PWPxBD Spherical −0.01 −0.08 1393 0.92 1.315E-06 0.99
FC: field capacity; PWP: permanent wilting point; C0: nugget variance; C1: structural variance; C0+C1: sill; A: range of influence
(m), RSS: residual sums of squares for the theoretical semivariogram and crossvariogram models.
Table 2.Correlations among soil characteristics.
S Si C BD FC PWP pH Lime OM EC S 1 Si −0.54** 1 C −0.80** −0.06 1 BD 0.91** −0.19** −0.95** 1 FC −0.66** 0.29** 0.57** −0.63** 1 PWP −0.70** 0.21** 0.68** −0.71** 0.68** 1 pH 0.19** 0. 23** −0.40** 0. 31** −0.23** −0.17* 1 Lime 0. 02 0. 31** −0.24** 0. 14* 0.16* −0.01 0.25** 1 OM 0. 06 0. 12* −0. 15** 0. 11* 0. 091 −0. 03 0. 03 0. 17** 1 EC −0. 14* −0. 06 0. 21** −0. 19** 0.03 0. 11 −0. 24** −0. 05 −0. 02 1 ** Correlation is significant at the 0.01 level (2-tailed).
* Correlation is significant at the 0.05 level (2-tailed).
Table 3.The regression model and coefficient for FC and PWP.
R2 R Regression model
FC (%, 0.033 MPa) 0.44 0.66 43.55–0.37 Sand
iogram and cross-semivariogram data very well in this study. Similarly, the RSS (residual sums of squares) values were extremely small.
Spatial distribution offield capacity and permanent wilting point
Ordinary kriging is one of the common methods used to determine spatial variability. It makes estimations by using kriging weightings computed by the semivariogram model – without any auxiliary variable – to minimize variance. Thus, the accuracy of the spatial distribution maps is lower than with other estimation methods (cokriging and regression kriging) having multivariate error rates and higher accuracy values. The spatial distribution maps produced with OK for the FC are similar to the spatial distribution maps obtained with CK and RK. However, estimated MSE value (2.84%, for FC) was higher than for the other methods, as expected. The spatial distribution maps produced for the PWP were partially different, but the estimated MSE value (2.03%) for PWP performed with OK was still higher (Table 5andFigure 4).
In the CK method, which includes correlation analysis results to improve estimation quality, sand was used as an auxiliary variable for FC, and BD was used as an auxiliary variable for PWP; spatial distribution maps were produced by using the cross semivariogram models. MAE (0.23%) and MSE (0.11%) values calculated for PWP with the CK method were lower than those calculated with the OK and RK methods. Similarly, Ahmed and De Marsily (1987) indicated that it is important for selecting between interpolation methods by taking into consideration of auxiliary variables to
Table 5.Statistics of the interpolation models for FC and PWP.
OK CK RK
Variables Sample Mean OK mean CK mean RK mean MAE MSE MAE MSE MAE
FC (%) 27.8 27.8 27.8 27.9 2.8 13.9 1.1 2.1 0.5
PWP (%) 16.0 16.0 16.0 16.0 2.0 7.9 0.2 0.1 0.5
Figure 3.Semivariogram and cross variogram models of FC and PWP in the study area (a, b and c are from OK, RK and CK for FC, respectively and d, e and f are from OK, RK and CK for PWP, respectively).
know the strength of the relation between the target variable and the auxiliary variable and the absence or presence of a spatial structure of the residuals of this relation. They also expressed that if the residuals are spatially correlated and the correlation coefficients between the variables are high then CK performs better than kriging combined with regression. Eldeiry and Garcia (2009) compared the capability of RK and CK to estimate soil salinity with Landsat Images. They reported that CK is more precise when the primary variable of interest is less intensively sampled than the other variables. In the current study, the negative correlation between PWP and BD (−0.719) was higher than that between FC and S (−0.661), which might have affected the CK estimation, leading to lower error values. A high correlation coefficient affected the CK estimation values more than the RK estimation values. The auxiliary variable reflected the spatial structure well and was highly correlated with target variable. In addition, the CK technique produced smaller error values than the OK and RK methods. The low coefficients of determination in the regression equation lead to an increase in the error terms; the reflection of these terms in the estimation maps also reduced the quality of the RK-produced maps and caused higher values of MAE and MSE.
For the RK process, log transformations were performed on the residual values obtained from the regression model to achieve normal distribution prior to geostatistical analysis, and RK maps were generated to combine the regression results and the residual values of the regression. The minimum error in the estimation maps generated for FC was revealed with the RK technique (Figure 4).
When the results of prediction performance of RK and OK were compared (Table 5), RK produced smaller MAE and MSE values than OK. Based on MSE values, RK yielded more accurate results than OK. Moreover, Kalivas et al. (2002) stated that CK had smaller MAE and MSE values than residual regression for clay and sand content values at 0–25 cm in Western Greece. Both smoothing the kriging values, due to the nature of the technique, and reflecting the residuals (uncertainty) in the estimations performed with RK improved the estimation quality. Odeh et al. (1995) reported that the regression residuals represent uncertainty which needs to be incorporated into kriging
Figure 4.Interpolation maps of FC and PWP in the study area (a, b and c were generated by OK, RK and CK for FC, respectively, and d, e and f were created by OK, RK and CK for PWP, respectively).
systems. They further stated that in the RK model CK, the regressed values and kriged residuals were added together at unsampled points and used as the predicted values.
Overall results showed that RK produced smaller MAE and MSE values than OK and CK for FC, but CK was superior for the interpolation of PWP in terms of smaller MAE and MSE obtained at the sampling point.
Conclusions
Our research was performed to spatially predict and prepare distribution maps of soil water retention properties of a semi-arid agricultural area of Central Anatolia in Turkey using the regression kriging, cokriging and kriging models for optimal estimation and spatial interpolation of values at unsampled points. The coefficient of variation of the soil water retention characteristics and other auxiliary variables exhibited high variability across the study area. The resulting maps were examined carefully and compared with the real field assessment. Based on the expert knowledge and supported by the low error values, the following models were chosen as the most suitable. It can be concluded, based on the results and considering all factors that all tested models can be used for spatial interpolation of the field capacity, but cogriging model is recom-mended for permanent wilting point spatial interpolation, as it provides the most realistic outcomes.
Soil properties relationships are important parameters that affect the quality of CK and RK estimates. Accuracy of CK and RK estimates vary depending on correlation enhancement, sample size and locational structures. Consequently, thefinal maps and calculated results can also be used in the decision processes for land and water managements and soil conservation practices by authorities, as well as by farmers at irrigatedfields of terrestrial ecosystem of the semi-arid areas.
Acknowledgements
The authors thank TUBITAK in Ankara, Turkey forfinancial support for this project (TOVAG Project No: 110 O 729) and Gregory T. Sullivan of the School of Earth and Environmental Sciences at the University of Queensland in Brisbane, Australia and Marketa Mihalikova from Czech University of Life Sciences Prague, Czech Republic for editing the English and rereading this manuscript.
Disclosure statement
No potential conflict of interest was reported by the authors.
Funding
This work was supported by the TUBITAK in Ankara; [110 O 729].
References
Adhikary SK, Muttil N, Yilmaz AG.2017. Cokriging for enhanced spatial interpolation of rainfall in two Australian catchments. Hydrol Processes. 31:2143–2161.
Ahmed S, De Marsily G.1987. Comparison of geostatistical methods for estimating transmission data on transsivity and specific capacity. Water Resour Res. 23:1717–1737.
Anderson S2001. An evaluation of spatial interpolation methods on air temperature in Phoenix, AZ. accessed 2017 Apr 18.http://www.cobblestoneconcepts.com.
Bayraminİ, Kılıç Ş, Dengiz O, Başkan O, Tunçay T, Yıldırım A, Koç A, Öğütmen Ç.2013. Radar görüntülerinin toprak etüt ve haritalama çalışmalarında kullanımı. Ankara: TUBITAK Projesi Sonuç Raporu; p. 110O729.
Bishop TFA, McBratney AB.2001. A comparison of prediction methods for the creation offield-extent soil property maps. Geoderma. 103:149–160.
Blake GR, Hartge KH.1986. Methods of soil analysis, Part I, Physical and mineralogical Methods. In: Bulk density and particle density. 2nd ed. Madison: ASA and SSSA Agronomy Monograph No 9; p. 363–381.
Bouyoucos GJ.1951. A recalibration of the hydrometer method for making mechanical analysis of soils. Agron J. 43:435–438.
Driessen PM.1970. Soil salinity and alkalinity in the great Konya Basin, Turkey. Agric Res Rep. 743:99. Centre for Agricultural Publishing and Documentation, Wageningen.
Eldeiry A, Garcia LA.2009. Comparison of regression kriging and cokriging techniques to estimate soil salinity using landsat images. J Irrig Drain Eng. 136:27–38.
Goovaerts P.1997. Geostatistics in soil science: state of- the- art and perspectives. Geoderma. 89:1–45.
Goovaerts P.2000. Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall. J Hydrol. 228:113–129.
Hengl T, Heuvelink G, Stein A.2004. A generic framework for spatial prediction of soil variables based on regression-kriging. Geoderma. 122:75–93.
Hengl TA.2007. A pratical guide to geospatial mapping of environmental variables. Luxembourg: JRC, Scentific and Technical Research series, office for official Publications of European Communities; p. 143.
Iqbal J, Thomasson JA, Jenkins JN, Owens PR, Whisler FD.2005. Spatial variability analysis of soil physical properties of alluvial soils. Soil Sci Soc Am J. 69:1338–1350.
Isaaks HE, Srivastava RM.1989. An introduction to applied geostatistics. New York: Oxford University Press. Jabro JD, Stevens WB, Evans RG, Iversen WM.2010. Spatial variability and correlation of selected soil properties in the
Ap horizon of a CRP grassland. Appl Eng Agric. 26:419–428.
Kalivas DP, Triantakonstantis DP, Kollias VJ. 2002. Spatial prediction of two soil properties using topographic information. Global Nest: Int J. 4:41–49.
Khosla R, Fleming K, Delgado JA, Shaver T, Westfall DG.2002. Use of site-specific management zones to improve nitrogen management for precision agriculture. J Soil Water Conserv. 57:513–518.
Klute A.1986. Water retention: laboratory methods. In: Klute A, ed. methods of soil analysis. Part 1. Physical and mineralogical methods. 2nd ed. Madison (WisconsIn:) SSSAJ; p. 635–662.
Kravchenko A, Bullock DG.1999. A comparative study of interpolation methods for mapping soil properties. Agron J. 91:393–400.
Lark RM.2000. Estimating variograms of soil properties by the method of moments and maximum likelihood. Eur J Soil Sci. 51:717–728.
Liu D, Wang Z, Zhang B, Song K, Li X, Li J, Li F, Duan H.2006. Spatial distribution of soil organic carbon analysis of related factors in cropland of the black soil region, Northeast China. Agric Ecosystems Environ. 113:73–81. Lloyd CD.2005. Assessing the effect of integrating elevation data into the estimation of monthly precipitation in Great
Britain. Hydrol J. 308:128–150.
Matheron G.1969. Le krigeage universal. Vol 1. Cahiers du centre de morphologie mathematique. Fontainebleau: Ecole des Mines de Paris.
McBratney A, Odeh I, Bishop T, Dunbar M, Shatar T.2000. An overview of pedometric techniques of use in soil survey. Geoderma. 97:293–327.
Meng Q, Borders M, Madden M.2010. High resolution image fusion using regression kriging. Int J Remote Sens. 31:1857–1876.
Mesić KI.2016. Comparison of ordinary and universal Kriging interpolation techniques on a depth variable (a case of linear spatial trend), case study of theŠandrovac Field. The Mining-Geology-Petroleum Engineering Bulletin; MGPB, p. 41–58.
Miháliková M, Başkan O, Dengiz O.2015. Capability of different interpolation models and pedotransfer functions to estimate soil hydraulic properties in Büyükçay Watershed. Environ Earth Sci. 74:2425–2437.
Miháliková M, Matula S, Doležal F.2014. Application of k-Nearest code to the improvement of class pedotransfer functions and countrywide Field Capacity and Wilting Point maps. Soil Water Res. 9:1–8.
Miháliková M, Özyazici MA, Dengiz O.2016. Mapping soil water retention on agricultural lands in central and Eastern parts of the black sea region in Turkey. J Irrig Drain Eng. 142:05016008.
Minasny B, Hartemink AE.2011. Predicting soil properties in the tropics. Earth Science Reviews. 106:52–62. Mzuku M, Khosla R, Reich R, Inman D, Smith F, MacDonald L.2005. Spatial variability of measured soil properties
across site-specific management zones. Soil Sci Soc Am J. 69:1572–1579.
Odeh I, McBratney A, Chittleborough D.1995. Further results on prediction of soil properties from terrain attributes: heterotopic cokriging and regression-kriging. Geoderma. 67:215–226.
Pachepsky YA, Rawls WJ.2004. Development of pedotransfer functions in soil hydrology. Amsterdam: Elsevier; p. 512. Pachepsky YA, Timlin DJ, Rawls WJ.2001. Soil water retention as related to topographic variables. Soil Sci Soc Am J.
65:1787–1795.
Qingmin M, Zhijun L, Borders BE.2013. Assessment of regression kriging for spatial interpolation- comparisons of seven GIS interpolation methods. Cartography and Cartogr Geogr Inf Sc. 40:28–39.
Richards LA1954. Diagnosis and improvement of saline and alkali soils. U.S. Department of Agriculture, Handbook No 60 p.109,
Robinson TP, Metternicht G. 2006. Testing the performance of spatial interpolation techniques for mapping soil properties. Comput Electron Agr. 50:97–108.
Soil Surveystaff. 2015. Illustrated Guide to Soil Taxonomy. Version 1.1. U.S. Department of Agriculture. National Resources Conservation Service, Lincoln, Nebraska.
Sullivan DG, Shaw JN, Rickman D.2005. IKONOS imagery to estimate surface soil property variability in two Alabama physiographies. Soil Sci Soc Am J. 69:1789–1798.
TİGEM. 1997. Altınova tarım işletmesi topraklarının detaylı toprak etüt ve haritalaması”, Tarım İşletmeleri Genel Müdürlüğü Yayınları, No: 23, 100 sayfa, Ankara (in Turkish).
Vergouwen L.1981. Salt minerals and waters from soils in Konya and Kenya. PhD Thesis. Agricultural University of Wageningen, The Netherlands; p. 140.
Wackernagel H.2003. Multivariate geostatistics. New York: Springer.
Webster R, Oliver M.2001. Geostatistics for environmental scientists. Chichester: Statistics in Practice. John Wiley & Sons; p. 265.
Wösten JHM, Hoogland T, Massop H, De Vries F, Vroon HRJ, Wesseling JG.2012. New soil physical schematization of Netherlands. Wageningen (The Netherlands): Alterra Report.
Wösten JHM, Pachepsky YA, Rawls WJ.2001. Pedotransfer functions: bridging the gap between available basic soil data and missing soil hydraulic characteristics. Hydrol J. 25:123–150.
WRB.2014. World reference base for soil resources. international soil classification system for naming soils and creating legends for soil maps. Food and agriculture organization of the United Nations. World Soil Resources Reports, no: 106, 203p.
Zhu Q, Lin HS.2010. Comparing ordinary kriging and regression kriging for soil properties in contrasting landscapes. Pedosphere. 20:594–606.
