Temporal Downscaling of Crop Coefficients for Winter Wheat in the North China Plain: a Case Study at the Gucheng Agro-meteorological Experimental Station

The crop coefficient (K c) is widely used for operational estimation of actual evapotranspiration (ET a) and crop water requirements. The standard method for obtaining K c is via a lookup table from FAO-56 (Food and Agriculture Organization of the United Nations Irrigation and Drainage Paper No. 56), which broadly treats K c as a function of four crop-growing stages. However, the distinctive physiological characteristics of overwintering crops, such as winter wheat (Triticum aestivum L.), which is extensively planted in the North China Plain (NCP), are not addressed in this method. In this study, we propose a stage-wise method that accounts for K c variations for winter wheat at each critical phenological stage, thereby estimating K c at finer temporal scales. Compared with the conventional FAO method, the proposed stage-wise method successfully captures the bimodal pattern in K c time series for winter wheat, which is shown at both ten-day and phenological time scales. In addition, the accuracies of the proposed stage-wise K c method and the FAO method were evaluated using micro-meteorological measurements of ET a collected at the Gucheng agro-meteorological experimental station in the NCP. Using a leave-one-out strategy, the evaluation revealed that the stage-wise method significantly outperformed the FAO method at both daily and critical phenological time scales, with root-mean-square errors in ET a for the stage-wise method and the FAO method being 0.07 mm·day −1 and 0.16 mm·day −1 , respectively, at the daily time scale, and 0.01 mm·day −1 and 0.27 mm·day −1 at the critical phenological time scale. Generally, the FAO method underestimates ET a during the initial stage and overestimates ET a during both the development and mid-season stages. It is shown that the proposed stage-wise method is important for the water-stressed NCP where precision irrigation is highly desirable, especially during the critical phenological stages. Results from this study provide insight into accurate estimation of water requirements for winter wheat at phenological time scales.


Introduction
The crop coefficient (K c ) is a crucial parameter that is widely applied in agricultural water management and irrigation scheduling, as it reflects the impacts of inherent crop biological characteristics and planting conditions on water requirements.K c is defined as the ratio of crop evapotranspiration under well-watered conditions (ET c ) to a reference crop evapotranspiration (ET 0 ).ET c is the evapotranspiration expected from disease-free, well-fertilized crops grown in large fields, under optimum soil water conditions and achieving full production under the given climatic conditions [1].ET 0 is the evapotranspiration rate expected from a reference surface, which is hypothesized as a grass reference crop with specific characteristics limited only by atmospheric demand [1], considering only a constraint from the available energy.
Among these, the FAO method is most popular given its simplicity and limited demands on data and processing steps.The FAO approach considers the crop coefficient as a function of four crop-growing stages, namely the initial, development, mid-season and late-season stages.However, the inherent biological characteristics of crops in different climatic zones are not considered in the FAO method.Specifically, for overwintering crops with dual peaks in their growth curves, the rapid growth stage before the overwintering period is smoothed or missed in the "four-stage" FAO method.
This distinctive phenological feature of overwintering crops has been confirmed in recent studies on the temporal downscaling of crop coefficients for winter wheat and other related crops, using both in situ measurements and soil water balance models.For example, studies investigating monthly K c of winter wheat in Yangling of Shaanxi Province in China [3] and weekly K c of rice-wheat rotation crops grown in the Indo-Gangetic plains of India [29] have all demonstrated the above-mentioned bimodal features in winter wheat.Generally, investigations on crop coefficients for winter wheat have been focusing on a particular growth stage, such as the post-elongation stage [30], the elongation to heading stage [31] and the green-up to ripening stage [32].In contrast, research on K c estimation covering the entire crop growth period from seedling emergence to the mature stage, as well as at multiple temporal scales (i.e., from the phenological stage scale to the ten-day fine scale) is rarely seen.
The North China Plain (NCP) region, regarded as the 'Granary of China', is characterized by water-intensive agricultural production with dominant cropping systems of winter wheat and summer maize rotation.Since precipitation during the winter wheat growing period is highly insufficient to meet water requirements, agricultural production in the NCP heavily relies on effective irrigation schemes, which require accurate information on ET a and K c .
As previously stated, for FAO methods, the inaccuracy in determining K c of overwintering crops can cause non-negligible errors in estimated actual evapotranspiration (ET a ) during different phenological stages.To alleviate this problem, in this study, we propose a stage-wise method to temporally downscale K c .The proposed K c method makes full use of meteorological and eddy covariance data collected at the Gucheng agro-meteorological experimental station, maintained by the Chinese Academy of Meteorological Sciences (CAMS).The performances of the proposed K c method and the conventional FAO method at both daily and phenological scales are evaluated against eddy covariance measurements of latent heat flux.The paper is organized as follows.Section 2 describes the establishment of the proposed stage-wise K c model and the required datasets, including meteorological observations, eddy covariance measurements and soil moisture datasets.Section 3 reports the temporal dynamics of ET 0 and ET a utilized for K c estimation, and the K c dynamics at different timescales including: (i) the "four-stage" scale (the same temporal scale used in the FAO method); (ii) critical phenological stages and (iii) the ten-day scale.The performance of the proposed method and the standard FAO method is assessed in comparison with surface flux measurements using a leave-one-year-out cross-validation method.

Gucheng Agro-Meteorological Experimental Station
Gucheng station is located in Hebei Province of the NCP region (115.73• E, 39.15 • N, 15.2 m above mean sea level) and includes an area of approximately 150,000 m 2 (Figure 1).The station is in a temperate continental monsoon climate zone with mean annual precipitation of 474.9 mm, about 65%-85% of which is concentrated from June to September [33].The meteorological records from 1997 to 2013 show an annual mean temperature of 12.0 • C and sunshine duration of 2158.8 h•year −1 .The dominant soil type at Gucheng station is a neutral loam, with a gentle slope and a deep arable soil layer.The station represents a high-yielding agricultural zone for winter wheat production in the NCP region, with relevant agro-meteorological variables extensively monitored in this experimental station.covariance data collected at the Gucheng agro-meteorological experimental station, maintained by the Chinese Academy of Meteorological Sciences (CAMS).The performances of the proposed Kc method and the conventional FAO method at both daily and phenological scales are evaluated against eddy covariance measurements of latent heat flux.The paper is organized as follows.Section 2 describes the establishment of the proposed stage-wise Kc model and the required datasets, including meteorological observations, eddy covariance measurements and soil moisture datasets.Section 3 reports the temporal dynamics of ET0 and ETa utilized for Kc estimation, and the Kc dynamics at different timescales including: (i) the "four-stage" scale (the same temporal scale used in the FAO method); (ii) critical phenological stages and (iii) the ten-day scale.The performance of the proposed method and the standard FAO method is assessed in comparison with surface flux measurements using a leave-one-year-out cross-validation method.

Gucheng Agro-Meteorological Experimental Station
Gucheng station is located in Hebei Province of the NCP region (115.73°E, 39.15° N, 15.2 m above mean sea level) and includes an area of approximately 150,000 m 2 (Figure 1).The station is in a temperate continental monsoon climate zone with mean annual precipitation of 474.9 mm, about 65%-85% of which is concentrated from June to September [33].The meteorological records from 1997 to 2013 show an annual mean temperature of 12.0 °C and sunshine duration of 2158.8 h•year −1 .The dominant soil type at Gucheng station is a neutral loam, with a gentle slope and a deep arable soil layer.The station represents a high-yielding agricultural zone for winter wheat production in the NCP region, with relevant agro-meteorological variables extensively monitored in this experimental station.

Reference Crop Evapotranspiration (ET0) Estimation
The reference crop evapotranspiration, ET0 (mm•day −1 ), is estimated in this study using the FAO Penman-Monteith method [1]: where Rn is the net radiation (MJ•m −2 •day −1 ) and G is the soil heat flux (MJ•m −2 •day −1 ), which is assumed to be approximately 0 at the daily time scale.T is the average air temperature at a 2-m height (°C).U2 is the average wind velocity at a 2-m height (m•s −1 ).es and ea are the saturated water vapor pressure and actual vapor pressure (kPa), respectively.Δ is the slope of the saturation vapor pressure curves (kPa•°C −1 ).γ is the psychrometric constant (kPa•°C −1 ).

Reference Crop Evapotranspiration (ET 0 ) Estimation
The reference crop evapotranspiration, ET 0 (mm•day −1 ), is estimated in this study using the FAO Penman-Monteith method [1]: where R n is the net radiation (MJ•m −2 •day −1 ) and G is the soil heat flux (MJ•m −2 •day −1 ), which is assumed to be approximately 0 at the daily time scale.T is the average air temperature at a 2-m height ( • C).U 2 is the average wind velocity at a 2-m height (m•s −1 ).e s and e a are the saturated water vapor pressure and actual vapor pressure (kPa), respectively.∆ is the slope of the saturation vapor pressure curves (kPa Solar radiation was calculated with the Ångström-Prescott formula using extraterrestrial radiation and relative sunshine duration [1]: where R s is the solar or shortwave radiation (MJ•m −2 •day −1 ) and R a is the extraterrestrial radiation (MJ•m −2 •day −1 ).N is the maximum possible duration of sunshine or daylight hours (hour), n is the actual duration of sunshine (hour), and n/N is the relative sunshine duration.The parameters a s and b s are Ångström values (dimensionless), which vary with location and time of year.In this study, both a s and b s are adopted following Mao et al. [34].

Estimation of Crop Evapotranspiration under Standard Conditions (ET c )
With proper management within the footprint of the flux tower at Gucheng station, there was sufficient supply of water and nutrition for the winter wheat crop and no pest or disease damage.Therefore, the actual evapotranspiration ET a can be considered as well-watered ET c .In this study, the ET a was acquired by two approaches: (1) direct measurements from the eddy covariance (EC) system, and (2) indirect estimation using the water balance equation.Tower-based ET a measurements collected during the growing season in a two-year study period are used to establish a crop coefficient model for winter wheat, and the estimated K c is validated using a leave-one-year-out cross-validation (LOOCV) method.Meanwhile, the ET a derived from an independent approach based on water balance equation is also used to evaluate K c for winter wheat.The details for deriving ET a from the two approaches are described in the following two subsections.
ET a Calculated Using Eddy Flux Measurements Daily latent heat flux was accumulated from half-hour data collected with the EC system and was then converted to ET a using the latent heat of water vaporization: where ET a is the daily actual evapotranspiration (mm•day −1 ), LE i is the half-hour latent heat flux (W•m −2 ), t is the time interval of sampling (1800 s), λ is the latent heat of water vaporization (2.44 MJ•kg −1 ), ρ is the density of water (kg•m −3 ) and n is number of daily samples.The flux tower measurements and energy closure assessments are described further in Section 2.3.2.

ET a Calculated by Water Balance Equation
The ET a can also be calculated using a water balance equation in the following form: where P is the precipitation (mm), I is the irrigation (mm), Q g is the contribution from the water table (mm), ET a is the actual evapotranspiration (mm), R is the surface runoff (mm) and D is the deep drainage (mm).The left-hand side of Equation ( 4) is the water budget during the study period, while ∆W s and ∆W v are the water content variations (mm) during the same period for soil and canopy, respectively.As the topography of the experimental site is relatively flat and the precipitation is not intensive during the growing period of winter wheat (annual average precipitation is about 138.7 mm during the study period), the surface runoff and the deep drainage can be neglected in this case [8,11,35,36].Contributions from groundwater are also negligible as the water table is over 30 m deep.As the canopy equivalent water thickness (EWT) of winter wheat during the entire growing period is less than 1 mm [37], the water content variations of crop (∆W v ) can be neglected.Thus, the ET a based on the water balance equation can be expressed as: where θ v1 and θ v2 are the volumetric soil water contents measured with time domain reflectometry (TDR, TRIME ® -TDR-PICO-IPH/T3, Imko, Ettlingen, Germany) at the beginning and end of the study period, respectively.

Crop Coefficient (K c ) Estimation
In this paper, the effects of both crop transpiration and soil evaporation on ET a were integrated into the single parameter of K c .The K c of a specific crop reflects the crop characteristics and its climatic features, i.e., changes in vegetation are reflected in K c variations during the growing period.As stated in the Section 1, the standard FAO method is employed in the study as a reference to benchmark the proposed stage-wise K c method.Both methods are introduced in detail in the following subsections.

Proposed Stage-Wise Method
The crop coefficient of winter wheat during each phenological period (Equation ( 6)) was estimated using the temporally-averaged reference crop evapotranspiration ET 0 (Equation ( 1)) and the measured crop evapotranspiration under standard conditions ET c (represented by ET a ): where the subscript i denotes the different phenological periods of interest.K c,i > 1 indicates that ET c,i is higher than ET 0,i , which usually occurs during the rapid growing stage; while K c,i < 1 is generally observed during the initial growing stage, when the plant height and leaf area index are low.In this study, K c calculated by the water balance equation using TDR soil moisture measurement is referred to as K c _TDR, so as to distinguish it from the K c calculated using eddy fluxes, denoted as K c _Eddy.

The FAO Method
In the FAO method, K c is considered as a function of crop growth stages, namely the initial, development, mid-season and late-season stages.Only three values are required to construct the K c curve: K c during the initial stage (K c,ini ), the mid-season stage (K c,mid ) and the late-season stage (K c,end ).According to the lookup table for single crop coefficients under the FAO method [1], typical values of K c,ini , K c,mid and K c,end for winter wheat with frozen soils are taken as 0.40, 1.15 and 0.25-0.40.At Gucheng station in the NCP region, winter wheat is harvested shortly after ripening; therefore, a high value of 0.40 is used for K c,end .As for the K c,mid , the effect of aerodynamic properties is not only crop-specific, but also depends on the climatic conditions and crop height.Thus, K c,mid can be adjusted based on the meteorological observations as below [1]: where K c,mid(tab) is the typical value of K c,mid listed in the FAO handbook.U 2 , RH min and h are the mean values for daily wind speed at a 2-m height (m•s −1 ), daily minimum relative humidity (%) and plant height (m) during the mid-season stage.After adjustment, K c,mid is 1.15 at Gucheng station, the same as the typical value for crop coefficients of winter wheat during the mid-season stage in the FAO method.

Leave-One-Year-Out Cross-Validation
Cross-validation [38] is a model validation technique for assessing how the results of statistical analysis will generalize to an independent dataset.When there is not enough data to partition into separate training and validating sets without losing significant modeling or validation capability, cross-validation provides an effective means to evaluate model performance.
Leave-one-year-out cross-validation (LOOCV) [39] involves using one year of observations as the validation set and observations from the remaining years as the training set.This is repeated n times (n is the number of observation years in the original sample) to partition the original sample into training and validation sets.

Meteorological Observations
The meteorological data used in the computation of ET 0 , including daily average air temperature, maximum air temperature, minimum air temperature, vapor pressure, atmospheric pressure, relative humidity, wind speed and sunshine duration, were collected from an automatic weather station in Gucheng during 2009-2014.All observations were quality controlled and contained no missing records during the study period.

Eddy Covariance Measurements
The eddy covariance system was installed at a height of 4 m above the ground level on a flux tower in the eastern part of Gucheng station, and winter wheat was planted in the tower fetch of up to 500 m in the prevailing wind direction.The system consists of a sonic anemometer and CO 2 /H 2 O analyzer, sampling fluxes of sensible heat, latent heat and CO 2 at a high frequency of 10 Hz, with data stored in the data logger at both 10 Hz and half-hour intervals.Additional instruments include a four-component net radiometer and a soil heat flux plate buried under the soil surface at a 5-cm depth.
Previous studies showed that the performance of coordinate rotation can be omitted as the corresponding errors were less than 1% at Gucheng station due to flat terrain [40].Additionally, these studies demonstrated that half-hour is an acceptable interval for averaging the high-frequency flux observations obtained from the height below 30-m level [40][41][42][43].Meanwhile, the 30-min averaging period is commonly used by the micro-meteorological community for flux measurements over agricultural sites [14,15].Thus, in this study, the 30-min EC data we used were processed from the original high-frequency observations after a series of processing steps, including spike removal, WPL (Webb, Pearman, and Leuning) correction, sonic virtual temperature correction and frequency response correction.In addition, quality control was performed to minimize the effects of extreme weather conditions and equipment malfunction on the flux data.Specifically, the obvious outliers in the time series and invalid nocturnal measurements taken with friction velocity <0.08 m•s −1 were excluded [44].Following this, the missing data for small gaps (2-3 half-hourly missing) were replaced by linear interpolation using the adjacent value(s) [45], whereas large data gaps were filled using exponential regression equations [46] developed for each individual stage, namely the seedling growth stage, dormancy stage and rapid growth stage.To avoid inaccurate interpolation, days with continuous data gaps exceeding 1/3 of the daily record length, resulting from either missing measurements or quality control, were completely removed during model construction and validation.After these treatments, two years with reliable flux observations, i.e., the winter wheat growing periods for 2009-2010 and 2012-2013 (denoted as the study periods of 2009 and 2012 hereinafter), were included in the analysis.
The typical energy closure condition in the EC measurements was evaluated for each phenological stage of the two-year study period using the energy balance ratio (EBR) [47].EBR < 100% indicates that the measurements do not fully characterize the local flux conditions, with errors typically attributed to the underestimation of sensible and latent fluxes and/or the overestimation of available energy during nocturnal periods, especially weak turbulent periods [48].Closure ratios of 80% or higher are considered to be reasonable at most sites with climate zones ranging from Mediterranean to temperate and arctic [48].

Soil Water Content Measurements
Soil water content was measured by TDR once every 10 days (i.e., on the 8th, 18th and 28th days of each month) during the growing period of winter wheat for 2013 and 2014.According to ground measurements in the southern NCP region, the mean of total root length density (RLD) for winter wheat is 57.7% from surface down to a 0.5 m depth, and the value is 81.1% from surface down to 1.0 m [49].Therefore, in this study, soil moisture was measured at 0.1 m intervals from the surface down to a 1.0 m depth to account for the change of soil water content within the root zone.The soil moisture was recorded in percentage units (%) with an accuracy of ±4% according to instrument manual.

Remotely Sensed Leaf Area Index (LAI) Time Series
The Moderate Resolution Imaging Spectroradiometer (MODIS) 8-day LAI product (MOD15A2, https://reverb.echo.nasa.gov/reverb/)at 1-km spatial resolution during the winter wheat growing season from two years (2009 and 2012) was used to indirectly evaluate the K c estimation, as temporal dynamics of LAI and K c have proven to be highly correlated [3].The time series of MODIS LAI in the pixel collocated with Gucheng station during the two-year study period was extracted, and the QC flags were used to exclude unreliable observations, retaining data retrieved by the main retrieval method.Consequently, a linear interpolation method was used to construct a complete 8-day time series of LAI and, finally, re-sampled at a ten-day interval so as to be consistent with K c at the temporal scale.

Meteorological Conditions during the Winter Wheat Growing Period
The winter wheat used in this experiment was a medium-early maturing variety.The detailed growing conditions for each winter wheat variety, including growing dates and field management, during the two-year study period are listed in Table 1.For the two-year study period, the average temperature, average maximum temperature and average minimum temperature were 5.9 • C, 12.1 • C and 0.5 • C, respectively, with little inter-annual variation during the winter wheat growing period (Table 1).Total precipitations were 125.5 mm and 193.8 mm, which temporally coincide with the dry season.The percentages of rainy days occupied 17.4% and 12.8% for the two growing seasons respectively, and flood irrigations were applied 4-5 times for each growing season.On average, water supply throughout the growing periods was generally regulated by weekly rainfall and/or irrigation.Therefore, precipitation has little constraining impact on winter wheat growth and, thus, is excluded from the discussion in this paper.

Energy Balance Closure of Flux Measurements
A scatter plot comparing the surface fluxes (LE + H) and the available energy (R n −G) measured at the flux tower during all critical phenological stages is shown in Figure 2. The EBRs for the two growing periods range from 90.39% to 91.35%.The mean EBRs for the seedling, dormancy and rapid-growth periods are 94.14%, 84.83% and 91.58%, respectively, all higher than the mean EBR of 80% for a diverse range of vegetation covers and climate types [48].The good closure conditions observed in Figure 2 indicate that the ET a measurements used to compute K c _Eddy are reasonable.

ET0 and ETa Dynamics
Temporal dynamics of ET0 calculated using the FAO Penman-Monteith method and ETa obtained from eddy covariance observations at Gucheng station during growing seasons are shown in Figure 3.These curves are computed on 10-day intervals and represent averages over the two-year study period.The vertical lines in Figure 3 represent the nine generalized growth stages (S1-S9) for winter wheat, as tabulated in Table 2.

ET 0 and ET a Dynamics
Temporal dynamics of ET 0 calculated using the FAO Penman-Monteith method and ET a obtained from eddy covariance observations at Gucheng station during growing seasons are shown in Figure 3.These curves are computed on 10-day intervals and represent averages over the two-year study period.The vertical lines in Figure 3 represent the nine generalized growth stages (S1-S9) for winter wheat, as tabulated in Table 2.

ET0 and ETa Dynamics
Temporal dynamics of ET0 calculated using the FAO Penman-Monteith method and ETa obtained from eddy covariance observations at Gucheng station during growing seasons are shown in Figure 3.These curves are computed on 10-day intervals and represent averages over the two-year study period.The vertical lines in Figure 3 represent the nine generalized growth stages (S1-S9) for winter wheat, as tabulated in Table 2.In the early growth stages from seedling to green-up (S1-S4), both ET 0 and ET a are low due to low insolation levels and low air temperature during fall to winter.Starting from the S4 stage, the rising temperatures and solar radiation have resulted in increasing ET 0 and ET a .As the canopy starts to close during the erect-growing to blossom stages (S5-S8), ET a approaches the potential rate.Both ET (ET 0 and ET a ) curves show the highest degree of inter-annual variability during S5-S8.In particular, the inter-annual variations of ET a mainly occurred during the rapid growing periods, i.e., late March-early May, due to variability in crop canopy height, LAI and leaf morphology between years.ET a downturn occurs in the milking to ripening stage (S9), as the green canopy covers begin to decline.

Comparison with the Four-Stage FAO Method
In order to directly compare the capabilities of characterizing winter wheat K c between the proposed stage-wise methods and the "four-stage" FAO method, we regrouped the nine finer-scale phenological stages of winter wheat in our method to match the "four-stage" in the FAO method.Specifically, the initial stage in the FAO method corresponds to the seeding to green-up stages; the crop development stage and mid-season stage cover the erect to elongation stages and booting to blossom stages, respectively; and the late-season stage corresponds to the milk to harvest stages (Table 2).
Based on this regrouping scheme, K c time series derived from both the FAO method (referred to as K c _FAO) and the proposed stage-wise method (based on two years of quality-controlled eddy covariance fluxes (see Section 2.3.2),referred to as K c_ Eddy) are illustrated in Figure 4, showing noticeable differences between two sets of K c values.In particular, K c_ Eddy exceeds K c _FAO in the initial and late-season stage and is lower than K c _FAO during the mid-season stage.This is mainly attributed to the difference in climatic features between the sub-humid region (where the FAO method was originally developed) and this study region in the NCP.Specifically, the FAO method for winter wheat K c was originally developed for sub-humid regions with average daily minimum relative humidity (RH min ) around 45% and moderate wind speeds of 2 m•s −1 .However, meteorological conditions at Gucheng station are quite different, representing semi-arid and sub-humid climates with sparse rainfall during the winter wheat growing period.In addition, the average RH min for the two-year study period was about 31.6%, well below the expected 45% that generated the typical K c _FAO values.
As soil moisture is not a limiting factor for evapotranspiration with regard to the experimental design at the Gucheng site (which provides sufficient irrigation, as mentioned in Section 2.3.5), the relatively large atmospheric demand enhance the potential soil evaporation and crop transpiration, and therefore, increased ET a , resulting in higher K c_ Eddy than K c _FAO during the initial stage.Differences in K c values during the mid-season stage may be attributed to lower crop heights at Gucheng station (0.7 m) in comparison with the expected maximum height of winter wheat in the FAO method (1.0 m).Since crop transpiration accounted for considerable portion of ET a during the mid-season stage, the relatively shorter crop height at Gucheng may be responsible for the lower K c _Eddy than K c _FAO.As for the late-season stage, winter wheat at Gucheng station was harvested immediately after reaching mature stage, so as to prepare for the sowing of second crop-rotation (maize).Therefore, K c _Eddy of winter wheat at Gucheng station for the late-season stage was significantly higher than K c _FAO before ET a dropped to the minimum value.It is worth noting that due to the prevailing maize-wheat rotation system in NCP, this common agricultural practice of immediate harvesting of winter wheat after the mature stage is a human-induced factor that contributes to the unsuitability of the FAO method across this region.
Water 2017, 9, 155 10 of 15 Differences in Kc values during the mid-season stage may be attributed to lower crop heights at Gucheng station (0.7 m) in comparison with the expected maximum height of winter wheat in the FAO method (1.0 m).Since crop transpiration accounted for considerable portion of ETa during the mid-season stage, the relatively shorter crop height at Gucheng may be responsible for the lower Kc_Eddy than Kc_FAO.As for the late-season stage, winter wheat at Gucheng station was harvested immediately after reaching mature stage, so as to prepare for the sowing of second crop-rotation (maize).Therefore, Kc_Eddy of winter wheat at Gucheng station for the late-season stage was significantly higher than Kc_FAO before ETa dropped to the minimum value.It is worth noting that due to the prevailing maize-wheat rotation system in NCP, this common agricultural practice of immediate harvesting of winter wheat after the mature stage is a human-induced factor that contributes to the unsuitability of the FAO method across this region.

Kc Dynamics at the Critical Phenological Stages
To examine Kc dynamics at Gucheng station at finer temporal scales, Kc values for winter wheat developed using ETa estimates from eddy covariance (Kc_Eddy) and water balance (Kc_TDR) are compared at the time scale of critical phenological stages in Figure 5.

K c Dynamics at the Critical Phenological Stages
To examine K c dynamics at Gucheng station at finer temporal scales, K c values for winter wheat developed using ET a estimates from eddy covariance (K c _Eddy) and water balance (K c _TDR) are compared at the time scale of critical phenological stages in Figure 5.
Water 2017, 9, 155 10 of 15 Differences in Kc values during the mid-season stage may be attributed to lower crop heights at Gucheng station (0.7 m) in comparison with the expected maximum height of winter wheat in the FAO method (1.0 m).Since crop transpiration accounted for considerable portion of ETa during the mid-season stage, the relatively shorter crop height at Gucheng may be responsible for the lower Kc_Eddy than Kc_FAO.As for the late-season stage, winter wheat at Gucheng station was harvested immediately after reaching mature stage, so as to prepare for the sowing of second crop-rotation (maize).Therefore, Kc_Eddy of winter wheat at Gucheng station for the late-season stage was significantly higher than Kc_FAO before ETa dropped to the minimum value.It is worth noting that due to the prevailing maize-wheat rotation system in NCP, this common agricultural practice of immediate harvesting of winter wheat after the mature stage is a human-induced factor that contributes to the unsuitability of the FAO method across this region.

Kc Dynamics at the Critical Phenological Stages
To examine Kc dynamics at Gucheng station at finer temporal scales, Kc values for winter wheat developed using ETa estimates from eddy covariance (Kc_Eddy) and water balance (Kc_TDR) are compared at the time scale of critical phenological stages in Figure 5.As seen in Figure 5, the temporal dynamics of both K c _Eddy and K c _TDR exhibit similar bimodal features during the growing period.The winter wheat development cycle is characterized by two rapid growing periods, namely the three-leaf to dormancy stages and the erect-growing to booting stages, which occur before and after overwintering, respectively.This is similar to other studies of K c features of winter wheat using both in situ observations [3] and models [50,51].It is noted that due to the neglect of the rapid growing period before overwintering in the FAO method, the actual water consumption prior to the green-up stage was considerably underestimated by approximately 9 mm, accounting for 10.5% of actual water consumption during the initial stage.The average K c _Eddy dynamics at Gucheng station are compared with that of LAI observations at the 10-day scale in Figure 6.Overall, the temporal dynamics of K c and LAI at the 10-day scale were highly correlated as expected.Specifically, at this 10-day time scale, the two peaks in the K c _Eddy curve from late November to early December and from mid-April to late May are very evident, as well as a clear overwintering feature.During the overwintering to erect-growing stages, i.e., from mid-December to the next mid-March, K c _Eddy stayed at a low level and experienced little variability.Starting in March, K c _Eddy increases significantly as the winter wheat starts rapid growth, with the main peak around mid-April at the booting stage.K c _Eddy remains above 1.0, signifying high rates of water consumption until early June when winter wheat approached the mature stage.As seen in Figure 5, the temporal dynamics of both Kc_Eddy and Kc_TDR exhibit similar bimodal features during the growing period.The winter wheat development cycle is characterized by two rapid growing periods, namely the three-leaf to dormancy stages and the erect-growing to booting stages, which occur before and after overwintering, respectively.This is similar to other studies of Kc features of winter wheat using both in situ observations [3] and models [50,51].It is noted that due to the neglect of the rapid growing period before overwintering in the FAO method, the actual water consumption prior to the green-up stage was considerably underestimated by approximately 9 mm, accounting for 10.5% of actual water consumption during the initial stage.

Kc Dynamics at the Ten-Day Scale
The average Kc_Eddy dynamics at Gucheng station are compared with that of LAI observations at the 10-day scale in Figure 6.Overall, the temporal dynamics of Kc and LAI at the 10-day scale were highly correlated as expected.Specifically, at this 10-day time scale, the two peaks in the Kc_Eddy curve from late November to early December and from mid-April to late May are very evident, as well as a clear overwintering feature.During the overwintering to erect-growing stages, i.e., from mid-December to the next mid-March, Kc_Eddy stayed at a low level and experienced little variability.Starting in March, Kc_Eddy increases significantly as the winter wheat starts rapid growth, with the main peak around mid-April at the booting stage.Kc_Eddy remains above 1.0, signifying high rates of water consumption until early June when winter wheat approached the mature stage.

LOOCV of the Estimated ETa Using the Proposed Stage-Wise Method
LOOCV was used to evaluate the performance of the stage-wise Kc method proposed in this study.The daily ETa during the entire growing period estimated from the proposed stage-wise Kc method (Kc_Eddy) and the FAO method (Kc_FAO) is validated against measurements from the eddy covariance system at Gucheng station.The comparison statistics at both daily and critical phenological scales with the LOOCV technique are summarized in Tables 3 and 4.
Generally, the proposed stage-wise method outperforms the FAO method in terms of linear regression coefficients k and root-mean-square errors (RMSEs).It is shown that the slope of the linear regression line (against ground measurements) for the stage-wise method is more approximate to unity compared to the FAO method (0.97 versus 1.08 at the daily scale and 1.00 versus 1.09 at the critical phenological scale).In addition, the outperformance of the stage-wise method over the FAO

LOOCV of the Estimated ET a Using the Proposed Stage-Wise Method
LOOCV was used to evaluate the performance of the stage-wise K c method proposed in this study.The daily ET a during the entire growing period estimated from the proposed stage-wise K c method (K c _Eddy) and the FAO method (K c _FAO) is validated against measurements from the eddy covariance system at Gucheng station.The comparison statistics at both daily and critical phenological scales with the LOOCV technique are summarized in Tables 3 and 4.
Generally, the proposed stage-wise method outperforms the FAO method in terms of linear regression coefficients k and root-mean-square errors (RMSEs).It is shown that the slope of the linear regression line (against ground measurements) for the stage-wise method is more approximate to unity compared to the FAO method (0.97 versus 1.08 at the daily scale and 1.00 versus 1.09 at the critical phenological scale).In addition, the outperformance of the stage-wise method over the FAO method is also supported by the metrics of RMSEs (0.07 mm•day −1 versus 0.16 mm•day −1 at the daily scale and 0.01 mm•day −1 versus 0.27 mm•day −1 at the critical phenological scale).In summary, our method showed slightly better performance than the FAO method in ET a estimation during the growing period.The errors in estimated ET a for both the stage-wise and the FAO method at critical phenological stages are shown in Figure 7. From Figure 7, we can see that the uncertainties of ET a from the stage-wise method are much smaller than those from the FAO method.Generally, during the initial stages, from seeding to the erect-growing stage, the ET a derived from the FAO method is underestimated by −0.2 mm•day −1 , higher than that from the stage-wise method (−0.004 mm•day −1 ).For the development stage and mid-season stage, the ET a from the FAO method is notably overestimated by 0.63 mm•day −1 and 0.64 mm•day −1 , which are higher than the overestimation from the stage-wise method (−0.01 mm•day −1 and 0.04 mm•day −1 ).The underestimations of ET a in S9 are −0.02mm•day −1 and −0.43 mm•day −1 from the stage-wise and FAO methods, respectively.Overall, a slight superiority of the proposed stage-wise method over the FAO method is exhibited for estimating actual evapotranspiration at phenological time scales.method is also supported by the metrics of RMSEs (0.07 mm•day −1 versus 0.16 mm•day −1 at the daily scale and 0.01 mm•day −1 versus 0.27 mm•day −1 at the critical phenological scale).In summary, our method showed slightly better performance than the FAO method in ETa estimation during the growing period.The errors in estimated ETa for both the stage-wise and the FAO method at critical phenological stages are shown in Figure 7. From Figure 7, we can see that the uncertainties of ETa from the stage-wise method are much smaller than those from the FAO method.Generally, during the initial stages, from seeding to the erect-growing stage, the ETa derived from the FAO method is underestimated by −0.2 mm•day −1 , higher than that from the stage-wise method (−0.004 mm•day −1 ).For the development stage and mid-season stage, the ETa from the FAO method is notably overestimated by 0.63 mm•day −1 and 0.64 mm•day −1 , which are higher than the overestimation from the stage-wise method (−0.01 mm•day −1 and 0.04 mm•day −1 ).The underestimations of ETa in S9 are −0.02mm•day −1 and −0.43 mm•day −1 from the stage-wise and FAO methods, respectively.Overall, a slight superiority of the proposed stage-wise method over the FAO method is exhibited for estimating actual evapotranspiration at phenological time scales.

Conclusions
The crop coefficient (K c ) is a crucial parameter in the estimation of actual water consumption and water requirement for agricultural development.In this study, a method to estimate temporal fine-scale K c of winter wheat was developed using the meteorological and eddy covariance observations, as well as the winter wheat growing records during different phenological stages.This downscaling model is motivated by the significant bimodal feature in K c time series of winter wheat and its distinctive overwintering characteristic, which are not addressed in the conventional FAO method.
Results show that compared to the "four-stage" FAO method, our proposed stage-wise method accurately captures the K c peak before the overwintering stage and the trough during overwintering.In addition, our method can estimate K c at temporally fine resolution at the crop development and mid-season stages and, thus, can serve as the basis for the accurate estimation of actual water consumption during the critical phenological stages of winter wheat.
The uncertainties in ET a based on K c derived from our proposed stage-wise method and the conventional FAO method are estimated during all critical phenological stages.Compared to the stage-wise method, the overestimation of ET a from the FAO method is obvious with the maximum value of 0.77 mm•day −1 during all of the critical phenological stages.The RMSEs of ET a for the stage-wise and FAO methods are 0.07 mm•day −1 and 0.16 mm•day −1 at the daily scale and 0.01 mm•day −1 and 0.27 mm•day −1 at the critical phenological scale, respectively.
As Gucheng station represents a typical high-yield agricultural zone in the NCP region, the winter wheat K c at the critical phenological stages generated from our method can be used as a reference for the entire NCP region, where in situ measurements of K c are relatively scarce.However, accurate and spatially-distributed K c at the crop phenological stages is crucial for precisely estimating crop water requirements, improving crop water use efficiency and yield, and consequently achieving the goal of highly-promoted water-saving irrigation.Further investigation on the spatial characteristics of winter wheat K c at the phenological stages at the regional scale is needed.Meanwhile, the validation of the proposed framework in a larger observation sample is crucial in the future.

Figure 1 .
Figure 1.Location of the Gucheng agro-meteorological station in the North China Plain (NCP) region of China (the extent of the station is outlined in yellow line in the rightmost sub-figure).

Figure 1 .
Figure 1.Location of the Gucheng agro-meteorological station in the North China Plain (NCP) region of China (the extent of the station is outlined in yellow line in the rightmost sub-figure).

Figure 2 .
Figure 2. Energy balance closure during all critical phenological stages of winter wheat in two-year study period at Gucheng station.

Figure 3 .
Figure 3. ET0 and ETa for winter wheat at 10-day intervals during the growing period averaged over the two-year study period.The error bars show the standard deviation in ET0 and ETa for each interval over the two years.Vertical lines indicate nine generalized growth stages (S1-S9) for winter wheat.

Figure 2 .
Figure 2. Energy balance closure during all critical phenological stages of winter wheat in two-year study period at Gucheng station.

Figure 2 .
Figure 2. Energy balance closure during all critical phenological stages of winter wheat in two-year study period at Gucheng station.

Figure 3 .
Figure 3. ET0 and ETa for winter wheat at 10-day intervals during the growing period averaged over the two-year study period.The error bars show the standard deviation in ET0 and ETa for each interval over the two years.Vertical lines indicate nine generalized growth stages (S1-S9) for winter wheat.

Figure 3 .
Figure 3. ET 0 and ET a for winter wheat at 10-day intervals during the growing period averaged over the two-year study period.The error bars show the standard deviation in ET 0 and ET a for each interval over the two years.Vertical lines indicate nine generalized growth stages (S1-S9) for winter wheat.

Figure 4 .
Figure 4. Comparison of Kc between the proposed stage-wise method and the FAO method.

Figure 5 .
Figure 5.The proposed stage-wise Kc during the critical phenological stages.

Figure 4 .
Figure 4. Comparison of K c between the proposed stage-wise method and the FAO method.

Figure 4 .
Figure 4. Comparison of Kc between the proposed stage-wise method and the FAO method.

Figure 5 .
Figure 5.The proposed stage-wise Kc during the critical phenological stages.

Figure 5 .
Figure 5.The proposed stage-wise K c during the critical phenological stages.

3. 3 . 3 .
K c Dynamics at the Ten-Day Scale

Figure 7 .
Figure 7.Comparison of delta ETa for winter wheat between the stage-wise method and the FAO method during the critical phenological stages.

Figure 7 .
Figure 7.Comparison of delta ET a for winter wheat between the stage-wise method and the FAO method during the critical phenological stages.

Table 1 .
Winter wheat varieties and meteorological conditions at Gucheng station during the growing seasons of the two-year study period.

Table 2 .
The growing stages of winter wheat and its corresponding phenological periods according to the FAO method.

Table 3 .
Statistical metrics of the comparison between the stage-wise and the FAO estimations of ET a with the LOOCV method at the daily time scale.

Table 4 .
Statistical metrics of the comparison between the stage-wise and the FAO estimations of ET a with the LOOCV method at the critical phenological time scale.
Note: * k is the regression coefficient with the linear equation "y = kx".

Table 3 .
Statistical metrics of the comparison between the stage-wise and the FAO estimations of ETa with the LOOCV method at the daily time scale.

Table 4 .
Statistical metrics of the comparison between the stage-wise and the FAO estimations of ETa with the LOOCV method at the critical phenological time scale.
Note: * k is the regression coefficient with the linear equation "y = kx".