Identifying Within-Field Spatial and Temporal Crop Water Stress to Conserve Irrigation Resources with Variable-Rate Irrigation

: Addressing within-ﬁeld and within-season variability of crop water stress is critical for spatially variable irrigation. This study measures interactions between spatially variable soil properties and temporally variable crop water dynamics; and whether modelling soil water depletion is an effective approach to guide variable-rate irrigation (VRI). Energy and water balance equations were used to model crop water stress at 85 locations within a 22 ha ﬁeld of winter wheat ( Triticum aestivum L.) under uniform and spatially variable irrigation. Signiﬁcant within-ﬁeld variability of soil water holding capacity (SWHC; 145–360 mm 1.2 m − 1 ), soil electrical conductivity (0.22–49 mS m − 1 ), spring soil water (314–471 mm 1.2 m − 1 ), and the onset of crop water stress were observed. Topographic features and modelled onset of crop water stress were signiﬁcant predictors of crop yield while soil moisture at spring green-up, elevation, and soil electrical conductivity were signiﬁcant predictors of the onset of crop water stress. These results show that modelling soil water depletion can be an effective scheduling tool in VRI. Irrigation zones and scheduling efforts should consider expanding to include temporally dynamic factors, including spring soil water content and the onset of crop water stress.


Introduction
Worldwide demand for freshwater is increasing to accommodate growing populations, declining groundwater levels, deteriorating water quality, environmental regulation and other factors [1]. Irrigation in agricultural and urban systems currently consume 80-90% of the United States freshwater use, and the Food and Agriculture Organization of the United Nations (FAO) estimates irrigation withdrawals will need to increase 14% by 2030 to meet world food demand [2]. Climate conditions often do not allow for increasing freshwater demand, such as the Western United States which is in the midst of a "mega drought" [3]. With continued demand and scarcity of freshwater resources, irrigated agriculture must reduce irrigation resources while maintaining food production [1].
Variable-rate irrigation (VRI) is a promoted tool to improve water use efficiency by spatially matching irrigation rates to crop water demand. Field scale spatial variability of soil and crop water status is well documented, and utilizing VRI to match spatial crop water demand during a growing season could decrease irrigation inputs while maintaining crop yield [4][5][6]. This was observed in potato (Solanum tuberosum L.) production where the total irrigation consumption was not significantly greater, but water productivity improved under VRI management [7]. Delineating irrigation zones from apparent electrical conductivity and utilizing a water balance approach to inform irrigation rates, VRI reduced

Site Description, Field Observations, and Management Practices
The study was conducted for two growing seasons (2016 and 2017) in a commercial soft white winter wheat (Triticum aestivum L.) production field (22 ha) near Grace, ID, USA, (42.611113 N, 111.784244 W). The soil is a coarse-silty, mixed, superactive, frigid Calcic Haploxerolls with 1-4% slopes. The field contains 0.3 ha of unfarmed rock outcroppings ( Figure 1a) with slight topographic and top-soil textural variability with analysis at 42 sites each classified as silty clay loam ( Table 1). The mean annual precipitation is 392 mm with much of the precipitation falling as winter snowfall and spring rain. Irrigation is from a 380 m long center-pivot irrigation system equipped with VRI control (GrowSmart Precision VRI, Lindsay Zimmatic, Omaha, NE, USA).
In 2016, the irrigation management followed the grower's standard practice (GSP) with irrigation events occurring about every four days. The VRI system was utilized to eliminate irrigation to uncropped, rock outcroppings within the field, but otherwise irrigation rates were uniform over the cropped areas. Thus, the 2016 year was used to evaluate the ability to model soil water variation under uniform irrigation. A spatially variable irrigation pattern was used in 2017 to evaluate the model with variation in management. The spatially variable irrigation pattern in 2017 included two irrigation zones with superimposed GSP control strips. The two irrigation zones were a low irrigation treatment (Low = 0.7 * GSP) and a high irrigation treatment (High = 1.3 * GSP). Irrigation zones were delineated from three years of historic yield, and growers experience with inherent production and dynamic water properties. Locations where productivity was average or above received the high irrigation rate. The low irrigation rates were applied at locations with historic low yield potential [19]. Planting dates were 5 October 2015 and 10 October 2016 and harvest dates Agronomy 2021, 11, 1377 3 of 14 were 16 August 2016 and 30 August 2017. Yield data was processed prior to analysis by removing erroneous data points which were defined as outside the limits of ±75% of the median yield and where a full header width was not harvested [20]. Electrical conductivity (EC) was collected using a Veris EC 3100 (Veris Technologies, Salina, KS, USA) in parallel transects spaced at 20 m following field rows. Data was collected at two depths "deep" or "shallow" readings, however the two were highly correlated and for this manuscript, only the "shallow" (0-0.9 m) reading was included in analysis. = 0.7 * GSP) and a high irrigation treatment (High = 1.3 * GSP). Irrigation zones were delineated from three years of historic yield, and growers experience with inherent production and dynamic water properties. Locations where productivity was average or above received the high irrigation rate. The low irrigation rates were applied at locations with historic low yield potential [19]. Planting dates were October 5, 2015 and October 10, 2016 and harvest dates were August 16, 2016 andAugust 30, 2017. Yield data was processed prior to analysis by removing erroneous data points which were defined as outside the limits of ±75% of the median yield and where a full header width was not harvested [20]. Electrical conductivity (EC) was collected using a Veris EC 3100 (Veris Technologies, Salina, KS, USA) in parallel transects spaced at 20 m following field rows. Data was collected at two depths "deep" or "shallow" readings, however the two were highly correlated and for this manuscript, only the "shallow" (0-0.9 m) reading was included in analysis. Figure 1. (a) Bare soil image with superimposed elevation contours and soil sample locations. Dark colored areas are unfarmed exposed bedrock, and (b) irrigation treatments for the 2017 growing season. GSP = grower's standard practice, High = 1.3 * GSP rate, Low = 0.7 * GSP rate, and Zero = no irrigation. Small squares are the paired t-test plots used to compare yields across irrigation treatments.

Modelling Water Dynamics and Crop Stress
Modeling water dynamics for the 2016 and 2017 seasons followed the American Society of Civil Engineers FAO Penman-Monteith estimation of ET to adjust daily soil water levels [21]. This model uses the following equation to estimate ET: where ETc adj is the adjusted crop ET, Ks is the transpiration reduction factor reflecting conditions where soil water availability limits the rate of ET, Kc is the crop coefficient that reflects the crop development stage, and ETo is the reference crop ET [8,12,21,22]. Soil core samples were collected at 85 points within the field site for each season at spring crop green-up and post-harvest to measure soil water content and estimate plant available soil water (Figure 1a). Root zone soil profile water depletion was modeled with a daily time step at each sample point. Modelling began on the day of spring sampling and continued through post-harvest sampling. Spatially uniform non-stressed reference Dark colored areas are unfarmed exposed bedrock, and (b) irrigation treatments for the 2017 growing season. GSP = grower's standard practice, High = 1.3 * GSP rate, Low = 0.7 * GSP rate, and Zero = no irrigation. Small squares are the paired t-test plots used to compare yields across irrigation treatments.

Modelling Water Dynamics and Crop Stress
Modeling water dynamics for the 2016 and 2017 seasons followed the American Society of Civil Engineers FAO Penman-Monteith estimation of ET to adjust daily soil water levels [21]. This model uses the following equation to estimate ET: where ET c-adj is the adjusted crop ET, K s is the transpiration reduction factor reflecting conditions where soil water availability limits the rate of ET, K c is the crop coefficient that reflects the crop development stage, and ET o is the reference crop ET [8,12,21,22]. Soil core samples were collected at 85 points within the field site for each season at spring crop green-up and post-harvest to measure soil water content and estimate plant available soil water (Figure 1a). Root zone soil profile water depletion was modeled with a daily time step at each sample point. Modelling began on the day of spring sampling and continued through post-harvest sampling. Spatially uniform non-stressed reference evapotranspiration (ET o ) and a crop coefficient (K c ) were used to estimate the crop specific ET (ET c ) for wheat and the remaining soil water in the root zone. When depth of soil water fell below a site-specific readily available water content (RAW) threshold, ET c was adjusted with a site-specific crop stress coefficient (K s ) to estimate the adjusted crop water stressed ET (ET c-adj ). A two-fold approach was used to validate the model output. First, the season long sum of modeled daily ET c-adj was compared to the ET measured from the following season-long water balance equation: where P represents the season-long total precipitation, ∆S represents the measured difference in soil water between spring green-up and harvest, I is season-long total applied irrigation, and D is season-long total drainage of water below a soil depth of 1.  (Figure 1a). The second approach used to validate the model was to compare the model predicted soil water content at the end of each growing season with the measured soil water content. Statistical measures of model fit included: (1) relative error (RE), (2) normalized objective function (NOF) (3) root mean square error (RMSE), and (4) r-squared where RE represents model bias, NOF indicates the fit of the model (where NOF < 1 is less than one standard deviation from the mean), and RMSE is the difference between predicted and measured values [23].
Weather data were collected at a weather station~2 km from the field site (42.51496 N, −111.73606 W; Cooperative Agricultural Weather Network AgriMet). The weather data were used to calculate daily evapotranspiration (ET o ) for a reference crop using ASCE Standard Penman-Monteith Evapotranspiration model [21]. The crop coefficient approach was used to calculate daily ET c-adj using Equation (1), where K s was calculated daily for each sample point and the K c curve was estimated from table values for a similar environment [21] and validated with field observations. The K s was calculated as: where SWHC is the SWHC in the 1.2 m deep soil profile, p is the table value 0.55 that represents the average fraction of SWHC that can be depleted before crop stress for winter wheat, and D r is the current root zone water depletion [21]. Daily depth of soil water in the 1.2 m deep soil profile was modeled based on SWHC and readily available soil water (RAW) using: where SWHC is calculated using Equation (4) for the crop rooting depth, θ FC is the volumetric water content (m 3 m −3 ) at FC, θ WP is the volumetric water content at wilting point (m 3 m −3 ), and Z r is the estimated rooting depth of the crop and was assumed to be 1.2 m for this study. RAW is the readily available soil water content in the root zone (mm), and p is defined above. Field capacity measurements followed Martin et al. (1990) which states that a good method to estimate FC is sampling 1-3 d after a thorough soil-wetting event [24]. At this site, winter snowfall and spring thaw act as the wetting event and FC was estimated for each sampling point as the greatest observed depth of soil water for each depth increment as measured at that site over a two field and lab measurement error (n = 39, RMSE = 11 mm). Daily root zone depletion was estimated using the following equation: where D r,i is the root zone soil water depletion at the end of d i , D r,i−1 is the root zone depletion at the end of the previous day, P i is the precipitation on d i , I i is the irrigation depth on d i , ET c,i is the crop evapotranspiration on d i , and DP i is deep percolation out of the root zone on d i .

Statistical Analysis
For 2017, paired t-tests were performed using 128, 10 m 2 paired locations to evaluate responses in yield located parallel to the GSP control strips and spatially variable irrigation zones ( Figure 1b). A buffer distance of 10 m from treatment edges was applied to avoid irrigation drift between zones. A significance level of p < 0.05 was used. One of the physical joints in the center pivot irrigation system leaked during irrigation events. This created inaccurate data points, which were excluded from all analysis. JMP Pro 13 (SAS Institute, Cary, NC, USA) was used for paired t-test and spatial analysis. Multiple linear regression was performed using a best subset approach using the olsrr package in R studio (R Foundation for Statistical Computing, Vienna, Austria) [25]. Four regression models were fitted, with yield and day of onset of crop water stress for 2016 and 2017 being used as the dependent variables. The potential independent variables for onset of crop water stress were: depth of soil water in the 1.2 m deep soil profile at spring green-up, apparent electrical conductivity 0-0.9 m depth (ECa −0.9 ), slope, elevation, SWHC, and the spatially variable irrigation treatments in 2017. The same potential independent variables were used, in addition to the modelled onset of crop water stress, in regression models for predicting grain yield. Variograms were computed and all field and topographic attributes were kriged to a 5 m grid using SpaceStat (Jacquez et al., 2014). Elevation data was collected from multiple years of yield monitor data and topographic attributes computed at a 5 m grid using ArcMap (Environmental Systems Research Institute [ESRI], ArcGIS desktop: Release 10.5, Redlands, CA, USA).

Spatial Variation of Soil Properties and Topographic Features
Soil properties related to water dynamics varied throughout the field (  (Figure 2c). EC a can often reflect topsoil properties, but did not correlate well with elevation, slope or soil water holding properties (FC, WP, SWHC) at this field site ( Figure 2). Elevation was greatest (1712 m) at the north end of the field and lowest (1705 m) on the west half of the field (Figure 2e) and there was moderate variability in slope angle, with the greatest slope (4-7%) dissecting the field north to south (Figure 2d). Field variability in soil properties confirm other observations where significant variability in soil water dynamics are observed, even in seemingly uniform production fields [4,5,11,14,26].

Model Validation
Spatial and temporal variation of soil water content was modelled at the 85 locations throughout the growing season using a daily water balance (Equation (2)) and ASCE standardized equation to calculate daily ET (Equation (1)). Two approaches were used to validate model output. First, the modelled season-long total ET (sum of daily ET over the growing season) was compared to the ET measured from a season-long water balance for each sampling location (Equation (2)). Modelled seasonal ET values correlated well with measured ET with relative errors of 0.1% and 2.4% for uniform irrigation (2016) and spatially variable irrigation (2017)

Model Validation
Spatial and temporal variation of soil water content was modelled at the 85 locations throughout the growing season using a daily water balance (Equation (2)) and ASCE standardized equation to calculate daily ET (Equation (1)). Two approaches were used to validate model output. First, the modelled season-long total ET (sum of daily ET over the growing season) was compared to the ET measured from a season-long water balance for each sampling location (Equation (2) The two validation approaches support the modelling approach used to predict within field variation of soil water dynamics and crop water stress over time. The model performed better under uniform irrigation compared to spatially variable irrigation. There are two possible explanations for this. First, daily ET was calculated using spatially uniform crop coefficients, regardless of irrigation management. It is probable that the ±30% irrigation rate differences in 2017 altered crop canopy development enough to justify spatially variable Kc values. In practical settings, except where there is extreme textural variability, VRI rates will not require 30% differences in irrigation rates. This reduces the likelihood that variable Kc values are necessary in a production environment. However, further research into model improvement using spatially variable Kc is required. Second, the low irrigation rate associated with spatially variable irrigation led to significantly greater Ks values (reduction in ET from crop water stress), which likely introduced more uncertainty into the model. This error could be quantified and overcome by in-season soil moisture measurements to calibrate model estimated Ks values. The results support the use of this modelling approach to schedule variable irrigation rates, but the approach may be improved by using in-season dynamic factors such as soil moisture and crop canopy conditions. This idea is consistent with other VRI studies, that show improvement in irrigation management when dynamic crop conditions are used [27]. Currently, these dynamic measurements can require significant manual labor which can be cost and time prohibitive, however, remote sensing techniques to estimate Kc are currently being developed to overcome these restrictions [28]. The two validation approaches support the modelling approach used to predict within field variation of soil water dynamics and crop water stress over time. The model performed better under uniform irrigation compared to spatially variable irrigation. There are two possible explanations for this. First, daily ET was calculated using spatially uniform crop coefficients, regardless of irrigation management. It is probable that the ±30% irrigation rate differences in 2017 altered crop canopy development enough to justify spatially variable K c values. In practical settings, except where there is extreme textural variability, VRI rates will not require 30% differences in irrigation rates. This reduces the likelihood that variable K c values are necessary in a production environment. However, further research into model improvement using spatially variable K c is required. Second, the low irrigation rate associated with spatially variable irrigation led to significantly greater K s values (reduction in ET from crop water stress), which likely introduced more uncertainty into the model. This error could be quantified and overcome by in-season soil moisture measurements to calibrate model estimated K s values. The results support the use of this modelling approach to schedule variable irrigation rates, but the approach may be improved by using in-season dynamic factors such as soil moisture and crop canopy conditions. This idea is consistent with other VRI studies, that show improvement in irrigation management when dynamic crop conditions are used [27]. Currently, these dynamic measurements can require significant manual labor which can be cost and time prohibitive, however, remote sensing techniques to estimate K c are currently being developed to overcome these restrictions [28].

Uniform Irrigation Season (2016)
Significant spatial variation in measured soil water content at spring green-up was observed in 2016, with a field average initial soil profile water depletion of 12 mm below FC (Figure 4a). Spatial variability of initial water content values varied between some sites at the site-specific FC with other areas initial water depletions near 50 mm below FC (Figure 4a). Irrigation events began on day 151 (3 May) with a season-long irrigation total of 198 mm applied over 11 events and the last irrigation on day 193 (12 July). Crop physiological peak crop water demand occurred between days 170 and 191 (19 June-10 July), which corresponded with critical grain development stages. The onset of crop water stress was defined as the first day that model depletion of soil water in the 1.2 m profile to be greater than the site-specific RAW for a given field location; or simply the first day the crop experienced water stress. Crop water stress first occurred on day 175 (24 July) at field location 20, while the last field location to begin experiencing crop water stress began day 196 ( Table 2; 15 July) at field location 37 (Figure 1). This 19-day variation in the onset of crop water stress was evident in the modelled differences in daily ET. For example, on day 196 (15 July), ET c-adj values were 8.4 mm and 5.4 mm for field location 37 and 20, respectively. This 36% decrease in ET c-adj reflects the spatial differences in water depletion that are inherent at this field site under uniform irrigation. Addressing these spatial differences in water depletion illustrates the potential for VRI to improve irrigation management and yield by limiting time spent in water stress for various field locations. However, while significant variability exists, there is limited research identifying whether this degree of variability would justify VRI implementation at this site. Developing these benchmarks will guide determination of suitable field locations for VRI and act as a springboard to farmer adoption and implementation in water scarce regions. Table 2. Field average, and season totals in precipitation, irrigation, and evapotranspiration (ET) under uniform irrigation (2016), and partitioned between spatially variable irrigation treatments (2017). Included are the average and first and last ordinal day that locations experienced crop water stress; the earlier in the season, the more water stress experienced and lost yield. Average  First  Last  mm  Ordinal Day   2016  95  198  520  188  175  196  2017  90  188  155  211  2017 Irrigation Treatments  Low  158  437  174  155  187  GSP  225  509  187  158  195  High  291  573  203  186  211 Linear regression models were used to identify key factors affecting the onset of crop water stress. The best-fit model (p < 0.01; r 2 = 0.52) identified elevation, EC a , and soil water content at spring green-up as significant variables (Table 3). Elevation and EC a are already used to delineate VRI zones [8,12] while at this study site, spring soil moisture at green up was more significant in predicting the onset of crop water stress than the variability of model parameters such as FC and WP. This phenomenon is likely because the stored soil moisture at spring green-up represented 47% of the total seasonal water consumption of the grain crop. Under uniform irrigation, the spatial patterns of stored soil moisture appears to be a potentially important factor to inform spatially variable irrigation rates later in the season.

Season Precipitation Irrigation ET Onset of Crop Water Stress
Significant variation was observed in 2016 grain yield, which ranged from 2.3-10.6 Mg ha −1 and averaged 7.5 Mg ha −1 ( Figure 5). Multiple linear regression models were evaluated to identify key factors affecting grain yield. The best-fit model (p < 0.01; r 2 = 0.25) identified slope, elevation, and predicted onset of crop water stress as significant variables (Table 3). This observation supported the hypothesis that the onset of crop water stress would be significantly correlated to variability in grain yield. However, the poor r 2 value reflects the numerous variables that influence yield. Historic yield patterns (data not shown) for this field generally confirm that the areas with greatest yield are at the lowest elevations while the lowest yielding locations are areas with the greatest slopes. The model used in this study does not account for potential lateral flow of soil water, which may partially explain the significant relationships between yield and topography. Additionally, the soil properties that are commonly used to delineate management zones (e.g., texture, soil EC a ) were not as significant in yield production as topography or dynamic factors found in this field. Table 3. Regression model coefficients and performance for predicting yield and the onset of modeled crop water stress as identified through best subset model selection. For variables included in the regression a coefficient and an indication of statistical significance are given. Variables not included in the regression are marked with (-). Variables are grouped according to topographic features, dynamic factors, and soil properties. Each listed regression model was significant (p-value < 0.001).

Spatially Variable Irrigation (2017)
Spatial variation in measured soil water content at spring green-up was greater in 2017 than it was in 2016, with an average soil profile water depletion of 35 mm in the 1.2 m deep profile and some areas of the field having soil water depletions as high as 100 mm (Figure 3). At this field site, the variation of water content at spring green-up appears to be a function of snow accumulation and melting patterns over the winter and early spring. Variable irrigation treatments were applied at each irrigation event with season-long total irrigation applications of 158, 225, and 291 mm for the low, GSP, and high treatments, respectively (Figure 1b). Modeled water stress first occurred on day 155 (4 June) at point 29 in a zone with the low irrigation treatment (Figure 4c). The last modelled onset of crop water stress occurred on day 211 (30 July) at point 39 in a zone with the high irrigation treatment (Figure 4e). Peak ET demand occurred between days 166-206 (15 June-25 July). Variable irrigation treatments, together with higher spatial variation in soil water content at spring green-up, resulted in a span of 56 days between the first and last onset of crop water stress, compared to a span of 19 days under uniform irrigation in 2016. Irrigation treatment averages for the modelled day of onset of crop water stress were days 176, 184, and 205 (25 June, 3 July, and 24 July) for the low, GSP and high irrigation rate zones, respectively. Additionally, the within irrigation zone variability of the onset of crop water stress was 32, 37 and 25 days between the first and last point to experience crop water stress for the low, GSP, and high treatments, respectively (Table 2). ET c-adj was more spatially variable under variable irrigation treatments than uniform irrigation. For example, on day 211 (30 July) ET c-adj was 7.0 mm for field location 39 and was 2.9 mm for field location 29, a 59% decrease in daily ET c-adj . The winter wheat at field location 29 experienced water stress for the entire period of peak crop water uptake, while the winter wheat at field location 39 did not experience water stress once during that period (Figure 4). This further demonstrates the linkage between localized crop water demand and variation of onset of crop water stress.
The regression models for predicting the onset of crop water stress under spatially variable irrigation (p < 0.01; r 2 = 0.84) showed dependence on irrigation treatment, but also identified elevation, EC a , and soil water content at spring green-up as significant variables ( Table 3). All of these factors were identified in the best-fit models for both study years, demonstrating their influence on the onset of crop water stress for this field ( Table 3). The impact of soil water content at spring green-up on predicted late season crop water stress was confirmed in both study years, signifying that VRI management can be improved with consideration of dynamic factors such as the variation of soil water content at spring green-up and the onset of crop water stress.
Significant spatial variation was observed in 2017 grain yield, which ranged from 1.8-8.4 Mg ha −1 and averaged 5.8 Mg ha −1 ( Figure 5). The field average yield in 2017 decreased 23% compared to 2016 (Table 4). Declines in second year wheat yield are common occurrence in wheat-wheat-potato crop rotations [29]. For this farm site and operation, second year yield decreases of 27-33% are common and can reach up to 64% excluding disease or similar factors. The yield reduction of only 22% (less than the average at this site for second year wheat) suggests an overall yield benefit from spatially variable irrigation (Figure 5c). Yield averages for the three irrigation zones were 4.9, 6.0, and 6.4 Mg ha −1 for the low, GSP, and high irrigation treatments, respectively (Table 4). Paired t-test results showed significant differences in yield between the irrigation zones, with a p-value < 0.05 for both the high (6% yield increase) and low (11% yield decrease) irrigation treatments relative to the GSP. However, yield decreases between the GSP and low treatment zones are common as seen in Table 4 and suggests yield is limited by other factors than water availability. If not economically feasible to overcome (poor growing environment), perhaps decreasing irrigation applications at these locations will produce greater production per unit of applied water (crop water productivity) compared to the GSP [30]. Water scarce regions commonly recommend measuring crop water productivity to evaluate water management practices. Further work is necessary to investigate the spatial patterns of crop water productivity in production fields and how to VRI can optimize this metric to improve water conservation in irrigated agriculture.  Similar to observations in 2016, the best fit multiple linear regression model for predicting yield (p < 0.01; r 2 = 0.33) in 2017 identified elevation and the predicted onset of crop water stress as significant variables ( Table 3). The regression model showed that the later onset of crop water stress corresponded to greater grain yield. The variable irrigation treatments in 2017 created more extreme differences in the onset of crop water stress, thus increasing the importance of this factor on yield compared to uniform irrigation management in 2016. These observations further confirm that modelling of the onset of crop water stress can be useful in predicting within-field yield variation, or more importantly can be used to inform variable-rate irrigation zones and rates to minimize crop water stress. As observed in 2016, however, the relatively poor r 2 for the best-fit model to predict yield suggests that many other factors influence crop yield. In 2017, slope was not included as a significant topographic feature in the best-fit model but it was in 2016. However, elevation was significant in both years, demonstrating the importance of topography in the spatial variation of water dynamics and crop yield. Similar to observations in 2016, the best fit multiple linear regression model for predicting yield (p < 0.01; r 2 = 0.33) in 2017 identified elevation and the predicted onset of crop water stress as significant variables ( Table 3). The regression model showed that the later onset of crop water stress corresponded to greater grain yield. The variable irrigation treatments in 2017 created more extreme differences in the onset of crop water stress, thus increasing the importance of this factor on yield compared to uniform irrigation management in 2016. These observations further confirm that modelling of the onset of crop water stress can be useful in predicting within-field yield variation, or more importantly can be used to inform variable-rate irrigation zones and rates to minimize crop water stress. As observed in 2016, however, the relatively poor r 2 for the best-fit model to predict yield suggests that many other factors influence crop yield. In 2017, slope was not included as a significant topographic feature in the best-fit model but it was in 2016. However, elevation was significant in both years, demonstrating the importance of topography in the spatial variation of water dynamics and crop yield.
The variables included in predicting the onset of crop water stress and yield included topographic features (slope and/or elevation) and soil properties (soil EC a or SWHC). Many studies have shown use of soil factors for establishing static VRI zones [7,8,12]. Study results support these earlier findings and further suggest including topographic and dynamic soil/crop factors could improve VRI zone delineation and irrigation management. While there was significant variation in soil properties and their associated temporal variation in the onset of crop water stress, the topographic soil properties were more strongly related to grain yield. Soil profile water content at spring green-up was important for modelling the onset of crop water stress and the onset of stress was a key factor in predicting yield. Onset of crop water stress and soil profile water content at spring green-up are dynamic variables that could improve irrigation decisions for VRI systems. Others have proposed different approaches for incorporating variation in crop water stress into VRI management, such as use of plant feedback systems [27]. This compliments other recommendations that integrated decision support tools are required to optimize the allocation of limited water [1,10]. Incorporating these dynamic and topographic features into spatially variable irrigation could improve VRI decision support tools. The variables included in predicting the onset of crop water stress and yield included topographic features (slope and/or elevation) and soil properties (soil ECa or SWHC). Many studies have shown use of soil factors for establishing static VRI zones [7,8,12]. Study results support these earlier findings and further suggest including topographic and dynamic soil/crop factors could improve VRI zone delineation and irrigation management. While there was significant variation in soil properties and their associated temporal variation in the onset of crop water stress, the topographic soil properties were more strongly related to grain yield. Soil profile water content at spring green-up was important for modelling the onset of crop water stress and the onset of stress was a key factor in predicting yield. Onset of crop water stress and soil profile water content at spring greenup are dynamic variables that could improve irrigation decisions for VRI systems. Others have proposed different approaches for incorporating variation in crop water stress into VRI management, such as use of plant feedback systems [27]. This compliments other recommendations that integrated decision support tools are required to optimize the allocation of limited water [1,10]. Incorporating these dynamic and topographic features into spatially variable irrigation could improve VRI decision support tools.

Conclusions
This study characterized and modelled within-field variation of soil water relationships with crop water stress and observed yield of winter wheat. These results demonstrate the potential of using soil water depletion modeling to estimate spatial and temporal crop water demand in VRI. In this study, spatially variable irrigation treatments were created based on an analysis of historic yield. Modelled spatial variability in soil water depletion and the onset of crop water stress, however, was highly related to soil properties, topographic factors, and the variation in soil water content in early spring. Modelled within-field variation of the onset of crop water stress was strongly related to observed grain yield for both uniform and spatially variable irrigation treatments. Predicting the onset of crop water stress was highly dependent on soil properties, topography, and initial soil water content in early spring. Integrating soil properties, topographic factors, and soil water at planting, into decision support tools could improve water conservation under VRI management. Some dynamic factors can be universal [onset of crop

Conclusions
This study characterized and modelled within-field variation of soil water relationships with crop water stress and observed yield of winter wheat. These results demonstrate the potential of using soil water depletion modeling to estimate spatial and temporal crop water demand in VRI. In this study, spatially variable irrigation treatments were created based on an analysis of historic yield. Modelled spatial variability in soil water depletion and the onset of crop water stress, however, was highly related to soil properties, topographic factors, and the variation in soil water content in early spring. Modelled within-field variation of the onset of crop water stress was strongly related to observed grain yield for both uniform and spatially variable irrigation treatments. Predicting the onset of crop water stress was highly dependent on soil properties, topography, and initial soil water content in early spring. Integrating soil properties, topographic factors, and soil water at planting, into decision support tools could improve water conservation under VRI management. Some dynamic factors can be universal [onset of crop water stress], but others are likely field specific [spring soil moisture] and require further investigation to identify and employ these into decision support tools. Further work should focus on establishing thresholds in spatial or temporal factors that justify adopting spatially variable irrigation. Leveraging this information will allow farmers to identify environments where conditions justify VRI implementation, thereby maintaining or increasing crop yield while limiting freshwater consumption.