Estimating Soil Water Content and Evapotranspiration of Winter Wheat under Deficit Irrigation Based on SWAP Model

Deficit irrigation strategy is essential for sustainable agricultural development in arid regions. A two−year deficit irrigation field experiment was conducted to study the water dynamics of winter wheat under deficit irrigation in Guanzhong Plain in Northwest China. Three irrigation levels were implemented during four growth stages of winter wheat: 100%, 80% and 60% of actual evapotranspiration (ET) measured by the lysimeter with sufficient irrigation treatment (CK). The agro−hydrological model soil−water−atmosphere−plant (SWAP) was used to simulate the components of the farmland water budget. Sensitivity analysis for parameters of SWAP indicated that the saturated water content and water content shape factor n were more sensitive than the other parameters. The verification results showed that the SWAP model accurately simulated soil water content (average relative error (MRE) < 21.66%, root mean square error (RMSE) < 0.07 cm3 cm−3) and ET (R2 = 0.975, p < 0.01). Irrigation had an important impact on actual plant transpiration, but the actual soil evaporation had little change among different treatments. The average deep percolation was 14.54 mm and positively correlated with the total irrigation amount. The model established using path analysis and regression methods for estimating ET performed well (R2 = 0.727, p < 0.01). This study provided effective guidance for SWAP model parameter calibration and a convenient way to accurately estimate ET with fewer variables.


Introduction
The Guanzhong Plain is rich in fertile land and abundant in light and heat resources. The plain is an important winter wheat (Triticum aestivum L.) producing area in Northwest China [1]. The total planting area and yield of winter wheat accounts for 80% and 85% of Shaanxi Province, respectively [2]. Winter wheat production has great significance for ensuring food security, social stability and agricultural economic development [3]. Restricted water for irrigation is an important factor limiting the growth and production of winter wheat [1,[4][5][6]. In recent years, the rapid development of urbanization and industrialization has caused serious water pollution and a reduction of freshwater resources for irrigation [7][8][9]. Due to climate change and the uneven spatial and temporal distribution of water resources in the region, the drought risk in Northwest China has increased [10]. In addition, unrestrained planting and irrigation have caused a serious waste of water resources. In moderate

Experimental Area
The experiment was carried out from October 2012 to June 2014 at the Institute of Water-saving Agriculture in Arid Areas of China (108°04′ E, 34°17′ N, 521 m above sea level), Northwest A&F University, Yangling, Shaanxi Province, China. The area has a continental temperate monsoon climate, and the multi-year average temperature and precipitation are 12.5 °C and 609 mm, respectively. The average temperature and precipitation of the two growing seasons are 8.42 °C, 220.10 mm (2012-2013 growing season) and 8.35 °C, 290.30 mm (2013-2014 growing season), as shown in Figure 1. The local soil texture is silty clay loam [21]. The average field water capacity of the 1 m soil layer is 0.32 cm 3 cm −3 , and the average bulk density of the soil was 1.39 g cm −3 ( Table 1). The groundwater is below 50 m, and groundwater recharge can be ignored. The area of each plot is 6.6 m 2 (length: 3 m, width: 2.2 m). A concrete side wall was utilized to prevent water side seepage. The bottom of the plot was provided with a filter layer for free drainage. A mobile electric shelter was installed above the test area, which was closed (manually) during precipitation and opened when precipitation did not occur ( Figure 2).

Experimental Design
The winter wheat variety was "Xiaoyan 22" and the planting density was 357 seeds m −2 . The sowing times of the two−year experiment were 18 October 2012 and October 15, 2013. The harvest times were 3 June 2013 and 8 June 2014. Each treatment received N 244 kg ha −1 and P2O5 270 kg ha −1 at the sowing stage. The irrigation method was surface irrigation. The entire growth period of the winter wheat was divided into four stages: emergence-jointing stage, jointing-heading stage, heading-grouting stage and grouting-mature stage. Three irrigation levels were implemented during each stage: 100% (no water deficit), 80% (moderate deficit) and 60% (severe deficit) of actual ET measured by the large weighing lysimeter with sufficient irrigation treatment (CK, namely W1 and T1) ( Figure 3). The range of the device was 0 to 6 t, and the measurement accuracy was less than 150 g [35]. The size of the lysimeter is consistent with the plot. A 4-factor 3-level partial orthogonal test design was used with nine treatments that were replicated three times ( Table 2).

Experimental Design
The winter wheat variety was "Xiaoyan 22" and the planting density was 357 seeds m −2 . The sowing times of the two−year experiment were 18 October 2012 and October 15, 2013. The harvest times were 3 June 2013 and 8 June 2014. Each treatment received N 244 kg ha −1 and P 2 O 5 270 kg ha −1 at the sowing stage. The irrigation method was surface irrigation. The entire growth period of the winter wheat was divided into four stages: emergence-jointing stage, jointing-heading stage, heading-grouting stage and grouting-mature stage. Three irrigation levels were implemented during each stage: 100% (no water deficit), 80% (moderate deficit) and 60% (severe deficit) of actual ET measured by the large weighing lysimeter with sufficient irrigation treatment (CK, namely W1 and T1) ( Figure 3). The range of the device was 0 to 6 t, and the measurement accuracy was less than 150 g [35]. The size of the lysimeter is consistent with the plot. A 4-factor 3-level partial orthogonal test design was used with nine treatments that were replicated three times (Table 2). xperimental Design he winter wheat variety was "Xiaoyan 22" and the planting density was 357 seeds m −2 g times of the two−year experiment were 18 October 2012 and October 15, 2013. The ha were 3 June 2013 and 8 June 2014. Each treatment received N 244 kg ha −1 and P2O5 270 kg sowing stage. The irrigation method was surface irrigation. The entire growth period o r wheat was divided into four stages: emergence-jointing stage, jointing-heading s ng-grouting stage and grouting-mature stage. Three irrigation levels were impleme g each stage: 100% (no water deficit), 80% (moderate deficit) and 60% (severe deficit) of a easured by the large weighing lysimeter with sufficient irrigation treatment (CK, namely 1) (Figure 3). The range of the device was 0 to 6 t, and the measurement accuracy was less [35]. The size of the lysimeter is consistent with the plot. A 4-factor 3-level partial orthog esign was used with nine treatments that were replicated three times (Table 2).    T2  100  80  80  0  T3  100  60  60  0  T4  80  100  60  0  T5  80  80  100  0  T6  80  60  80  0  T7  60  100  80  0  T8  60  80  60  0  T9  60

Meteorological data
Daily weather data were obtained from the national general weather station near the plots (The distance between them is approximately 0.3 km), including maximum temperature (T max , • C), minimum temperature (T min , • C), average relative humidity (RH, %), wind speed (u 10 , m s −1 , the wind speed sensor was 10.2 m above the ground), sunshine duration (S d , h) and precipitation (P, mm).

2.
Soil data Soil profile data was obtained from field measurements. The detailed physical properties of each soil layer are shown in Table 1. The soil water content (SWC) of the 1 m soil layer was measured by the method of weighing after drying. Every 10 cm soil depth was a measurement point, and the specific sampling time is shown in Figure 4.

3.
Crop ET The crop ET of each treatment was calculated using the water balance method: where I is irrigation (mm); D is deep percolation (mm); ∆W is the change of soil water storage in the root zone (mm); Z r is root zone depth (m); θ(t 1 ) and θ(t 2 ) are the average soil volume water content of the root zone at t 1 and t 2 respectively.

Crop ET
The crop ET of each treatment was calculated using the water balance method: where I is irrigation (mm); D is deep percolation (mm); ∆W is the change of soil water storage in the root zone (mm); Zr is root zone depth (m); θ(t1) and θ(t2) are the average soil volume water content of the root zone at t1 and t2 respectively.

Crop height (CH)
Ten wheat plants were selected in each plot, and CH was measured using a tape measure with an accuracy of 1 mm every week.

Leaf area index (LAI)
The SunScan-SS1 canopy analyzer (Delta-T Devices Ltd., Cambridge, UK) was used at each growth stage of winter wheat to measure LAI. The SWAP model uses one-dimensional Richards equations to describe the water flow and has a strong physical foundation [30,36]:

SWAP Model Introduction and Sensitivity Analysis
where C(h) is the differential water capacity (cm −1 ), indicating the change in water content caused by the change of the unit matrix potential, and is numerically equal to the slope of the soil water

Crop height (CH)
Ten wheat plants were selected in each plot, and CH was measured using a tape measure with an accuracy of 1 mm every week.

5.
Leaf area index (LAI) The SunScan-SS1 canopy analyzer (Delta-T Devices Ltd., Cambridge, UK) was used at each growth stage of winter wheat to measure LAI. Water flow The SWAP model uses one-dimensional Richards equations to describe the water flow and has a strong physical foundation [30,36]: where C(h) is the differential water capacity (cm −1 ), indicating the change in water content caused by the change of the unit matrix potential, and is numerically equal to the slope of the soil water characteristic; t is time (d); K is the hydraulic conductivity (cm d −1 ); h is the soil water pressure head (cm); z is the vertical coordinate (cm) taken positively upward; and S(z) is the root water uptake at a certain depth (cm 3 cm −3 d −1 ).

Evapotranspiration
SWAP offers two methods for the distribution of plant transpiration and soil evaporation. One is based on crop and soil factors and the other is a direct application of the Penman-Monteith equation. The direct application of the Penman-Monteith equation was used in this study. This method does not require crop or soil factors to translate reference ET to crop ET and avoids the complicated calibration process [28]. The functions of potential plant transpiration (T p ) and soil evaporation (E p ) are as follows [21,28]: where W frac is the fraction of a day that the canopy is wet (-), since the shelter was closed during precipitation events, the rainfall interception was not considered in this study, and W frac was set to 0 [28]; V c is the vegetation cover (-) (function 6); ∆ v is the slope of the water vapor pressure curve (kPa • C −1 ); λ w is the latent heat of evaporation (J kg −1 ); R n is the net radiant flux at the canopy surface (J m −2 d −1 ); G is the heat flux of the soil (J m −2 d −1 ), which is negligible on a daily basis; p 1 is the unit conversion factor, 86400 (s d −1 ); ρ a is the air density (kg m −3 ); C a is heat capacity of moist air (J kg −1 • C −1 ); e s is the saturation vapor pressure (kPa); e a is the actual vapor pressure (kPa); r a c is the aerodynamic resistance of the canopy surface corrected according to V c ; γ is the wet and dry thermometer constant (kPa • C −1 ); r c, min is the minimum canopy resistance (s m −1 ), with a value of 20 selected according to Feng et al. [37]; r a s is the aerodynamic resistance of the soil surface corrected according to V c ; r s, min is the minimum soil resistance, i.e., the wet soil resistance (s m −1 ), and a value of zero was selected according to the model value range.
where k dir and k dif are the direct and diffuse optical extinction coefficients (-), respectively, taking the model default values of 0.60 and 0.75 [38]. The potential root water uptake (S p ) is equal to T p [28], and the actual root water uptake (S a ) should consider the effects of water, salt and low temperature stress: S a (z) = α rw α rd α rs α r f S p (z) (7) where α rw , α rd , α rs , and α rf (-) are the wet, dry, salt and low temperature stress coefficients, respectively. This test used local groundwater for surface irrigation, regardless of salt and low temperature stress. The water stress in SWAP is a piecewise linear function proposed by Feddes et al. 1978 [39]. The parameters of the water stress function are based on previous research [37,40,41]. The actual plant transpiration (T a ) is obtained by integrating S a . The actual soil evaporation (E a ) is the smaller value between E p and maximum soil evaporation [28,29]: where E max is maximum soil evaporation (cm d −1 ); K 1/2 is the average hydraulic conductivity (cm d −1 ) between the soil surface and the first node; h atm is the soil water pressure (cm) balanced with the relative humidity of the atmosphere; h 1 is the first soil water pressure (cm) of the node; and z 1 is the soil depth (cm) of the first node.

Model Input
The data required by the meteorological module were obtained from the weather station near the test area. The soil module was needed to discretize the vertical soil profile. In our study, the 1 m soil depth was divided into five layers. Based on the recommended value of the model and previous research, each soil layer was further divided into compartments of different thicknesses [28]. The crop module needed to input the growth stage to establish a function of LAI, CH and root length with the development stage (DVS). Due to the lack of root length data, the longest root length of winter wheat was set to 1 m following Jha et al. [42]. The root length of all other stages was obtained based on interpolation using CH [37]. There are two types of simulations for the crop growth process: a simple module and a detailed module [28]. This paper focuses on field water movement and therefore selects the simple crop module. The upper boundary of the simulation was located above the crop canopy, and the lower boundary was set to 1 m below the soil surface. The bottom drained freely, and lateral percolation and groundwater recharge were not considered. Other model theory and calculation methods can refer to Kroes et al. [28].

Sensitivity Analysis
Sensitivity analysis is a prerequisite for model calibration [33]. The calibration of the simple module in the SWAP model is the process of continuously adjusting the van Genuchten-Mualem model (VGM) parameters to minimize the deviation between the simulation SWC and the observed values [31,35,37]. Therefore, we mainly analyzed the sensitivity of VGM parameters. This study used the parameter change method for sensitivity analysis [43,44]. The sensitivity index (SI) was used to evaluate the influence of parameters on the model results [45,46]. The specific steps were as follows:

3.
Determine the initial value of the VGM parameters The measured soil physical properties (Table 1) were fitted to the VGM equations with the RETC software [47]. The fitted values were considered as the initial soil parameters in model calibration [31,48].

4.
Determine the change range The range of local VGM parameters was measured or simulated according to previous studies [49,50]. Based on the initial values obtained by the first step above and the local range, changes of ±10%, ±20%, ±30%, and ±50% were set as model inputs.

5.
Calculate the SI The SI was calculated by the following function [43,46]: where Y 0 is the model output (SWC or ET in this study) under the condition that the parameters are the initial value and Y i represents the output after changing the parameters.

Sensitivity ranking consistency test
The consistency of the sensitivity ranking was evaluated by the TDCC coefficient (top-down coefficient of concordance) [32,51,52]. The specific calculation process was as follows: First, the absolute values of the SI were sorted to get the following matrix: where O ij is the outcome (i.e., the absolute SI value) for the i th parameter in the j th sensitivity analysis test; r(O ij ) is the order of the i th parameter in the j th sensitivity analysis test (the absolute SI value ranges from high to low); Then, the Savage score matrix was calculated [53], and the following matrix was obtained: where ss(O ij ) is the Savage score of the i th parameter in the j th sensitivity analysis test; and n is the number of model parameters (the value is 5 in this study, which represents the five soil VGM parameters: θ r , θ s , K s , α and n).
Finally, the TDCC coefficient is: where m is the number of sensitivity analysis tests (the value is 8 in this study, which represents eight different changes of ±10%, ±20%, ±30%, and ±50%). The significant P value of the TDCC coefficient is calculated by the statistic T: where the statistic T in the function approximates the χ 2 distribution with a degree of freedom of n-1. Generally, the p value is less than 0.05 at a larger TDCC value. This means that the parameter sensitivity ranking of each sensitivity analysis test has significant consistency [32,52,54].

Model Calibration and Verification
Based on the sensitivity analysis results, we adopt different parameter adjustment strategies to calibrate the SWAP model. This study selected the W4, W5 and T4, T5 treatments for model calibration because their irrigation amounts were at a medium level compared with the other treatments. We avoided using extreme treatments (irrigated too much or too little) for model calibration to ensure that the calibration results also apply to general irrigation. The degree of agreement between the simulated and measured SWC was evaluated by the average relative error (MRE) and the root mean square error (RMSE): where n is the number of measured values; P i is the i th predicted value; and O i is the i th observed value. We adjusted the corresponding VGM parameters in SWAP until MRE ≤ 20% and RMSE ≤ 0.05 cm 3 cm −3 [30,35,37,55]. The adjusted parameters were used for model verification for all other treatments. The calibrated VGM parameters are shown in Table 3.

Path Analysis and Establishment of the ET Estimation Model
Path analysis is a multivariate statistical analysis technique that studies the relationship between independent variables and dependent variables [56]. The reliability of this method has been verified in many studies [57][58][59][60]. This method can help to analyze the direct and indirect effects of the independent variable on the dependent variable, and get the decision coefficient (DC) of each factor: where DC(i) is the decision coefficient of the independent variable, i; p i is the direct path coefficient of the independent variable, i; and r iy is the correlation coefficient between the independent variable i and the dependent variable y. DC(i) > 0 indicates that the independent variable promotes the dependent variable, and DC(i) < 0 indicates that the independent variable has an inhibitory effect on the dependent variable.
Based on ET simulated by the SWAP model, correlation analysis and path analysis were carried out with meteorological factors (radiation, wind speed at 2 m, RH and average temperature), crop factors (CH and LAI) and soil factors (SWC of different layers). We selected the most relevant factors according to the ranking of DC, and established the ET estimation model under deficit irrigation conditions by multiple nonlinear regression method.

Data Processing and Statistical Analysis
Path analysis, significance analysis, and regression model construction were performed using SPSS 21 (IBM, Armonk, NY, USA).

Overall Parameter Sensitivity
Taking the ±30% change of the VGM parameters as an example, the sensitivity of SWC and ET to the parameters was explained. In terms of total sensitivity, SWC was more sensitive to the VGM parameters than ET. The average SI of the VGM parameter for SWC was 0.045 and −0.017 for ET.
When the parameter values were changed by +30%, θ r and θ s of the 0-20 cm layer were positively correlated with SWC. The SI values of θ r and θ s were 0.058 and 0.349, respectively (Figure 5a). K s , α and n were negatively correlated with SWC, and their SI were −0.025, −0.033 and −0.262, respectively. The absolute SI values of θ s and n were much larger than those of θ r , α and K s . For other soil profiles, the correlation between the VGM parameters and SWC was consistent with the 0-20 cm soil layer (Figure 5b-e). For the 20-100 cm soil layer, the average absolute SI value of each parameter was θ s (0.284) > n (0.203) > θ r (0.049) > K s (0.028) > α (0.005). When the parameter values changed by -30%, the positive and negative correlation between the VGM parameters and the SWC of the root zone (0-100 cm soil layer) was unchanged. The absolute SI values were slightly different from the +30% change: n (0.633) > θ s (0.284) > θ r (0.050) > K s (0.037) > α (0.015). However, the absolute SI values of the first two parameters (θ s and n) were still much larger than those of three other parameters.
Sustainability 2020, 12, x FOR PEER REVIEW 12 of 32 consistent with results with the parameter value change of +30%. The two most sensitive parameters during this growth stage were also θs (-0.141) and n (-0.144), but the corresponding absolute SI values were generally higher than those observed in response to the change of +30% (0.072 and -0.080, respectively).  For ET (Figure 6a-d), when the parameter values were changed by +30%, the average SI values of the emergence-heading stages were negative (θ r , K s , α and n were −0.014, −0.006, −0.011 and −0.061, respectively) except for θ s (0.084). θ s and α were positively correlated with ET during the heading-maturity stages, and the other parameters were negatively correlated. The absolute SI values of θ s and n during the heading-grouting stage (0.012 and 0.014) were much lower than in the other growth stages. When the parameter values were changed by -30%, θ s and n were positively correlated with ET at the emergence-heading stage. The other parameters were negatively correlated. θ r and K s were negatively correlated with ET at the heading-maturity stage, and all other parameters were positively correlated with ET. The ranking of the absolute SI values in the emergence-jointing stage was consistent with results with the parameter value change of +30%. The two most sensitive parameters during this growth stage were also θ s (−0.141) and n (−0.144), but the corresponding absolute SI values were generally higher than those observed in response to the change of +30% (0.072 and −0.080, respectively).

Influence of Different Change Ranges on Sensitivity Analysis
Taking the high sensitivity parameters θs and n affecting SWC and ET as examples, the influences of different change ranges on the sensitivity analysis were explained.
For SWC in the root zone, the average value and standard deviation of the SI for the two parameters increased with the parameter change range. When the range of change increased from ±10% to ±50%, the SI of θs increased by 441% (from +10% to +50%) and 396% (from -10% to -50%), respectively. The absolute SI values of n increased by 203% (from +10% to +50%) and 624% (from -10% to -50%), respectively.
For ET throughout the growth period, the average value and the standard deviation of the SI of θs increased with the parameter change range. When the range of change increased from ±10% to ±50%, the SI of θs increased by 306% (from +10% to +50%) and 946% (from -10% to -50%), respectively. The SI of n varied with the positive and negative parameter change range during the emergence-heading stage. When the positive change gradually increased from +10% to +50%, the absolute SI values increased gradually. When the negative change gradually increased from -10% to -20%, the SI of n gradually decreased from a positive value (0.016 and 0.018, respectively) to -0.144 and -0.083 (when the change was -30%). As the negative change increased, the SI of n further decreased to -0.916 and -0.817 (when the change was -50%). When the negative change gradually increased from -10% to -50%, the SI of n in the heading-maturity stage was always negative, and the absolute SI values increased gradually.
The Savage scores of the VGM parameter were calculated using the sensitivity ranking results under different change ranges in different soil layers (Table 4). Furthermore, we obtained the TDCC coefficient to evaluate the consistency of sensitivity ranking. The results of the sensitivity consistency test for the five parameters under the eight changed conditions are shown in Table 5. The TDCC

Influence of Different Change Ranges on Sensitivity Analysis
Taking the high sensitivity parameters θ s and n affecting SWC and ET as examples, the influences of different change ranges on the sensitivity analysis were explained.
For SWC in the root zone, the average value and standard deviation of the SI for the two parameters increased with the parameter change range. When the range of change increased from ±10% to ±50%, the SI of θ s increased by 441% (from +10% to +50%) and 396% (from -10% to -50%), respectively. The absolute SI values of n increased by 203% (from +10% to +50%) and 624% (from -10% to -50%), respectively.
For ET throughout the growth period, the average value and the standard deviation of the SI of θ s increased with the parameter change range. When the range of change increased from ±10% to ±50%, the SI of θ s increased by 306% (from +10% to +50%) and 946% (from -10% to -50%), respectively. The SI of n varied with the positive and negative parameter change range during the emergence-heading stage. When the positive change gradually increased from +10% to +50%, the absolute SI values increased gradually. When the negative change gradually increased from -10% to -20%, the SI of n gradually decreased from a positive value (0.016 and 0.018, respectively) to −0.144 and −0.083 (when the change was -30%). As the negative change increased, the SI of n further decreased to −0.916 and −0.817 (when the change was -50%). When the negative change gradually increased from -10% to -50%, the SI of n in the heading-maturity stage was always negative, and the absolute SI values increased gradually.
The Savage scores of the VGM parameter were calculated using the sensitivity ranking results under different change ranges in different soil layers (Table 4). Furthermore, we obtained the TDCC coefficient to evaluate the consistency of sensitivity ranking. The results of the sensitivity consistency test for the five parameters under the eight changed conditions are shown in Table 5. The TDCC coefficients of different soil layers exceeded 0.778 (with an average of 0.789) and the P values were less than 0.01. This indicated that the ranking of parameter sensitivity for SWC had extremely significant (p < 0.01) consistency under different ranges in the root zone. The Savage scores of the VGM parameter were calculated by the sensitivity ranking results under different change ranges in different growth stages ( Table 6). The TDCC coefficients of different growth stages (Table 5) exceeded 0.622 (the average was 0.695) and the P values were all less than 0.01. This indicated that the ranking of parameter sensitivity for ET had extremely significant (p < 0.01) consistency under different ranges throughout the growth period.  Note: TDCC (top-down coefficient of concordance) is the consistency test coefficient; T value is the statistic calculated according to TDCC and approximates the χ 2 distribution with a degree of freedom of n-1 (n is the number of parameters, the value is 5 in this study, which represents 5 VGM parameters: residual water content θ r , saturated water content θ s , saturated hydraulic conductivity K s , and water content shape factor α and n); p value is the result of the saliency analysis. When p < 0.01, the parameter sensitivity ranking results have extremely significant consistency.

SWAP Model Calibration and Verification
The SWC simulated by the SWAP model were in good agreement with the measured values ( Figure 7). The model could clearly characterize the change of SWC caused by irrigation. The MREs of each treatment were 4.71%-21.66%, and the RMSEs were 0.01−0.07 cm 3 cm −3 (Table 7). These values were all within an acceptable range [37]. This meant that the verified SWAP model can be used to accurately predict SWC.

Pattern of Soil Water
During the period of no irrigation, the water flux of the surface soil was mainly positive due to soil evaporation ( Figure 8). As the experiment progressed, the positive water flux gradually extended to the deep soil. The soil water flux was negative under irrigation, and the greater the irrigation amount, the larger the negative water flux in the lower soil layers. At the heading-grouting stage, the irrigation amount of the T1 treatment was 40% higher than that of T8. The average water flux of T1 in the root zone was negative. The T8 treatment only maintained the water flux of −0.035 mm d −1 in the 60 cm soil layer; however, the flux in the deeper soil had become positive. In addition, the experiment excluded the influences of precipitation, making irrigation and soil storage water the two main sources of water consumption in the winter wheat field. Since the deficit irrigation did not meet the demand of crop ET, all treatments needed to consume 15.99-35.58% of water stored in the soil ( Table 8). The deep percolations of all treatments were positive (the average is 14.54 mm), and the total deep percolation generally increased with irrigation. However, some treatments did not conform to this pattern. For example, irrigation of W3 was 24 mm less than W9, and the deep percolation of W3 was 16.72 mm higher than the latter. This was also observed in many treatments during the 2013-2014 growing season.

Pattern of ET
The simulated ET in the entire growth period was 257.05-389.63 mm, which had good agreement (R 2 = 0.975, p < 0.01) with the calculated value ( Figure 9). In general, ET was positively correlated with the amount of irrigation in the same year. The difference in ET between two growing seasons was large. For example, the irrigation amount of T5 was 10.42% lower than that of W9, but the ET was 7.97 mm higher. The daily ET of winter wheat increased first and then decreased with growth ( Figure 10   T a of the two−year test accounted for 88.40% of the total water consumption, which was the main consumption for the entire growth stage of winter wheat. T a increased gradually with the growth period, while E a was the opposite (Table 9). Irrigation had a great influence on T a but limited effects on E a . For instance, the difference in T a between W3 and W8 reached 30.47% at the grouting-mature stage. However, E a was stable at 6.96 mm. The effects of water stress on ET of the rewetting stage were different. For example, W1 and W7 received the same irrigation amount at the jointing-heading stage, and the difference in ET was only 0.68 mm. The irrigation amount of T5 at the jointing-heading stage was 20% lower than that of T7. The irrigation in T5 at the heading-grouting stage was 20% higher than that in T7, but the ET at the heading-grouting stage was 6.43% lower than that in T7. The situation was similar for W4 and W8, and W7 and W9. In the case of water stress during the heading-grouting stage and rewetting during the grouting-mature period, the above phenomenon still existed but was weakened. For example, irrigation of W9 during the heading-grouting stage was 20% lower than that of W6. W9 was irrigated more during the grouting-mature stage, while ET of this stage was 1.89% lower than that of W6. W5 and W6, W7 and W8 showed similar results. Note: the unit of T a and E a is mm.

Correlation and Path Analysis of ET
In order to screen out the main influential factors of ET, the meteorological factors: radiation (10_RAD), wind speed at 2 m (X11_U2), relative humidity (X9_RH), average temperature (X8_T); crop factors: crop height (X2_CH) and leaf area index (X1_LAI); soil factors: average volume SWC of each soil layer (X3_SWC20, X4_SWC40, X5_SWC60, X6_SWC80, X7_SWC100) were taken as independent variables. The simulated daily ET was taken as a dependent variable (Y_ET). Subsequently, a stepwise linear regression was performed. The multivariate linear regression estimated the daily ET well (R 2 = 0.815, p < 0.01) after removing the X4_SWC40 and X7_SWC100. However, the estimation model requires too many parameters. By running the path analysis, the number of parameters needed for estimation model can be reduced which would lead to model simplification. The results of the correlation (Table 10) and path analysis (Table 11) showed that, except for X9_RH and X6_SWC80, the direct effect coefficients (DE) of other factors for daily ET were positive. The absolute DE of X2_CH, X10_RAD, X5_SWC60, X9_RH and X6_SWC80 for daily ET were all greater than 0.100. They had large direct influences on Y_ET, and the ranking of influence was: X2_CH (0.545) > X5_SWC60 (0.420) > X10_RAD (0.340) > X6_SWC80 (−0.185) > X9_RH (−0.117). From the perspective of the total indirect effect coefficients (IE), except for X5_SWC60 and X9_RH, the IE of other factors on daily ET was positive. The simple correlation coefficient (CC) with the dependent variable was equal to the sum of the DE and the total IE. Except for X5_SWC60 and X6_SWC80, the influence of other factors on the Y_ET was enhanced by the IE (compared with the absolute DE, the absolute CC increased). The strong positive DE of X5_SWC60 on ET was mainly obscured by the negative IE of X2_CH (−0.213) and X6_SWC80 (−0.181). This resulted in a reduced CC between X5_SWC60 and Y_ET. The CC between X1_LAI, X8_T and Y_ET was enhanced by the IE of X2_CH (0.448 and 0.431, respectively) and X10_RAD (0.146 and 0.235, respectively). The total DC reflects the comprehensive effect of the independent variable on the dependent variable (Table 12). The DC of the error term was small (0.029), indicating that the variables affecting Y_ET were considered comprehensively in the path analysis. The absolute DC of meteorological factors was X10_RAD (0.359) > X8_T (0.089) > X9_RH (0.043) > X11_U2 (0.032); the crop factors was X2_CH (0.550) > X1_LAI (0.084); and the soil factors was X5_SWC60 (−0.235) > X6_SWC80 (0.029) > X3_SWC20 (0.021).

Establishment and Verification of Crop and Water Stress Coefficient Estimation Model
According to the DC ranking mentioned above, the first factor was selected from the crop and soil factors, namely, X2_CH and X5_SWC60 were used to construct the respective crop and water stress coefficient estimation models.
The reference crop ET (ET 0 ) for the 2012-2013 and 2013-2014 growing seasons was calculated using the Penman-Monteith method [28]. ET of the contrasting treatment was divided by ET 0 to obtain the crop coefficient (K c ). We obtained the relationship between CH and K c through curve estimation (R 2 = 0.425, p < 0.01): K c = 0.250CH 0.311 (19) where K c is the crop coefficient (-) and CH is the crop height of winter wheat (cm). The ET of deficit irrigation was divided by ET of sufficient irrigation (W1 and T1) to obtain a water stress coefficient (K θ ). In order to avoid the influence of soil texture, and increase the applicability of the fitting equation, we took the relative water content (the ratio of SWC to field water capacity, RWC) of the 0-60cm soil layer as the independent variable and K θ as the dependent variable, and the following regression model was obtained by curve estimation analysis of all deficit irrigation treatments (R 2 = 0.357, p < 0.01): K θ = −0.983RWC 2 +2.44RWC (20) where K θ is the water stress coefficient (-) and RWC is the relative water content of the 0-60 cm soil layer (cm 3 cm −3 ). Based on Equations (12) and (13), the following ET estimation model was constructed: where a, b, c and d are empirical coefficients. Note: r ij represents the correlation coefficient between the i th variable and the j th variable; Dependent variable: Y_ET; Independent variables: X2_CH represents crop height (cm); X10_RAD represents radiation (MJ m −2 ); X11_U2 represents wind speed at a height of 2 m (m s −1 ); X1_LAI represents leaf area index (−); X9_RH represents relative humidity (%); X6_SWC80 represents the average water content of the 0-80 cm soil layer (cm 3 cm −3 ); X5_SWC60 represents the average water content of the 0-60 cm soil layer (cm 3 cm −3 ); X8_T represents the average temperature ( • C); X3_SWC20 represents the average water content of the 0-20 cm soil layer (cm 3 cm −3 ); ** represents significant difference at p < 0.01; * represents significant difference at p < 0.05. The footnotes of this table also apply to Tables 11 and 12.  Note: DC represents the decision coefficient, which is an evaluation index of the comprehensive effect of independent variables on dependent variables. Total DC represents the sum of the decision coefficients of different independent variables on the dependent variable; e represents the error term in the path analysis, and the smaller the decision coefficient, the more comprehensive the factors affecting evapotranspiration.

Parameter Estimation and Verification of the ET Estimation Model
The coefficients of the independent variables in Equations (12) and (13) were taken as initial values, and the parameters of Equation (14) were estimated by multivariate nonlinear regression analysis of SPSS. After several iterations, the obtained parameter values were a = 0.791, b = −0.002, c = −2.463, d = 3.610 (R 2 = 0.996). The parameter estimation results were substituted into Equation (14) to calculate the daily ET. Figure 11 compares the calculated values and outputs of the SWAP model, showing good agreement (R 2 = 0.727, p < 0.01). This estimation model provides an effective method for estimating actual ET under deficit irrigation through a few easy-to-measure indicators.

Discussion
From the sensitivity analysis of VGM parameters, SWC and ET were more sensitive to θs and n, while θr, Ks and α were relatively less sensitive. This is consistent with the results of Sun [45] and slightly different from those of Singh et al. [61]. This may be because n is the exponential term of the VGM function [28], and θs is the main influential factor of the exponential product term. Thus, small changes of θs and n have a great impact on SWC. The standard deviation of SI for θs and n show an increasing trend with an increase in the range of the parameter changes. Therefore, we should control the parameter adjustment within a reasonable range based on the results of previous studies [49,50]. This avoids large deviations between the model output and the observed values. In particular, when the range of parameter changes reaches -50%, n is already close to the minimum threshold of 1.001. Therefore, it is not recommended to take n around this value. According to previous research results, θr, Ks and α have a large range of reference values. Therefore, θs and n should be adjusted preferentially when the model parameters are calibrated, and then θr, Ks and α are adjusted appropriately based on the initial range of values. This can effectively shorten the time of model calibration and also further improve the simulation accuracy.
The simulated SWC in this study was generally accurate. However, the simulation accuracy of deep SWC was better than that of shallow soil. The average MRE of 0-60cm and 60-100cm soil layers were 13.94% and 11.81%, respectively. This may be related to the variability of water content in different soil layers. The SWC of the 0-60 cm soil layer varied greatly, which was consistent with Ma et al. [31], who believed that the combined effect of irrigation and ET was more likely to affect the water content of the upper soil. In addition, the surface soil was in close contact with the atmosphere. According to research results of Wang et al. [23], soil evaporation was greatly affected by temperature, radiation and ground cover, which was likely to cause changes in SWC. Furthermore, the root activity area of winter wheat was concentrated in the 0-40 cm soil layer [62]. After rewetting, the root activity was stimulated, and the water absorption intensity was enhanced [55,63,64], resulting in a reduction in SWC. However, the 60-100 cm soil layer had a small change throughout

Discussion
From the sensitivity analysis of VGM parameters, SWC and ET were more sensitive to θ s and n, while θ r , K s and α were relatively less sensitive. This is consistent with the results of Sun [45] and slightly different from those of Singh et al. [61]. This may be because n is the exponential term of the VGM function [28], and θ s is the main influential factor of the exponential product term. Thus, small changes of θ s and n have a great impact on SWC. The standard deviation of SI for θ s and n show an increasing trend with an increase in the range of the parameter changes. Therefore, we should control the parameter adjustment within a reasonable range based on the results of previous studies [49,50]. This avoids large deviations between the model output and the observed values. In particular, when the range of parameter changes reaches -50%, n is already close to the minimum threshold of 1.001. Therefore, it is not recommended to take n around this value. According to previous research results, θ r , K s and α have a large range of reference values. Therefore, θ s and n should be adjusted preferentially when the model parameters are calibrated, and then θ r , K s and α are adjusted appropriately based on the initial range of values. This can effectively shorten the time of model calibration and also further improve the simulation accuracy.
The simulated SWC in this study was generally accurate. However, the simulation accuracy of deep SWC was better than that of shallow soil. The average MRE of 0-60cm and 60-100cm soil layers were 13.94% and 11.81%, respectively. This may be related to the variability of water content in different soil layers. The SWC of the 0-60 cm soil layer varied greatly, which was consistent with Ma et al. [31], who believed that the combined effect of irrigation and ET was more likely to affect the water content of the upper soil. In addition, the surface soil was in close contact with the atmosphere. According to research results of Wang et al. [23], soil evaporation was greatly affected by temperature, radiation and ground cover, which was likely to cause changes in SWC. Furthermore, the root activity area of winter wheat was concentrated in the 0-40 cm soil layer [62]. After rewetting, the root activity was stimulated, and the water absorption intensity was enhanced [55,63,64], resulting in a reduction in SWC. However, the 60-100 cm soil layer had a small change throughout the growth stage. The irrigation amount of the T2 and T7 treatments was similar, but the deep percolation of the former was much higher than that of T7. The reason is that the distribution of irrigation was different during the growth stage. T2 was fully irrigated at the emergence-jointing stage. However, the root system was still underdeveloped during the early stage of growth. Most of the irrigation water was not utilized by plants but directly infiltrated into the deep soil. Therefore, the irrigation at the early growth period can be reduced moderately, and more water should be distributed during the middle and late stages. This irrigation system can effectively reduce deep percolation and increase the utilization of soil storage water.
The simulated actual ET is in good agreement with the observed values. However, the simulated E a is low compared with the study by Yu et al. [21]. The main reasons are as follows: (1) The test was carried out under a shelter, and water can only enter the farmland water circulation system by irrigation at one time. Consequently, the surface soil was dry for a long time. (2) The SWAP model introduces the vegetation cover parameter V c when calculating E a and applies it to the correction of soil surface resistance r a s [28]. As a result, the radiation term and aerodynamic term used to calculate E p are reduced. The water requirement of winter wheat in the reproductive growth stage (emergence-heading stage) was obviously higher than that of the vegetative growth stage (heading-maturation stage), and the average ET was 84.96 mm and 72.77 mm, respectively. This indicated that ensuring adequate water supply in the later growth stages of winter wheat was essential for its healthy growth. The change of ET in different growth stages was closely related to the ratio of E a and T a , especially in the early growth period of winter wheat [23]. The plants are young at the emergence-regreening stage, and transpiration is not high. At this stage, the soil evaporation is dominant and does not change obviously. During the wintering stage, the temperature drops and the SWC decreases result in ET gradually declining. The daily ET remains at a low level (less than 2 mm d −1 ). The winter wheat grows rapidly during the heading-grouting stage, the root system develops, and the leaf area increases. At the same time, the temperature gradually rises. Therefore, the daily ET increases rapidly. The plants withered at the end of the growth stage, the root activity decreased, leaves wilted, and ET decreased. Irrigation had a stronger effect on T a than on E a . This was because T a was closely related to root water uptake. Different irrigation amounts caused variations in soil water status, and these effects were reflected in the changes in root morphology (root length and density distribution) and physiological functions (root vitality and water uptake) [35,[65][66][67]. Under water stress, roots were vulnerable to drought and died. The decrease in the number of living roots affected the ability of the roots to absorb water, resulting in large differences in T a . However, in addition to weather parameters such as temperature and radiation, ground cover was the most important factor affecting E a [23]. Comparing the measured data, the LAI had a slight difference among treatments in the same stage. Thus, the E a fluctuation was small. The water stress at the emergence-jointing stage had little effect on winter wheat, and ET returned to a normal level after rewetting. Water stress during the heading-grouting stage had a great impact on ET, and ET barely recovered after rewetting. This was because the rapid growth of winter wheat needed to consume a large amount of water at this stage. Water shortage inhibited plant growth, root development and assimilation. This once again emphasizes the importance of irrigation for the later growth stages of winter wheat. Therefore, when designing a deficit irrigation schedule for winter wheat, it can be subjected to a certain degree of water stress before heading. This drought exercise can enhance the drought resistance of winter wheat and improve water use efficiency [4]. However, due to the large water requirement of winter wheat in the later growth stage, it is necessary to provide sufficient irrigation to minimize adverse effects [20,23].

Conclusions
This study conducts a sensitivity analysis of the SWAP model through the parameter change method, and concludes that θ s and n are high sensitive parameters that affect the simulation performance of SWC and ET. This clarifies the direction of SWAP model calibration and greatly improves the efficiency of parameter estimation. The method used in this study provides a reference of parameter calibration for crop models without parameter automatic optimization modules. Based on the measured data of the deficit irrigation experiment of winter wheat in Guanzhong Plain, as well as the sensitivity analysis, the SWAP model is calibrated and validated. The calibrated SWAP model has good performance in predicting SWC. The MRE of each treatment is lower than 21.66%, and the RMSE is lower than 0.07 cm 3 cm −3 . The simulated ET of winter wheat during the entire growth stage is in good agreement (R 2 = 0.975, p < 0.01) with the measured values. This makes SWAP a powerful tool to quantitatively analyze the impacts of irrigation on the growth process of winter wheat. The transformation of the farmland water system under different irrigation conditions is analyzed. By investigating the response pattern of winter wheat to different water stress conditions, this study points out the importance of irrigation in the later growth stage of winter wheat, which provides guidance for the formulation of scientific and reasonable irrigation schedules in the actual agricultural production process. Path analysis showed that CH and SWC of the 0-60 cm layer are strongly correlated with daily ET. The ET calculation model constructed by the above key factors has good performance (R 2 = 0.727, p < 0.01) and can be used as an effective method to quickly estimate ET under limited parameters.
However, this study paid less attention to the growth process and dry matter accumulation of winter wheat; therefore, further research should be carried out using detailed module of the SWAP model. It is necessary to directly observe the root growth in future studies to improve the simulation accuracy. Finally, the test plots were arranged under a shelter. When precipitation occurred, the air humidity, light and ventilation conditions inside the shelter were slightly different from those outside. The growth of winter wheat was different from that in the larger field. Therefore, additional field trials are needed to enhance the applicability of the research.