Site-Specific Leaching Map of a Salt Affected Soil in Egypt

Soil salinity is a crucial problem facing countries in arid and semi arid regions which negatively affects sustainable agriculture Eilers et al. [1]. This problem needs to be accurately quantified and mapped in order to mitigate and control its negative effects Herrero and Pérez-Coveta [2]; Benyamini et al. [3]; Wang [4] especially in irrigated soils Amezketa [5]. Soil salinity cab be solved by leaching soluble salts from the root zone by adding water with good drainage conditions. In Egypt, salt-affected soils are located in the NorthernCentral part of the Nile Delta and on its Eastern and Western sides. Also, in Wadi El-Natroun, El-Kebeir, the Oases, many parts of the Nile Delta and Valley and El-Fayoum province are characterized as salt affected soils. In Egyptian irrigated lands, about nine hundred thousand hectares suffer from salinization. These are distributed as follows: 60 % is in Northern Delta, 20 % in Southern Delta and Middle Egypt and 25 % in Upper Egypt FAO [6]. Conventionally, data of soil properties are summarized by averages of collected soil samples with no consideration of the spatial variation either at macro scale or at micro scale (within-field) (Webster and Oliver [7] and Navarro-Pedreño et al. [8].


Introduction
Soil salinity is a crucial problem facing countries in arid and semi arid regions which negatively affects sustainable agriculture Eilers et al. [1]. This problem needs to be accurately quantified and mapped in order to mitigate and control its negative effects Herrero and Pérez-Coveta [2]; Benyamini et al. [3]; Wang [4] especially in irrigated soils Amezketa [5]. Soil salinity cab be solved by leaching soluble salts from the root zone by adding water with good drainage conditions. In Egypt, salt-affected soils are located in the Northern-Central part of the Nile Delta and on its Eastern and Western sides. Also, in Wadi El-Natroun, El-Kebeir, the Oases, many parts of the Nile Delta and Valley and El-Fayoum province are characterized as salt affected soils. In Egyptian irrigated lands, about nine hundred thousand hectares suffer from salinization. These are distributed as follows: 60 % is in Northern Delta, 20 % in Southern Delta and Middle Egypt and 25 % in Upper Egypt FAO [6]. Conventionally, data of soil properties are summarized by averages of collected soil samples with no consideration of the spatial variation either at macro scale or at micro scale (within-field) (Webster and Oliver [7] and Navarro-Pedreño et al. [8].
Geostatistics provides good tools in capturing within field spatial variation of soil variables that is used to delineate management zones and consequently know where and how much farm inputs should be added Oliver and Webster [7]. Geostatistical methods are non-destructive, time saving and cost effective as comparing with traditional methods and provide finescale information about soil variables. Spatial distribution of soil properties has been recognized by various researchers Burgess and Webster [9]; Warrick et al. [10]; Odeh et al. [11]; Juang & Lee [12]. Oridinary kriging as a geostatistical interpolation method was the most common used to predict and map soil parameters at unsampled locations Lopez-Granados et al. [13]; Meul & Van Meirvenne [14]; Sumfleth & Duttmann [15]. This study aimed at assessing the spatial distribution and predicting soil salinity over a field in Sharkia governorate using ordinary kriging for the purpose of preparing a prescription map of leaching requirements.

Site Description, Sampling and Laboratory Analysis
This study was conducted in 4.6 ha field (30° 51' 53.31'' N, 32° 02' 50.45'' E), located in Bahrelbaqar, Sharkia governorate. A total of 100 soil samples based on an almost regular 25x25 m grid at 0.20 m depth were collected and transported to the soil laboratory, where they were dried, ground and sieved with a 2 mm sieve and electrical conductivity (EC) was measured. EC was measured in 1:5 soil water extract according to Jackson [16]; after that the values were converted to soil paste extract by using multiplicative factors for converting EC 1:5 to EC of soil paste extract (ECe) according to Slavich and Petterson [17]. Soil texture fractions (sand, silt and clay) were determined according to the international pipette method. The soil texture down to 0.20 m includes two textures of light clay and heavy clay according to the Moeys [18].

Data Processing
The whole data set was divided randomly into two subsets: the calibration data set containing 2/3 of soil samples and the validation data set containing 1/3 of soil samples for the property under study ( Table 1). For the validation data set being representative of the calibration data set, the means of the variable (ECe) in calibration and validation subsets should not be statistically different. So a t test with unequal variances was carried out to test if there is a significant difference between means.

Interpolation Procedure
Ordinary kriging based on 1m x1 m grid was applied on the calibration data set using ArcGIS 10.1 (ESRI, USA) mapping software. The interpolated map was classified into 3 geometric intervals classes (low, medium, high) to be easy for managing. Validation of prediction was implemented using cross-validation and validation consisted of 63 and 28 soil samples respectively. The mean error as a measure of bias and root mean square error as a measure of inaccuracy of prediction were calculated for both.

An overview of kriging
Kriging is a type of interpolation technique. The procedure is similar to averaging techniques of interpolation but the weights are derived from the spatial arrangement as well as from the distance between nearby points, i.e. from the variogram. The fitted variogram or the directional variogram (for anisotropic variation) is/are used to determine the weights λi needed for local interpolation. The weights are chosen so that the estimate is unbiased, and that the estimation variance is less than for any other linear combination of the observed values. Mapping the spatial distribution of a soil property involves interpolation or spatial prediction. Geostatistical interpolation uses the variogram to optimize prediction by kriging. The most basic form of kriging is ordinary punctual kriging in which the unknown value z(x0) of a given realization of Z(x0) in an unsampled point x0 is predicted from the known values z(xi) i= 1, 2, . . ., N, at the support points xi of the same realization using a so-called "best linear unbiased estimator" (BLUE). The best linear unbiased predictor z*(x0) of Z(x0) is given by a linear combination of the bservations: Where i are weights. The weights are chosen in such a way that the estimator is unbiased: And the estimation is minimized.
Using a Lagrangian multiplier µ, minimization of the estimation variance under the constraint of unbiasedness yields a set of N+1 linear equations: From which the and µ can be calculated. The estimation variance is then given by: Or equivalently: And represents the uncertainty in the prediction in x0. ( ) , g xi xj and ( ) , 0 g xi x are the semi-variances between the observed locations xi and xj and between the observed location xi and the interpolated location x0, respectively. In the case of spatial dependence, semivariance tends to increase with the distance between observations; therefore errors decrease with the density of data and not just with their total number, as is the case with traditional statistical models. It needs to be pointed out that kriging is optimal and unbiased only on the condition that the model is correct. However, predictions are only slightly affected by the choice of the model, provided it is reasonable of course. This is one of the strengths of kriging that is robust enough in this sense; however, error variances can be seriously affected by the model Castrignanò [19]. neighbourhood are all zero and the estimation variance equals the observation. e) Only the nearest few samples are spatially correlated to the kriged location and therefore they are the most weighted. Two important advantages become clear: firstly, the variogram needs to be accurate only in the first few lags; secondly whatever is gained from introducing distant points is limited also because sample locations interposed between the kriged point and more distant samples screen the distant ones by reducing their weights Castrignanò [19].

Validation Test
In geostatistical practice, the usual method of testing is crossvalidation. However, its results are actually biased and somewhat too optimistic Creutin and Obled [20], because it retains the same variogram, whereas the variogram should be recomputed and fitted every time that an observation is removed Laslett et al. [21]. Moreover, cross-validation is not true validation, because the same sample data set is used for both estimation and validation. All these shortcomings can be avoided by using a separate independent set of data for validation. So a validation data set consisting of 28 soil samples was used.

Prediction Assessment
The performance of the prediction was assessed using cross validation Isaaks & Srivastava [22], whereby one observation (z) at a time is temporally removed from the data set and re-estimated (z*) from the remaining data. To assess the precision and accuracy of estimation of the soil variable under study, two statistical criteria were considered: mean error (ME), as an indicator of bias, and mean standardized squared error (MSSE) (scaled by the predicted standard deviation of estimation), as an accuracy measure: Where N is the number of active observations and σ the cokriging standard deviation, If the estimation is unbiased and accurate, the first statistic should be close to zero, where as the second one should be close to one, because the latter corresponds to the ratio between an experimental variance and a theoretical one Carroll & Cressie [23].
Calculations of leaching requirements for different delineated zones, Reclamation requirement was calculated using Reeve equation [24]. Where D iw is the depth of leaching water (cm), Ds is the depth of soil (cm), EC ei and EC ef are soil salinity in dS m -1 before and after leaching respectively. Leaching requirement was calculated to reduce soil salinity to be 2 dS m -1 for only zones having a mean value of ECe higher than 4 dS.m -1 for 0.20m depth. The mean value and area of each delineated zones were calculated as follows: a) Reclassification of the interpolated raster map using "reclassify" tool, available in spatial analyst tools. b) Applying of "raster to polygon" tool, available in conversion tools to convert a raster dataset to polygon features c) Use of "dissolve" tool , available in data management tools to unify polygons of the same class based on grid code number d) Use of "zonal statistics as table" tool, available in Spatial analyst tools to obtain the mean value of each class e) Use of "calculate areas", available in spatial statistics tools to calculate the area of each class.
The leaching requirements was also calculated based on the average value of all samples (91) to quantify the water saved and make a good comparison between the two geostatistical and traditional methods. Table 2 shows that the calibration and validation data sets were extracted from the same population since the p-value (0.06) is greater than α value (0.05) in the achieved t test. This is an important criterion to use different data set in validation (Table 2).   Table 3 and Figure 1 show variogram model parameters where a spherical model with nugget effect equal to 1.18 was fitted to the experimental variogram. The range of autocorrelation was 57.39 m and the partial sill was 1.60 (Table 3) & (Figure 1). Table 4 shows the assessment statistics in cross-validation and validation data sets. It can be seen that results of cross-validation are better than those of validation. However results of validation are more reliable. The mean error and MSSE in cross-validation and in validation were 0.09 and -0.95 and 1.07 and 1.27 respectively. Although the MSSE values is somewhat far from 1 but still within the tolerance interval (±3√2/N, N is number of observations) Chiles & Delfiner [25], which means that the model is accurate (Table 4). Table 5 shows the different quantities of water to be added to the soil to reduce soil salinity to 2 dS.m -1 in the different classes ( Figure  2). There is no need to leach the first class, named as low, whereas in the medium class, 3625m 3 of water should be added to reduce the salinity in that area to reach 2 dS.m -1 . To reduce the salinity of the high class, 1149m 3 of water should be added to reduce its salinity to 2 dS.m -1 . It is worth mentioning that assuming the ECe mean value for 91 soil samples (5.12 dS.m -1 ) -without geostatistical interpolation -a quantity of 6203 m 3 of water should be added to reduce salinity of the whole field to 2 dS.m -1 . So 1429 m 3 can be saved and then used for irrigation [26,27]. These results emphasize the importance of geostatistical techniques in detecting within field variability and hence applying site-specific management (Table 5) & ( Figure 2).     The low ECe values found in the south western area of the field whereas high and medium ECe values were observed in the remaining areas of the field, In order to interpret the reason beyond the existence of different areas with different ECe values, a quick check of soil texture (Figure 3) was performed to interpret the low, medium and high ECe values in the delineated zones. The results of soil texture proved that the low values coincided with the coarser soil texture (light clay) and then increased infiltration rate. The high and medium values detected in the remaining parts of the field coincided with heavier soil texture (heavy clay) (Figure 3).

Conclusion
Soil properties vary greatly over space and/or time. So it is necessary to recognize the spatial and temporal variation in soil properties. This cannot be achieved by the traditional methods based on averages of laboratory measured data and hence the essence of the usage of the developed tools such as geostatistics. Geostatistics allows predicting soil properties in unsampled locations and consequently to delineate site-specific management zones which are used for variable rate applications. In this study, a prescription map of leaching requirements was developed using ordinary kriging. The prediction map was classified into 3 management zones with different mean values of salinity. One of the three management zones had a mean value of ECe lower than 4 whereas the remaining two zones had ECe values higher than 4. So the leaching requirements were calculated only for zones having ECe higher than 4. A comparison between traditional method of calculating leaching requirements based on average of collected samples and geostatistical method was achieved and a quantity of 1429 m3 water was saved by applying the geostatistical method. These results encourage the use of geostatistics to develop sitespecific management zones for leaching requirements especially in arid and semi arid regions and under water scarcity conditions.