Winter Wheat Water Productivity Evaluated by the Developed Remote Sensing Evapotranspiration Model in Hebei Plain, China

Agricultural water is the main reason for the rapid decline of the NCP groundwater levels. It is of vital importance for the NCP sustainable agricultural development to master the ETa and its CWP. In this paper, the EBEM model was developed according to the theory of energy balance. From 2001 to 2006, the winter wheat ETa and CWP were estimated, and the spatial and temporal variations and their influencing factors were studied in the Hebei Plain. The results indicate that the EBEM model performed well by applying MODIS data to estimate the daily net radiation and ETa. For the daytime net radiation, the relative error between the estimation and the measurement amounted to 8.2% and the SEE was 0.82 MJ m−2/day. The average ETa deviation between the estimates and the measures amounted to 0.86 mm daily, and the SEE was 1.2 mm. The spatial variations indicated that the major distribution of ETa ranged from 350 to 450 mm, which trended downward within the study area from west to east. In the study period, the winter wheat CWP was mainly distributed between 0.29 and 1.67 kg/m3. In space, the CWP was higher in the west than in the east.


Introduction
With the worldwide population growth and economic development, the global water demand and the degree of water stress are both increasing. Thus, determining how to acquire a better management of water resources and to achieve sustainable development becomes an important issue that every country's government must consider and solve. Nearly a decade ago, 87% of global water extraction (extraction minus return) was consumed by irrigation. Meanwhile, worldwide food production of those irrigated fields accounted for 40-45% [1]. Therefore, reasonable estimates and scientific management of irrigation water are not only important for the sustainable use of water resources but also related to the issue of food security [2]. However, the scientific management of irrigation depends on the accurate estimation of both crop water productivity (CWP) and crop water consumption (CWC) and can be expressed as active evapotranspiration (ETa) [3,4]. In the regional perspective, remote sensing refers to one of the suitable means for reasonable estimation of the ETa. In the past few decades, remote sensing-based regional ETa estimation has been used to plan and manage water irrigation with high accuracy and precision [5,6]. Among all of the remote sensing data being used to estimate ETa, the Moderate Resolution Imaging Spectroradiometer (MODIS) is more suitable for medium-to-large regional-scale application because of its advantage of high time resolution, which means medium spatial resolution, high spectral resolution, and free access to the information [7][8][9], even in a data scarce and heterogeneous landscape [10].
The North China Plain (NCP), an area of 320,000 km 2 with a population of over 300 million, is one of the most  important grain-producing areas in China. The arable land area of the NCP amounts to 1.2 million hectares, accounting for 18.3% of the national total, and its food production provides 21.6% of China's national total. In this region, the main planting agency is the winter wheat-summer corn rotation each year, while the growing season of the winter wheat is from early October to the end of May in the following year. While the NCP is in the monsoon climate with very little rainfall, the multiyear average rainfall  is only 119 mm. The amount of precipitation is even smaller than one-third of the requirement for achieving the average yield [11].
Therefore, irrigation is necessary to ensure adequate food production. Generally, there are three to five time periods of irrigation for a total of 300-400 mm of water per year. As a result, irrigation water accounts for approximately 70% of the total water consumption of the local area, while winter wheat consumes almost 70% of the total amount of irrigation water. Therefore, approximately 50% of the water is consumed by winter wheat [12].
The relationship between ETa and CWP was primarily studied by adjusting the irrigation methods and the irrigation amount in the NCP [13,14]. Although these experiments achieved meaningful results, the scale of the field experiments cannot reflect the spatial variability and is not sufficiently representative. The GIS-based model was also used to research the CWP variation in the NCP [15]. However, the crop water consumption in that research study was obtained by indirect calculation of electricity consumption. Mo et al. simulated the ETa and CWP temporal and spatial variations in the SVAT model of the NCP. This model requires a number of input data, such as DEM, soil texture, and meteorological data, and, thus, the additional required data and parameters are very inconvenient for utilization of the model [16]. Li et al. estimated the distribution of the ETa and CWP by applying AVHRR data input to the SEBAL model in the year of 2003 [12]. However, the model did not adequately describe either the space representation or the interannual variation. Tian et al. estimated the ETa from 2001 to 2009 but the CWP was not considered [17]. Thus, it is necessary to perform further study on the temporal variation and the spatial distribution of ETa and CWP to provide a reference for agricultural water management and sustainable regional development.
Although the statistical data at the spatial resolution is less than the yield estimation based on the remote sensing data through the model to determine production, the county scale represents the largest scale of production data available in the NCP [15] and is also the basic unit of the agricultural water use and the water resources management. Therefore, in this research, the county scale was selected, and the ETa and CWP of winter wheat were estimated based on the remote sensing data in the NCP. The ETa of the winter wheat planting area was estimated using an energy balance-ET model called the energy balance and evapotranspiration model (EBEM). Through the application of the remote sensing data, the winter wheat areas were classified to determine the counties' average ETa, and the CWP was defined as the ETa divided by the yield. Therefore, the factors influencing the CWP and ETa of winter wheat were analyzed and their temporal and spatial variations were investigated.

Study Area.
In this study, the Hebei Plain (HP), located in the northern part of the NCP, was chosen to represent the NCP. A total of 80 counties with an elevation of below 100 m were included in this area ( Figure 1). The local climate is a typical temperate semiarid one with a strong summer monsoon season during a hot and rainy summer and a dry and cold winter. The average monthly temperatures range from −2.5 ∘ C to 26 ∘ C, and the average annual rainfall is approximately 400 to 500 mm, with approximately 70% of this rainfall occurring Regarding the scope of this study, a total of four groups of data were required to cover the entire area, namely, "h26v04, " "h26v05, " "h27v04, " and "h27v05. " Such data were processed by splicing, projection conversion, and denoising through the use of the quality control data and the BRDF correction using the parameters of the BRDF.

Meteorological and Ground Observation Data.
In this paper, the meteorological data include the following values: daily maximum and minimum temperatures ( ∘ C), wind speed (m s −1 ), relative humidity (%), and number of sunshine hours (h), which are obtained from the China Meteorological Data Sharing System (http://cdc.cma.gov.cn/), with the site distribution shown in Figure 1. The data were interpolated to the entire study area by using the Kriging interpolation methods. In addition, the data provided by the radiation observing system and large weighing Lysimeter (Lysimeter) of the Luancheng Agro-Ecosystem Experimental Station of the Chinese Academy of Sciences (Luancheng Station) were provided for the model in this paper. The Lysimeter, with an area of 1.5 m × 2 m, a depth of 2.5 m, an empty weight of two tons and a weight of 14 tons when filled with soil, has an accuracy of up to 0.02 mm [18] where the crops, irrigation, and fertilization are the same as the surrounding farmland. In this region, the main crop type is the winter wheat-summer corn rotation each year.
In this research, the production statistics, including the acreage and winter wheat production from 2000 to 2006 (Hebei Provincial People's Government Office and the Bureau of Hebei Province, 2001-2007), were obtained from the Hebei Rural Statistical Yearbook.

EBEM Theory.
The EBEM calculates the latent heat flux by taking advantage of the energy balance principle, as shown in the following equation: where LE represents the latent heat flux, refers to net radiation reaching the earth's surface, represents the soil heat flux, and refers to the sensible heat flux, with the units in W m −2 .

Instantaneous Net Radiation.
Latent heat flux can be obtained by determining each component of the right-hand side of (1). However, we can only obtain the instantaneous surface information using remote sensing data, as shown in the following equation: where ↓ represents the incident shortwave radiation, refers to the surface albedo, ↓ represents the incident longwave radiation, ↑ refers to the outgoing longwave radiation, and denotes the wide band surface thermal emission rate. The shortwave radiation reaching the surface is calculated as follows: where sc represents the solar constant (1367 W m −2 ), rel refers to the sun incident angle, sw represents the atmospheric transmittance, and 2 refers to an average distance squared for one day [19]. The atmospheric transmittance sw is obtained using the air pressure (kPa), (mm) of the water vapor content in the atmosphere, the solar elevation angle, and turbidity coefficient , as shown in (4), where 0 ≤ ≤ 1.0; when air quality is better, is equal to 1; and when the air is extremely turbid or seriously polluted, is equal to 0.5 [19]: Allen et al. [19] calculated by using the near-surface water vapor pressure. However, acquiring the near-surface water vapor pressure was limited to the ground station. Therefore, we obtained the water vapor pressure by using the MOD08 E3 of MODIS. The surface albedo was calculated using the narrowband reflectance estimation method, as shown in the following equation: where refers to each band narrowband reflectance [20]. The outing longwave radiation was calculated using (6), where represents the surface thermal emission rate and ( ) refers to the surface temperature. Both and are obtained by the MOD11 product, and is the Stefan-Boltzmann constant (5.67 * 10 −8 W m −2 K −4 ): 4 The Scientific World Journal The incident longwave radiation was calculated using (7), where is the air thermal emission rate, which can be obtained using (8). refers to the air temperature: Air temperature can be obtained by the surface temperature of the vegetation coverage concept space, as shown in Figure 2. The space evolved by the surface temperature-vegetation index ( -VI) space vegetation fraction ( ) can better reflect the relationship between vegetation and soil [21] from the pixel scale because the NDVI or the EVI only represent the green vegetation distribution of the surface.

Calculating the Instantaneous Soil Heat Flux and the Sensible Heat Flux.
After obtaining the instantaneous , both and must be obtained. To achieve instantaneous LE, was obtained by taking advantage of its relationship with the ratio with the other surface parameters. By using the Luancheng Station observation data and the remote sensing observations of surface temperature, was calculated according to the proposed equation where ( ) refers to the surface temperature of the MODIS observations.
The transient thermal flux and the estimation of the relevant parameters, using the procedures of Bastiaanssen et al. [22] and Allen et al. [19], are obtained using an iterative process involving the aerodynamic equation and the Monin-Obukhov similarity theory. was calculated according to the following equation: where air represents the air density (kg m −3 ); refers to the specific heat of air (J kg −1 K −1 ); ah represents the aerodynamic drag (s m −1 ) near the ground between the heights 1 and 2 ; 1 represents 0.1 m and 2 refers to 2 m; and denotes the temperature gradient between 1 and 2 . Bastiaanssen et al. have proposed the use of a linear relationship between and , as described in the following equation, to calculate [22]: where and refer to empirical parameters, which can be obtained by estimating the area under each ETa regression of the remote sensing images, which can be obtained from at least two points (hot and cold spots). Under normal circumstances, by artificially selecting good vegetation growth, the adequate water supply point is selected as the cold spot, while the area without vegetation and relatively dry is selected as the hot spot. However, such a method adds much uncertainty. In this paper, to avoid this uncertainty, the -space automatically obtains the cold and hot spots, which can make the model more appropriate for automated batch computing. As shown in Figure 2, the hot point defined as is very small (close to 0) and has a high surface temperature, while is close to 1 and the surface temperature is the lowest pixel point, which is defined as the cold spot.

Whole Day ETa.
Because the value of all-day evapotranspiration is more meaningful than the instantaneous value [19], the daily ET (ET day ) value was obtained using (12)- (14): where EF represents the evaporation ratio and (J kg −1 ) refers to the latent heat of water vaporization at the specific temperature, which can be obtained from the following equations, respectively: Samani et al. assumed that the instantaneous net radiation ( ) and the daytime net radiation (positive net radiation day ) with instantaneous shortwave radiation ( ) are proportional to the daily shortwave radiation ( 24 ) [23] and proposed the following equation to calculate the all-day net radiation day : The Scientific World Journal 5 where day ( ) represents the average temperature for the entire day; ( ) denotes the instantaneous temperature; refers to the instantaneous shortwave radiation by using (3); and day represents the shortwave radiation for the entire day by the following equation: where max represents the maximum temperature ( ∘ C), min represents the minimum temperature ( ∘ C), and refers to the adjustment factor. For the large coastal areas of air masses approaching water influence, was set to 0.19. For the air mass in the inland areas, which are relatively less affected by water, was set to 0.16 [24]. refers to the atmospheric outer radiation (extraterrestrial radiation) in the following equation: where sc is the solar constant (0.0820 MJ m −2 min −1 ); denotes an average distance countdown for a day; refers to the sunset hour angle (rad); marks the latitude (rad); and represents the real declination (rad). Given the fluctuations in soil heat flux ( day ) within a full day, Bastiaanssen et al. considered that it can be ignored [22]:

ETa Calculation in a Certain
Period. The seasonal or annual total ETa is significant for irrigation and water resources management, which can be calculated using the following equation: where ET period refers to the cumulative evapotranspiration in a certain period, which may be from the beginning of one day to the end of the other day. ET represents the transpiration ratio for this period; ET day denotes the daily reference evapotranspiration, which can be calculated using the meteorological data. The data input in the remote sensing model represents a synthetic product of 8 days and represents the best weather conditions of the 8 days. Therefore, ET was calculated according to the daily evapotranspiration obtained using the remote sensing data and the highest reference evapotranspiration in this period, as shown in the following equation: where ET max period refers to the maximum reference evapotranspiration of the synthesis period of 8 days.

ETa Calculation at the County Level.
Because there are also many other crops, natural vegetation, and nonvegetation in the study area, a classification method based on time series NDVI data was used for determining the actual winter wheat area so that the water consumption by winter wheat can be accurately estimated [25].

CWP Calculation.
In the study area, the crop water productivity of the winter wheat was calculated according to [26] where CWP (kg m −3 ) refers to the crop water productivity; GY (kg ha −1 ) represents the crop yield; and ET (mm) denotes the sum of the entire growing season of ETa from October 1 to June 15 in the following year.

Verification of Daily Net
Radiation. LE refers to the energy, as described by the energy balance equation, which constitutes the source of energy occurring during evapotranspiration. Accordingly, it is very important to accurately estimate the daily net radiation ( day ) for obtaining accurate LE values. In this paper, based on (15), the all-day net radiation was calculated by the instantaneous net radiation. The method was proposed by [23] and verified. However, the method has differences among the study area, vegetation type, and the focus of the research, and it is necessary to evaluate the applicability of the method in the region. In Figure 3, the 91 observed data ranges from 2.37 to 18.19 MJ m −2 day −1 were derived from the Luancheng Station and compared with the estimated value based on the remote sensing data. The average ratio was 0.976, the SEE amounted to 0.82 MJ m −2 day −1 , and the relative error was 6.7%, which was between the estimates and the measured values.

ETa Validation.
To evaluate the ETa using the models, the observed evapotranspiration value, which was obtained by the Lysimeter located in the Luancheng Station, was compared with the ETa calculated using the MODIS data. The ground observation data have been called into question because of the differences among the wide ranges of the existing space. However, in previous studies, the data derived from the Luancheng Station Lysimeter have been repeatedly successfully applied to the remote sensing data to estimate the ET verification work (see [12,27]). Therefore, we also used the previous data as evidence. The observation data provided by the Lysimeter were compared with the average ETa provided by the Luancheng Station 3 × 3 around a 1 km pixel to reduce the interference of the ambient uncertainties. Figure 4 shows the comparison between the estimated ETa and the measured value from 2005 to 2006. It can be seen that the estimated ETa using the model is in a good agreement with the measured values. Fluctuations in the state of the two are almost the same. The data show two distinct peak times, one in 4-5 months and one in 7-9 months, because of the winter wheat-summer corn rotation. Based on the statistical analysis, the daily average deviation is 0.86 mm and the SEE is 1.2 mm between the estimation and the observed values. Figure 5 shows the ETa spatial distribution of the Hebei Plain 2001-2006 winter wheat growing season (October 5 to June 10 of the following year). The figure shows that the annual ETa is greater in the western than the eastern regions of the study area. Because the western region is one of the major winter wheat-producing areas, its soil and climate are more appropriate; in addition, the developed irrigation systems ensure the moisture requirement for the growth of winter wheat. In contrast, in the eastern region, particularly in the northeast (top right), the winter wheat growing fields are mostly dependent on the rain. Thus, the ETa is much lower. The reason why the ETa is relatively low in the southeast area (lower right) is also the same as that of the northeast area. In Table 1 However, in 2003, the rainfall was up to 248.7 mm because the rainfall was above normal in the fall. During the period from October 10 to October 12, the strongest rainfall occurred which represents the highest amount in nearly 50 years of the same period in the history. However, ETa was not the highest For the 480 points of county-level data (80 counties and 6 years), the ETa accounted for 38.3% (up to 184 points) of the total range from 400 mm to 450 mm and for 73.1% of the total range from 350 mm to 450 mm. The result is consistent with the result obtained using the models and the remote sensing data [12,[16][17][18][19][20]. The result is also in accordance with the ETa range of the winter wheat growing season under the conditions of full irrigation determined by field experiments [11,18,[28][29][30][31][32][33][34]. Figure 6 represents the average CWP distribution from 2001 to 2006. According to Figure 6, the CWP in the western and central areas was higher than that in the other regions, ranging from 1.2 kg/m 3 to 2.2 kg/m 3 . In particular, in the central west region, the CWP is the highest, ranging from 1.6 kg/m 3 to 2.2 kg/m 3 . In the eastern region, the CWP is mostly lower than 1 kg/m 3 . An average of 1.32 kg/m 3 , a maximum of 2.23 kg/m 3 , and a minimum value of 0.05 kg/m 3 were found in the 480 CWP points of the research area (80 counties, 6 years). After removal of less than 5% and greater than 95% of the cumulative frequency distribution, the frequency distribution ranges from 0.29 to 1.67 kg/m 3 . The other findings of the region were 0.12-2.15 kg/m 3 [15], 0.05-1.58 kg/m 3 [16], 0.2-1.63 kg/m 3 [29], and 0.66-1.06 [30][31][32][33][34]. Through studying the relationship between the yield and the ETa with the average CWP in many years, it can be concluded that the CWP was positively correlated with the yield, with an 2 of 0.9, as shown in Figure 7. Although the relationship was more complex between the CWP and the ETa, it was nearly a cosine curve. In the eastern region, the ETa was low, while the decrease of yields of the regions was lower than the decrease of ETa. Thus, CWP was maintained at a relatively high level. In contrast, both the ETa and yield increased, whereas the CWP declined. In the middle region, the magnitude of the increasing yield was greater than the ETa and the CWP also increased. When the range of the ETa was