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

Abstract

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.

1. 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 areas [1]. Irrigation water accounts for 90% of the total agriculture water use [2]. With the increasing domestic, industrial and environmental water uses, the amount of water available for agriculture has been decreasing in recent years [3,4]. As a result, agricultural production in arid and semi-arid areas is facing severer challenge. Especially for water shortage areas, the uneven spatiotemporal distribution of water resources aggravates the challenge of agricultural production. Therefore, quantifying agriculture water consumption for different crops over large area is important for agricultural water management in arid and semi-arid areas [5].

Crop evapotranspiration (ET), the total amount of water evaporation from soil (E) and canopy interception and transpiration by vegetation (T) [6,7], is the largest sink of water for agriculture [8]. Accurate estimation of the spatial distribution and temporal variation of crop growing season ET is the basis for agricultural water management. Conventional crop ET estimation methods can generally be categorized into three groups. The first group is ground-based measurement techniques, such as Bowen ratio (BR) [9,10], eddy covariance (EC) [11,12], and weighing lysimeter systems [13,14]. These methods are generally expensive and time-consuming [6], and applicable only for small scales. Consequently, their accuracies degrade when extrapolating to regional scale [15]. The second group includes weather-based calculation formulas, such as the Penman–Monteith equation, recommended by the Food and Agriculture Organization [16,17], which estimate crop ET by multiplying the reference evapotranspiration (ET₀) calculated from meteorological data by crop coefficients (K_c) based on crop types and growth stages. However, it is generally difficult to predict the idealized K_c over large regions [8]. The third group is physically based simulation models that estimate crop ET from simulation models for water flow in the soil-plant-atmosphere continuum (SPAC) and/or crop growth [18,19,20,21]. However, these models require a large number of in situ measured data as parameters and inputs [22]. Therefore, the above-stated methods are generally applicable to small-scale crop ET estimation and difficult to be extended to large areas.

Fortunately, the application of remote sensing (RS) technology in earth observation provides an alternative approach for mapping ET at regional and even global scales [5,8,23,24]. Accordingly, various RS-based ET models based on surface energy balance equation have been developed to adapt to different underlying surface characteristics [25], such as the Surface Energy Balance Algorithm for Land (SEBAL) [26], the Two-Source Energy Balance Model (TSEB) [27], and the Hybrid dual-source scheme and Trapezoid framework based ET Model (HTEM) [28]. SEBAL is a single-source model, which cannot distinguish E and T from ET. While TSEB and HTEM are dual-source models, which can distinguish E and T by decomposing component (soil and vegetation) temperatures from the radiative land surface temperature (LST) observed by satellite RS. TSEB determines component temperatures by using an iterative approach [27], while HTEM uses the trapezoidal vegetation index (VI) and LST space [28]. Compared with the measured data of EC flux tower, HTEM performed better than TSEB in arid area [29]. However, these RS-based ET models usually give instantaneous ET estimates at the satellite overpassing time or ET on specified periods after temporal upscaling. To study crop ET in the whole crop growing season or different growing stages using these ET models, knowledge on crop distribution [30] and crop growing season [22] is necessary.

The fragmented and inter-annual varied crop planting pattern is very common in China and many other parts of the world [19,31,32], which restricts the availability of timely crop planting map. Hetao Irrigation District (HID) in Inner Mongolia is a typical example of such case, which has a fragmented crop planting pattern due to the special ‘‘farmland use rights” [18,33]. The most frequently used RS-based crop identification method is supervised classification with the basic unit of a pixel, which classifies each pixel into a specified type of crop according to the training data based on the sampling points [34,35]. Based on supervised classification, Yu and Shang (2017) [36] mapped multi-year distribution of major crops in HID at 30 m resolution accurately using ellipse space of Normalized Difference Vegetation Index (NDVI) and phenological characteristics derived from HJ-1A/1B charge coupled device (CCD) data.

The determination of crop growing season is another major factor in estimating growing season ET. Unlike crop ET simulated by crop growth model [18], crop ET estimation based on RS during growing season is usually extracted on a monthly scale [30,37]. Crop growing season is related not 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.

2. 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 (ET₀).

2.1. 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].

2.2. 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).

2.3. 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

$$\mathrm{NDVI}=a+(b/k)\cdot {(1+n)}^{-(k+1)/k}\cdot n\cdot {(k+1)}^{(k+1)/k},$$

(1)

$$n=\exp \left[\left(t+d\cdot \ln (k)-c\right)/d\right],$$

(2)

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:

$${\mathrm{SOS}}_{\mathrm{m}}=t\_\inf \_1-60,$$

(3)

$${\mathrm{EOS}}_{\mathrm{m}}=t\_\inf \_2+15,$$

(4)

$${\mathrm{SOS}}_{\mathrm{sf}}=t\_\inf \_1-45,$$

(5)

$${\mathrm{EOS}}_{\mathrm{sf}}=t\_\inf \_2+25,$$

(6)

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).

2.4. A Brief Description of HTEM Model

The HTEM model uses a layer approach to allocate total net radiation (R_n) (W·m⁻²) for canopy (R_nc) and soil (R_ns) components (Equations (7,8)), and adopts a patch approach to partition available energy into the latent heat (LE), sensible heat (H), and ground heat fluxes (G) (Equations (9) and (10)) [28].

$$R_{\mathrm{nc}}=R_{\mathrm{n}}(1-\exp (-k_{\mathrm{c}}\mathrm{LAI})),$$

(7)

$$R_{\mathrm{ns}}=R_{\mathrm{n}}-R_{\mathrm{nc}},$$

(8)

$$R_{\mathrm{nc}}=F_{\mathrm{c}}(L{E}_{\mathrm{c}}+H_{\mathrm{c}}),$$

(9)

$$R_{\mathrm{ns}}-G=(1-F_{\mathrm{c}})(L{E}_s+H_s),$$

(10)

where substitutes c and s represent the canopy and soil components, respectively; k_c is the extinction coefficient of radiation attenuation within the canopy (dimensionless); LAI is the leaf area index (m²·m⁻²); and F_c is the fractional vegetation cover that can be computed from remotely sensed NDVI [49].

$$F_c=1-{\left(\frac{{\mathrm{NDVI}}_{\max }-\mathrm{NDVI}}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}\right)}^{0.696},$$

(11)

where NDVI_max is the NDVI for complete canopy cover, and NDVI_min is the NDVI for bare soil.

Ground heat flux in HTEM is assumed to be a constant fraction of the soil net radiation [50]:

$$G=c{R}_{\mathrm{ns}},$$

(12)

where c is the ratio of G over R_ns, 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):

$$F_{\mathrm{c}}L{E}_{\mathrm{c}}=R_{\mathrm{nc}}-F_{\mathrm{c}}{H}_{\mathrm{c}}=R_{\mathrm{nc}}-F_{\mathrm{c}}\rho C_{\mathrm{p}}\frac{T_{\mathrm{c}}-T_{\mathrm{a}}}{r_{\mathrm{a}}^{\mathrm{h}}},$$

(13)

$$(1-F_{\mathrm{c}})L{E}_{\mathrm{s}}=R_{\mathrm{ns}}-G-(1-F_{\mathrm{c}})H_{\mathrm{s}}=R_{\mathrm{ns}}-G-(1-F_{\mathrm{c}})\rho C_{\mathrm{p}}\frac{T_{\mathrm{s}}-T_{\mathrm{a}}}{r_{\mathrm{a}}^{\mathrm{a}}+r_{\mathrm{a}}^{\mathrm{s}}},$$

(14)

where ρ is air density (kg·m⁻³); C_p is specific heat capacity of air at constant pressure (J·kg⁻¹ K⁻¹); T_c and T_s are canopy and soil surface temperatures (K), respectively; T_a is air temperature (K) at the reference height; r_a^h, r_a^a and r_a^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 F_c/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 F_c/LST space via soil-vegetation-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

$$T_{\mathrm{s}}=F_{\mathrm{c}}\frac{b}{a+b}(T_{\mathrm{s}\_\max }-T_{\mathrm{c}\_\max })+\mathrm{LST}$$

(15)

According to Sánchez et al. (2008) [51], the canopy temperature can be estimated from

$$T_{\mathrm{c}}=\frac{1}{F_{\mathrm{c}}}\left[\mathrm{LST}-(1-F_{\mathrm{c}})T_{\mathrm{s}}\right]$$

(16)

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 ET₀, 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 ET₀ were calculated using Penman–Monteith method, recommended by the Food and Agriculture Organization (FAO) [17].

2.5. 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) in 2009 [41,57] and modelled ET and T at Yangchang canal command area (YCA) observation site. Specifically, during the growing season of 2009, soil moisture at 0–20, 20–40, 40–60, 60–80, and 80–100 cm soil layers were measured by both oven drying method and Time domain Reflectometry with Intelligent MicroElements (TRIME-IPH) methods with the sampling interval of 5–7 days. ET between sampling days were obtained through the calculation of soil water balance. Meanwhile, the estimated ET and T in 2012 were compared with the simulation results of Ren et al. (2016) [18] during the same period at YCA observation site (Figure 2) using the HYDRUS-1D Software package for simulating the one-dimensional 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].

$$\mathrm{ET}\_\mathrm{Water} \mathrm{Balance}=P+I-D-\mathsf{\Delta }S,$$

(17)

$$\mathsf{\Delta }S=\mathsf{\Delta }G\times A\times \mu ,$$

(18)

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]:

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left(S_i-O_i\right)}^2}},$$

(19)

$$\mathrm{MRE}=\frac{1}{N}{\displaystyle \sum _{i=1}^N\frac{\left(S_i-O_i\right)\times 100}{O_i}},$$

(20)

where N is the number of data pairs to evaluate the model performance; S_i and O_i, i = 1, 2, …, N, are the ith model simulated and observed (or other model simulated) values of ET, respectively.

3. Results and Discussion

3.1. 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. Validation of estimated start (SOS) and end (EOS) of growing seasons of maize and sunflower at Line (LH) Meteorology Station and Yangchang canal command area (YCA).

Crop Type	Estimation/Observation	SOS at LH (DOY)				SOS and EOS at YCA (DOY)
		2010	2011	2012	2013	SOS (2013)	EOS (2013)
Maize	Estimation	120	129	125	124	120	268
	Observation	127	118	117	118	122	265
Sunflower	Estimation	157	155	161	158	151	262
	Observation	---	155	---	157	154	265

). 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.

Furthermore, we extracted the SOS and EOS of maize and sunflower pixels from 2009 to 2015 (Figure 5 and Figure 6) based on crop planting map [36].

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].

3.2. 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.

3.3. 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 (Figure 5 and Figure 6), we extracted the ET and T distribution of maize and sunflower during the growing season, as shown in Figure 8 and Figure 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), ET and T both reach their minimum values in Dengkou with the multi-year average values of 473 and 313 mm, respectively, while reach their maximum values in Hangjinhouqi with the multi-year average values of 537 and 364 mm, respectively. Similarly, for sunflower (Figure 9), ET and T are also minimum in Dengkou with the multi-year average values of 397 and 270 mm, respectively, while maximum in Hangjinhouqi with the multi-year average values of 460 and 326 mm, respectively. In the study area, the planting density of crops in Dengkou is the smallest, while that in 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 Yang et al. (2012) [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.

3.4. Temporal Variations of Crop Growing Season ET and T

Figure 10 and Figure 11 show the temporal variations of crop growing season ET and T during growing season of different counties from 2009 to 2015.

For maize, the multi-year average ET values (Figure 10a1–d1) range between 383 (10th percentile) and 558 mm (90th percentile) in Dengkou where crop ET is lower, while between 411 (10th percentile) and 608 mm (90th percentile) in Hangjinhouqi where crop ET is higher. The multi-year average T values (Figure 10a2–d2) range between 198 (10th percentile) and 424 mm (90th percentile) in Dengkou, while between 262 (10th percentile) and 461 mm (90th percentile) in Hangjinhouqi.

In terms of the inter-annual variations of maize ET and T from 2009 to 2015, the average maize ET over the whole study area range from 456 (2014) to 591 mm (2009), while average maize T range from 282 (2015) to 394 mm (2012). In addition, both maize ET and T show decreasing trends due to the annual decrease of irrigation water diversion [6], while T/ET of maize show an increasing trend.

For sunflower, the multi-year average ET values (Figure 10a2–d2) range between 315 (10th percentile) and 480 mm (90th percentile) in Dengkou where crop ET is lower, while between 395 (10th percentile) and 525 mm (90th percentile) in Hangjinhouqi where crop ET is higher. The multi-year average T values (Figure 11a2–d2) are between 159 (10th percentile) and 380 mm (90th percentile) in Dengkou, while between 227 (10th percentile) and 419 mm (90th percentile) in Hangjinhouqi.

In terms of the inter-annual variations of sunflower ET and T from 2009 to 2015, the average sunflower ET over the whole study area range from 391 (2014) to 492 mm (2012), while average T from 260 (2015) to 351 mm (2012). Similar to maize, sunflower growing season ET and T both show decreasing trends due to the annual decrease of irrigation water diversion [6], while T/ET of sunflower 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.

3.5. Crop Coefficients

Based on the calculation results of daily crop ET and T and daily ET₀ obtained by Penman–Monteith method [17], we calculated daily crop evapotranspiration coefficient (K_c, the ratio of ET over ET₀) and transpiration coefficient (K_cb, the ratio of T over ET₀) per crop pixel in the study area. We extracted the average values of maximum and mean K_c and K_cb of maize and sunflower in each county of the study area during growing season (Figure 12 and Figure 13).

As shown in Figure 12, the multi-year average values of maximum K_c (K_{c, max}) of maize and sunflower over the whole study area are 1.003 and 0.996, respectively, while those of maximum K_cb (K_{cb, max}) of maize and sunflower are 0.894 and 0.898, respectively. The maximum K_c of maize is slightly greater than that of sunflower, while the maximum K_cb of maize is slightly smaller than that of sunflower. In addition, K_c 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 K_{c, max} and K_{cb, 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 K_{c, max} and K_{cb, 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 K_c (K_{c, mean}) of maize and sunflower over the whole study area are 0.594 and 0.621, respectively, while those of mean K_cb (K_{cb, mean}) of maize and sunflower are 0.396 and 0.407, respectively (Figure 13). Both the K_{c, mean} and K_{cb, mean} of sunflower are slightly higher than those of maize. Similar to K_{c, max} and K_{cb, max}, the K_{c, mean} and K_{cb, mean} of maize and sunflower both reach the minimum values in Dengkou; however, the differences of K_{c, mean} and K_{cb, mean} between different counties are less than those of K_{c, max} and K_{cb, max}.

4. 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.