Optimal Irrigation Mode and Spatio-Temporal Variability Characteristics of Soil Moisture Content in Different Growth Stages of Winter Wheat

To advance site-specific management of soil volumetric moisture content (VMC), this study analyzed the spatio-temporal evolution characteristics of soil VMC using the method of sequential Gaussian simulation (SGS) during the different growth stages of winter wheat. This was compared with data measured by time domain reflectometry (TDR) which is a well-established electromagnetic technique to measure soil VMC. The spatial autocorrelation coefficient of VMC indicated the strongest clustering of VMC in the tillering stage, and the least clustering of VMC in the harvest stage. A threshold of VMC in topsoil in the jointing stage of winter wheat was put forward. This threshold is 26, signifying that at a lower value, irrigation should be performed and irrigation efficiencies can be improved. Stable and sub-stable areas in the spatial variability maps of VMC were identified in the winter wheat jointing stage. Furthermore, the optimal irrigation stage was the early jointing stage, and irrigation was performed once as a guide. A loose-couple spatial model was constructed using the VMC in topsoil and the volume of water for irrigation. The VMC in the jointing stage of winter wheat was linked with efficient and water-saving irrigation.


Introduction
Soil moisture is a fundamental property affecting plant growth, transport, and transformation of soil nutrients, water, and energy budgets in the soil-crop system.Soil moisture content is used for determination of the soil-water balance and precision irrigation management.It is not only an important variable that affects atmospheric and hydrologic process by influencing the partitioning of incoming radiation into sensible and latent heat flues and by separating precipitation into infiltration and surface runoff but is also a state variable controlling a wide array of ecological, hydrological, geotechnical, and meteorological processes.With the rapid development of China's economy and increasing population, the water available for agriculture is becoming increasingly scarce [1,2].Limited water resources often result in reduced crop yield.In the Yellow-Huai River Plain, winter wheat (Triticum aestivum L.)/summer maize (Zea mays L.) rotation is a staple planting system and crop production is generally dependent on rainfall and irrigation water.Therefore, improving water use efficiency, the balance between soil moisture in topsoil and the volume of water for irrigation, and narrowing the gap between winter wheat water demands and soil water supply are of great significance in improving agricultural production [3][4][5].Soil moisture, especially soil water in the topsoil, is one of the controlling factors for hydrological, meteorological, and ecological processes [6].The stability of soil moisture content, irrigation volumes, and the appropriate timing of irrigation are critical for developing effective agricultural water management practices and improving agricultural water use efficiency [7].For winter wheat growth stages, there are different levels of sensitivity to water demands in the tillering, jointing, heading, milky dough, and harvest stages across the entire winter wheat growth season.Obviously, large-scale in situ monitoring of soil moisture is both time-consuming and costly across the winter wheat growth stage.Accurately predicting the regional distribution of soil moisture content based on point-scale data is of important interest for the guidance for soil and water conservation and decreases costs and required time.Spatio-temporal characteristics of soil moisture content can be analyzed using these measured data [8,9].The spatial variation of soil moisture contains critical implications for the processes as soil moisture varies greatly across space and time based on different soil properties and field topography.An understanding of the spatio-temporal variability of soil moisture is instrumental in quantifying the linkage between regional hydrology, ecology, and physiography.This spatio-temporal variability model of soil moisture content is more applicable to different growth stages of soil moisture and irrigation water management with winter wheat cropping systems.Besides, this model is also readily applicable to the endogenous and exogenous factors influences on winter wheat growth and the optimal irrigation mode [10].
Spatial analysis is an effective and worthwhile tool for the estimation of the spatial pattern of the regional ecological variables [11][12][13].Spatial interpolation and estimation methods are commonly used for estimating soil characteristics of sampled specific-site data [14].A water-saving planting scheme with leaf area index (LAI), dry matter accumulation, and grain yield of winter wheat was developed in the Yellow-Huai River Plain in winter wheat growing seasons [15].Spatial distribution of light interception (LI) of the canopy was explored using a geostatistical sampling-based method [16].Researchers developed time-sequenced risk assessments of soil moisture scarcity in an oasis, and further reported that stochastic simulation was an effective a technique by exhibiting the possible spatial variations of soil moisture [17].
Many of soil properties that are relevant to precision agriculture vary both in space and time.Understanding spatial soil variability is therefore a crucial component for precision irrigation [18].The conventional approach using soil survey and dense sampling is the most accurate but analyzing a large number of samples is time-consuming and a major financial and resource constraint.Precision irrigation can aid in improving crop productivity and increasing water and energy efficiency in irrigated agriculture has been studied for more effective use of rainfall and irrigation water to help farmers with the capability to demonstrate environmentally sustainable practice [1,19].Precision irrigation offers the potential to eliminate over-irrigation and apply water in a deliberate non-uniform or variable manner in response to the specific irrigation requirements of different discrete management units, hence maximizing crop response and minimizing any adverse environmental impacts [20].The primary objective of optimizing the spatial scale and timing of irrigation applications is therefore to increase the crop yield, and simultaneously reduce fertilizer losses.The AquaCrop model, as water-driven cultivation model, simulates the yield response of herbaceous crops to water and is particularly well-suited to conditions in which water is a key limiting factor in crop production [21,22].The AquaCrop model of winter wheat area should be more significant in the future work.
It is impractical to collect soil samples for adequate assessment of soil moisture that represents both spatial distributions at the field-scale and temporal dynamics across the growing season because it is laborious and time-consuming.Consistency and early assessment of soil moisture is vital for crop production.An optimal spatial sampling scheme to reduce sampling density and estimation of unsampled values can save time and money significantly [23].However, the effectiveness of sampling relies on the accuracy of the spatial interpolation used to determine the spatial variability.Geostatistics provide descriptive tools such as kriging to directly implement property prediction at an unsampled location according to known data points within a local neighborhood; they can be used as an alternative to predict regional variables [24,25].Grid-based sampling increases with the variogram range and spatially structured variance [26], and the prediction accuracy of soil moisture may improve with advances in computerization.Summarily, the method has the drawback of smoothing effects for the prediction of regional variables.Use of geostatistical methods can optimize irrigation management [1].
Soil moisture content is not only the main factor for crop production in the soil-crop system, but also an important factor in the evolution of the soil physical, chemical, and biological processes.In particular, topsoil moisture content is a key factor in crop establishment and subsequent production, soil detachment, runoff generation, and other soil physical, chemical, and biological processes.The spatio-temporal evolution of soil properties is an important topic in research.In prior research, the variability of soil moisture content was studied, a point-scale distributed parameter model was constructed, and the temporal stability of the spatial pattern of soil moisture content from the moisture zones was investigated [27,28].For soil nutrients, the spatial distribution of total potassium was explored, and the results showed that available potassium and total potassium were distinctly influenced by different agricultural practices [29].Using 58,559 soil samples, the spatial variability of soil available potassium (K) was studied and the variability was confirmed by the related yield results.Moreover, the anisotropy of soil available phosphorus (P) was analyzed as well [30,31].Soil organic carbon (SOC) plays an important role in soil quality and its conservation.The spatial variability of soil carbon was studied by spatial downscaling from the observed point data in Australia.The researchers mapped a farm field with a 100-m grid resolution to handle the prediction uncertainties of this spatial variability prediction map, which improved soil carbon dynamic precision monitoring at the farm scale [32,33].Soil nutrient supply capacity is important for assessing crop yield, and the relationship between soil nutrients and crop yield was also studied [34].Time domain reflectometry (TDR) has been widely used by the scientific community as a reliable method to indirectly measure the volumetric water content of soils; it has become an important tool for measuring soil moisture [35,36].
In contrast to the method of kriging (the geostatistical method) or Gaussian process regression, the conditional simulation emphasizes the accuracy of the simulated variable, reproducing its sample statistics (mean, histogram, covariance).Its pattern of spatial continuity, conditional simulation, also provides measures of local and global joint accuracy.The key difference between kriging and sequential Gaussian simulation (SGS) is in addressing the variance; that is to say, kriging loses variance.SGS is the most straightforward algorithm for the simulation of continuous variables, and it is also a process using a predefined probability distribution.In other words, a probability distribution function (PDF).This kind of simulation is regarded as the part of the Monte Carlo approach, which aims to solve mathematical problems by modelling random variables based on the generation of a large set of random numbers.Simulation accuracy is improved by sampling data with SGS [37][38][39].As an important conditional simulation method, SGS is a well-known approach based on the decomposition of the multivariate probability density function of a stationary and ergodic random process to the product of univariate posterior distribution functions [40,41].It has also been used to advance interpolation techniques [17,42].As a very powerful simulation technique, SGS can handle many different simulation problems, such as the simulation of a single variable, or the joint simulation of multiple correlated variables.Spatial heterogeneity is a natural feature of ecosystems.SGS has the advantages of representing high spatial entropy and producing more realistic output with spatial variability, and can be conducted with a joint probability model using the following formula: where Z 1 , Z 2 , . . ., Z N are the observation values at the sequence points U 1 , U 2 , . . ., U N ; P is the probability function; Z(u i ) is the simulation value of the variable Z at u i ; and SGS is a probabilistic approach that provides a joint distribution of multiple variables with equally probable realizations of the joint distribution of one or more property values spatially.SGS generates a model of uncertainty based on the Bayesian theory [43] describing the conditional or posterior probability of two events: This equation can be read as the probability that event A 1 will occur given event A 2 is equivalent to the ratio of the joint probability between the two events and the individual probability of an event.It can be expanded into multiple events.The events can be assigned to the nodes or cells to simulate on the map.The conditional probability distribution function (CPDF) is an important component in the SGS process.During the process, a one-point CPDF is modeled and sampled at each of the sampling points, which are named after nodes visited along a random sequence.To ensure the reproduction of the variable covariance and mean model, each one-node CPDF is made conditional not only to the original observed data but also to all values simulated at previously visited locations.In SGS the interpolations are based on a form of stochastic simulation in which data values are honored at their locations.This means that local details are not obscured by smoothing as they are in kriging.Therefore, we chose the SGS method to conduct the volumetric moisture content (VMC) prediction and simulation instead of kriging.
The previous study of soil digital mapping could be enhanced by the consideration of predictive variables with multivariate analysis [44].The related soil moisture measurement demonstrated that the accuracy and utility of soil maps for site-specific management can be improved using functional soil mapping approaches [45].Soil moisture measurements with increased temporal resolution have also been studied [46].
As summarized above, few studies have found spatial relationships between the ecological variables, and there are none on the variability of soil moisture across the whole winter wheat growth stages in the area of the north China plain.No study of a spatially-associated model in the soil-crop system was found with respect to soil characteristics such as soil physical/chemical properties, mineral characteristics, and crop growth indicators such as crop height, LAI, and the normalized difference vegetative index (NDVI) derived from remote sensing data.In this study, the spatial variability of soil volumetric moisture content (VMC) over a region in Shandong Province in the north China plain was investigated using SGS to quantify the spatial distribution of VMC in the different growth stages of winter wheat.To provide new methods and information for preserving water resources with the strengthened management of soil moisture in the farm, the data of VMC have been collected for a case study to fulfil the goal of the study to evaluate the spatial variation of VMC using the method of stochastic system simulation.The specific objectives are as follows.Firstly, we aim to investigate the spatial dependence of VMC in growth stages of winter wheat with the sampled soil data.The second objective is to analyze spatially varied characteristics and determine the storage amount of VMC using the method of SGS with the sampling data measured by time domain reflectometry (TDR) at different wheat growth stages.Thirdly, based on SGS simulation, we aim to explore the relationship between VMC in topsoil and the volume of water for irrigation, distinguishing stable and sub-stable irrigation areas for winter wheat, and combining the growth stage with the spatial variability maps of VMC for these areas.We will then discuss the optimization irrigation growth stage and the optimal irrigation water times.

Study Area and TDR Field Sampling
The study area (117 • 04.130 E, 36 • 42.979 N), covering approximately three hectares, is located in the north side of the Xiaoqinghe River in Shandong province, China.The core area is located in the Yellow-Huai River Plain.The area is in the warm temperate zone in the half moist monsoon climate region.In this area the terrain slope is high to low from south to north and the soil is sandy clay loam and the wheat-maize rotation was conducted yearly in the crop fields in the area.
For soil moisture content measurement, in each sampling site TDR probes were placed vertically on the sites of the sampling grid to monitor the VMC in topsoil.In our experiment, soil VMC was the ratio of the volume of moisture in a given volume of soil to the total soil volume, expressed as either a decimal or a percent.Specifically, the data logger's LCD (liquid crystal display) screen from the TDR instrument displayed the soil VMC data in "standard mode" when the probes were done, and the meter's built-in data logger could record data from different observation sites and eliminate the need to record data manually.At each sampling point, five VMC measurements were made within a 1-m diameter circle.The average reading for each point was used as a VMC datum point and displayed on an LCD screen.Through the software (a TDR-specific program), we were able to transform the soil VMC data from internal data logger measurements into a display on a personal computer screen.Based on TDR application instruction, the file which included VMC data was also saved; generally, data files are stored as a comma-delimited text files and may be viewed in text editor or spreadsheet software.In the study, we chose to save the soil VMC data using the Excel formation file to further prepare for the statistics characteristics of the soil VMC.Each VMC measurement was geo-referenced using a Differential Global Positioning System (DGPS).The sampling date was selected in the wheat cultivation season, and 499 VMC samples were collected from the surface layer (0-20 cm) in the tillering stage, jointing stage, heading stage, milky dough stage, and harvest stage during winter wheat growth in 2014.Detailed records were made on the coordinates of the sampling sites.Geological positions of sampling sites were transformed into a point-type shapefile (.shp), which was compatible with ArcGIS (ESRI, Redlands, CA, USA), and projected into plane coordinates in meters.The field sampling map was created using the method of SGS with spatial autocorrelation combined with GIS data.All statistical analyses, including descriptive statistics and validation, were calculated using SPSS 18.0 statistical software (IBM, Armonk, NY, USA).

Spatial Dependence of VMC and Sequential Gaussian Simulation
The correlation structure of a measurement in spatial autocorrelation provides insight into its spatial structure and describes how it changes in space in relation to the sampling frequency.Spatial autocorrelation is defined as the association of the values of one variable with the values of the same variable at all other localities, and spatial autocorrelation reflects the degree of spatial clustering and spatial heterogeneity, which depends on the degree to which values at one spatial locality are determined in part by the values at neighboring spatial locations.Moran's I can be used in a wide variety of circumstances and is one of the most commonly used measures for the analysis of spatial autocorrelation for regional variables in soil science and landscape ecology [47][48][49][50].Spatial autocorrelation reflects the degree of spatial clustering.Moran's I indicates not only the existence of spatial autocorrelation (positive or negative) but also the degree of spatial autocorrelation.So far, the most used global index for spatial autocorrelation is the Moran's I index [51], which is calculated using the following formula: where n is the total number of area units over the entire sampling area observed; X i and X j denote the data values at locations i and j, respectively; and W is a n × n spatial weight matrix where W ij takes 1 if the two locations are adjacent each other, otherwise it is 0. I ∈ [−1, 1].I = 0 stands for all independent variables, I > 0 stands for that they are positive correlation, and I < 0 stands for a negative correlation.Moran's I coefficient reflecting spatial dependency embodied with the similarity of nearby observations.Spatial autocorrelation can be defined as the coincidence of value similarity with site similarity.The Moran's I autocorrelation coefficient was used to describe the spatial dependency and the degree of correlation of the location-related attributes.In Moran's I autocorrelation coefficient, as the name suggests, autocorrelation is the correlation of a variable with itself.To assess the spatial dependence between VMC of samples, in the study we used Moran's I spatial autocorrelation coefficient.Moran's I autocorrelation coefficient of VMC was calculated by ArcGIS (ESRI, Redlands, CA, USA).The spatial dependency of sampling points separated by lag distance/separation distance h is given by the experimental variogram.It is expressed by the semivariance γ(h) at a given lag h.The semivariance model is shown by Formula (4): where γ(h) denotes semivariance for interval distance class h; Z (x i ) denotes measured sample value at point i; Z (x i +h) denotes measured sample value at point i + h; and N(h) denotes total number of sample couples for the lag interval h.According to the semivariance formula, the maximum semivariance is called the sill (C 0 + C); the sill variance comprises any nugget variance (C 0 ) and the spatially correlated variance (C).The finite distance at which some variograms reach their sill is the range (A 0 ).The semivariance model and SGS of VMC was conducted using geostatistics for the environmental science (GS+) and ArcGIS (ESRI, Redlands, CA, USA).A GS+ file on SGS results with respect to VMC was converted to an ArcGIS grid file using ArcToolbox.
SGS is the most straightforward algorithm for the simulation of continuous variables; the method can produce large fields of values because it is not constrained to keeping the covariance matrix of the entire domain in memory.In fact, it reduces the demand on memory by limiting the neighborhood using a specified search radius and/or a limited number of points to be included in the kriging system [47].SGS helps to optimize site-specific management, and the outcome can have immediate practical implications.In the case study, the specific steps for the node random sequences are as follows:

•
The first step, based on all available samples, involves computing the representative CPDF of the entire study area, and a normal-score transformation of the original data into Gaussian-restricted data.

•
In the next step, a random path is defined to visit each node in the grid once, with a specified number of items of neighboring conditioning data including both initially transformed data and previously simulated grid node values (the conditional restriction).Randomly drawing a value between 0 and 1 from the cumulative Gaussian distribution, the simulated value is added to the current dataset.

•
In the next step, simulation at different sampling locations, simulation of the next node occurs, looping around until all nodes are simulated (the sequential restriction).The simulated normal values are transformed back to the original variable for each location [43,48,52,53].

Results
In this study, we computed the Moran's I spatial autocorrelation coefficients of VMC for the stages of tillering, jointing, heading, milky dough, and harvest across the entire winter wheat growth season.The Moran's I coefficient describes the spatial dependency characteristics of VMC through spatial variability.The results indicated that the spatial dependency characteristics of VMC decreased step by step from the tillering stage to the harvest stage in the range of 0.7 to 0.2 on the whole.The Moran's I coefficients of VMC changed from the tillering stage through harvest stage for winter wheat growth stages, suggesting the more complicated tendency of Moran's I for VMC is profoundly affected by the random factors such as mean annual warming and rainfall.On the other hand, a different spatial structure and spatial dependence of VMC are presented, reflected by the different Moran's I value from the tillering stage to harvest stage of winter wheat.The larger the value of Moran's I, the greater the degree of spatial clustering degree.The strongest clustering of VMC was thus in the tillering stage, whereas the least clustering of VMC was in the harvest stage.Data normality analysis was necessary before semivariance modeling.After the normality analysis of VMC data, the parameters of semivariance models of VMC were estimated (Table 1).All the nugget values are greater than 0, which implies that the positive effect of VMC resulted from sampling error and short distance variability.The largest range value (75.9 m) indicated that the large block scope variability was mainly distributed in the study area in the jointing stage, while the lowest range value (6.9 m) indicated the low scope variability in the milky dough stage.The measured variables had a non-zero nugget effect, which was caused by the measurement error or micro-variability that can be detected at the sampling area.To compare the effect of different nugget sizes among variables and models, the effect can be expressed as a ratio of nugget to sill, and it is possible to divide the spatial dependence pattern into three groups based on the ratio between the nugget and sill.Generally, when the ratio is less than 0.25 the spatial correlation is strong; a ratio between 0.25 and 0.75 indicates a moderate spatial correlation; and a ratio is greater than 0.75 and less than 1 indicates a weak spatial correlation.The ratio of 1 exhibits a pure nugget effect, indicating no spatial correlation at the sampling area.In establishing the spherical and exponential semivariance, our study showed that the ratios of nugget and sill of VMC in the tillering stage and harvest stage were 0.44 and 0.50, respectively, which suggested a moderate spatial autocorrelation.The ratios of nugget and sill of VMC in jointing stage, heading stage, and milky dough stage were 0.24, 0.13 and 0.001, respectively.The results suggested the strong spatial autocorrelation of soil moisture.The coefficients of determination of the models were 0.97 (jointing stage), 0.96 (harvest stage), 0.91 (heading stage), 0.82 (tillering stage), and 0.35 (milky dough stage), indicating that the models at the joint, harvest stage, heading and tillering stage had goodness of fit, while the model of the milky dough stage had a poor fit.The results illustrated that irrigation, rainfall, and other random factors had significant influences on the spatial variability of VMC.The semivariance models of VMC at different winter wheat growth stages are plotted in Figures 1-5.
Water 2018, 10, x FOR PEER REVIEW 7 of 17 All the nugget values are greater than 0, which implies that the positive effect of VMC resulted from sampling error and short distance variability.The largest range value (75.9 m) indicated that the large block scope variability was mainly distributed in the study area in the jointing stage, while the lowest range value (6.9 m) indicated the low scope variability in the milky dough stage.The measured variables had a non-zero nugget effect, which was caused by the measurement error or micro-variability that can be detected at the sampling area.To compare the effect of different nugget sizes among variables and models, the effect can be expressed as a ratio of nugget to sill, and it is possible to divide the spatial dependence pattern into three groups based on the ratio between the nugget and sill.Generally, when the ratio is less than 0.25 the spatial correlation is strong; a ratio between 0.25 and 0.75 indicates a moderate spatial correlation; and a ratio is greater than 0.75 and less than 1 indicates a weak spatial correlation.The ratio of 1 exhibits a pure nugget effect, indicating no spatial correlation at the sampling area.In establishing the spherical and exponential semivariance, our study showed that the ratios of nugget and sill of VMC in the tillering stage and harvest stage were 0.44 and 0.50, respectively, which suggested a moderate spatial autocorrelation.The ratios of nugget and sill of VMC in jointing stage, heading stage, and milky dough stage were 0.24, 0.13 and 0.001, respectively.The results suggested the strong spatial autocorrelation of soil moisture.The coefficients of determination of the models were 0.97 (jointing stage), 0.96 (harvest stage), 0.91 (heading stage), 0.82 (tillering stage), and 0.35 (milky dough stage), indicating that the models at the joint, harvest stage, heading and tillering stage had goodness of fit, while the model of the milky dough stage had a poor fit.The results illustrated that irrigation, rainfall, and other random factors had significant influences on the spatial variability of VMC.The semivariance models of VMC at different winter wheat growth stages are plotted in Figures 1-5.SGS was conducted on the basis of the semivariance models of VMC.In the process of SGS, simulation times must be sufficient because few simulations cannot represent the characteristics of the simulated variables.Also, the local change pattern needs to show in-sequence simulations using different simulation realizations.We compared 100, 300, 500, 700, and 1000 conditional simulations.We used sufficient arguments and used trial and error on numerous occasions to improve the accuracy of VMC estimation.Finally, we chose the results of 500 realizations with a 1-m interval of VMC for SGS in all the growth stages.The simulation results are shown in Figures 6-10 and these prediction results showed the spatial variability of VMC at the field scale.SGS was conducted on the basis of the semivariance models of VMC.In the process of SGS, simulation times must be sufficient because few simulations cannot represent the characteristics of the simulated variables.Also, the local change pattern needs to show in-sequence simulations using different simulation realizations.We compared 100, 300, 500, 700, and 1000 conditional simulations.We used sufficient arguments and used trial and error on numerous occasions to improve the accuracy of VMC estimation.Finally, we chose the results of 500 realizations with a 1-m interval of VMC for SGS in all the growth stages.The simulation results are shown in Figures 6-10 and these prediction results showed the spatial variability of VMC at the field scale.SGS was conducted on the basis of the semivariance models of VMC.In the process of SGS, simulation times must be sufficient because few simulations cannot represent the characteristics of the simulated variables.Also, the local change pattern needs to show in-sequence simulations using different simulation realizations.We compared 100, 300, 500, 700, and 1000 conditional simulations.We used sufficient arguments and used trial and error on numerous occasions to improve the accuracy of VMC estimation.Finally, we chose the results of 500 realizations with a 1-m interval of VMC for SGS in all the growth stages.The simulation results are shown in Figures 6-10 and these prediction results showed the spatial variability of VMC at the field scale.As mentioned above, the spatial variability pattern of VMC was clarified by SGS from the tillering stage to the harvest stage.The high value area and the low value area of VMC were presented in spatial variability map of VMC for the growth stages, which provided the prescribed reference for the variable rate irrigation and other operations.Irrigation is an increasingly important practice for sustainable agriculture.The volume of water for irrigation should be considered in the spatial variability of VMC.With it, we can decrease the water volume for irrigation in the high-VMC storage area, and increase the water volume for irrigation in the low-VMC storage area.With the support of GIS map algebra and overlay analysis, we further explored the spatial relationship between the volume of water for irrigation and the spatial variability of VMC at different wheat growth stages.With the balance of the volume of water for irrigation and VMC in topsoil, we divided the field into stable and sub-stable VMC areas.The stable area of VMC has a lower volume of water for irrigation than the sub-stable area.As is known, the jointing stage is more important than other growth stages and needs water to promote growth.Thus, it is valuable to understand the stable and sub-stable areas with respect to VMC in the jointing stage.
On basis of GIS technology, including map algebra and overlay computerization, a comprehensive map was created to analyze the spatio-temporal characteristics of VMC.As illustrated in Figure 11, denser contours of VMC reflected higher VMC areas, and the sparse contours presented lower VMC areas.The deep color (the specific color is deep green) denotes the stable area of VMC and the light color (the specific color is Reseda green) denotes the sub-stable VMC area.In optimized areas, the VMC in the topsoil is more and uniform than in the sub-stable area before irrigation in the jointing growth stage.From another perspective, less irrigation is needed in the stable area than in the sub-stable area of VMC in the winter wheat jointing stage.As mentioned above, the spatial variability pattern of VMC was clarified by SGS from the tillering stage to the harvest stage.The high value area and the low value area of VMC were presented in spatial variability map of VMC for the growth stages, which provided the prescribed reference for the variable rate irrigation and other operations.Irrigation is an increasingly important practice for sustainable agriculture.The volume of water for irrigation should be considered in the spatial variability of VMC.With it, we can decrease the water volume for irrigation in the high-VMC storage area, and increase the water volume for irrigation in the low-VMC storage area.With the support of GIS map algebra and overlay analysis, we further explored the spatial relationship between the volume of water for irrigation and the spatial variability of VMC at different wheat growth stages.With the balance of the volume of water for irrigation and VMC in topsoil, we divided the field into stable and sub-stable VMC areas.The stable area of VMC has a lower volume of water for irrigation than the sub-stable area.As is known, the jointing stage is more important than other growth stages and needs water to promote growth.Thus, it is valuable to understand the stable and sub-stable areas with respect to VMC in the jointing stage.
On basis of GIS technology, including map algebra and overlay computerization, a comprehensive map was created to analyze the spatio-temporal characteristics of VMC.As illustrated in Figure 11, denser contours of VMC reflected higher VMC areas, and the sparse contours presented lower VMC areas.The deep color (the specific color is deep green) denotes the stable area of VMC and the light color (the specific color is Reseda green) denotes the sub-stable VMC area.In optimized areas, the VMC in the topsoil is more and uniform than in the sub-stable area before irrigation in the jointing growth stage.From another perspective, less irrigation is needed in the stable area than in the sub-stable area of VMC in the winter wheat jointing stage.Recent study results showed that irrigation should be limited to once or twice for optimum winter wheat growth on the premise of full soil moisture content before sowing or planting [54,55].Qiang Yu put forward that the tillering-jointing and heading stages of winter wheat were the key growth stages for irrigation.In terms of water, irrigation was performed twice during the tilleringjoint stage and heading stage, as otherwise a decreasing yield resulted.Together with the results of the study on integration technology of water and fertilizers based on grain crop, it is practical to irrigate depending on the different growth stages of crops, with irrigation using inorganic fertilizer and formula fertilizer on crops during the growth stages [55,56].As Figure 11 shows, independent of stable or sub-stable areas, the threshold of VMC is 26 in the topsoil in the jointing stage of winter wheat, which means if VMC value is less than 26, irrigation can be used to promote the winter wheat growth.If the VMC is greater than or equal to 26, including all of the green color areas and the Reseda green areas, then the chances of highly efficient irrigation are low.Spatial maps combined with the irrigation theory and practices, we put forward that irrigation in the early jointing stage in winter wheat growth is regarded as having a high chance of being efficient and water-saving, and irrigation promotes winter wheat growth in the jointing stage.Moreover, using combined spatial analysis with 500 SGS realizations, a loose-couple spatial model was constructed between the VMC in topsoil and the volume of water for irrigation.The VMC in topsoil during the winter wheat jointing stage is a certain indicator for the correct volume of water for irrigation.

Discussion
The loose-couple model, which was constructed using 500 SGS realizations integrated with spatial analysis technology, was formed to illustrate the relationship between VMC in topsoil and the volume of water for irrigation.Once irrigation was put forward in the early jointing stage in our study, with respect to the volume of water for irrigation, the threshold of VMC in the topsoil was 26.At a value lower than 26, we should consider irrigation for the region.Recent study results showed that irrigation should be limited to once or twice for optimum winter wheat growth on the premise of full soil moisture content before sowing or planting [54,55].Qiang Yu put forward that the tillering-jointing and heading stages of winter wheat were the key growth stages for irrigation.In terms of water, irrigation was performed twice during the tillering-joint stage and heading stage, as otherwise a decreasing yield resulted.Together with the results of the study on integration technology of water and fertilizers based on grain crop, it is practical to irrigate depending on the different growth stages of crops, with irrigation using inorganic fertilizer and formula fertilizer on crops during the growth stages [55,56].As Figure 11 shows, independent of stable or sub-stable areas, the threshold of VMC is 26 in the topsoil in the jointing stage of winter wheat, which means if VMC value is less than 26, irrigation can be used to promote the winter wheat growth.If the VMC is greater than or equal to 26, including all of the green color areas and the Reseda green areas, then the chances of highly efficient irrigation are low.Spatial maps combined with the irrigation theory and practices, we put forward that irrigation in the early jointing stage in winter wheat growth is regarded as having a high chance of being efficient and water-saving, and irrigation promotes winter wheat growth in the jointing stage.Moreover, using combined spatial analysis with 500 SGS realizations, a loose-couple spatial model was constructed between the VMC in topsoil and the volume of water for irrigation.The VMC in topsoil during the winter wheat jointing stage is a certain indicator for the correct volume of water for irrigation.

Discussion
The loose-couple model, which was constructed using 500 SGS realizations integrated with spatial analysis technology, was formed to illustrate the relationship between VMC in topsoil and the volume of water for irrigation.Once irrigation was put forward in the early jointing stage in our study, with respect to the volume of water for irrigation, the threshold of VMC in the topsoil was 26.At a value lower than 26, we should consider irrigation for the region.
From spatio-temporal of VMC perspective, we explored the spatial pattern of VMC ranging from the tillering, jointing, heading, milky dough, to harvest stages.In particular, taking into account the jointing stage importance for winter growth, the spatial relationship was explored between the volume of water for irrigation and the spatial variability of VMC in the topsoil during the jointing stage.The studied field was divided into stable and sub-stable VMC areas.The areas were visualized through VMC data mapping of the field.On integrating the spatial analysis technology with 500 SGS realizations, the threshold of VMC in topsoil in the jointing stage of winter wheat was put forward.This threshold is 26, which means that if the VMC is lower than this value, irrigation should be conducted.If the VMC is greater than or equal to 26, then irrigation is not recommended as a form of efficient irrigation.Based on the threshold, the regions for the winter wheat jointing stage with the spatial variability maps of VMC were divided into stable and sub-stable areas in the case study.The maps show that less irrigation is needed in the stable area as compared to the sub-stable VMC area in the jointing growth stage of winter wheat.Furthermore, the stable and sub-stable area contour map of VMC for the winter wheat jointing stage provided good indications for an optimized irrigation strategy on the regional scale.Correspondingly, the optimal irrigation stage was the early jointing stage of winter wheat growth, and irrigation was performed once to guide the irrigation.Moreover, as a quantitative relationship was expressed between VMC in topsoil and water-saving irrigation times, a loose-couple spatial model was constructed between VMC in topsoil and the volume of water for irrigation.VMC in the winter wheat jointing stage was a certain indicator for the volume of water for irrigation, and VMC was linked with the efficient and water-saving irrigation.A new mode of precision irrigation for winter wheat was created.
The winter wheat areas have a shortage of water resources in the Yellow-Huai River Plain.There are large volumes of water and fertilizer inputs with low utilization efficiency in the region, and the irrigation and fertilization regimes have a great influence on the utilization rate of water and fertilizer.The spatial distribution of VMC in different winter wheat growth stages not only shows the high-value area and low-value areas, but also provides foundational references for highly effective irrigation and scientific decision-making with respect to irrigation times of winter wheat at different growth stages.
Many factors affect winter wheat growth, including exogenous factors (temperature, precipitation and other climatic conditions, micro-topography, environmental conditions, water and fertilizer cultivation conditions, etc.) and endogenous factors (for example, the characteristics of varieties).A loose couple-model was constructed between VMC in topsoil and volume of water for irrigation in the jointing stage of winter wheat.The relationship characteristics between them were constructed using the method.A three-dimensional model in winter wheat growth stages will be created in further study, considering the other endogenous and exogenous factors with respect to winter wheat growth.

Conclusions
In conclusion, the Moran's I of VMC variability presented a decreasing tendency through the different wheat growth stages (from tillering, jointing, heading, milky dough to harvest), showing that irrigation and precipitation have a greater effect on VMC than other factors.The spatial autocorrelation coefficient indicated the strongest clustering of VMC in the tillering stage, whereas the lowest clustering of VMC occurred in the harvest stage.This study is particularly important for the analysis of spatio-temporal patterns of VMC with the strong and consistent relationship between VMC and winter wheat growth status.Through 500 SGS realizations, results were achieved with visualization of VMC to show the areas with high and low VMC values.The spatial distribution of SGS results agreed well with the original data, overcoming the smoothing issue of the kriging interpolation method.Correspondingly, a loose-couple model was constructed between VMC in topsoil and the volume of water for irrigation in the jointing stage of winter wheat, based on 500 simulations of VMC.
Irrigation was conducted to promote the winter wheat growth in the early jointing stage, and irrigation was put forward once in this stage in our study.As to the volume of water for irrigation, the threshold of VMC in topsoil was put forward.This value was 26; at a lower value we should consider irrigation for the case region from the perspective of highly efficient irrigation.

Figure 1 .
Figure 1.Semivariogram of VMC in the tillering stage with a spherical model.Figure 1. Semivariogram of VMC in the tillering stage with a spherical model.

Figure 1 .
Figure 1.Semivariogram of VMC in the tillering stage with a spherical model.Figure 1. Semivariogram of VMC in the tillering stage with a spherical model.

Figure 2 .
Figure 2. Semivariogram of VMC in the jointing stage with a spherical model.

Figure 3 .
Figure 3. Semivariogram of VMC in the heading stage with an exponential model.

Figure 2 .
Figure 2. Semivariogram of VMC in the jointing stage with a spherical model.

Figure 2 .
Figure 2. Semivariogram of VMC in the jointing stage with a spherical model.

Figure 3 .
Figure 3. Semivariogram of VMC in the heading stage with an exponential model.Figure 3. Semivariogram of VMC in the heading stage with an exponential model.

Figure 3 .
Figure 3. Semivariogram of VMC in the heading stage with an exponential model.Figure 3. Semivariogram of VMC in the heading stage with an exponential model.

Figure 4 .
Figure 4. Semivariogram of VMC in the milky dough stage with a spherical model.

Figure 5 .
Figure 5. Semivariogram of VMC in the harvest stage with an exponential model.

Figure 4 .
Figure 4. Semivariogram of VMC in the milky dough stage with a spherical model.

Figure 4 .
Figure 4. Semivariogram of VMC in the milky dough stage with a spherical model.

Figure 5 .
Figure 5. Semivariogram of VMC in the harvest stage with an exponential model.

Figure 5 .
Figure 5. Semivariogram of VMC in the harvest stage with an exponential model.

Figure 6 .
Figure 6.Grid map of VMC from 500 sequential Gaussian simulations (SGSs) in the tillering stage.

Figure 7 .
Figure 7. Grid map of VMC from 500 SGSs in the jointing stage.

Figure 6 .
Figure 6.Grid map of VMC from 500 sequential Gaussian simulations (SGSs) in the tillering stage.

Figure 6 .
Figure 6.Grid map of VMC from 500 sequential Gaussian simulations (SGSs) in the tillering stage.

Figure 7 .
Figure 7. Grid map of VMC from 500 SGSs in the jointing stage.Figure 7. Grid map of VMC from 500 SGSs in the jointing stage.

Figure 7 .
Figure 7. Grid map of VMC from 500 SGSs in the jointing stage.Figure 7. Grid map of VMC from 500 SGSs in the jointing stage.

Figure 8 .
Figure 8. Grid map of VMC from 500 SGSs in the heading stage.

Figure 9 .
Figure 9. Grid map of VMC from 500 SGSs in the milky dough stage.

Figure 8 .
Figure 8. Grid map of VMC from 500 SGSs in the heading stage.

Figure 8 .
Figure 8. Grid map of VMC from 500 SGSs in the heading stage.

Figure 9 .
Figure 9. Grid map of VMC from 500 SGSs in the milky dough stage.Figure 9. Grid map of VMC from 500 SGSs in the milky dough stage.

Figure 9 .
Figure 9. Grid map of VMC from 500 SGSs in the milky dough stage.Figure 9. Grid map of VMC from 500 SGSs in the milky dough stage.

Figure 10 .
Figure 10.Grid map of VMC from 500 SGSs in the harvest stage.

Figure 10 .
Figure 10.Grid map of VMC from 500 SGSs in the harvest stage.

Figure 11 .
Figure 11.Overlay map from 500 simulation contours of VMC with SGS in the jointing stage with a stable area grid of VMC.

Figure 11 .
Figure 11.Overlay map from 500 simulation contours of VMC with SGS in the jointing stage with a stable area grid of VMC.

Table 1 .
Parameters of the semivariance model of volumetric moisture content (VMC) during the winter wheat growth stages.

Table 1 .
Parameters of the semivariance model of volumetric moisture content (VMC) during the winter wheat growth stages.