Estimating Growing Season Evapotranspiration and Transpiration of Major Crops over a Large Irrigation District from HJ-1A/1B Data Using a Remote Sensing-Based Dual Source Evapotranspiration Model

: Crop evapotranspiration (ET) is the largest water consumer of agriculture water in an irrigation district. Remote sensing (RS) technique has provided an effective way to map regional ET using various RS-based ET models over the past several decades. To map growing season ET of different crops and partition ET into evaporation (E) and transpiration (T) at regional scale, appropriate ET models should be further integrated with crop distribution maps in different years and crop growing seasons determined for each crop pixel. In this study, a hybrid dual-source scheme and trapezoid framework-based ET Model (HTEM) fed with HJ-1A/1B data was applied in Hetao Irrigation District (HID) of China from 2009 to 2015 to map crop growing season ET and T at 30 m resolution. The HTEM model with HJ-1A/1B data performed well in estimating ET in HID, and the finer spatial resolution of model input data can improve the estimation accuracy of ET. Combined with the annual crop planting map identified in previous study, and crop growing seasons determined from fitted Normalized Difference Vegetation Index (NDVI) curves for crop pixels, the spatial and temporal variations of growing season ET and T of major crops (maize and sunflower) were examined. The results indicate that ET and T of maize and sunflower reach their minimum values in the southwest HID with smaller crop planting density, and reach their maximum values in northwest HID with higher crop planting density. Over the study period with a decreasing trend of available irrigation water, ET and T in maize and sunflower growing seasons show decreasing trends, while ratios of T/ET show increasing trends, which implies that the adverse effect of decreased irrigation water diversion on crop growth is diminished due to the favorable portioning of E and T in cropland of HID. In addition, the calculation results of crop coefficients show that there is water stress to crop growth in the study area. The present results are helpful to better understand the spatial pattern of crop water consumption and water stress of different crops during crop growing season , and provide the basis for optimizing the spatial distribution of crop planting with less water consumption and more crop yield.


Introduction
Agriculture is a major water consumer in China and around the world, and agricultural water use accounts for about 70% of the total water uses in the world. In arid and semi-arid areas where agriculture relies heavily on irrigation, the proportion of agriculture water use is higher than other only to crop type, but also to agricultural practice and crop growing environment. As a result, the growing season of the same crop may be different at different locations in an irrigation district [36]. Therefore, before extracting ET from crop growing season in large irrigation district, the crop growing season should be determined in advance. At small spatial scale, the crop growing season can be obtained through field survey [18]. While at regional scale, the crop growing season can be estimated by the Polyfit-Maximum method, which considers the period of greatest decrease in the seasonal NDVI time series as the beginning of vegetation dormancy [38][39][40].
The main objective of this study was to study the growing season evapotranspiration and transpiration of major crops at irrigation district scale by combining a RS-based dual source ET model (HTEM), annual crop distribution map, and growing season of maize and sunflower in HID determined by a simple algorithm from NDVI series. Compared with previous studies on RS-based ET estimates [6,30,37], this study aimed to estimate pixel-scale ET and T for different crops in their growing seasons, instead of pixel-scale ET and T in a whole region or in all croplands for a specified period of time. These results can provide the basis for optimizing the spatial distribution of crop planting with less water consumption and more crop yield.

Materials and Methods
The framework of the study is shown in Figure 1. Firstly, we used HJ-1A/1B satellite data as the input of HTEM model to estimate the daily ET and T of the whole study area. Secondly, the annual ET and T at crop pixels were extracted by combining the crop planting map provided by Yu and Shang (2017) [36]. Thirdly, we determined the crop growing season based on fitted NDVI curve [36], and finally we mapped crop ET and T during the growing season. Finally, the crop coefficients were calculated by combining the calculation results of ET and reference evapotranspiration (ET0).

Study Area
Hetao Irrigation District (HID) is a typical irrigation district of arid region located in the western Inner Mongolia Autonomous Region in North China. The HID is the third largest irrigation district and an important food production base in China, which has a good representativeness in arid area as the largest irrigation district in arid area of China.
In this study, four counties (Dengkou, Hangjinhouqi, Linhe, and Wuyuan) in western and middle HID were selected as the study area (Figure 2), which spans from 40.1°N to 41.4°N and from 106.1°E to 109.4°E, and covers an area of 910,000 ha with 44% being agricultural land ( Figure 2). The study area is characterized by arid continental climate with an annual precipitation of approximately 160 mm and pan evaporation (20 cm evaporation pan) of approximately 2240 mm during the study period from 2009 to 2015 [36]. The annual precipitation in the study area is much less than potential evaporation, which results in the high dependence of the agriculture on water diversion from the Yellow River for irrigation. However, a water-saving rehabilitation program has been carrying out since 1998 to reduce water diversion from the Yellow River to alleviate the rising of groundwater table and slow down the loss of river runoff [37,41,42]. Consequently, the reduction of irrigation water may have a negative impact on crop growth in the study area, and it is necessary to study whether the impact is beneficial or not.
Summer maize and sunflower are dominant crops in the study area, and the proportion of maize and sunflower planting has increased significantly in recent years to over 80% of the total agricultural land after 2015 [36,43]. Accordingly, this study focused on maize and sunflower. From 2009 to 2015, the annual average yields of maize and sunflower in the study area ranged from 10.61 t/ha to 11.41 t/ha and from 3.54 t/ha to 3.87 t/ha, respectively [43].

Satellite Data and Auxiliary Data
China launched the HJ-1A/1B satellite constellation for environment and disaster monitoring on 6 September, 2008, the satellite image of HJ-1A/1B has the same band division and spatial resolution (30 m) as that of Landsat, while the revisit time of HJ-1A/1B (2-4 days) is much shorter than that of Landsat (16 days). In this study, RS-based ET estimation requires four visible bands (B1-B4) data of the charge coupled device (CCD) carried by HJ-1A and HJ-1B satellites and thermal infrared band (B8) data of the infrared spectrograph (IRS) carried by HJ-1B [44,45]. Specifically, B1-B4 are blue (0.43-0.52 µm), green (0.52-0.60 µm), red (0.63-0.69 µm), and near-infrared (0.76-0.90 µm) bands, respectively, with a spatial resolution of 30 m and a temporal resolution of 2 days, and B8 is thermal infrared (10.5-12.5 µm) band with a spatial resolution of 300 m and a temporal resolution of 4 days [44]. The original images were downloaded from the China Centre for Resources Satellite Data and Application (CRESDA) (http://cresda.spacechina.com), and then were all re-projected into Universal Transverse Mercator (UTM) projection and resampled with a spatial resolution of 30 m. The satellite images used in this study consisted of 145 HJ-1A/1B CCD images and 145 HJ-1B IRS images of Level 2A for Hetao Irrigation District covering the years from 2009 to 2015, with the cloud cover less than 5%. Consequently, there was at least one date with the above satellite images per week during the crop growing season.
The preprocessing procedures for HJ-1A/1B images mainly include radiometric calibration and atmospheric and geometric corrections, and more details can refer to Li (2014) [45] and Yu and Shang (2017) [36]. Then, the parameters of HTEM model can be retrieved from the pre-processed satellite images, in which CCD image was used to retrieve vegetation index and surface albedo, and IRS image was used to retrieve land surface temperature, and the detailed calculation equations can be found in Li (2014) [45].
Land use map (1:100,000) of the study area ( Figure 3) for the year 2000 was provided by Environmental and Ecological Science Data Center for West China, National Natural Science Foundation of China (http://westdc.westgis.ac.cn), which was also resampled to a spatial resolution of 30 m.
Daily meteorological data at the Linhe (LH) (40.73°N and 107.37°E, and 1041.1 m above sea level) Meteorology Station (Figure 2) are available from China Meteorology Data Sharing Service System (http://data.cma.cn).
Statistics on water balance components of the study area were provided by the General Administration of Hetao Irrigation District (http://www.zghtgq.com).

Determination of Crop Growing Season
Many studies have shown that Normalized Difference Vegetation Index (NDVI) sequences of crops are closely related to crop phenological processes [46,47]. Due to the limitation of satellite revisit time and the effect of cloud cover, daily NDVI retrieved from RS is discontinuous, which needs further processing. Available studies have shown that an asymmetric logistic curve can be used to fit the NDVI time series [48], which has also been successfully applied in our study area [36,41]. Following Royo et al. (2004) [48], the fitting curve equation is where t is the day of year (DOY), a, b, c, d, and k are curve parameters ( Figure 3) and can be calibrated from NDVI series of each crop pixel using least-squares method [36]. As can be seen from Figure 3, the fitted NDVI curve has three characteristic points, that is, the left inflection point (t_inf_1, NDVI_inf_1) with the maximum growth rate, the peak point (t_max, NDVI_max) with the maximum NDVI, and the right inflection point (t_inf_2, NDVI_inf_2) with the maximum withering rate. The equations for calculating the above characteristic values can be obtained from Yu and Shang (2017) [36].
Based on the crop identification results in the study area [36], we extracted multi-year characteristic values of NDVI curves for maize and sunflower pixels and analyzed their distributions. It was found that t_inf_1 and t_inf_2 were closely related to the sowing and harvest dates of crops, respectively. In our study area, the multi-year average values of t_inf_1 and t_inf_2 of maize pixels are the 183rd day and 253rd day, respectively, while those for sunflower pixels are the 191st day and 243rd day, respectively. According to field survey, the actual sowing dates of maize and sunflower are early May and early June, respectively, and the actual harvest dates of both crops are late September [36]. Considering the relationship between the actual crop growing season and the characteristic values of NDVI curve, the estimation equations of sowing and harvest dates of maize and sunflower in HID are as follows: where SOS and EOS represent the start and end of growing seasons, respectively, and subscripts m and sf represent maize and sunflower, respectively. In addition, we compared the growing seasons of maize and sunflower determined by the above equations with the field measured data in Yangchang canal command area (YCA) observation site and the statistical data of Linhe (LH) Meteorology Station ( Figure 2) to validate the accuracy of the estimation Equations (3)-(6).
ns n nc where substitutes c and s represent the canopy and soil components, respectively; kc is the extinction coefficient of radiation attenuation within the canopy (dimensionless); LAI is the leaf area index (m 2 ·m −2 ); and Fc is the fractional vegetation cover that can be computed from remotely sensed NDVI [49].
where NDVImax is the NDVI for complete canopy cover, and NDVImin is the NDVI for bare soil. Ground heat flux in HTEM is assumed to be a constant fraction of the soil net radiation [50]: where c is the ratio of G over Rns, and is taken to be 0.30 here. For each component, the sensible heat flux is computed from Norman et al. (1995) [27]. Based on the surface energy balance, the latent heat flux for each component is computed as a residual in equations (13) and (14): where ρ is air density (kg·m −3 ); Cp is specific heat capacity of air at constant pressure (J·kg -1 K -1 ); Tc and Ts are canopy and soil surface temperatures (K), respectively; Ta is air temperature (K) at the reference height; ra h , ra a and ra s are the aerodynamic resistances that can be estimated following Sánchez et al. (2008) [51]. According to the above equations, the prerequisite for calculating the latent heat flux is to know surface temperatures for soil and canopy, which can be obtained from a trapezoidal Fc/LST space [52] in HTEM ( Figure 4). As shown in Figure 4, there are four critical points (A, B, C, and D) corresponding to four extreme conditions in the trapezoidal space. Specifically, A and B represent completely dry points in the bare soil and the fully vegetated areas, respectively, while C and D represent non-water stress points in the fully vegetated and the bare soil areas, respectively. Furthermore, it is assumed that the evapotranspiration rate on the warm edge AB is zero, and that on the cold edge CD is equal to its potential rate. Some existing soil wetness isolines were found in the Fc/LST space via soilvegetation-atmosphere transfer modeling [53]. It is assumed that soil with a uniform texture and sharing the same moisture content is also isothermal. Points E and F represent any two pixels on the same soil wetness isoline in the study area, which have the same soil wetness and temperature. The slope of line EF (any isoline) can be retrieved by interpolating the slopes of the warm edge (a) and cold edge (b) according to the temperature differences between the pixel and the warm and clod edges. Then, the soil temperature can be computed from After calculating the soil and canopy temperatures, the sensible and latent heat fluxes can be calculated from equations (13) and (14). Then the instantaneous ET can be calculated.
Assuming that the reference evapotranspiration fraction (EF), the ratio of ET over ET0, is constant throughout a day, the instantaneous ET can be scaled up to daily ET [54]. The instantaneous EF on days without RS image can be obtained by linearly interpolating the EF values between two adjacent days with RS images [55]. Specifically, the instantaneous and daily ET0 were calculated using Penman-Monteith method, recommended by the Food and Agriculture Organization (FAO) [17].

Evaluation of HTEM Performance
The HTEM model has been validated with ground measurements of ET for different crops in previous studies, including the Soil Moisture-Atmosphere Coupling Experiment (SMACEX) campaign conducted in the summer of 2002 in Iowa, USA, with the main crops of corn and soybean [28], the Weishan flux site in the North China Plain during the growing season of winter wheat and summer corn in 2007 [28], and the MUlti-Scale Observation EXperiment on Evapotranspiration over heterogeneous land surfaces 2012 (MUSOEXE-12) campaign in the Heihe River Basin in Northwest China [29]. Meanwhile, the model was also validated with ET and T measurements over a semiarid shrub ecosystem at the Lucky Hills study site of the U.S. Department of Agriculture-Agricultural Research Service (USDA-ARS) Walnut Gulch Experimental Watershed in southeastern Arizona, U.S. [56].
However, there was no ground measurement of T during the study period in the study region, and we validated the model with crop ET measurements in one experiment study, annual regional ET estimates from water balance analysis, and modelled ET and T in one experiment site. Similar validations were also used in previous studies of ET and T estimations in the same study region using RS-based ET models [6,37].
At field scale, the daily ET estimated by HTEM model was compared with average ET of 5-7 days measured at Dengkou Agricultural Experiment Station (Figure 2 [18] during the same period at YCA observation site (Figure 2) using the HYDRUS-1D Software package for simulating the onedimensional movement of water, heat, and multiple solutes in variably-saturated media [58].
At regional scale, Annual ET during growing seasons of 2009-2015 produced by HTEM model was compared with regional water balance analysis [6].
where ET_Water Balance is regional annual ET calculated from regional water balance equation (mm); P is annual precipitation; I is irrigation water; D is drainage water; ΔS is annual variation of total water storage in the unsaturated zone and unconfined aquifer; ΔG is annual change in groundwater table depth; A is total area in the study area; and µ is the specific yield taken as 0.046 [37].
To evaluate the performance of HTEM model, two indicators, including the root mean square error (RMSE) and mean relative error (MRE), were defined as follows [59]: where N is the number of data pairs to evaluate the model performance; Si and Oi, i = 1, 2, …, N, are the ith model simulated and observed (or other model simulated) values of ET, respectively.

Crop Growing Season
The average estimated values of SOS and EOS of maize and sunflower pixels were compared with the observed sowing and harvest dates of crops at YCA observation site [18] (Figure 2) in 2013 and the statistical sowing date of maize and sunflower from 2010 to 2013 provided by LH Meteorology Station (Figure 2), respectively ( Table 1). As shown in Table 1, the absolute errors between estimated and observed SOS of maize range from 2 days (YCA in 2013) to 11 days (LH in 2011), while those of sunflower range from 0 (LH in 2011) to 3 days (YCA in 2013). The estimated EOS of maize in 2013 at YCA is 3 days earlier than the observation, while that of sunflower is 3 days later than the observation. In general, the estimation accuracy of sunflower growing season is higher than that of maize. The main reason is that the SOS of sunflower changed little between years, while that of maize varied greatly, for example, there was a 10-day difference between 2010 and 2012. In general, the accuracy of estimated growing seasons of maize and sunflower is acceptable. In addition, crop ET is small in the early stage of growing season; therefore, the effect of estimation error of growing season on the total ET during crop growing season can be neglected.   Through comparison, it was found that SOS of crops was generally earlier in the western part, which is consistent with the agricultural practices in HID. For maize, the multi-year average value of SOS in Wuyuan is the 127th day, which is five days later than that in Dengkou. Similarly for sunflower, the multi-year average value of SOS in Wuyuan (the 159th day) is seven days later than that in Dengkou. Compared with SOS of crops in HID, the EOS values of maize and sunflower in different counties had little difference. The multi-year average values of SOS and EOS of maize in the whole study area are the 124th day and 268th day, respectively, while those of sunflower are the 157th day and 268th day, which are in good agreement with those of the field survey in 2012 [36].

Validation of HTEM
At field scale, the HTEM estimated daily ET closely matches the field observed values [57] (Figure 7a), with the RMSE of 0.52 mm/day and the MRE of 11.9%. Compared with the ET estimation results of HTEM-ABL (Atmospheric Boundary Layer) (0.63 mm/day for RMSE and 13.0% for MRE) [6] and SEBAL (0.53 mm/day for RMSE and 14.6% for MRE) [37] fed with Moderate Resolution Imaging Spectroradiometer (MODIS) data at 250 m resolution in the same study area, the model performance in this study at 30 m resolution is slightly better. Therefore, in HID, where the crop planting pattern is complex [19], the finer spatial resolution of satellite images may be favorable to improve the estimation accuracy of ET. However, the model precision was only assessed with limited field experiment data. More experiment data should be collected in the future to examine whether the spatial resolution of RS data influences the precision of ET models.
At regional scale, the comparison between ET during growing season estimated by HTEM model and the output of water balance model shows good agreement, with the RMSE of 48.9 mm and the MRE of 9.1%, respectively (Figure 7b). Compared with other similar studies [6,30], the ET estimation accuracy of this study is slightly lower at the regional scale (especially in 2011 with fewer satellite images). This is because the temporal resolution of HJ-1B IRS images is lower than that of MODIS and other data fusion products, which will result in slightly greater errors due to temporal interpolation errors. In areas with no ground measurements of ET and T, it is common to evaluate the model performance with ET estimates from regional water balance, which has been widely used in previous studies [6,30,37,41]. According to equations (17) and (18), precipitation (P), irrigation water (I), drainage water (D), and the change in groundwater table depth (ΔG) are all measured values, only the specific yield (µ) is an empirical value. Here the specific yield took the value of 0.046 that has been used in previous studies [6,37]. Moreover, the variation of total water storage in the unsaturated zone and unconfined aquifer is not significant at the annual scale. As a result, the ET estimates from regional water balance has relatively high accuracy, which can be taken as reference data to validate model outputs.
At the YCA observation site, relative differences of maize ET and T estimated in this study and those modelled by Ren et al. (2016) [18] are only 1% and 3%, respectively, while the differences for sunflower ET and T are both 11%, which also showed that the accuracy of crop growing season ET and T estimates in this study are satisfactory.

Spatial Patterns of Crop Growing Season ET and T
Compared with RS-based single-source ET model [30,37,57], the prominent advantage of the dual-source model is that it can obtain the spatial distribution of both ET and T [6,27]. According to the above calculation results of crop growing season (Figures 5 and 6), we extracted the ET and T distribution of maize and sunflower during the growing season, as shown in Figures 8 and 9.
Overall, no matter for maize or sunflower, the spatial patterns of ET and T in different years are similar, and the spatial variation of T is consistent with ET for each crop. Specifically, for maize (Figure 8 Hangjinhouqi is the largest [36]. It can be seen that the higher the planting density, the higher ET and T. Moreover, the main land use type in Dengkou is sandy land (Figure 2), which has poorer water holding capacity, which leads to poorer crop growth. The spatial patterns of crop ET are consistent with the results obtained by   [37] and Chen et al. (2018) [41] using the SEBAL model. Moreover, the spatial variation of crop T is greater than that of crop ET.    show an increasing trend. Our results imply that the adverse effect of the reduced irrigation water diversion on crop growth is narrowed due to new portioning of evaporation and transpiration in cropland.

Crop Coefficients
Based on the calculation results of daily crop ET and T and daily ET0 obtained by Penman-Monteith method [17], we calculated daily crop evapotranspiration coefficient (Kc, the ratio of ET over ET0) and transpiration coefficient (Kcb, the ratio of T over ET0) per crop pixel in the study area. We extracted the average values of maximum and mean Kc and Kcb of maize and sunflower in each county of the study area during growing season (Figures 12 and 13).  As shown in Figure 12, the multi-year average values of maximum Kc (Kc, max) of maize and sunflower over the whole study area are 1.003 and 0.996, respectively, while those of maximum Kcb (Kcb, max) of maize and sunflower are 0.894 and 0.898, respectively. The maximum Kc of maize is slightly greater than that of sunflower, while the maximum Kcb of maize is slightly smaller than that of sunflower. In addition, Kc of maize and sunflower both reach their maximum values in 2012, which are 1.101 and 1.087, respectively. The results of this study are slightly less than the standard crop coefficients recommended by FAO [17]. According to the values of crop coefficients, smaller Kc, max and Kcb, max in Dengkou indicate that maize and sunflower in Dengkou are likely to suffer from severer water or other stresses in the study region, while smaller Kc, max and Kcb, max in 2014 indicate severer stresses to crop growth in 2014 during the study years.
As shown in Figure 12, the multi-year average values of mean Kc (Kc, mean) of maize and sunflower over the whole study area are 0.594 and 0.621, respectively, while those of mean Kcb (Kcb, mean) of maize and sunflower are 0.396 and 0.407, respectively ( Figure 13). Both the Kc, mean and Kcb, mean of sunflower are slightly higher than those of maize. Similar to Kc, max and Kcb, max, the Kc, mean and Kcb, mean of maize and sunflower both reach the minimum values in Dengkou; however, the differences of Kc, mean and Kcb, mean between different counties are less than those of Kc, max and Kcb, max.

Conclusions
In this study, the HTEM model fed with HJ-1A/1B data with the spatial resolution of 30 m performed well in estimating crop ET and T during the growing seasons from 2009 to 2015 in Hetao Irrigation District. Compared with the ET estimates based on MODIS data with 250 m spatial resolution, the finer spatial resolution of model input data here may be favorable to improve the estimation accuracy of ET at both field and regional scales. We also developed a simple NDVI curve fitting method to estimate the growing seasons of maize and sunflower in the study area with relative higher accuracy. By integrating the estimated daily ET and T series, and maize and sunflower planting maps and growing seasons, we extracted the spatial and temporal distribution of growing season ET and T of maize and sunflower in the study area.
Our results showed that the total ET and T of maize and sunflower in the growing season both reach their minimum values in Dengkou, while reach their maximum values in Hangjinhouqi. The higher the crop planting density, the higher the crop ET. From 2009 to 2015, ET and T of both maize and sunflower show decreasing trends, while T/ET of both maize and sunflower show increasing trends, which imply that there is no obvious adverse effect on crop growth due to the reduction of irrigation water diversion in HID. The calculated crop coefficients showed that the crop growth in the study area is under certain water or other stresses. Crops in Dengkou suffered severer stress in the study region, while crops in 2014 suffered severer stress during the study period.
In future studies, crop yield [43] and crop water consumption estimated in this study can be integrated to assess water productivity and planting suitability of major crops in the study region, so as to find optimal crop planting pattern with less water consumption and more crop yield/benefit. Author Contributions: S.S. outlined the research topic, assisted with manuscript writing, and coordinated the revision activities. B.Y. performed data collection and processing, data analysis, model building, the interpretation of results, manuscript writing, and coordinated the revision activities.