1 Introduction

Crop phenology strongly influences crop yield, biomass accumulation, and other plant processes, such as the magnitude and timing of evapotranspiration (ET), in agricultural ecosystems. Crop yield is a key element of farm profit. ET is the primary hydrologic outflow from croplands and controls the fraction of precipitation and irrigation inflows that is available to contribute to soil water storage, streamflow, aquifer recharge, and ecosystem function [1, 2]. Crop growth and development interacts with management, water availability, and local conditions to influence yields and hydrologic responses [3]. Accurate representation of crop processes is essential to hydrologic simulation in watershed models. Yet watershed models typically use simplistic assumptions about crop growth and development and their influence on crop yield, biomass accumulation, and other plant process, such as ET.

The Environmental Policy Integrated Climate (EPIC [4, 5]) includes a crop growth submodel [6] that has been applied widely and modified in other models to meet specific modeling objectives including Soil Water Assessment Tool (SWAT [7, 8]), Wind Erosion Prediction System (WEPS [9]), and Unified Plant Growth Model (UPGM [10]). Important processes of the EPIC crop growth submodel are described in the Background Information.

Explicit representation of phenological timing in crop growth models is essential for simulating genotype (G) x environment (E) x management (M) crop responses. Phenological timing is strongly influenced by environment, primarily by the effects of temperature [11]. Other “secondary” E x M factors (light, nutrients, salinity, and CO2 [11]) and particularly the influence of water stress on phenology [12,13,14] are rarely considered in phenological simulation [15]. Liu et al. [16] simulated crop stress from waterlogging as a function of phenology and demonstrated the important role of phenological timing in understanding global response of barley to future climate change scenarios. Here, we assess the importance of explicitly simulating phenological timing (WEPS and UPGM) and the effects of water-deficit stress on phenological timing (UPGM only) in simulation of crop yield and biomass growth.

McMaster et al. [15] assessed performance of the three crop models in simulating the phenological timing of five winter-wheat developmental stages under various levels of water-deficit stress. Start of senescence (S) and physiological maturity (M) were compared among SWAT, WEPS and UPGM. Jointing (J) was compared between WEPS and UPGM. Heading (H) and start of anthesis (AS) were assessed for UPGM alone. UPGM was the only model to explicitly simulate all five developmental stages (J, S, H, AS, M). They found that UPGM simulated J, S, and M as well or better than WEPS and always better than SWAT. Unanswered by McMaster et al. [15] is the question of whether the demonstrated superiority in simulating phenological timing translates into improved simulation of crop yield and biomass growth components across varying environments.

The objective of the present study is to evaluate whether the improved phenology model in UPGM also improves the simulation of other important winter wheat growth and development processes (i.e., yield, aboveground biomass, canopy height, harvest index, and leaf area index [LAI]) under various levels of water-deficit stress compared to SWAT and WEPS. Importantly, all three crop models are executed within the Agricultural Ecosystems Services (AgES [17]) hydrology model, thereby using a common platform for other system processes that drive crop growth and development.

2 Background Information

The EPIC crop growth submodel is the basis for the three models studied here. In EPIC, biomass is accumulated daily as a product of radiation use efficiency (RUE) and intercepted photosynthetically active radiation, which is calculated based on simulated LAI. Daily biomass is partitioned among above- and below-ground fractions based on a ratio that changes from 40% (at germination) to 20% (at maturity) [18]. Canopy height is also calculated following a crop-specific LAI-based curve. LAI is accumulated daily according to heat units (HU) and is attenuated by crop stress. Biomass is translated to crop yield as a fraction (harvest index, HI) of aboveground biomass. HI is assumed to be maximized by breeding for a given cultivar and relatively stable across a range of environments [19]. A water-stress factor (Wstrs) is calculated daily based on a ratio of actual water transpiration (Ta) to potential water transpiration (Tp) as,

$$\mathrm{Wstrs}={\mathrm T}_{\mathrm a}/{\mathrm T}_{\mathrm p}$$
(1)

In this formulation, there is no water stress for Wstrs = 1, and the crop is fully water-stressed when Wstrs = 0. Subsequently, Wstrs affects both daily biomass accumulation and HI, where HI is influenced by water stress during crop stages between 0.3 and 0.9 of maturity, as measured by a heat-unit index (HUI), with maximum sensitivity at HUI of 0.6, near crop anthesis. The HUI on day i after planting (Eq. 2) is calculated as a fraction of accumulated heat units (HU, Eq. 3) to theoretical HUs from planting to maturity (THUM):

$$\mathrm{HUI_{i} }=\frac{\sum_{j=1}^{i}{HU}_{i}}{THUM}$$
(2)
$$\mathrm{HU_{i} = T_{avg,i}- T_{base}},\; \mathrm{for\; T_{avg,i} \ge T_{base}}$$
$$\mathrm{HU_{i} = 0,} \;\mathrm{for \;T_{avg,i} < T_{base}}$$
(3)

where HUI is heat unit index (fraction: 0 at planting, 1 at maturity) on day i after planting, HUi is heat units accumulated on day i (°C day, or growing degree-day [GDD]), THUM is theoretical HUs accumulated at physiological maturity (°C day), Tavg is average daily air temperature (°C), and Tbase is crop-specific base temperature (°C). Phenological development in EPIC follows the common approach of accumulating daily HUs (Eq. 3). Only a few developmental events are specified. The crop-specific input value of HU required to reach maturity (THUM) is fixed for each crop. The HUI is used to specify dates of leaf-area growth (different rate equations before and after HUI of crop-specific start of leaf decline), root depth (maximum when HUI reaches 0.4), economic yield (nonlinear HI increase from HUI of 0.5 to 0.9), and physiological maturity (HUI of 1). A couple potential problems with this formulation in Eq. 2 are that delayed seedling emergence and/or vernalization can result in significant errors in estimating maturity and are not reflected in the pre-determined HUI curves. Further, use of these fixed values restricts model ability to simulate the influence of water stress on phenology. Although water stress in EPIC influences biomass accumulation, crop yield, crop height, and LAI accumulation, which in turn affects water use, the HUs required to progress between developmental stages are fixed and independent of water stress. Therefore, while water stress affects the magnitude of crop growth processes, it does not reflect the important changes in timing of crop development.

Of the three models (SWAT, WEPS, and UPGM) evaluated in this study, the SWAT version 2009 [19] adopted the EPIC crop growth submodel with few modifications of the above processes. Each model provides a different approach to simulating phenology, as described by McMaster et al. [10, 15], based on the thermal-time approach to simulating phenological development (Eqs. 2 and 3). SWAT uses the approach described above, and the phenological approaches used in WEPS and UPGM will be described briefly herein, along with other modifications.

WEPS made some significant modifications and additions to the EPIC plant growth component to increase the accuracy and robustness of simulating the canopy important to simulating wind erosion. WEPS also uses the EPIC crop growth submodel with static values of phenological stages that do not respond to water stress. WEPS also defines the HUI beginning at planting (P) but defines a HU-based curve for reproductive partitioning. We associate the double-ridge (DR) event with the beginning of reproductive partitioning. Jointing (J) is assumed to occur 180 GDDs after DR [11, 15, 20]. Unlike EPIC, in which LAI is increased on a HU-based curve, WEPS accumulates and partitions biomass to various plant parts and then calculates the resulting LAI from leaf biomass using a crop-specific leaf area conversion factor [9]. This allows WEPS to develop LAI uniquely each growing season based on actual conditions, including water stress, rather than an assumed curve. WEPS uses an S-shaped curve to partition aboveground biomass accumulation to leaf, stem, and reproductive components based on HU (Eqs. 2 and 3). As such, WEPS uses an HU-based approach, rather than EPIC’s HI-based approach, to simulate yield [10]. HI is then calculated from total aboveground biomass and yield. The curve used for crop height increase is independent of LAI. It is a separately specified curve tied to GDD (stage of growth) and the incremental potential increase applied on each day is reduced by stress level (maximum of water and temperature stress) on that day. These WEPS modifications to crop growth make it more mechanistic in many processes, including determining plant height, partitioning biomass to plant components rather than dividing them up after the fact, and yield.

The underlying rationale for developing UPGM was to combine the various modifications to the EPIC plant growth submodel into one unified crop model that would facilitate parameterization and future model development, while maintaining important differences in addressing slightly differing model objectives. The WEPS plant growth submodel was used as the foundation for UPGM because the underlying structure of the various models was the most modified and therefore made it easier to incorporate other model improvements. A more mechanistic simulation approach was desired for UPGM model development, beginning with a more detailed phenology submodel that responded to varying water deficits that could then be used in defining the trajectory of biomass, yield, canopy height, etc. It also would allow for adjusting the timing of the aboveground biomass curves. UPGM integrated the WEPS plant growth component [10] with a phenological development component from PhenologyMMS [21]. UPGM explicitly defines the thermal time in HUs (GDDs) required to reach many important developmental events for both non-stressed (GN) and water-stressed (GS) growth conditions (Table 1). GN and GS parameters for each crop and developmental stage were taken from PhenologyMMS [20, 21]. UPGM adjusts the HUs required to reach the next developmental event using a piecewise linear interpolation between GN and GS within a range of 0.4 < Wstrs < 0.8 (see Fig. 1a for a schematic representation). Water stress can accelerate development (when GS < GN), slow development (when GS > GN), or have no affect (when GS = GN), depending on the crop and growth stage. Figure 1b shows quantitative functions of thermal time versus Wstrs for eight different developmental stages of winter wheat using the values given in Table 1. Knowing the HUs when a growth stage was reached would then allow this timing to be used rather than a set input parameter that did not respond to water stress.

Table 1 UPGM thermal times (growing degree-days, GDD, 0°C Tbase) between developmental events for winter wheat (Triticum aestivum L.) under non-stressed (GN) and water-stressed (GS) conditions (adapted from 15, 20)
Fig. 1
figure 1

UPGM water stress factor. a Conceptual diagram of thermal time to complete a developmental event with linear interpolation from non-stressed (GN, Wstr = 1) to water-stressed (GS, Wstr = 0) phenological parameters, where water stress can accelerate (decrease thermal time) or delay (increase thermal time) development. b Example of actual developmental stages for winter wheat showing 8 of 10 stages from Table 1 (E →TI (GS = GN) and AS→M (GS < GN) are not shown for graphical purposes only)

3 Materials and Methods

SWAT, WEPS, and UPGM were incorporated into the USDA-ARS Agricultural Ecosystems Services (AgES) hydrologic model [17, 22]. AgES is a spatially distributed hydrologic/water quality model that adapted hydrologic components from the European J2K model [23], nutrient components from the SWAT model [19], and crop growth components from either SWAT, WEPS, or UPGM in an object modeling framework [24]. AgES hydrologic components have been tested extensively [25,26,27] and have provided a consistent basis for applying climate, water stress, and other site factors to the three models under comparison.

For each test site below, parameters in AgES that control soil hydrology and ET were calibrated using “global” and soil layer parameters [15, 17]. Parameters specific to crop-growth simulation in each model were not calibrated; default model parameter values were used. Calibration of the phenology component of crop models is influenced by both model structure and decisions made in the calibration process [28]. In this study, the emphasis was on comparison of model structures rather than modeler decisions. Calibration of crop parameters would, no doubt, improve simulation of all models, but introducing variability from crop-model calibration decisions could also confound the results. Use of default crop parameters aligns with the common practice when applying these three crop models within a watershed model with many fields and management units; even though modelers calibrate other hydrologic parameters, they typically use the default crop parameters without any site-specific calibration. For these reasons, we did not optimize crop parameters for this study.

3.1 Experimental Comparison Data

Two experimental datasets that included phenological data for winter wheat (Triticum aestivum L.) in semi-arid northeastern Colorado were used. Site details are fully described elsewhere as noted below, and only key elements briefly summarized herein. Both data sets measured five developmental events: jointing (J; when the first node is above the soil surface), start of senescence (S; flag leaf blade growth completed as measured by formation of the ligule on the flag leaf, which is the end of leaf appearance on a shoot), heading (H; appearance of spike emerging from the flag leaf, not including awns), start of anthesis (AS; when the first anthers are visible in the spike), and physiological maturity (M; denoted by the absence of all green color in the spike and subtending peduncle) ([29] as described in [15]).

The Greeley dataset [15] was a 3-year study (Sep 2008 to Aug 2011) of various height and maturity classes of hard red and white winter wheat at the USDA-ARS Limited Irrigation Research Farm (LIRF) near Greeley, CO (40° 26′ 57.13″ N; 104° 38′ 12.04″ W; 1427 masl). Soils were deep and well drained fine sandy loam to clay loam. Weather data were available on-site (CoAgMet Station GLY04; Table 2). The full study was comprised of a split-plot design: 24 wheat varieties randomized within 5 irrigation treatments plots with 3 replications (360 plots, 7 m2 each). The mean of the 18 medium-maturity varieties (Above, Ankor, Avalanche, Baca, Bill Brown, BondCL, CO940610, Danby, Goodstreak, Hatcher, Jagalene, Keota, NuDakota, Platte, Prairie Red, RonL, Sandy, Yuma) and two of the irrigation treatments (fully irrigated [Wet], dryland [Dry]) were used in this modeling study. On fully irrigated plots (Wet), surface drip irrigation was applied to minimize water stress (maintaining soil water > 75%); dryland plots (Dry) received no irrigation. Fall planting operations consisted of disking (10-cm depth), fertilizer application (per soil test), roller harrowing, planting (2 Oct 2008, 6 Oct 2009, 8 Oct 2010), and portable sprinkler irrigation (24, 25, and 43 mm of water, respectively) to ensure emergence. Phenology data were typically collected 3 d per week. Dates on which the first shoot and half the shoots within a plot reached a developmental stage were recorded. Other plant measurements also were combined for the 18 medium-maturity varieties. Total aboveground biomass (and leaf blades, stem + leaf sheaths, inflorescence fractions if present) was measured each year at pre-jointing (early spring when still in the rosette form), jointing, anthesis, mid-grain filling (not in 2008–09), and maturity. All samplings before maturity harvested 0.5 m row lengths whereas maturity harvested 1 m row lengths. A plot combine was used to harvest the entire plot for yield and determine percent moisture. Harvest index (HI) was determined from the maturity aboveground biomass sampling. Canopy height was also generally measured each year at pre-jointing (early spring when still in the rosette form), jointing, anthesis, mid-grain filling, and maturity, with 4 to 9 randomly selected points in the plot.

Table 2 Weather and irrigation summary on experimental sites [15]

The Drake site [15, 30] was a 4-year dryland study (harvest years 2004 to 2007) on a 109-ha farm field (Drake Farm) east of Fort Collins, CO near Ault, CO (40.61°N, 104.84°W, 1559–1585 masl). Soils were Wagonwheel coarse silty-loam, Colby fine silty-loam, and Kim fine sandy-loam with 0.65 to 8.4% slopes. Weather data were available onsite ([30]; Table 2). Alternating strips (~ 120 m wide) were in 2-year wheat-fallow rotation. A medium-maturity winter wheat variety (Above) was planted in alternating strips (1 Oct 2003, 23–27 Sep 2004, 1 Oct 2005, 3 Nov 2006). Phenology data were typically collected 3 d per week as a repeated measurement of tagged plants from ten 30 × 30 m sites split between two of the strips (2 subsamples per site, 10 plants per subsample). Dates on which each shoot reached a developmental stage were recorded. Canopy height was measured for each tagged plant. Total aboveground biomass was sampled by hand within the phenological plot areas. Final grain yield was based on yield monitor data at 10-m resolution within each sample area.

3.2 Model Runs and Evaluation

AgES models were developed for each site using onsite weather data. The Greeley model was set up to run the two treatments (Wet, Dry) as independent hydrologic response units (HRUs) over the 3-year study (2009, 2010, 2011 named by year of harvest). The Drake model was simulated as a 56-ha watershed with 27 HRUs; two summit locations, HRU 23 (2004, 2006) and HRU 24 (2005, 2007), were used for this study to minimize run-on effects and help ensure water limitation. AgES was run using each of the three crop models with default crop parameters to focus on the primary objective of comparing performance of three similar (EPIC-based) crop models with enhanced phenology (WEPS, UPGM) and water-stress response (UPGM) without the influence of calibration skill. Models were assessed using several statistical performance measures: root-mean square error (RMSE; [31]), relative error (RE = [predicted mean – observed mean] / observed mean, %), index of agreement (d; [31]), and normalized objective function (NOF = RMSE / observed mean).

4 Results and Discussion

All three crop models performed well in simulating winter wheat yield (Table 3, Fig. 2). Across all sites, years, and soil water conditions, index of agreement (d) exceeded 0.9 for all models (d = 0.97 [SWAT], 0.96 [WEPS], 0.93 [UPGM]). Normalized objective function (NOF) demonstrated slightly decreasing performance from SWAT (0.16) to WEPS (0.23) to UPGM (0.32). WEPS and UPGM generally overpredicted yield (positive RE = 10.5% [WEPS], 19.6% [UPGM]) whereas SWAT slightly underpredicted yield (negative RE = -3.4%), particularly for Greeley Wet. Observed yield for Greeley Dry was similar over the 3 years (3270–3762 kg/ha), but all models simulated greater yield differences (SWAT: 3073–4046; WEPS: 2883–5159; UPGM: 3271–5631 kg/ha). Even without calibration, all models demonstrated skill in simulating wheat yield across the full range of water conditions spanning 10 crop-years. However, no clear benefit was evident for the HU approach to simulating yield used by WEPS and UPGM over the HI approach used by SWAT.

Table 3 Model performance statistics
Fig. 2
figure 2

Yield performance

All three crop models performed well in simulating final aboveground biomass (Table 3, Fig. 3). Performance in simulating final biomass declined slightly compared to yield performance for all models, but d exceeded 0.86 and NOF was less than 0.35 for all models. WEPS (RE = -11.6%, d = 0.91, NOF = 0.30) performed best, followed by UPGM (-13.7%, 0.87, 0.33) and SWAT (-20.2%, 0.86, 0.35). For biomass simulation, WEPS and UPGM, which used an HU partitioning approach outperformed SWAT with its standardized LAI-curve approach. This may point to a benefit in accumulating HUs as attenuated by water stress in simulating crop biomass.

Fig. 3
figure 3

Final aboveground (AG) biomass performance

Multiple in-season observations at the Greeley site demonstrated that WEPS and UPGM closely simulated aboveground biomass throughout the growing season, whereas SWAT results were consistently lower and underpredicted than observed values (Fig. 4). In cases where errors developed in biomass simulation, they generally occurred late in the season. Final biomass was substantially underpredicted by all models for two site years: Greeley Wet 2009 (Figs. 3 and 4) and Drake 2006 (Figs. 3 and 5). For Greeley Wet 2009, in-season data shows that all models performed reasonably well up to day 149, but simulated biomass plateaued after about day 180 while observed biomass increased up to day 199 (Fig. 4a). Other treatment-years at Greeley showed an end-of-season plateau in observed biomass that was captured by simulated results. This may indicate that the models are not capturing a periodic climate or management influence on biomass accumulation. For example, the underpredicted Greeley 2009 and Drake 2006 site years were the warmest and driest years during the grain-filling period (1 May to 30 June) at their respective sites during the study period (Table 2). In some cases, UPGM better tracked early-season while WEPS better captured end-of-season biomass accumulation (Greeley Wet 2011, Drake 2004) or vice versa, with WEPS better tracking early-season biomass accumulation (Greeley Wet 2010, Greeley Dry 2011), which suggests the need to calibrate timing of phenological development stages with biomass accumulation (Figs. 4 and 5).

Fig. 4
figure 4

Aboveground biomass simulated time series and observed data for the Greeley site

Fig. 5
figure 5

Aboveground biomass simulated time series and observed data for the Drake site

SWAT simulates daily accumulated total biomass; WEPS and UPGM also simulate daily biomass partitioned to leaves, stems, straw (leaves + stems), and seed heads (Fig. 6) as well as roots (not assessed in this study). SWAT consistently underpredicted total biomass compared to WEPS and UPGM except for three occasions: UPGM Greeley Wet 2009 (Fig. 6a) and WEPS Greeley Dry 2011 (Fig. 6f) at the very end of the season, and WEPS Drake 2004 HRU 23 (Fig. 6g) for most of the season.

Fig. 6
figure 6

UPGM and WEPS simulated biomass partitioning (total aboveground biomass, leaves, stems, straw [leaves + stems], seed heads) compared to SWAT (total only) and observed data for Greeley Wet (a: 2009, b: 2010, c: 2011), Greeley Dry (d: 2009, e: 2010, f: 2011), Drake HRU 23 (g: 2004, h: 2006), and Drake HRU 24 (i: 2005, j: 2007). Also shown are vertical lines for five observed phenological stages: jointing (J), start of senescence (S), heading (H), anthesis starts (AS), maturity (M)

The considerable effect of water stress on simulated timing of developmental stages by WEPS and UPGM is also evident in Fig. 6, as represented by the five vertical lines shown in each diagram for the simulated date of J, S, H, AS, and M stages. WEPS and UPGM differ in the magnitude and timing of biomass partitioning leading to the differences between the biomass accumulation curves. Biomass data collected near day 150 for 8 of the 10 site years often were well simulated for seed heads (WEPS, except Fig. 6c; and UPGM, except Fig. 6c, i), straw (WEPS, except Fig. 6a, b, g; and UPGM, except Fig. 6a, b, e), and total biomass (WEPS, except Fig. 6g; and UPGM, except Fig. 6a–e, i).

The wheat canopy was consistently simulated by SWAT to achieve maximum height regardless of site and climate factors, whereas UPGM and WEPS better reflected observed heights (Fig. 7). Simulated final canopy height increased in rough agreement with observed heights (d = 0.60 [WEPS], 0.68 [UPGM]), although height was underpredicted by both WEPS (RE = -42%) and, to a lesser extent, UPGM (-34%) (Table 3, Fig. 7). The SWAT method using a HU-adjusted fraction of maximum canopy height was accurate only for non-limiting conditions, such as the fully irrigated Greeley Wet conditions of 2009 and 2011 (Figs. 7 and 8). Overall, SWAT significantly overpredicted height (RE = 49%) (Table 3, Fig. 7c). Multiple in-season observations at the Greeley site demonstrated that all models simulated a mid-season increase followed by late-season plateau (Fig. 8). Although WEPS and UPGM allow simulation of season-specific canopy height development, both models have difficulty capturing both timing and magnitude of the growth without further adjustments to parameters.

Fig. 7
figure 7

Final canopy height (ht.) performance

Fig. 8
figure 8

Canopy height simulated time series and observed data for the Greeley site

By design, HI in SWAT is relatively stable (here, about 0.50) across environmental conditions, while HI in WEPS and UPGM better predicts the observed values (Fig. 9). HI ranged from 0.23 to 0.60 across sites and years, so WEPS and UPGM, which can account for this variability, have a clear advantage. All models generally overpredicted HI. WEPS (RE = 25%, d = 0.61, NOF = 0.36) performed best, followed by UPGM (35%, 0.54, 0.47) and SWAT (38%, 0.43, 0.48) (Table 3, Fig. 9). Interestingly, SWAT had the greatest overprediction of HI (RE = 38%), which combined with the greatest underprediction of aboveground biomass (-20%) and resulted in the most accurate prediction of yield (-3.4%). WEPS and UPGM do not use HI in calculation of yield, so a similar comparison is not relevant.

Fig. 9
figure 9

Harvest index performance

No LAI observations were made, so model performance can only be assessed by comparison among models (Fig. 10). For SWAT, the timing varied slightly among plots and years, but the pattern was similar with a consistent peak LAI of 3.0, because SWAT was not responding to variable water conditions. The magnitude and timing of LAI varied only slightly among site years for SWAT but varied considerably for WEPS and UPGM in response to site conditions. In many cases, WEPS and UPGM had greater peak LAI than SWAT, sometimes exceeding a reasonable maximum LAI of 5. Examining leaf senescence after peak LAI, all three models overpredicted final LAI.

Fig. 10
figure 10

Leaf area index (LAI) simulated time series comparison for each site and year

5 Conclusions

Three crop models with lineage to the EPIC crop growth submodel were compared to experimental winter wheat (Triticum aestivum L.) growth and yield data. Complexity in representing crop phenology and response to water stress increased from SWAT to WEPS to UPGM, and our hypothesis was that this increasing complexity would result in increasingly accurate simulation of all crop growth variables across a wide range of environmental conditions. This hypothesis was generally confirmed, particularly when considering most of the simulated processes were very similar in WEPS and UPGM. All crop models performed well in simulating yield (index of agreement [d] ≥ 0.93). WEPS and UPGM outperformed SWAT for other growth variables, including biomass, canopy height, and harvest index.

These results show benefits from more explicit simulation of phenology and response to water stress in crop models. We expect that the framework provided by explicit phenological simulation will allow further improvements in model performance that are not possible without the ability to simulate development-specific differences in biomass partitioning, leaf-area expansion, and grain filling. Addition of improved water-stress response will provide better simulation of future, non-stationary climate impacts on crop growth and development, which will inform farm to regional crop management and planning decisions. Such enhancements are also expected to improve simulation of other crop processes, such as evapotranspiration and biomass residue production and management.

This study focused on winter wheat, an important crop regionally and globally in semi-arid climates, where irrigation supply is limited and water stress is common during the growing season. Further development is needed to expand phenological parameterization to other crops, particularly those grown in water-limited production systems.