GUO Yan , ZHOU Yin, ZHOU Lian-qing, , LIU Ting WANG Lai-gang CHENG Yong-zheng HE JiaZHENG Guo-qing

1 Institute of Agricultural Economics and Information, Henan Academy of Agricultural Sciences, Zhengzhou 450002, P.R.China

2 Department of Resource Science, College of Environmental and Resource Sciences, Zhejiang University, Hangzhou 310058,P.R.China

3 Key Laboratory of Spectroscopy Sensing, Ministry of Agriculture, Hangzhou 310058, P.R.China

Abstract Over the past five decades, increased pressure caused by the rapidly growing population has resulted in a reclamation of agricultural and urban buffer zones along China’s coastline. However, information about the spatio–temporal variation of soil salinity in these reclaimed regions is limited. As such, obtaining this information is crucial for mapping the variation in saline areas and to identify suitable salinity management strategies. In this study, we employed EM38 data to conduct digital soil mapping of spatio–temporal variation and map these variations of different site-specific zones. The results indicated that the distribution of soil salinity was heterogeneous in the middle of, and that the leaching of salts was significant at the edges of, the study field. Afterwards, fuzzy-k means algorithm was used to divide the site-specific management zones within the time series apparent soil electrical conductivity (ECa) data and the spatial correlations of variation. We concluded that two management zones are optimal to guide precision management. Zone A had an average salinity level of about 165 mS m–1, in which salt-tolerant crops, such as cotton and barley can grow normally, while crops such as soybean and cowpeas may be planted using leaching and increasing the mulching film methods to reduce the accumulation of salt in surface soil. In Zone B, there was a low salinity level with a mean of 89 mS m–1 for ECa, which allows for rice, wheat, and a wide range of vegetables to be grown normally. In such situations, measures such as an optimized combination of irrigation and drainage, as well as soil amendment can be taken to adjust and control the salt content. Particularly, flattening the land with a large-scale machine was used to improve the ability of micro-topography to in fluence salt migration; rice and other dry, land crops were planted in rotation in combination with utilizing salt-leaching multiple times to speed up desalinization.

Keywords: apparent soil electrical conductivity (ECa), soil salinity, EM38, spatial variation, management zone

1. Introduction

Proximal soil sensors can collect soil information to within 2 m, or implanted in, a soil body. The sensors may actively or passively measure the soil properties with direct or indirect and invasive or non-invasive methods (Viscarra Rosselet al.2010; Guoet al.2013b). In this regard, proximal sensors have two notable advantages over traditional methods.First, proximal soil sensors can cheaply and rapidly acquire soil information. Secondly, they can easily collect enough soil data to conduct meaningful geo(statistical) analyses in order to con firm spatial soil attributes. In recent years,highly detailed soil sensor data have become an interesting source of auxiliary data to detect and characterize withinif eld soil variation. This primarily includes electrical and electromagnetic sensors (EM38; Myerset al.2010; Guoet al.2015), optical and radiometric sensors (ASD, Li Set al.2015; Shiet al.2015), mechanical sensors (soil strength pro file sensor; Chunget al.2006), and electrochemical sensors (ISFETs, Loreto and Morgan 1996). Of these, EM38 has the unique ability to detect soil properties (e.g., texture,water content, and salinity) by measurement of apparent soil electrical conductivity (ECa) with on-the-go or stop-and-go style measurements taken on the basis of electromagnetic induction. ECa mainly responds to soil salinity under saline conditions (Trianta filiset al.2000; Leschet al.2005; Li H Yet al.2015). For instance, a great deal of scientific literatures have well documented the relationship between ECa and soil salinity (Trianta filiset al.2000; Moralet al.2010; Eldeiry and Garcia 2011; Guoet al.2013a). Where yield correlates with ECa, maps of ECa are helpful for devising soil management or agricultural practices. Other important arguments in favor of the use of EM38 sensors include the robust construction of the equipment, their ease of use, and most significantly,the nondestructive nature of their measurements (Eigenberget al.2002; Corwin and Lesch 2003).

In reclaimed coastal areas, soil salinity is the main limiting factor for agricultural productivity. In such situations,farmers and land managers need to be able to identify the spatio–temporal variation of saline areas for the sake of understanding the agricultural and economic bene fits.With regard to identifying and mapping the spatio–temporal variation of soil salinity, geostatistical tools with spatial georeferences are universally used. For example, Douaiket al.(2005) used spatio–temporal kriging and Bayesian maximum entropy to predict soil salinity at unobserved spatial locations and time instances; Ali Aldabaaet al.(2015)quanti fied soil salinity at two playas in West Texas, USA with a combination of remote and proximal sensor data.

Reliable, timely, and cost-effective soil salinity monitoring and assessments are needed for an understanding of the dynamics of soil salinity spatio–temporal changes (Guoet al.2013a). However, there are very few publications on the spatio–temporal variation of soil properties in reclaimed coastal regions with continuous monitoring over multiple years with proximally sensed data. Hence, the main objectives were to: (1) characterize the spatio–temporal distribution of soil salinity by ECa in a reclaimed field; (2)analyze the spatio–temporal stability of soil salinity and map the temporal stability for better soil salinity management; and(3) use the fuzzy cluster means to divide the management zones and map the spatial distribution.

2. Materials and methods

2.1. Study area

The study area is located in the northern region of Shangyu City, Zhejiang Province which is southeast of China’s Hangzhou Gulf of the Yangtze delta. Approximately 17 000 ha of coastal land has been reclaimed around Shangyu City since the 1960s (Fig. 1). The study site is characterized by a subtropical climate, with a mean annual temperature of 16.5°C and annual precipitation of 1 300 mm. The soil is derived from recent marine and fluvial deposits. At present, the dominant land use is for rice farming, vegetable cultivation, and growing nursery plants. The experiment was conducted on a 4.25-ha paddy rice field which was reclaimed in 1996.

Fig. 1 Location of the study field in Hangzhou Gulf and the apparent soil electrical conductivity (ECa) sampling points.

2.2. Data collection and processing

In this study, we used the program “EM38 proto record data”with the Geonics EM38 Data Logging System consisting of a field computer, Allegro CX. Data files created with this program were used to position a survey according to the locations recorded separately by a Global Positioning System (GPS). These results can be combined with EM38 records through NMEA-0183 compatible data (McNeill 1980).

ECa measurements were acquired along an approximate 20-m grid after the rice was harvested and the field was drained in 2009, 2010, and 2011. There were 251, 256, and 339 ECa measurements with GPS gathered, respectively.In order to calculate the coef ficient of variation over time,ECa measurements in 2010 and 2011 were harmonized onto a common grid based on the 251 ECa measurement locations in 2009 (Fig. 1) by the nearest neighbor algorithm.

2.3. Spatio–temporal variation analysis of ECa

The geostatistical method is a subset of statistics that specializes in analysis and interpretation of spatially (and temporally) referenced data, often employed to estimate the values of variables between data points and quantify the variance structure, to characterize the spatial distribution,and to track trend changes of soil attributes (Maliva 2016;Silaet al.2017; Zhanget al. 2017). This approach is used to first identify and analyze the spatial structure of the variables of concern (i.e., ECa), and then to interpolate the values of variables from neighboring values, taking into consideration their spatial structure. In the present study, a semi-variogram was applied to quantify the spatial variation of soil salinity. It can be defined by eq. (1) (Webster and Oliver 2007):

Where,γ(h) is empirical semi-variogram, and estimates the averaged variation between two ECa pointsxandx+h;Z(xi) andZ(xi+h) denotes the values of the variableZ(ECa)at locationsxiandxi+h;N(h) is the set of pairs of ECa observations separated by the lag distanceh.

The predication ECa value is interpolated by ordinary kriging (OK) which can be expressed by eq. (2) (Webster and Oliver 2007):

Where,Z＊(x0) is interpreted as a random ECa value at locationx0, as well as the values of neighboring samplesZ(xi);λiis the weight assigned to the observation Z(xi);nis the number of ECa observations.

Temporal stability analysis is a common method employed to identify spatial variation (Wanget al.2015).In this paper, the temporal stability was expressed with the coef ficient of variation (CVti) calculated over time at each measurement. The method was used by Guoet al.(2015)to characterize the soil properties in coastal areas.

Where,CVtiis the coef ficient of variation over time at theith observation of ECa in thetth time;nis the ECa numbers.

2.4. Fuzzy k-means (FKM) algorithm

The FKM algorithm is an unsupervised clustering method.During the process, each training variable can be associated to more than one quantization point cell, with some degree of pertinence to each point cell. It has been widely used in climatology, remote sensing, geology, environmental studies, computer science, and soil landscape mapping(Gorsevskilet al.2003; Liet al.2008). In this study, we used the FKM algorithm to con firm the number of classes within ECa and calculated the coef ficient of variation.Trianta filiset al.(2003) gave a detailed description of the FKM approach. The aim is to minimize a set of data points into self-similar groups such that the points that belong to the same group are more similar than the points belonging to different groups. In brief, it calculates a measure of similarity between an individualiand a clustercin multivariable space (Bezdek 1981). When the objective functionJ(M, C) is minimized, the best outcome is:

Where, M=micis ann×kmatrix representing the degree of membership values (nis the number of objects); C==(ccv)is ank×pmatrix representing the existing centroids (pis the number of variables;ccvis the value of the center of classcfor variablev;xi==(xi1,…,xip)Tis a vector representing the individuali;cc==(cc1,…,ccp)Tis the vector representing the center of classis the square distance betweenxiandccaccording to a chosen definition of distance.

The distance function (dic) is applied to measure similarity or dissimilarity between two individual variables, and then later between two clusters. In this study, Mahalanobis distance was used as the metric. It is defined as follows:

The two above formulas form the mathematical basis for the FKM algorithmJ(M,C) is minimized iteratively. The given solution can be considered as an optimally stable way for the clusters when it achieves a minimum value (Bezdek 1981). During the processing, the fuzziness performance index (FPI) and the modi fied partition entropy (MPE) were used to determine a suitable value forkas a function of the fuzziness exponent (ϕ). In the FKM algorithm, the exponentϕ(1＜ϕ＜∞)determines the degree of fuzziness, withϕas a real number. The closerϕis to in finity (∞), the greater the fuzziness of the solution and the closerϕisto 1, then the solution becomes increasingly similar to the clustering of binaryk-means. The goal ofJ(M, C) is to find an optimal balance between structure and continuity; choosing a value ofϕthat maximizes the objective function generates the‘hardest’ fuzzy clustering solution. McBratney and Moore(1985) devised the measure of fuzziness for determining the objective function –dJ/dϕ.Therefore, choosing the optimal combination of classes and exponent with the least maximum of –dJ/dϕis established.Detailed information of the two indices was discussed by McBratney and Moore (1985) and Odehet al.(1992). The FKM analysis was completed using the FuzME Program (Minasny and McBratney 2002).

3. Results and discussion

3.1. Summary statistics of EM38 data

Table 1 shows summary statistics of ECa measured by EM38 sensor in 2009, 2010, and 2011. The averaged ECa values show a substantial decreasing trend from 166.19 mS m–1(2009) to 134.02 mS m–1(2010), with a smaller decrease from 2010 to 2011 (20.73 mS m–1). The quartile estimates also show a decreasing trend from 2009 to 2011,with a smaller difference between 2010 and 2011 than from 2009 to 2010. Field-scale variation is characterized by the coef ficient of variation (CV). In general, ECa from 2009 to 2010 shows increasing CVs, with CVs of 26.44 and 49.06%,respectively. This is most likely due to the presence of rice cultivation with a continuous circle of irrigation and drainage operations, which can leach the salts into a deeper soil pro file or into ditches. Given the negative skew in the ECa,and that the means are smaller than the median values, this suggests that the landscape is the predominant variable.

In order to characterize these differences, we chose the Tukey-Kramer multiple comparison procedure. Fig. 2 shows the plot of the calculated Tukey-Kramer means and the analysis of mean for ECa (α=0.05). The means for ECa in 2009 and the ECa in 2011 are statistically different from the overall mean. The center line indicates the overall mean,and the lines above and below the mean are the upper decision limits (UDL) and lower decision limits (LDL). If a group mean falls outside of the decision limits, then that mean is significantly different from the overall mean. If a group standard deviation falls outside of the decision limits, then that standard deviation is significantly different from the root mean square error. Thus, the ECa mean in 2009 is lower and the mean in 2011 is higher. That is to say, ECa was significantly different from 2009 to 2011, and 2009 to 2010.

3.2. Digital soil mapping of spatio–temporal variation of soil salinity

Semi-variance analyses were carried out to express the spatial dependence of the three ECa datasets. These processes were carried out using geostatistical software,GS+ ver. 9.0 for Windows (Gamma Design Software,Plainwell, MI). Based on the regression coef ficient of determination (r2), the model with the least sum of squares was chosen as the “best fit.” The semi-variogram models for the ECa across three years demonstrated that the spatial behavior had good continuity in space, and the three datasets can be modeled quite well with exponential models. Fig. 3 shows the plot of experimental Semivariances and the fitted semi-variogram models for the ECa between the years 2009–2011. The parameters for the three models are shown in Table 2. The “nugget” (C0) was the random variation that usually results from the inaccuracy of measurements or variations of the properties that cannot be detected in the sample range; the “sill” (C0+C) was the total ECa variation at which the semi-variogram levels for the patterned data; the structure variance (C) was the difference between the sill and nugget; the “range” was the maximum distance over which the measured ECa exhibited significant spatial autocorrelation. Generally, the ratio relatingC0to sill(C+C0), which is an especially important parameter, can be used to characterize the spatial dependency. If the ratio is less than 25%, this implies a strong spatial dependency of the observations; if the ratio is between 25 and 75%, the ratio indicated moderate spatial dependency; if the percentage is greater than 75%, it indicated a weak spatial dependency in the present sampling resolution (Šamonilet al.2016). As presented in Fig. 3 and Table 2,C0decreased from 2009 to 2011, indicating that the variations for ECa became smaller.Moreover, the ratio ofC0to (C+C0) also declined sharply,from 17.07% in 2009 to 0.26% in 2011. Following this, we concluded that the autocorrelation of ECa was becoming stronger, which may have been induced by rice cultivation with alternating irrigation and drainage farming.

Table 1 Descriptive statistics of soil electrical conductivity (ECa, mS m–1) in 2009, 2010 and 2011

The spatial distribution characteristic of ECa contributed to assessing the extent and level of soil salinity. We used kriging interpolation to map the variations of ECa across these three years. Several homogenous zones were divided according to similar values of ECa. The smoothed contour maps obtained for the three years are presented in Fig. 4-A–C. The three smoothed contour maps of ECa express quite similar patterns, which demonstrate that the salinity distribution has a similar trend and represents consistently high- and low-salinity areas of the field. In particular, the kriging maps of each ECa indicates a high ECa level (＞175 mS m–1) in the middle and a low ECa levels (＜100 mS m–1) in the surrounding regions. ECa is decreasing within the operational years for all ECa levels.

Fig. 2 One-way ANONA analysis of soil electrical conductivity(ECa, mS m–1), in 2009, 2010, and 2011.

For example, the maximum ECa value is 181.8 mS m–1in 2011 versus 226.7 mS m–1in 2009, with the minimum ECa value also decreasing from 2009 (51.3 mS m–1) to 2011(10.5 mS m–1). The high ECa level is distributed in the middle, while the low ECa level is found in the surrounding fields. These are mainly caused by anthropogenic activities like tillage by large-sized tractors which tended to result in uniform field surface topography, irrigation, and drainage for rice cultivation, which leach the surface salt into deeper soil pro files and ditches which channel around the field.

In order to assess the temporal stability, the coef ficient of variation (CVti) was calculated over time at each ECa survey location (Fig. 4-D). Based on this, we consider the variation to be stable whenCVti＜10%, moderately stable whenCVtiis between 10 and 25%, and unstable whenCVtiis between 25 and 100%, according to Shiet al.(2005). Interestingly, we found that the middle region with a higher ECa value displayed temporal stability, in contrast to the surrounding region with a lower ECa value that displayed temporal instability. This is consistent with the reports by Shiet al.(2005). It is not readily apparent why these performances occurred, and this matter requires further investigation and analysis.

3.3. FKM algorithm with sequential ECa and CVti

Fig. 3 Variograms for soil electrical conductivity (ECa, mS m–1)in 2009, 2010, and 2011.

Table 2 Models and parameters of semi-variogram for soil electrical conductivity (ECa, mS m–1) in 2009, 2010 and 20111)

Fig. 4 Maps of soil electrical conductivity (ECa), in 2009 (A), 2010 (B) and 2011(C) with the plot (D) of coef ficient of variation(CVti) over time.

The detection of variation in salinity can help farmers be aware of the regions within their fields where they may be experiencing poor yields due to high salinity or management conditions. The time series ECa data and the spatio–temporal correlations of variations were imported into FuzME software for management zone analysis to classify the four datasets into MZs, where FKM algorithm was performed. The values of FPI and MPE are plotted against the number ofkandϕin Fig. 5.

The values of FPI and MPE againstk(Fig. 5-A and B)showed that the gradient increased first askincreased from 2 to 3 (depending onϕ), and decreased askincreased from 3 to 4. After this, the gradient slowly increased askcontinued to increase up to 8 whileϕincreased from 1.2 to 1.8. However, all values of FPI and MPE are higher thank=2. This suggests that the sensitivity of FPI and MPE toϕwas the highest with aϕof approximately 1.8. Given the ambiguous nature of the results, the relationship ofϕversus-dJ/dϕwas plotted in Fig. 5-C. The function’s sensitivity toϕwas related not only tokbut also toϕ. When choosingϕ, the value considered most suitable is when –dJ/dϕis at a maximum. The value of –dJ/dϕshows local minima fork=2 for both the FPI and MPE withϕ=1.8.

3.4. Statistical comparison of FKM classes and site-specific management of soil salinity

In classifying the serial ECa measurements andCVtiintok=2 classes, we have more or less imposed a “blocking” structure whereby discrete, compartmentalized units are evident.To make ensure the classification result could effectively denote spatial variation of soil salinity, we first conducted statistical comparisons of the two classes andCVti,with the resultsshown in Table 3. There was greater standard deviation and standard error in class B than in class A. The 95% con fidence interval in class B was far more distant to the mean than in class A. This indicated that the variations in class B are larger than in class A. Afterwards,t-test was used to validate the classes, demonstrating that significant differences existed between the two classes (Fig. 6-A). This indicates that the proposed solution is optimal in terms of the organization of the classes with respect tok=2 withϕ=1.8.

Table 3 Statistics of the two classes for soil electrical conductivity (ECa, mS m–1) in 2009, 2010, 2011 and the coef ficient of variation CVti (%)

Moreover, the significance of class centers does not prove in itself the validity of fuzzy continuous classification from the point of view of salinity-forming processes, which need to be supported by a spatial distribution of memberships consistent with the soil salinity-forming processes of the surveyed environment. Confusion index (CI) was used to display the spatial structure of membership variation.The graphical representation in Fig. 6-B indicated the complex variation. Here, a membership value (more than 0.5) was employed to qualify each of the ECa locations to membership of classes A or B. Class A had the larger number of members (198) and class B had the smaller(53). Class A consisted of approximately four times as many members as B. The means of the auxiliary data within the two classes show that class A has the larger ECa(e.g., 184.35 mS m–1for 2009) and smallerCVti(15.68%).Conversely, class B has the smaller ECa (e.g., 98.35 mS m–1for 2009) and largerCVti(36.59%).

As above, two site-specific management zones are applicable in this study field. Digital mapping of the two management zones can be seen in Fig. 6-C. In zone A, with an average salinity level of about 165 mS m–1,salt-tolerant crops, such as cotton and barley, can grow normally, while soybean and cowpeas can be planted with additional treatments including leaching and increasing the mulching film methods which reduce the accumulation of salt in the surface soil. By contrast, in zone B, the salinity level was low, with an ECa mean of 89 mS m–1, which allows for rice, wheat, and a wide range of vegetables to be grown normally. In such situations, measures such as combinatorial optimization of irrigation and drainage and soil amendment can be taken to adjust and control the salt content. Flattening the land with a large-scale machine was particularly useful when employed to improve the in fluence of micro-topography on salt migration; the rotation of rice and dry, land crops with multiple applications of salt-leaching can speed up desalinization.

Fig. 6 Significance testing of the two classes, and the spatial distribution of membership and management zones. A, t-test for the two classes. ＊ indicates the range from 1 to 99% of ECa and CVti values. Different letters indicate significant differences in level of 0.05. B, confusion index (CI) of the spatial structure of membership variation. Letters A and B indicate ECa locations belong to classes A or B. C, spatial distribution map of soil salinity management zones.

4. Conclusion

In this study, geostatistical tools combined with a fuzzy clustering algorithm were used to identify spatio–temporal features of soil salinity by proximally sensed EM38 data in a coastal field to assess the spatio–temporal variation of soil salinity. We found that soil salinity decreased with the passage of time. This conclusion is of importance for farmers who can take measures to decrease soil salinity content to achieve high crop productivity and sustain economic returns. We mapped the temporally stable middle region with a high salinity content, as opposed to the surrounding region with a lower salinity level that showed temporal instability. Afterwards, fuzzy-kmeans were used to divide the field into site-specific management zones,finding that two management zones are optimal to guide precision management. In this situation, differentiation into two management zones may be very helpful for farmers to adopt site-specific management, which satis fies the criteria that management zones be simple, functional, easy to understand, and economically feasible. To reiterate,different crops and measures can be used in the two types of management zones.

These case study results con firm that soil salinity substantially changes with the tillage years and the particular agricultural practices used in rice paddy fields.This is likely due to rice planting with continuous irrigation and drainage. We will conduct further research in other crop types. In addition, our findings help increase the understanding of how salinity changes with time, allowing us to better understand the in fluence of these changes on crops in various conditions, and can help direct the design of soil sampling.

Acknowledgements

This material is based upon work funded by the National Natural Science Foundation of China (41601213), the National Key Research and Development Program of China(2017YFD0700501) and the Major Science and Technology Projects of Henan, China (171100110600).