Abstract

Knowledge of spatial variability is important for management of land affected by various anthropogenic activities. This study was conducted at West Mesa land application site to determine the spatial variability of electrical conductivity (EC1:1) and suggest suitable management strategy. Study area was divided into five classes with EC increasing from class I to V. According to the coefficient of variation (CV), during 2009 and 2010, EC1:1 values for different classes were low to moderately variable at each depth. Semivariogram analysis showed that EC1:1 displayed both short and long range variability. Area coverage of classes I and II were much higher than classes III, IV, and V during 2009. However, during 2010 area coverage decreased from 26% to 14.91% for class II, increased from 12.11% to 22.97%, and 10.95% to 20.55 for classes IV and V, respectively. Overall area under EC1:1≥ 4 dS/m increased during 2009. Soil EC map showed EC classes IV (4.1–5 dS/m) and V (>5.1 dS/m) were concentrated at northwest and southeast and classes I and II were at the middle of the study plot. Thus, higher wastewater should be applied in the center and lower in the northwest and southwest part of the field.

1. Introduction

Land application of treated wastewater is a cost-effective alternative to reduce the pressure on freshwater resources especially in the arid and semiarid regions of the world [1]. Application of treated wastewater often raises a question—“Is it environmentally sound to apply?” One of the possible risks for soil and plant due to the land application is salt accumulation in soil, which is further accelerated by the lower precipitation and higher evapotranspiration in the arid regions. Electrical conductivity (EC) is an important chemical parameter to describe the soluble salt and soil salinity and can be used for monitoring soil quality changes due to wastewater application. Information on soil EC is relatively inexpensive to determine and can be gathered in great intensity from the field.

Soil EC maps were used to study the soil condition and apply the field management strategies in the Malaysian paddy fields, where higher EC levels were determined in the southern and central regions of the study area [2]. ECa surface map produced using 400 square-foot grid cells showed that the highest salinity levels were located in the eastern portion of the study plot dominated by the denser clay soil [2]. Coefficients of variations (CVs) and semivariogram analysis were used to characterize the variability and spatial structure of EC in a 33 ha drainage experimental plot in Cauto River Valley, Cuba [3]. The variogram analysis revealed that EC and soil physical properties displayed similar spatial scales for the arable field of Southern England [4]. Knowledge of spatial variability of soil properties is useful for site-specific management and applying inputs in accordance with the specific requirements of soil and plant [57].

Most of the EC researches were conducted in the agricultural fields located in the humid areas with high rates of precipitation. To the best of our knowledge very little information is available for land application sites located in semiarid regions. Semiarid and arid environments also have heterogeneous soils and vegetation characteristics [8]. Along with soil physical and hydraulic properties, nonuniform applications of treated wastewater can further exacerbate the spatial variability of soil EC. The primary objective of this study was to determine the spatial variability of EC, due to the nonuniform application of treated wastewater in the West Mesa land application site. In addition this research is aimed to identify the hot and cold EC spots in the land applications site and suggest appropriate management strategy for wastewater application and soil environment management.

2. Materials and Methods

2.1. Experimental Site, Soil Sampling, and Soil Analysis

Study was conducted at West Mesa land application site (longitude W 106° 54.408 latitude N 32° 15.99, altitude 1298 m) located west of Las Cruces, New Mexico, USA. Treated municipal and industrial wastewater has been applied to the 36 ha mesquite (Prosopis glandulosa Torr. var glandulosa) and creosote (Larrea tridentata, (DC) Co7v.) dominated shrub land since 2002 by fixed head sprinkler irrigation system. The dominant soil series in the study area are Onite-Pajarito association and Bluepoint loamy sand. Land application site is vertically dissected by an arroyo approximately into 17.5 ha and 18.5 ha plots [9]. The 17.5 ha plot was selected for the study because the sprinklers on this plot were properly working and landscape is level compared to the 18.5 ha plot. To make the area perfectly rectangular for the ease of constructing transit lines, some area beyond the land application was also included. Six horizontal and nine vertical transect lines were drawn, which resulted in 54 grids for the entire 17.5 ha plot (Figure 1). Soil sampling was performed at the center of each gird ( m). Five grids were chosen randomly for additional lags, and soil samples were taken at the distances of 2 m, 7 m, and 22 m from the center of the ( m) grid sampling point in two directions (Figure 1).

About 168 bulk and core soil samples were collected at 0–20 cm depths during March 2009. Again same numbers of bulk samples were collected during March 2010. Four random bulk soil samples were also selected, and EC was determined in 1 : 1, 1 : 2, 1 : 3, and 1 : 4 soil : water. A linear relationship was obtained between EC and soil-water ratio with an of 0.94 () (Figure 2). Remaining soil samples were arid dried, passed through 2 mm sieve, and analyzed for EC in 1 : 2 ratio of soil : water using EC meter. The EC values were converted to 1 : 1 soil-water ratio using the linear relationship obtained in Figure 2. In addition, EC on soil saturation paste extract (ECe) for 12 randomly selected samples was also determined. Soil particle size distribution (sand, silt, and clay contents) was analyzed using hydrometer method [10], bulk density (BD) using core method [11], and OM using loss on ignition method at Harris Lab, Columbus, Nebraska.

2.2. Data Analysis

The EC1:1 values were tested for normality using Statistical Analysis System (SAS) [12]. A histogram, box plot, and normal plot were constructed for the EC1:1 at 0–20 and 20–40 cm depths during 2009 and 2010. There were two outliers for EC1:1> 10 dS/m at 0–20 cm depth during 2009 which quantitatively affected the normality of distribution. Outliers were also observed in the EC data set at 20–40 cm depth during 2009 and 0–20 and 20–40 cm depths during 2010. However, these outliers were important to the analysis of salinity, and hence, they were kept in the dataset. The data were also log and square root transformed; however, no improvements were observed. Therefore, original dataset was used for further analysis.

EC1:1data were classified with smart quintile, equal interval, and natural breaks technique introduced by ArcGIS software. All the three techniques did not visualize much variability; most of the data fell into the same class. So a boundary was set where there was a relatively big jump in the data values, and attempts were made to keep number of counts similar in each class. Using these two criteria, EC values (dS/m) were divided into five classes; class I contained EC1:1≤ 1, class II 1 < EC1:1≤ 2, class III 2 < EC1:1≤ 4, class IV 4 < EC1:1≤ 5, and class V EC1:1> 5.

Descriptive statistics including mean, number of counts, standard error (SE), maximum, minimum, skewness, and kurtosis were determined using SAS for each class, separately. The coefficient of variation (CV) was also calculated for each classes, as a general index of variability among the samples taken in each site and was classified as least (CV < 0.15), moderate (0.15 < CV < 0.35), or most (CV > 0.35) variable [13]. GS+ software [14] was used to obtain semivariograms for the EC1:1separately by depth and year of sampling. The spherical, linear, exponential, and Gaussian models shown below were fitted to the semivariogram.

Spherical model: Exponential model: Linear: Gaussian: where is the lag distance, is the nugget effect, which is the local variation occurring at scales finer than the sampling interval or fine scale variability, measurement, or sampling error; is the sill or total variance, and a is the range of spatial dependence at which the semivariogram levels off to reach the sill value and beyond which sampling variables are not correlated [1517].

Semivariograms were calculated both isotropically and anisotropically. The anisotropic calculations were performed in four directions (0°, 45°, 90°, and 135°) with a tolerance of 22.5° to determine whether semivariogram functions depended on sampling orientation and direction. The parameters of the model including nugget, range, and total sill, and residual sum of square (RSS) and regression (R) were determined for each semivariogram. Nugget semivariance is the variance at zero distance; sill is the lag distance between measurements at which one value for a variable does not influence neighboring values; range is the distance at which values of one variable become spatially independent from others; RSS provides an exact measure of how well the model fits the variogram data, and provides an indication of how well the model fits the variogram data. Model selection for semivariograms was done on the basis of RSS, , and visual fitting. Nugget to total sill ratio (NSR) was expressed as the percentage of total semivariance and was used to define for spatial dependency of EC1:1. NSR < 0.25 indicated strong spatial dependence, 0.25 < NSR < 0.75 indicated moderate spatial dependence, and NSR > 0.75 indicated weak spatial dependence [18]. Several other authors have also mentioned NSR-based classification when discussing spatial dependence [1925]; however none of the papers mentioned NSR as well as distribution map patterns together to discuss spatial dependence. EC1:1spatial distributions maps were produced by using Surfer version 8 [26] to determine the spatial coverage distribution of EC1:1 or hot and cold spots in the land application site.

Spatial autocorrelation was used to express spatial changes in the soil EC1:1and degree of dependence among neighboring sampling locations. Information on autocorrelation helps to identify the soil sampling interval for which observations remains spatially correlated and can be used for designing soil-sampling schemes [27]. Moran’s statistics [28] was used to calculate the coefficient at selected lag distance. The Moran’s is the conventional method used to measure autocorrelation, similar interpretation to the Pearson’s correlation statistics for independent samples. In both Moran’s and Pearson’s correlation, statistics range from +1.0 to −1 with +1 indicating strong positive spatial autocorrelation, 0 a random pattern, and −1.0 strong negative autocorrelation. The statistic is defined as where is the number of points, the variables of interest, the mean of and the spatial weight describing the adjacency or distance between the th and th point.

3. Results and Discussion

3.1. Soil Texture and Wastewater Quality

Soil textural analysis revealed that average sand, silt, and clay contents of the land application site were , , and , respectively, at 0−20 cm depth. Similarly average OM content of the land application site was 0.56 ± 0.02% (Table 1). Sodium adsorption ratio (SAR), pH, chloride (Cl), nitrate (NO3  −), and EC are important criteria for irrigation quality particularly at the west mesa land application facility. The amount of wastewater applied and chemical properties of the wastewater were provided by the City of Las Cruces. Total of 49 cm3/cm2 of wastewater was applied during 2008 and 16 cm3/cm2 till November, 2009 (Table 2). Wastewater analysis showed that two-year average concentration of Cl was  mg/kg, NO3  − mg/kg, EC  mg/kg,, SAR was , and pH was (Table 2).

3.2. Descriptive Statistics

Each EC1:1 class covered variable percentages of areas and numbers of counts at 0–20 cm depth. About 6.89 ha (39.37%) of area was occupied by class I, 4.69 ha (26%) by class II, 1.9 ha (10. 85%) by class III, 2.2 ha (12.11%) by class IV, and 1.9 ha (10. 85%) by class V during 2009 at 0–20 cm depth (Table 3). The EC for classes II and V were moderately variable with CV between 0.15 and 0.35, and EC for classes I, III, and, IV were least variable with CV < 0.15. EC distribution map at 0–20 cm showed that classes III, IV, and V were concentrated in northwest direction and southeast direction of the plot, while the classes I and II were concentrated in the middle of the plot (Figure 3(a)). Mean EC1:1 for class IV was 4.70 dS/m for class V was 8.12 dS/m which was >4 dS/m, and for classes I, II, and III was <4 dS/m. Average mean EC1:1 of class III, IV, and V was lower at 20-40 cm depth than at 0–20 cm depth during 2009 (Table 4). CV was least for class IV, and moderate for classes I, II, III, and V. This indicated that class IV was least variable and classes I, II, III, and V were moderately variable.

The EC1:1 at 20–40 cm depth during 2009 ranged from 0.66 to 6.51 dS/m with the standard error of 0.03 dS/m to 0.63 dS/m, respectively, (Table 4). The total area covered by EC1:1 classes I, IV, and V was similar to 0–20 cm depth with the similar number of counts (Table 4). For Class II (20.55%) area coverage was less at 20–40 cm than 0–20 cm depth (26%). In contrast, area coverage for class III was higher at 20–40 cm (17.07%) than 0–20 cm depth (10.85%). The occurrence of hot and cold spots (high and low) indicated that wastewater application was not uniform, although inherent variability of soil hydraulic properties could have caused variable lateral and horizontal movement of wastewater in the soil profile. Classes I, II, III, and V were moderately variable with CV > 0.15, and class IV was least variable with CV < 0.15. EC1:1 map showed similar trends at 20–40 cm depth, where higher EC1:1 values were concentrated in northwest and southeast direction of the plot, while the class I and II were concentrated in the middle of the plot (Figure 1(b)).

During 2010 at 0–20 cm depth, coverage area of class I was 39.37%, which was similar to 2009. Area coverage for class II was 14.91%, which was 11% less than during 2009. However, area under classes III, IV, and V increased by 5.27%, 1.37%, and 5.22%, respectively (Table 5). This indicated that areas under higher EC1:1 classes III, IV, and V increased in the experimental site and area under lower EC1:1 class II was decreased. EC1:1 classes I, III, and IV were least variable with CV < 0.15, and classes II and V were moderately variable with CV between 0.15 and 0.35. At 20–40 cm depth area coverage under class I decreased than in 0–20 cm depth from 39.37% to 28.25%, and classes IV and V increased from 13.48% to 22.97% and 16.12% to 20.55%, respectively (Tables 5 and 6). This indicated that leaching patterns were slightly different between 0–20 and 20–40 cm depth. Similarly, EC spatial distribution map showed less area under EC class II at 0–20 cm depth during 2009 than at the same depth during 2010. In general higher EC classes were encroaching on the lower EC class areas in the southeast and northwest direction of the study plot (Figure 4(b)).

3.3. Spatial Variability

Four semivariogram models, linear, spherical, exponential, and Gaussian, were fitted to the data for 0–20 and 20–40 cm depths for 2009 and 2010 (Tables 710). Anisotropic semivariograms did not show any differences in spatial dependence based on direction, indicating that the semivariance was only a function of the separation distance between pairs at each depth and time of sampling. Therefore, only isotropic semivariograms were chosen. The geostatistical analysis indicated different spatial distribution models and spatial dependence levels for the soil EC1:1 at 0–20 and 20–40 cm depth during 2009 and 2010 (Tables 710). Among four models fitted to the EC1:1 data at 0–20 cm depth, exponential model was the best fit (Figure 4(a)) with lowest RSS and highest R and a range of 306 m (Table 7). The NSR was 0.28 indicating that EC values were moderately spatially dependent. EC distribution map (Figure 2(a)) with many smaller numbers of patches in southwest and northeast directions and few bigger patches in the center agrees with the NSR of moderate spatial dependence. Range of 306 m indicated that soil samples ≥306 m apart were independent. In a separate study in Centralia, Montana, spherical model was obtained as best fit for ECa with NRS of 0.74 and a range of 162 m [29]. Gaussian model was the best fit at 20–40 cm depth during 2009 (Figure 4(b)). NSR was 0.31, indicating that EC1:1 values were moderately spatially dependent (Table 8). Similar to our study Gaussian model was obtained as the best fit for EC1:5 at 20–40 cm depth with a range of 128 m in a study conducted at Cauto River Valley, Cuba [3]. During 2010 at 0–20 cm depth, exponential model was the best fit (Figure 4(c)) with NSR of 0.17, indicating that EC1:1 values were strongly spatially dependent (Table 9). Exponential model was the best fit for the EC1:1 at 20–40 cm depth during 2010, with NSR 0.10 indicating strong spatial dependency (Table 10). EC distribution maps (Figures 3(a) and 3(b)) showed few bigger patches in the center and many smaller patches in the southeast and northwest direction as in 2009, however did not agree with the NSR of strong spatial dependency.

3.4. Spatial Autocorrelation

Soil EC1:1 correlograms showed positive as well as negative spatial autocorrelation for both 0–20 and 20–40 cm depth during 2009 and 2010 (Figure 5). Positive autocorrelation at shorter lag distances and negative autocorrelation at longer lag distance were observed at both depths and years. Moran’s statistics was 1 at zero lag distances and starts to decline as the lag distances increase. The autocorrelation of EC1:1 at both 0–20 and 20–40 cm depths during 2009 and 2010 became negative as the lag distances approached 125 m. Initial rapid decline of Moran’s at the shorter lag distances indicated that data were tending to a random distribution at the sampling interval of 25 m at 0–20 cm depth and 50 m at 20–40 cm depth. Correlograms of soil EC1:1 for both years and depths showed significant () autocorrelation at lag distances of 25 m and started to become insignificant at a separation distance of 125 m. Thus the sampling interval beyond 25 m is suitable for collecting independent soil samples at 0–20 m depth. In contrast, study conducted in a cotton field at Perthshire, Minnesota, showed that correlograms of soil bulk density showed significant Moran’s up to the lag distance of 100 m for surface soil [23]. Correlograms at 20–40 cm depth showed significant () autocorrelations at lag distances of 50 m. Therefore, independent soil samples could be collected for separation distances of ≥50 m at 20–40 cm depth. Few oscillations in the correlograms curve were also observed at larger lag distances and may be due to one or several periodic components [30] and some random variation in the data due to the nature of the soil properties or fewer pairs available at larger lag distances.

3.5. Management and Reclamation Strategy

EC classes IV and V covered nearly 23% of the areas with EC > 4.1 dS/m during 2009 at 0–20 cm depth. Similarly during 2010 area under classes IV and V was nearly 30%, showing an increase of 7% in a year. In contrast, EC class II decreased by 11% from 2009 to 2010 at 0–20 cm depth. Generally for an agricultural field, EC1:1 above 0.8–1.0 dS/m is considered as the threshold level, beyond which the growth and activity of several plants and microorganism can be altered [31]. Similarly, ECe> 4 dS/m is one of the indicators of saline soil that can affect the growth and development of plants. The predominant vegetation of the study site is creosote (8.7%), mesquite (14.4%) [32], perennial and annual weeds (0–80%). The tolerance limit of EC e for creosote is up to 7.51 dS/m [33], and EC e for mesquite is up to 9.36 dS/m [34]. Analysis of EC e of randomly selected samples () showed that EC e was approximately 2.13 times greater than EC1:1in the study site. Thus, the actual tolerance limit based on is 3.52 dS/m for creosote and 4.39 dS/m for mesquite. Thus, soil EC has exceeded the critical level in several parts of study site. Rooting depths of mesquite and creosote are about 12 m and 3 m, respectively. Majority of mesquite roots are distributed within 0–100 cm depth [35], creosote within 0–25 cm depth [36], and perennial and annual weeds within 0–20 cm depth. Visual observation during 2009, which is also one of the classical methods of observing EC [37], showed some leaf burn symptoms and dead bushes of annual and perennial grasses in the southeast and northwest direction of the study area. This is true because the EC1:1 values for these areas were >4.1 dS/m, which might have killed some of the annual and perennial grasses and caused leaf burn symptoms on creosote. Symptoms of high EC were not seen in the mesquite probably due to its deep rooting system and higher EC tolerance.

The EC class I, which has the EC level <1 dS/m, is not a problem for any vegetation including annual and perennial weeds; however this area was decreasing due to the increases in EC in general in the experimental site. EC class II with EC level between 1.1 to 2 dS/m poses no problems for mesquite and creosote but is not favorable for annual and perennial grasses. The EC levels of classes IV, V and some areas of class III are problem for the growth and development of creosote. Similarly, class V and some areas of class IV are not favorable even for mesquite along with other vegetations where mean EC1:1was >4.50 dS/m. One of the strategies to reclaim these soils could be by increasing the amounts of wastewater application in area under class I and no wastewater application for certain period of time in areas under class V and IV. Natural precipitation that usually occurs in the form of thunderstorms could leach the salts in class V and IV. In addition soil EC is directly related to the EC level of wastewater. It has been observed that EC of water in the holding ponds was as high as 5.64 dS/m, so it is necessary to take steps to decrease the evaporation from the holding ponds. Since noticeable changes in the areas under different EC classes were seen within a year, it is necessary to analyze the soil EC after every quarter to adjust the irrigation strategies accordingly. The low variability [38] of sand and clay contents and bulk density of soil indicated that treated wastewater application did not impact these properties. Usually soil textural properties depend on pedogenesis and do not significantly change for a long duration. Some of hydraulic, water retention, and transport properties of soil [39] can change however change at much smaller timescale due to the application of wastewater. This study demonstrated that although chemical properties of the applied wastewater for both years remained similar, the areas under EC classes II decreased and those under classes IV and V increased indicating a cumulative effect of wastewater chemical concentration. Further research should be conducted to determine the spatial correlation of EC with other soil physical and chemical properties.

4. Conclusions

In general, soil EC classes in the West Mesa land application site were moderately variable as well as moderately spatially dependent. Correlogram analyses showed significant autocorrelation for about 25 m lag distances at 0–20 cm and 50 m at 20–40 cm depths for both years. EC values in the northwest and southeast sides have exceeded the critical levels of tolerance of the existing vegetation. Areas under classes I and II were higher than areas under classes III, IV, and V during 2009. During 2010 area covered by class II decreased from 26% to 14.91%, and areas under classes III, IV, and V increased from 33.91% to 45.72%. Immediate steps are necessary to control the increase of higher EC area. If the irrigation pattern stays as is and EC keeps on increasing, the likelihood of creosote and mesquite being displaced by other high salt tolerant desert shrubs like salt bush may also increase. It is necessary to control evaporation in holding ponds and change the wastewater application pattern in such a way that more wastewater is applied in the middle and less in the northwest and southeast portion of the study site while simultaneously monitoring soil EC at least once every four months.

Acknowledgments

The authors are thankful to New Mexico State University Agricultural Experiment Station for support, New Mexico Water Resources Research Institute for funding, and City of Las Cruces for providing help and access to the experimental site.