A Modified Multi-Source Parallel Model for Estimating Urban Surface Evapotranspiration Based on ASTER Thermal Infrared Data

To date, little attention has been given to remote sensing-based algorithms for inferring urban surface evapotranspiration. A multi-source parallel model based on ASTER data was one of the first examples, but its accuracy can be improved. We therefore present a modified multi-source parallel model in this study, which has made improvements in parameterization and model accuracy. The new features of our modified model are: (1) a characterization of spectrally heterogeneous urban impervious surfaces using two endmembers (highand low-albedo urban impervious surface), instead of a single endmember, in linear spectral mixture analysis; (2) inclusion of an algorithm for deriving roughness length for each land surface component in order to better approximate to the actual land surface characteristic; and (3) a novel algorithm for calculating the component net radiant flux with a full consideration of the fraction and the characteristics of each land surface component. HJ-1 and ASTER data from the Chinese city of Hefei were used to test our model’s result with the China–ASEAN ET product. The sensitivity of the model to vegetation and soil fractions was analyzed and the applicability of the model was tested in another built-up area in the central Chinese city of Wuhan. We conclude that our modified model outperforms the initial multi-source parallel model in accuracy. It can obtain the highest accuracy when applied to vegetation-dominated (vegetation proportion > 50%) areas. Sensitivity analysis shows that vegetation and soil fractions are two important parameters that can affect the ET estimation. Our model is applicable to estimate evapotranspiration in other urban areas.


Introduction
Land surface evapotranspiration (ET) is the largest component of the land surface water and energy balance [1].It consists of soil evaporation and vegetation transpiration from land surfaces into the atmosphere, as well as the evaporation of canopy interception and surface water [2].ET is fundamental to water balance in natural environments, particularly in semi-arid and arid regions where ecosystem processes are often impacted by the limited availability of water and where water loss is dominated by ET [3,4].In urban environments, ET plays an essential role in alleviating urban heat island effect and regulating urban climate [5].Nevertheless, urban land surface evapotranspiration remains little studied due to the complexity of urban land covers and of urban ET processes.
A variety of techniques were used to assess evapotranspiration on different scales, which have been summarized by Verstraeten et al. [6].On point/leaf/plant and field scales, which rarely exceed 1-2 km [4,7], ET observation methods including lysimetry (LM), water balance (WB), Bowen ratio-energy balance (BR), Scintillometer (SCM) Sap-flow (SF), and the Penman-Monteith (PM) equation have been widely applied.On a landscape scale, methods such as Eddy covariance (EC), WB, and PROMET remain dominant.On large scales (i.e., regional and continental), the spatial and topographic heterogeneity of land surfaces are more complex [4].Therefore, thermal remote sensing, which can record heat information of large-scale land surfaces in the thermal infrared region of electromagnetic spectrum, has been combined with mathematical models and meteorological data to estimate ET [4,[8][9][10][11][12][13][14].Based on the turbulent heat flux exchange mode between land surfaces and near-ground atmosphere, remote sensing-based ET estimation models fall into two main categories: single-source and dual-source.
In single-source models, a land surface is considered homogeneous without distinguishing between vegetation and soil, i.e., the energy exchange interface is assumed to be a big leaf with a uniform component.A number of single-source models have been developed, such as the surface energy balance algorithm for land (SEBAL) model [15], mapping evapotranspiration with internalized calibration (METRIC) model [16,17], the surface energy balance system (SEBS) model [18], and the simplified surface energy balance index (S-SEBI) model [19,20].These models are simple and extensively used as their parameters are easy to obtain.However, they have large errors because of heterogeneous land cover, particularly in urban areas [4].
Under most conditions, vegetation canopy is not compact and is therefore mixed with vegetation and soil.Shuttleworth and Wallance [21] proposed a dual-source model in which evapotranspiration is considered consisting of soil evaporation and vegetation transpiration [22].In addition, according to the interaction mechanism and resistance connection model between soil and vegetation, dual-source models are either serial or parallel models.A serial model takes soil and vegetation canopy as a continuum where vapor and heat of the canopy converge at a hypothetical height inside the vegetation canopy from which it is then exchanged to the atmosphere, i.e., the water and energy flux of the entire land surface is equal to the sum of each layer's flux.Due to the difficulty of obtaining the required parameters, application of serial models is restricted.
In arid and semi-arid regions, vegetation distribution is sparse and uneven, and the coupling between soil and vegetation canopy is weak and can be ignored.Based on this, parallel models such as the two-source model (TSM) [23][24][25], the atmosphere-land exchange inverse model (ALEXI) [26][27][28] and the two-source trapezoid model for evapotranspiration (TTME) [29] consider vegetation and soil as two independent components without interaction, i.e., the turbulent exchanges with atmosphere for vegetation and soil are two parallel processes.Parallel models are in fact simplified versions of dual-source models, and thus easier to use.
As urban land surface energy balance has four main components, i.e., vegetation, soil, impervious surfaces, and water [30], dual-source models only take vegetation and soil into account, and the proportions of vegetation and soil in a mixed pixel are always ignored, and thus are not suitable for urban areas.In order to solve the problem, Zheng [5] first developed a multi-source parallel model in which evapotranspiration of each pixel in a satellite image of urban land surface was inferred based on an energy residual algorithm through combining the transpiration and evaporation of the vegetation component and the soil component that are extracted by linear spectral mixture analysis.The multi-source parallel model (hereinafter referred to as Zheng's model) has improved evapotranspiration inversion accuracy for urban areas [5].However, Zheng's model ignored the spectral heterogeneity of urban impervious surfaces in linear spectral mixture analysis, and used general empirical algorithms or coefficients for calculating momentum roughness length, heat roughness length, and component net radiant flux, and failed to consider the characteristic of each component in the study area.As such, Zheng's model can be improved by optimizing its parameterization and simplifying part of its algorithm.Based on Zheng's model, this study presents a modified multi-source parallel model (hereinafter referred to as the modified model or our model).We aim to improve the accuracy of estimating urban surface evapotranspiration through our model by: (i) Characterizing spectrally heterogeneous urban impervious surfaces with two spectral endmembers (high-and low-albedo); (ii) Revising the methods of deriving roughness length for each land surface component; (iii) Recalculating the component net radiant flux with a full consideration of the fraction and the characteristic of each land surface component.

Study Area
Hefei (31 • 40 -31 • 59 N, 117 • 4 -117 • 26 E), the capital city of the eastern Chinese province of Anhui, was selected for this study.Hefei is located on the western side of the Yangtze River Delta between the Yangtze River and Huaihe River.It covers an area of 11,445.1 km 2 , including the 770 km 2 area of Chaohu Lake, one of the five largest freshwater lakes in China.Situated at mid-latitudes, Hefei has a subtropical humid monsoon climate with four distinguishable seasons and large precipitation and insolation, with an average annual temperature of 15.7 • C, an average annual rainfall of 1000 mm, an average annual sunshine period of 2000 h and an average annual frost-free period of 228 days [31].Hefei is generally flat and characterized by low elevations ranging from 20 m to 40 m.
By the end of 2015, the urban population of Hefei reached 5.484 million, with an urbanization rate of 70.4%.This suggests that Hefei has become a highly urbanized city.In order to focus on urban land surface evapotranspiration, we extracted the urban built-up area from Hefei's HJ-1A satellite imagery (presented in Section 2.2) as the study area (Figure 1).The selected area highlights the four main land cover types of Hefei, i.e., impervious surfaces, vegetation, bare soil, and water.
Remote Sens. 2017, 9, x FOR PEER REVIEW 3 of 32 Based on Zheng's model, this study presents a modified multi-source parallel model (hereinafter referred to as the modified model or our model).We aim to improve the accuracy of estimating urban surface evapotranspiration through our model by: (i) Characterizing spectrally heterogeneous urban impervious surfaces with two spectral endmembers (high-and low-albedo); (ii) Revising the methods of deriving roughness length for each land surface component; (iii) Recalculating the component net radiant flux with a full consideration of the fraction and the characteristic of each land surface component.

Study Area
Hefei (31°40′-31°59′N, 117°4′-117°26′E), the capital city of the eastern Chinese province of Anhui, was selected for this study.Hefei is located on the western side of the Yangtze River Delta between the Yangtze River and Huaihe River.It covers an area of 11,445.1 km 2 , including the 770 km 2 area of Chaohu Lake, one of the five largest freshwater lakes in China.Situated at mid-latitudes, Hefei has a subtropical humid monsoon climate with four distinguishable seasons and large precipitation and insolation, with an average annual temperature of 15.7 °C, an average annual rainfall of 1000 mm, an average annual sunshine period of 2000 hours and an average annual frost-free period of 228 days [31].Hefei is generally flat and characterized by low elevations ranging from 20 m to 40 m.
By the end of 2015, the urban population of Hefei reached 5.484 million, with an urbanization rate of 70.4%.This suggests that Hefei has become a highly urbanized city.In order to focus on urban land surface evapotranspiration, we extracted the urban built-up area from Hefei's HJ-1A satellite imagery (presented in Section 2.2) as the study area (Figure 1).The selected area highlights the four main land cover types of Hefei, i.e., impervious surfaces, vegetation, bare soil, and water.

Data and Preprocessing
In order to inverse urban land surface ET through a multi-source model, satellite image data must have at least three different thermal infrared (TIR) bands in the wavelength range of 8-13 μm-

Data and Preprocessing
In order to inverse urban land surface ET through a multi-source model, satellite image data must have at least three different thermal infrared (TIR) bands in the wavelength range of 8-13 µm-ASTER (carried by the satellite of Terra) and MODIS (carried by the two satellites of Terra and Aqua) are the only two sensors producing such freely available data (note that ASTER data are free as of 1 April 2016) [32].ASTER has a lower time resolution but a higher spatial resolution than MODIS (16 days vs. four days; 90 m vs. 1 km).Since our focus was on an urban built-up area with heterogeneous land cover types and a higher spatial resolution can result in a finer ET map, we decided to select the 90-meter-resolution ASTER for inferring component temperatures.An ASTER L1B level image of Hefei acquired on 13 October 2013 was downloaded from MADAS (METI AIST Data Archive System) of the Geological Survey of Japan for the study.
ASTER SWIR detectors malfunctioned as of April 2008, and thus no usable ASTER SWIR data have been provided since then [32].Although the remaining three VNIR bands can be used to decompose mixed pixels through spectral unmixing, they do not suffice for more than three endmembers due to inter-band correlations.As a result, HJ-1A CCD2 data at 30 m spatial resolution, acquired on 13 October 2013, were used to provide visible and near-infrared (VNIR) bands for decomposing mixed pixels.Technical specifications of ASTER TIR and HJ-1A CCD2 VNIR bands are shown in Table 1.The ASTER image was georeferenced to the HJ-1A image with an overall error of less than half a pixel.An FLAASH based atmospheric correction was applied to the HJ-1A CCD bands as they were used for linear spectral mixture analysis.All ASTER thermal bands were then spatially resampled to 30 m to match the HJ-1A CCD bands.In this study, the four VNIR bands of HJ-1A CCD (band 1 to 4) and the four TIR bands of ASTER (band 11-14) were required for decomposing mixed pixels and inferring component temperatures respectively.
In the absence of field-based latent heat radiation data of the study area, the 1-km resolution MODIS global evapotranspiration product (hereinafter referred to as the MOD16 ET product) [33] was firstly considered as the validation data for assessing the modeled ET because it has a good correlation with flux tower data (r = 0.86) [34].However, the data of the study area on 13 October 2013 (i.e., the 286th day of 2013) is missing from the MOD16 ET product due to its coarse temporal resolution (eight days).Another freely available ET dataset is the 1 km/daily evapotranspiration product over China and the Association of Southeast Asian Nations (ASEAN) for 2013 (hereinafter referred to as the China-ASEAN ET product), provided by the Chinese Academy of Sciences [35].The China-ASEAN ET product is distributed in 40 adjacent non-overlapping sinusoidal tiles that are approximately 10 • × 10 • (at the equator).It has a higher temporal resolution (one day) but an equal spatial resolution (1 km) in relation to the MOD16 ET product [35].
Because the accuracy of the China-ASEAN ET product has not been reported by its provider, we here assessed it by comparison with the MOD16 ET product of October 8, 2013 (i.e., the 281st day of 2013).As the ET data of all urban area were excluded in the MOD16 ET product, a rectangular test site Remote Sens. 2017, 9, 1029 5 of 33 consisting of 81250 pixels (~1 km × 1 km) near the study area was used to perform the comparative analysis.A total of 1000 sample points were randomly generated, using ArcGIS, in the rectangular test site with a minimum allowed distance of 1 km.Values of the pixels containing the 1000 sample points were extracted from the two ET datasets and compared in Figure 2. The observations are nearly evenly distributed around the 1:1 line with a correlation coefficient of 0.7253.This strong correlation suggests [36,37] that the China-ASEAN ET product has a reliable quality and thus can be used as alternative validation data for assessing our model's accuracy.
Remote Sens. 2017, 9, x FOR PEER REVIEW 5 of 32 nearly evenly distributed around the 1:1 line with a correlation coefficient of 0.7253.This strong correlation suggests [36,37] that the China-ASEAN ET product has a reliable quality and thus can be used as alternative validation data for assessing our model's accuracy.

Methodology
In contrast to natural surfaces, urban built-up areas have uneven and sparse vegetation distribution.As a result, the coupling between soil evaporation and vegetation transpiration is weak and ignorable, which enables parallel models to estimate surface evapotranspiration [23,38].However, traditional parallel models only take two components (vegetation and soil) into account and neglect the important role that impervious surface plays in surface energy balance [5].In lowresolution satellite imagery, a pixel is usually a mixture of different land cover types in urban areas.It is, therefore, necessary to consider different components in the pixel when estimating evapotranspiration in urban areas.Multi-source models, however, take vegetation, soil, impervious surfaces (water can be excluded and its evapotranspiration can be calculated separately) as the three components for urban mixed pixels.The urban land surface energy balance relationship of the net radiant flux, the sensible heat flux, the latent heat flux and the internal heat flux for each pixel can be considered as a combination of the energy balance of the three components, which could be expressed as follows [5]:

Methodology
In contrast to natural surfaces, urban built-up areas have uneven and sparse vegetation distribution.As a result, the coupling between soil evaporation and vegetation transpiration is weak and ignorable, which enables parallel models to estimate surface evapotranspiration [23,38].However, traditional parallel models only take two components (vegetation and soil) into account and neglect the important role that impervious surface plays in surface energy balance [5].In low-resolution satellite imagery, a pixel is usually a mixture of different land cover types in urban areas.It is, therefore, necessary to consider different components in the pixel when estimating evapotranspiration in urban areas.Multi-source models, however, take vegetation, soil, impervious surfaces (water can be excluded and its evapotranspiration can be calculated separately) as the three components for urban mixed pixels.The urban land surface energy balance relationship of the net radiant flux, the sensible heat flux, the latent heat flux and the internal heat flux for each pixel can be considered as a combination of the energy balance of the three components, which could be expressed as follows [5]: where R n,v , R n,s , and R n,imp refer respectively to the net radiant flux of vegetation, soil and impervious surface; H v , H s and H imp refer respectively to the sensible heat flux of vegetation, soil and impervious surface; LE v and LE s refer respectively to the latent heat flux of vegetation and soil (evapotranspiration does not apply to impervious surfaces except in rainy weathers); G s refers to soil heat flux, and G imp refers to the heat flux inside impervious surface.
Based on Equations ( 1)-( 3), the total latent heat flux (LE) (i.e., the sum of vegetation latent heat flux LE v and soil latent heat flux LE s ) can be calculated if the net radiant flux (R n ), the sensible heat flux (H) and the internal heat flux (G) of each component are calculated.But prior to that, two important intermediate parameters be, namely pixel component fraction ( f k , k refers to component type, k = 1, 2, 3, 4 refer respectively to vegetation, soil, high-and low-albedo impervious surfaces) and pixel component temperature (T k ), should be acquired first.
In order to better illustrate the methodology of this study, three flowcharts detailing the input data, parameters, and equations are given in different sections of the paper.The flowchart for the first step of our modified model (Sections 3.1 and 3.2) is shown in Figure 3.This figure demonstrates estimation of component temperatures from remotely sensed imagery.

Linear Spectral Mixture Analysis
In a linear spectral mixture analysis (LSMA) of an urban area, the spectrum of a pixel is considered as a linear combination of spectra of pure endmembers within the pixel weighted by their component abundance [30,39].The fully constrained linear spectral mixture analysis where the sum of all components is one and no component is negative [39,40] was applied to inverse pixel component fraction using a least squares method.The fully constrained LSMA is expressed by the following equations: where is the reflectance for each band in a pixel; , is the reflectance of endmember in band for that pixel; is the fraction of endmember in a pixel; is the residual; and is the number of components.
As it adversely affects impervious surface estimation [30,41] and, in the form of open water bodies, is hardly mixed with other components, water was masked through land cover classification before endmember selection [30].Endmember selection is a primary step of LSMA.Brightness  A1 in the Appendix A.

Linear Spectral Mixture Analysis
In a linear spectral mixture analysis (LSMA) of an urban area, the spectrum of a pixel is considered as a linear combination of spectra of pure endmembers within the pixel weighted by their component abundance [30,39].The fully constrained linear spectral mixture analysis where the sum of all components is one and no component is negative [39,40] was applied to inverse pixel component fraction using a least squares method.The fully constrained LSMA is expressed by the following equations: where R b is the reflectance for each band b in a pixel; R k,b is the reflectance of endmember k in band b for that pixel; f k is the fraction of endmember k in a pixel; e b is the residual; and N is the number of components.
As it adversely affects impervious surface estimation [30,41] and, in the form of open water bodies, is hardly mixed with other components, water was masked through land cover classification before endmember selection [30].Endmember selection is a primary step of LSMA.Brightness normalization [42][43][44] was first applied to the HJ-1A CCD bands to decrease intra-class spectral variability, followed by a minimum noise fraction (MNF) transform to reduce data dimensionality [45].After that, the pixel purity index (PPI) was calculated to extract the pure pixels [46,47], then the pure pixels were loaded into the interactive N-Dimensional visualizer of software package ENVI.By trial and error, four different endmembers were determined in the pixel clouds: high-albedo feature (e.g., cement concrete, glass, and composite materials), low-albedo feature (e.g., asphalt road and asphalt roof), vegetation (e.g., grassland and tree crown), soil (e.g., bare soil in construction site and unused wasteland).Figure 4 shows the normalized spectral reflectance for HJ-1 VNIR bands of the four selected endmembers.The fraction of each endmember in a pixel was finally estimated through LSMA.
Remote Sens. 2017, 9, x FOR PEER REVIEW 7 of 32 pure pixels were loaded into the interactive N-Dimensional visualizer of software package ENVI.By trial and error, four different endmembers were determined in the pixel clouds: high-albedo feature (e.g., cement concrete, glass, and composite materials), low-albedo feature (e.g., asphalt road and asphalt roof), vegetation (e.g., grassland and tree crown), soil (e.g., bare soil in construction site and unused wasteland).Figure 4 shows the normalized spectral reflectance for HJ-1 VNIR bands of the four selected endmembers.The fraction of each endmember in a pixel was finally estimated through LSMA.Surface component fractions are fundamental model parameters of our model and should be assessed before they are used for calculating other parameters.Statistical analysis of the RMSE image (Figure 5) resulting from the fully constrained LSMA showed that the average RMSE was 0.006.Despite a maximum value of 0.12, approximately 99% of pixels on the RSME image had values ranging from 0 to 0.03.This suggests that the accuracy of the FLCS result was high.Surface component fractions are fundamental model parameters of our model and should be assessed before they are used for calculating other parameters.Statistical analysis of the RMSE image (Figure 5) resulting from the fully constrained LSMA showed that the average RMSE was 0.006.Despite a maximum value of 0.12, approximately 99% of pixels on the RSME image had values ranging from 0 to 0.03.This suggests that the accuracy of the FLCS result was high.
To further validate the accuracy of surface component fractions, a 2 m resolution satellite image on 2 October 2013 of the study area acquired from Google Earth as validation data of the spectral unmixing.A total of 100 validation samples, each consisting of 3 × 3 pixels (90 × 90 m 2 ) of the HJ-1A image, were randomly generated on the Google Earth image (shown in Figure 6), and four surface component fractions of each validation sample were extracted by manual interpretation.
Surface component fractions are fundamental model parameters of our model and should be assessed before they are used for calculating other parameters.Statistical analysis of the RMSE image (Figure 5) resulting from the fully constrained LSMA showed that the average RMSE was 0.006.Despite a maximum value of 0.12, approximately 99% of pixels on the RSME image had values ranging from 0 to 0.03.This suggests that the accuracy of the FLCS result was high.Despite four surface component fractions, we here chose only vegetation and soil fractions for accuracy assessment for the ease of assessment and, more importantly, because only they were used for calculating other parameters.Scatter plots in Figure 7 shows high correlation coefficients (both approximately 0.92) and low RMSEs (both only slightly higher than 0.10) for the two different datasets of surface component fractions.Figure 7 indicates that LSMA produced good surface component fractions, which are acceptable and allow the following calculations.Despite four surface component fractions, we here chose only vegetation and soil fractions for accuracy assessment for the ease of assessment and, more importantly, because only they were used for calculating other parameters.Scatter plots in Figure 7 shows high correlation coefficients (both approximately 0.92) and low RMSEs (both only slightly higher than 0.10) for the two different datasets of surface component fractions.Figure 7 indicates that LSMA produced good surface component fractions, which are acceptable and allow the following calculations.
Despite four surface component fractions, we here chose only vegetation and soil fractions for accuracy assessment for the ease of assessment and, more importantly, because only they were used for calculating other parameters.Scatter plots in Figure 7 shows high correlation coefficients (both approximately 0.92) and low RMSEs (both only slightly higher than 0.10) for the two different datasets of surface component fractions.Figure 7 indicates that LSMA produced good surface component fractions, which are acceptable and allow the following calculations.

Retrieval of Land Surface Component Temperature
According to Planck's law, the thermal radiation of isothermal surface is related to the blackbody radiation, which can be expressed as [48]: where L λ is the thermal radiation of land surface at wavelength λ; ε λ is the land surface emissivity at wavelength λ; B λ (T sur ) is the Planck blackbody radiation when the land surface temperature is T sur .Qin et al. [49] have found the following relationship between B λ (T sur ) and at-sensor thermal radiation B λ,sens : where τ λ is atmospheric transmittance at wavelength λ, T air is near-surface average atmospheric temperature.The thermal radiation of each mixture pixel can be considered as a linear combination of the thermal radiation of each component [50]: where ε λk is the emissivity of component k at wavelength λ; B λ (T k ) is the Planck blackbody radiation of component k at wavelength λ when the temperature of component k is T k .The Planck blackbody radiation shows an approximately linear relationship with temperature in the range between 273.15 K and 330.15K [51] (Figure 8 and Table 2).
where is the emissivity of component at wavelength ; ( ) is the Planck blackbody radiation of component at wavelength when the temperature of component is .The Planck blackbody radiation shows an approximately linear relationship with temperature in the range between 273.15 K and 330.15K [51] (Figure 8 and Table 2).Based on Table 2 and Equations ( 6)-( 8), the following equation was obtained: Based on Table 2 and Equations ( 6)-( 8), the following equation was obtained: where f k have already been obtained in Section 3.1, when ε kλ , ε λ , B λ,sens , τ λ and B λ (T air ) have been obtained in the next calculations, the land surface component temperatures T k would be the four unknowns.Based on the four equations established with the four ASTER thermal infrared bands, optimal surface component temperatures T k were solved using a least squares method.

Average Land Surface Emissivity
For a pixel mixed only with vegetation and soil, Mao et al. [52] have proposed that average pixel emissivity can be obtained based on the relationship between emission and vegetation fraction.As two different impervious surface endmembers (high-and low-albedo), in addition to vegetation and soil, were identified in this study, we assume that average surface emissivity can be given by: where ε λ is average land surface emissivity; f v , f s , f imp_h and f imp_l are the fraction of vegetation, soil, high-and low-albedo imperious surfaces respectively; ε λv , ε λs , ε λi_h and ε λi_l are the emissivity of each component at wavelength λ respectively; R v , R s and R m are the temperature ratio [53] of vegetation, soil and man-made construction with R k = (T k /T sur ) 4 respectively.As it has a good linear relationship with vegetation coverage P v , which is calculated through NDV I, the temperature ratio for each of the components is given by: R v = 0.9332 where NDV I v and NDV I s , which are the NDVIs for dense vegetation and bare soil, respectively.
In the context of our study area, NDV I v and NDV I s were given approximate values of 0.65 and 0.05, respectively.When radiation makes its way to a surface, the sum of reflectivity, transmissivity, and emissivity for the surface is 100% [54].The transmissivity is 0% for an opaque surface-a surface component's emissivity can be calculated if its reflectivity is given.We investigated the typical urban ground objects in Hefei through a field survey conducted on 2 October 2016 and extracted their emissivity from the ASTER Spectral Library (Version 2.0) [55]: three types of vegetation (conifer, deciduous and grass), two types kinds of soil (light yellowish brown clay and light yellowish brown loam), four types of construction concrete (concrete paving solid 0092uuu, 0397uuu, 0424uuu and construction cement solid 0432uuu) and two types of construction asphalt (concrete paving solid 0095uuu and 0096uuu).The extracted emissivity was averaged to approximate the emissivity of the four endmembers used in the study, i.e., vegetation, soil, high-and low-albedo impervious surfaces (Figure 9 and Table 3).Since the geographical size of the study area is small and the image used was extracted from a high-quality ASTER scene, the atmospheric transmittance was assumed to be homogenous across the study area.As such, the atmospheric transmittance can be taken as a constant for the study area.Mao et al. [56] have noted that there is a good linear relationship between atmospheric transmittance and atmospheric vapor content (Table 4).Since the geographical size of the study area is small and the image used was extracted from a high-quality ASTER scene, the atmospheric transmittance was assumed to be homogenous across the study area.As such, the atmospheric transmittance can be taken as a constant for the study area.Mao et al. [56] have noted that there is a good linear relationship between atmospheric transmittance and atmospheric vapor content (Table 4).Since atmospheric vapor content of the study area was not directly observed by meteorological stations, Yang and Qiu [57] have proposed an algorithm for calculating atmospheric vapor content W of Chinese areas by using average atmospheric water vapor pressure data: where e is the average atmospheric water vapor, which can be obtained through the Daily Dataset of China Surface Climate [58] (detailed in Table 5); ϕ = 31.04• N, is the latitude of the center of the study area; and E 0 = 200 m, the average elevation of the study area.This equation is applicable throughout China, except for part of south China [57].

Retrieval of Land Surface Component Sensible Heat Flux
The flowchart for the second step of our modified model is shown in Figure 10.This figure demonstrates the estimation of component sensible heat flux followed by Step 1 in Figure 3.
In our modified multi-source model, land surface sensible heat flux is considered as the sum of the component sensible heat flux for mixed pixels.The component sensible heat flux can be described with a bulk transfer equation [59]: where H k is the component sensible heat flux; ρ is the air density (kg•m −3 ); C p is the heat capacity of air at constant pressure (J•kg −1

Retrieval of Land Surface Component Sensible Heat Flux
The flowchart for the second step of our modified model is shown in Figure 10.This figure demonstrates the estimation of component sensible heat flux followed by Step 1 in Figure 3.In our modified multi-source model, land surface sensible heat flux is considered as the sum of the component sensible heat flux for mixed pixels.The component sensible heat flux can be described with a bulk transfer equation [59]: where is the component sensible heat flux; is the air density (kg•m −3 ); is the heat capacity of air at constant pressure (J•kg −1 •K −1 ); is the air temperature (K), which can be obtained from a

Aerodynamic Resistance to Heat Transfer
The aerodynamic resistance to heat transfer of each component can be estimated separately as follows [60][61][62]: For stable conditions: for unstable conditions: where Z is the elevation at which wind speed is observed (the elevation of Hefei meteorological observatory is 27 m); d is the zero displacement height; Z 0m and Z 0h are the roughness lengths for momentum and for heat, respectively; k c is von Karman's constant with k c = 0.41; u z is the wind speed (m/s); ψ m and ψ h are stability correction functions for momentum and heat respectively; and L M is the Monin-Obukhov length with Z/L M > 0 being for stable conditions and Z/L M < 0 being for unstable conditions.L M can be estimated as follows [63]: where g is the gravitational acceleration with g = 9.8 m•s −2 ; T * is temperature scale with T * = T sur − T air ; u * is friction velocity estimated by [62]: 3.3.2.Momentum Roughness Length, Heat Roughness Length, and Zero Displacement Height Roughness length and zero displacement height are the most important parameters for characterizing the aerodynamics of land surface.Compared with natural land surfaces, the structure of urban surfaces is more complex.For a mixed pixel consisting of several land cover types, it is therefore less accurate to use a single value to represent roughness length or zero displacement height.As such, vegetation, soil, and impervious surfaces (high-and low-albedo) were calculated separately for roughness length and zero displacement in this study.

• Roughness length and zero displacement of vegetation
For vegetated surfaces, the relationship between average vegetation height h and Z 0m and d is given by [64]: where the average vegetation height h was previously evaluated at 5 m for urban vegetation [65], and Z 0h is given by [64]: • Roughness length and zero displacement of impervious surfaces For impervious surfaces, Grimmond et al. [66] have proposed empirical equations for calculating Z 0m and d: where the average building height h was evaluated at 20 m [67].The thermal roughness length was defined in previous studies [68,69] for bluff-rough situations: where Re * is the roughness Reynolds number, with a kinematic molecular viscosity v = 1.461 × 10 −5 m•s −1 .
• Roughness length and zero displacement of soil surfaces For bare soil surfaces, the values of momentum roughness length and zero displacement were obtained from Liu's measurements [70]: Stewart et al. [71] have found a mean value of 4.5 for KB −1 (i.e., ln(Z 0m /Z 0h ) [64]) for bare soil, so the thermal roughness length can be calculated as: of the land surface.The effective atmospheric emissivity and the average atmospheric temperature are only affected by the atmosphere.This indicates that only land surface albedo , land surface emissivity and land surface temperature are related to the characteristics of the land surface.Therefore, for one single pixel, its surface net radiant flux for four surface components is given by:   Land surface net radiant flux is the difference between incoming solar radiation to a land surface and outgoing radiation from the land surface [72].The incoming radiation is mainly from solar radiation at shortwave lengths and atmospheric reflection received by the land surface, the outgoing radiation is the energy emitted into the atmosphere from the land surface itself.Therefore, land surface net radiant flux R n is given by [73]: where R s ↓ and R L ↓ are solar shortwave radiation (i.e., incident solar radiation on ground Q [74]) and longwave radiation emitted from the atmosphere to the land surface respectively; R s ↑ and R L ↑ are shortwave radiation reflected by the land surface and longwave radiation emitted by the land surface.In Equation (37), the difference between shortwave radiations can be calculated as (1 − a)Q, and the difference between longwave radiations can be calculated as ε air σT 4 air − ε λ σT 4 sur [74].Therefore, the equation of calculating land surface net radiant flux can be also expressed as [18]: where a is the albedo of the land surface; ε air is the effective emissivity of the atmosphere; and σ is the Stefan Boltzmann constant with σ = 5.67 In Equation (38), the incident solar radiation on ground Q is not affected by the characteristics of the land surface.The effective atmospheric emissivity ε air and the average atmospheric temperature T air are only affected by the atmosphere.This indicates that only land surface albedo a, land surface emissivity ε λ and land surface temperature T sur are related to the characteristics of the land surface.Therefore, for one single pixel, its surface net radiant flux for four surface components is given by: where R n,v , R n,s , R n,imp_h , and R n,imp_l are component surface net radiant flux of vegetation, soil, high-and low-albedo impervious surfaces; a v , a s , a imp_h , and a imp_l are the surface albedo of these components.

Incident Solar Radiation on Ground
Since all information extracted from satellite image are instantaneous data, the instantaneous incident solar radiation on the ground should be used in the calculation as follows [75][76][77]: where Q 0 is instantaneous astronomical solar radiation; P m is atmospheric transparency; and AM is the air mass that defines the direct optical path length through the Earth's atmosphere.For the mid-latitude area, AM can be estimated as 1.5 [78,79].

• Instantaneous astronomical solar radiation
Because the height of the sun varies with time, astronomical solar radiation radiated on the horizontal plane of upper atmosphere can be computed as [80]: where I sc is the solar constant with I sc = 1367 w•m −2 ; d m is the correction coefficient of sun-earth distance; δ is the solar declination (rad); ω is solar hour angle (rad); d m , ω and δ are given by: where θ is day angle (rad); ST is real solar time (hour).θ and ST are given by [81]: where D n is the number of days in the year; h b is Beijing time (UTC+8); λ * s is the longitude of local standard time (the longitude of Beijing is 120 • E); λ * is local longitude (the central longitude of the study area is 117.3 • E); η is the time lag (rad), it can be calculated as: η = 0.000043 + 0.002061 cos θ − 0.032040 sin θ − 0.014974 cos 2θ − 0.040685 sin 2θ. (50)

• Atmospheric Transparency
Atmospheric transparency P m is the ratio of incident solar radiation on ground Q to astronomical solar radiation Q 0 [75].Since it is challenging to acquire instantaneous atmospheric transparency data, the daily radiation data from ground radiation observation station in the study area was used to estimate P m by the following equation [75][76][77]: where Q d is the daily incident solar radiation on ground, the value of Q d from Hefei radiation observation station on 13 October 2013 is 1887 MJ•m −2 •d −1 [58]; Q 0d is daily astronomical solar radiation, which can be calculated as: where T d is the time in one day in seconds, i.e., 86400 seconds; ω 0 and −ω 0 are the solar hour angles at the time of sunrise and sunset, respectively; and ω 0 can be calculated as:

Atmosphere Effective Emissivity
Atmosphere effective emissivity ε air is the function of average atmospheric water vapor pressure e and average atmospheric temperature T air , under the cloudless condition, ε air can be estimated as [82,83]: (54)

Land Surface Albedo
Liu et al. [84] have proposed an algorithm of land surface albedo for HJ-1 A/B CCD data, which transforms narrow band albedo to broad band albedo through regression analysis.It is described as follows: where a HJ is the broad band albedo; band i (i = 1, 2, 3, 4) is the reflectance of each CCD band.Since a HJ is average albedo of one pixel, for obtaining the albedo of each typical component, endmember pixels (i.e., pure pixels) of each component was used to calculate the average albedo of the four typical land cover types (Table 6).

Retrieval of Land Surface Component Internal Heat Flux
For soil component, soil heat flux is the heat exchange between soil and the atmosphere.Friedl [85] established a proportional relationship between soil heat flux G s and soil surface net radiant flux: So far there is no widely-used empirical equation for estimating the internal heat flux for impervious surfaces (high-and low-albedo) but nevertheless, an algorithm based on surface energy balance for each component can be used to approximate internal heat flux of impervious surfaces: where G imp_h , R n,imp_h and H imp_h are the internal heat flux, net radiant flux and sensible heat flux of high-albedo impervious surfaces, respectively; G imp_l , R n,imp_l and H imp_l are the internal heat flux, net radiant flux and sensible heat flux of low-albedo impervious surfaces, respectively.The results of sensible heat flux and net radiant flux calculation have already been presented in Sections 3.3 and 3.4.

Retrieval of Daily Evapotranspiration
After all parameters of the modified model were determined, the latent heat flux LE was obtained through adding together the component latent heat flux of vegetation and soil.Latent heat flux LE can be expressed as a product of evapotranspiration (ET, kg•m −2 •s −1 or mm•s −1 ) and latent heat of vaporization (L, L = (2.501− 0.00237 × T air ) × 10 6 , J•kg −1 ) [86].For transforming instantaneous evapotranspiration (ET i ) into daily evapotranspiration (ET d ), Xie's [87] research has established the relationship between ET d and ET i of any moment, which is expressed as: where t is the time interval from sunrise to the acquisition time of satellite imagery (the local sunrise time on 13 October 2013 in Hefei was 06:12 AM Beijing Time (UTC+08:00) and N E is the evapotranspiration duration-which is two hours shorter than the insolation duration.

Results
Following the steps in Section 3 and the flowcharts (Figures 3, 10 and 11), we acquired the daily evapotranspiration of the study area on 13 October 2013 (Figure 12).Compared with the 2 m resolution satellite images of the study area, it was observed that the high ET values are mostly located in densely vegetated areas (e.g., green belts, city parks) and areas currently being developed, while low ET values are mainly distributed in intensive building areas, such as commercial and residential areas in downtowns, industrial areas in suburbs, and roads of contrasting levels.
Following the steps in Section 3 and the flowcharts (Figures 3, 10 and 11), we acquired the daily evapotranspiration of the study area on 13 October 2013 (Figure 12).Compared with the 2 m resolution satellite images of the study area, it was observed that the high ET values are mostly located in densely vegetated areas (e.g., green belts, city parks) and areas currently being developed, while low ET values are mainly distributed in intensive building areas, such as commercial and residential areas in downtowns, industrial areas in suburbs, and roads of contrasting levels.In order to make the result comparable with the China-ASEAN ET product, the image of daily evapotranspiration (Figure 12) was re-sampled to 1 km using the bilinear interpolation method.Scatter plots were used to compare the results and validation data (the China-ASEAN ET product) based on 552 sample points (Figure 13).The sample points could be linearly fitted with an R 2 (R 2 ) value of 0.5806 (Figure 13a).Such linear fitting varies with land cover type (vegetation, soil and impervious surfaces) (Figure 13b).The fraction image of each endmember was re-sampled to 1 km with pixels, with fractions >50% being selected as validation samples for each land cover type.Since the focus here was to explore the accuracy of the modeled ET for basic urban land cover types, highand low-albedo impervious surfaces were combined as impervious surfaces in the following analysis.In order to make the result comparable with the China-ASEAN ET product, the image of daily evapotranspiration (Figure 12) was re-sampled to 1 km using the bilinear interpolation method.Scatter plots were used to compare the results and validation data (the China-ASEAN ET product) based on 552 sample points (Figure 13).The sample points could be linearly fitted with an R 2 (R 2 ) value of 0.5806 (Figure 13a).Such linear fitting varies with land cover type (vegetation, soil and impervious surfaces) (Figure 13b).The fraction image of each endmember was re-sampled to 1 km with pixels, with fractions >50% being selected as validation samples for each land cover type.Since the focus here was to explore the accuracy of the modeled ET for basic urban land cover types, high-and low-albedo impervious surfaces were combined as impervious surfaces in the following analysis.The relationship between the two ET datasets for each land cover type was found to be linear, with each correlation coefficient >0.7.The linear fitting models are shown in Table 7.The R 2 value for the linear relationship between our modeled ET and the validation ET in vegetation-dominated areas was 0.7495 (124 sampled points), obviously higher than that in soil-dominated areas (130 sampled points, R 2 = 0.5349) and impervious surface-dominated areas (102 sampled points, R 2 = 0.5200).However, the lowest MRE and RMSE values were obtained for the impervious surfaces.
The relationship between the two ET datasets for each land cover type was found to be linear, with each correlation coefficient >0.7.The linear fitting models are shown in Table 7.The R 2 value for the linear relationship between our modeled ET and the validation ET in vegetation-dominated areas was 0.7495 (124 sampled points), obviously higher than that in soil-dominated areas (130 sampled points, R = 0.5349) and impervious surface-dominated areas (102 sampled points, R = 0.5200).However, the lowest MRE and RMSE values were obtained for the impervious surfaces.

Discussion
Estimating urban land surface evapotranspiration did not receive sufficient attention in the community of quantitative remote sensing until the emergence of the multi-source parallel model proposed by Zheng [5].Zheng's model allows evapotranspiration to be estimated but its accuracy remains to be further improved.Based on his model, we therefore present a modified multi-source parallel model using ASTER thermal infrared bands.
Important improvements are discussed below.

Using HJ-1A CCD VNIR Data Alternative to ASTER SWIR Data
Component temperature is an important intermediate parameter of evapotranspiration estimation in the multi-source parallel model.In order to infer component temperatures of the four endmembers (vegetation, soil and high-and low-albedo impervious surface) in urban surface, remotely sensed data with at least four thermal infrared bands are required.As such, ASTER data should have been considered as most appropriate for such purpose-but its SWIR detectors have failed since April 2008 [32].Instead, Zheng [5] made use of Landsat TM data to estimate endmember fraction through linear spectral mixture analysis.Despite the same revisit cycle of Terra and Landsat satellites (16 days), it is difficult to acquire satellite images of the same location at the same overpass time.The CCD sensor mounted on the HJ-1A satellite that was launched together with the HJ-1B on

Discussion
Estimating urban land surface evapotranspiration did not receive sufficient attention in the community of quantitative remote sensing until the emergence of the multi-source parallel model proposed by Zheng [5].Zheng's model allows evapotranspiration to be estimated but its accuracy remains to be further improved.Based on his model, we therefore present a modified multi-source parallel model using ASTER thermal infrared bands.
Important improvements are discussed below.

Using HJ-1A CCD VNIR Data Alternative to ASTER SWIR Data
Component temperature is an important intermediate parameter of evapotranspiration estimation in the multi-source parallel model.In order to infer component temperatures of the four endmembers (vegetation, soil and high-and low-albedo impervious surface) in urban surface, remotely sensed data with at least four thermal infrared bands are required.As such, ASTER data should have been considered as most appropriate for such purpose-but its SWIR detectors have failed since April 2008 [32].Instead, Zheng [5] made use of Landsat TM data to estimate endmember fraction through linear spectral mixture analysis.Despite the same revisit cycle of Terra and Landsat satellites (16 days), it is difficult to acquire satellite images of the same location at the same overpass time.The CCD sensor mounted on the HJ-1A satellite that was launched together with the HJ-1B on 6 September 2008 has, however, offered a solution to such a problem by imaging land surfaces with its four spectral bands (blue, green, red, and near-infrared) at 30 m resolution with a two-day revisit cycle [88].Therefore, when the multi-source parallel model is applied for time series analysis, HJ-1A CCD data can be used as preferable data for spectral unmixing.

Optimized Endmembers of Urban Land Surface
Zheng's model has extracted vegetation, soil and impervious surfaces for urban land surface, among which impervious surfaces were considered as a single endmember.This might be not an appropriate practice for urban areas as impervious surfaces are often spectrally heterogeneous [30].If the heterogeneity of impervious surfaces was ignored, the accuracy of mixed pixel decomposition would be reduced, thus adversely affecting the estimation of component energy flux of urban land surfaces.In our modified model, pure impervious surface pixels were therefore divided into highand low-albedo endmembers [30,89].Impervious surfaces' fraction in each pixel was the sum of the fractions of high-and low-albedo endmembers [30].

Improved Parameterization of Z om and Z oh
For the calculation of Z 0m and Z 0h , Zheng [5] used typical land surface roughness for vegetation, soil, and impervious surfaces, that is, constant values for each component: Z 0m,veg = 100Z 0h,veg = 0.1 (vegetation), Z 0m,soil = 50Z 0h,soil = 0.001(soil), and Z 0m,imp = 1000Z 0h,imp = 1.5 (impervious surfaces).Our modified model, however, used spatially specific Z 0m and Z 0h values for each component through the equations proposed by Kustas et al. [64,65] for Z 0m,veg and Z 0h,veg , the equations proposed by Grimmond et al. [66] and Brutsaert et al. [65,66] for Z 0m,soil and Z 0h,soil , the equations proposed by Liu et al. [67] and Stewart et al. [71] for Z 0m,imp and Z 0h,imp .As these equations take the spectral characteristics of different land surfaces into account, the results of would consequently be more reliable.

Algorithm Improvement of Component Net Radiant Flux
Regarding Equation ( 46), Zheng [5] assumed that it was challenging to obtain each component's albedo.Instead, he estimated vegetation net radiant flux R n,v and soil net radiant flux R n,s through the empirical relationships between total net radiant flux R n , leaf area index (LAI) and day angle θ [23], and obtained the net radiant flux of impervious surface R n,imp by R n − R n,v − R n,s .As a general empirical algorithm, Zheng's component net radiant flux calculation has not considered the actual characteristics of each component in the study area.Therefore, the result might be less accurate.We have proposed a new algorithm for calculating component net radiant flux.Thanks to the ease of deriving the broad albedo of each pixel, the average albedo of pure pixels of each component can be considered as the component albedo, which was consistent with the surface parameter values proposed by De Ridder et al. [90].After obtaining component albedos, we developed an algorithm for estimating component net radiant flux (Equations ( 39)-( 42)) based on Equation (38).Compared with Zheng's universal empirical algorithm of component net radiant flux, our improved algorithm, which considers the actual surface parameters (e.g., component fraction, component albedo, component emissivity, and component temperature), has the potential to result in higher accuracy.

Sensitivity Analysis
As shown in the three flowcharts (Figure 3, Figure 10, and Figure 11), the most important feature of the modified multi-source parallel model is introducing surface component fractions into ET estimation.Surface component fractions are fundamental model parameters that play an important role in the calculations of temperature, sensible heat flux, and net radiant flux.It is therefore necessary to assess the model's sensitivity to the change of surface component fractions.Since only vegetation and soil components ( f v and f s ) were directly involved in ET estimation, the effect of variation in vegetation and soil fractions on ET estimation are discussed here.
We here assumed that all other surface component fractions are stable when the controlled variable changes.Although it is inconsistent with the sum-to-unity nature of the fully constrained LSMA (i.e., the sum of different component's fraction is [39]) used in this study, this assumption can better help to explore the effect of a single surface component fraction ( f v or f s ) on ET estimation.Based on the initial results of linear spectral mixture analysis, a change ranging from -0.3 to 0.3 with a f v (or f s ) increase of 0.1 has been set, and four types of sample areas with the initial f v (or f s ) value ranging from 0.3 to 0.4, from 0.4 to 0.5, from 0.5 to 0.6, and from 0.7 to 0.8 were chosen as the testing areas for sensitivity analysis.The average value of the corresponding ET of each testing area for each change of f v (or f s ) is shown in Figure 14.It is revealed that the variation of vegetation fraction or soil fraction shows a positive correlation with the corresponding ET (i.e., ET resulted from vegetation/soil fraction changes in the sensitivity analysis), and introduced obvious error into the ET estimation result.The difference between the initial ET and the corresponding ET for every increase of 0.1 in f v (or f s ) has been calculated and shown in Table 8.
Remote Sens. 2017, 9, x FOR PEER REVIEW 21 of 32 We here assumed that all other surface component fractions are stable when the controlled variable changes.Although it is inconsistent with the sum-to-unity nature of the fully constrained LSMA (i.e., the sum of different component's fraction is [39]) used in this study, this assumption can better help to explore the effect of a single surface component fraction ( or ) on ET estimation.Based on the initial results of linear spectral mixture analysis, a change ranging from -0.3 to 0.3 with a (or ) increase of 0.1 has been set, and four types of sample areas with the initial (or ) value ranging from 0.3 to 0.4, from 0.4 to 0.5, from 0.5 to 0.6, and from 0.7 to 0.8 were chosen as the testing areas for sensitivity analysis.The average value of the corresponding ET of each testing area for each change of (or ) is shown in Figure 14.It is revealed that the variation of vegetation fraction or soil fraction shows a positive correlation with the corresponding ET (i.e., ET resulted from vegetation/soil fraction changes in the sensitivity analysis), and introduced obvious error into the ET estimation result.The difference between the initial ET and the corresponding ET for every increase of 0.1 in (or ) has been calculated and shown in Table 8.Table 8.ET error between the initial ET and the corresponding ET with the variation of (or ).It is observed from Table 8 that the ET error increases with the change of vegetation/soil fraction (i.e., ∆ and ∆ ).The maximal ET error has reached 1.717 mm•day −1 and 0.965 mm•day −1 , respectively, when or has increased by 0.3.It also indicates that: (1) with ∆ ranging from −0.3 to 0.3, the average ET error rate of all vegetation fraction testing area (testing area 1 to 4) for each ∆ are 45.63%, 23.86%, 6.68%, 23.33%, 43.36% and 67.95%; and (2) with ∆ variated from −0.3 to 0.3, the average ET error rate of all soil fraction testing area (testing area 5 to 8) for each ∆ are 30.62%,17.71%, 12.64%, 14.65%, 22.60% and 38.34%.Sensitivity analysis has revealed that vegetation fraction and soil fraction are two factors sensitive to ET estimation, and that the sensitivity of ET to the vegetation fraction is higher than to the soil fraction.It is observed from Table 8 that the ET error increases with the change of vegetation/soil fraction (i.e., ∆ f v and ∆ f s ).The maximal ET error has reached 1.717 mm•day −1 and 0.965 mm•day −1 , respectively, when f v or f s has increased by 0.3.It also indicates that: (1) with ∆ f v ranging from −0.3 to 0.3, the average ET error rate of all vegetation fraction testing area (testing area 1 to 4) for each ∆ f v are 45.63%, 23.86%, 6.68%, 23.33%, 43.36% and 67.95%; and (2) with ∆ f s variated from −0.3 to 0.3, the average ET error rate of all soil fraction testing area (testing area 5 to 8) for each ∆ f s are 30.62%,17.71%, 12.64%, 14.65%, 22.60% and 38.34%.Sensitivity analysis has revealed that vegetation fraction and soil fraction are two factors sensitive to ET estimation, and that the sensitivity of ET to the vegetation fraction is higher than to the soil fraction.

Comparison with Other Models
On an urban scale, remote sensing-based algorithms are considered most appropriate.Although our model is designed for an urban environment, it is interesting, in terms of estimation accuracy, to compare it with other remote sensing-based models that are often applied to natural or agricultural surfaces, e.g., SEBS, SEBI, SEBAL, S-SEBI, METRIC, and TSM [91].Su et al. [92] used SEBS modeled surface heat fluxes of regional scale by Landsat data and obtained a result with an error of approximately 5%.Roerink et al. [19] compared measured evaporative fraction values with the result from S-SEBI, and the maximum relative difference between them was 8%.SEBAL is a widely tested model whose accuracy has been reported as 85% for 1 day compared to in situ data, and can increase to 95% on a seasonal basis [93].Bastiaanssen et al. [17] indicated the overall error of the result from remote sensing-based SEBAL does not exceed 15%.Based on SEBAL, Allen et al. [94] proposed the METRIC model, whose error, reported by Gowda et al. [95], was 12.7±8.1% (mean bias error ± RMSE) on DOY 178 and −4.7 ± 9.4% on DOY 210.Zheng's model and our modified model were developed from TSM and consider impervious surfaces for urban areas in order to remove the interference components from mixed pixels.Long and Singh [29] developed a TSM model (TTME) and concluded that TTME is capable of reproducing latent heat flux, the mean absolute percentage difference is about 10%.Compared with the abovementioned models, our model's accuracy is acceptable with mean relative errors of 6.08%, 6.12%, and 5.90% for vegetation-, soil-, and impervious surface-dominated areas respectively.

Comparison with Zheng's Model
In order to demonstrate the improvement of our model through comparative analysis, we also inferred the daily evapotranspiration of the study area on 13 October 2013 using Zheng's model.The ET inferred from Zheng's model and the ET difference map (i.e., the ET derived by our model minus the ET derived by Zheng' model) are shown in Figure 15.The difference map (Figure 15b) shows the difference value are mostly positive, and the ET difference values in vegetated and bare soil areas, indicating that the ET distributions inferred from our and Zheng's models are almost consistent with each other in city parks and undeveloped zones in suburbs.However, in the downtown or in intensively built areas, the differences between the two ET results are obvious.This suggests that there is a more significant difference in the ET estimation of the two models in impervious surface-dominated areas.

Comparison with Other Models
On an urban scale, remote sensing-based algorithms are considered most appropriate.Although our model is designed for an urban environment, it is interesting, in terms of estimation accuracy, to compare it with other remote sensing-based models that are often applied to natural or agricultural surfaces, e.g., SEBS, SEBI, SEBAL, S-SEBI, METRIC, and TSM [91].Su et al. [92] used SEBS modeled surface heat fluxes of regional scale by Landsat data and obtained a result with an error of approximately 5%.Roerink et al. [19] compared measured evaporative fraction values with the result from S-SEBI, and the maximum relative difference between them was 8%.SEBAL is a widely tested model whose accuracy has been reported as 85% for 1 day compared to in situ data, and can increase to 95% on a seasonal basis [93].Bastiaanssen et al. [17] indicated the overall error of the result from remote sensing-based SEBAL does not exceed 15%.Based on SEBAL, Allen et al. [94] proposed the METRIC model, whose error, reported by Gowda et al. [95], was 12.7±8.1% (mean bias error ± RMSE) on DOY 178 and −4.7 ± 9.4% on DOY 210.Zheng's model and our modified model were developed from TSM and consider impervious surfaces for urban areas in order to remove the interference components from mixed pixels.Long and Singh [29] developed a TSM model (TTME) and concluded that TTME is capable of reproducing latent heat flux, the mean absolute percentage difference is about 10%.Compared with the abovementioned models, our model's accuracy is acceptable with mean relative errors of 6.08%, 6.12%, and 5.90% for vegetation-, soil-, and impervious surfacedominated areas respectively.

Comparison with Zheng's Model
In order to demonstrate the improvement of our model through comparative analysis, we also inferred the daily evapotranspiration of the study area on 13 October 2013 using Zheng's model.The ET inferred from Zheng's model and the ET difference map (i.e., the ET derived by our model minus the ET derived by Zheng' model) are shown in Figure 15.The difference map (Figure 15b) shows the difference value are mostly positive, and the ET difference values in vegetated and bare soil areas, indicating that the ET distributions inferred from our and Zheng's models are almost consistent with each other in city parks and undeveloped zones in suburbs.However, in the downtown or in intensively built areas, the differences between the two ET results are obvious.This suggests that there is a more significant difference in the ET estimation of the two models in impervious surfacedominated areas.The ET result inferred from Zheng's model is also validated by the China-ASEAN ET product, as shown in Figure 16 and Table 9 (based on Figure 16b).Compared with Figure 13a, the distribution of the observations is more discrete, and the R 2 value for the linear relationship between the ET estimated from Zheng's model and the China-ASEAN ET product is lower in Figure 16a.
The ET result inferred from Zheng's model is also validated by the China-ASEAN ET product, as shown in Figure 16 and Table 9 (based on Figure 16b).Compared with Figure 13a, the distribution of the observations is more discrete, and the R 2 value for the linear relationship between the ET estimated from Zheng's model and the China-ASEAN ET product is lower in Figure 16a.Table 9 shows that the accuracy of Zheng's model for impervious surface-dominated areas is obviously lower than that for vegetation-dominated and soil-dominated areas.Tables 7 and 9 suggest that the two models both performed best for vegetation-dominated areas.It is also noted that the correlation, fitting goodness, and RMSE between our result and the validation data were all improved.
The Taylor diagram [96] (Figure 17) reveals the accuracy improvement of the modified model.In the Taylor diagram, a single point indicates the correlation degree ( = ) and the ratio of the standard deviations ( ) between the modelled result ( ) and the reference (true) data ( ) ( = / ).An ideal model has a standard deviation ratio of 1.0 and a correlation degree of 1.0, i.e., the reference point (REF) on the x-axis [80].
Taylor skill ( ) [96][97][98] is defined by Taylor as a single value summary of a Taylor diagram in which unity indicates perfect agreement with inversion values.Traditionally, skill scores vary from 0 (least skillful) to 1 (most skillful) and each point for any arbitrary data group can be scored as follows [96] The calculated Taylor skill values are given in Table 10.The results suggest that our model outperforms Zheng's model.Table 9 shows that the accuracy of Zheng's model for impervious surface-dominated areas is obviously lower than that for vegetation-dominated and soil-dominated areas.Tables 7 and 9 suggest that the two models both performed best for vegetation-dominated areas.It is also noted that the correlation, fitting goodness, and RMSE between our result and the validation data were all improved.
The Taylor diagram [96] (Figure 17) reveals the accuracy improvement of the modified model.In the Taylor diagram, a single point indicates the correlation degree (R d = r) and the ratio of the standard deviations (σ) between the modelled result (σ i ) and the reference (true) data (σ 0 ) (σ norm =σ i / σ 0 ).An ideal model has a standard deviation ratio of 1.0 and a correlation degree of 1.0, i.e., the reference point (REF) on the x-axis [80].
Taylor skill (S) [96][97][98] is defined by Taylor as a single value summary of a Taylor diagram in which unity indicates perfect agreement with inversion values.Traditionally, skill scores vary from 0 (least skillful) to 1 (most skillful) and each point for any arbitrary data group can be scored as follows [96]: The calculated Taylor skill values are given in Table 10.The results suggest that our model outperforms Zheng's model.

Applicability of Our Modified Model
To further explore its applicability and robustness, our modified model is also tested to the urban built-up area of Wuhan, the capital city of Hubei Province, central China.The ASTER L1B level data and the HJ-1B CCD1 data of Wuhan acquired on 8 August 2013, were used as testing data.Detailed information of the satellite images and meteorological data are shown in Appendix Table A2.The result of daily evapotranspiration of Wuhan's urban built-up area on 8 August 2013 is obtained through our modified model and given in Figure 18.
The Wuhan ET estimation is also validated by the China-ASEAN ET product, as illustrated in Figure 19 and Table 11.

Applicability of Our Modified Model
To further explore its applicability and robustness, our modified model is also tested to the urban built-up area of Wuhan, the capital city of Hubei Province, central China.The ASTER L1B level data and the HJ-1B CCD1 data of Wuhan acquired on 8 August 2013, were used as testing data.Detailed information of the satellite images and meteorological data are shown in Appendix A Table A2.The result of daily evapotranspiration of Wuhan's urban built-up area on 8 August 2013 is obtained through our modified model and given in Figure 18.For both the entire testing area and component-dominated areas, the ET estimated by our modified model shows a good agreement with the China-ASEAN ET product.In particular, the highest accuracy is also obtained in the vegetation-dominated area of Wuhan, which is similar to the case study of Hefei.The test on Wuhan has confirmed that our modified model is applicable for different urban areas.

Conclusions
This study has presented a modified multi-source parallel model to improve the accuracy of estimating urban land surface evapotranspiration.Key modifications include: The Wuhan ET estimation is also validated by the China-ASEAN ET product, as illustrated in Figure 19 and Table 11.For both the entire testing area and component-dominated areas, the ET estimated by our modified model shows a good agreement with the China-ASEAN ET product.In particular, the highest accuracy is also obtained in the vegetation-dominated area of Wuhan, which is similar to the case study of Hefei.The test on Wuhan has confirmed that our modified model is applicable for different urban areas.

Conclusions
This study has presented a modified multi-source parallel model to improve the accuracy of estimating urban land surface evapotranspiration.Key modifications include: For both the entire testing area and component-dominated areas, the ET estimated by our modified model shows a good agreement with the China-ASEAN ET product.In particular, the highest accuracy is also obtained in the vegetation-dominated area of Wuhan, which is similar to the case study of Hefei.The test on Wuhan has confirmed that our modified model is applicable for different urban areas.

Conclusions
This study has presented a modified multi-source parallel model to improve the accuracy of estimating urban land surface evapotranspiration.Key modifications include:

•
Instead of a single endmember, impervious surfaces are characterized by two different ones (highand low-albedo) in linear spectral mixture analysis.

•
Rather than constant empirical values, the roughness length for each land surface component (vegetation, soil, high-and low-albedo impervious surfaces) is calculated specifically to better approximate the real conditions of those land surfaces.

•
Our model includes a new algorithm for estimating component net radiant flux by taking the fraction and characteristics of each land surface component into account.
By comparing the evapotranspiration of an urban built-up area, inferred from our modified model and from Zheng's model with the China-ASEAN ET product, it is concluded that both models can produce good results for urban land surface, but our modified model outperforms Zheng's-which can be explained by the abovementioned modifications.It is also noted that both models can result in more accurate evapotranspiration for vegetated areas.
As fundamental parameters in the ET estimation model, surface component fractions are used multiple times in the calculations of other parameters.Therefore, surface component fractions, particularly vegetation and soil fractions, are sensitive factors for the final ET result.The accuracy of linear spectral unmixing has an important impact on the accuracy of ET estimation.
Our modified model is practical and easy to apply, as it allows for estimation of daily evapotranspiration for urban areas using only freely available optical and thermal image data.This has been demonstrated by an additional test of our model in Wuhan's built-up area.Due to the absence of field-based latent heat radiation data from the study area, however, this study had to use the China-ASEAN ET product as validation data.Despite having the same acquisition time as the image data used in the study, this product has a low spatial resolution, thus making it less accurate than field-based measurements.Future work should investigate whether the algorithms of sub-pixel ET estimation can improve other remote sensing-based ET models.

Figure 1 .
Figure 1.Study area: (a) Anhui province in China; (b) study area in Anhui's capital city of Hefei; (c) study area highlighted in red in the true color composite of HJ-1A satellite imagery (13 October 2013).

Figure 1 .
Figure 1.Study area: (a) Anhui province in China; (b) study area in Anhui's capital city of Hefei; (c) study area highlighted in red in the true color composite of HJ-1A satellite imagery (13 October 2013).

Figure 1 .
Figure 1.Scatter plot of the China-ASEAN ET product against the MOD16 ET product on the 281st day of 2013 for a rectangular site near the study area.

Figure 2 .
Figure 2. Scatter plot of the China-ASEAN ET product against the MOD16 ET product on the 281st day of 2013 for a rectangular site near the study area.
Remote Sens. 2017, 9, x FOR PEER REVIEW 6 of 32In order to better illustrate the methodology of this study, three flowcharts detailing the input data, parameters, and equations are given in different sections of the paper.The flowchart for the first step of our modified model (Sections 3.1 and 3.2) is shown in Figure3.This figure demonstrates estimation of component temperatures from remotely sensed imagery.

Figure 3 .
Figure 3. Flowchart of the first step (Step 1) of our model, as detailed in Sections 3.1 and 3.2.The gray rectangles represent the input data, and the other color rectangles indicate the parameters that were used more than once in the next steps.Equations used in each of the processes are highlighted in blue.Full forms of the variables in the flowcharts are listed in Table A1 in the Appendix.

Figure 3 .
Figure 3. Flowchart of the first step (Step 1) of our model, as detailed in Sections 3.1 and 3.2.The gray rectangles represent the input data, and the other color rectangles indicate the parameters that were used more than once in the next steps.Equations used in each of the processes are highlighted in blue.Full forms of the variables in the flowcharts are listed in Table A1 in the Appendix A.

Figure 4 .
Figure 4. Normalized spectral reflectance for HJ-1 VNIR bands of the four selected endmembers.

Figure 5 .
Figure 5. RMSE image resulting from the fully constrained LSMA.

Figure 4 .
Figure 4. Normalized spectral reflectance for HJ-1 VNIR bands of the four selected endmembers.

Figure 5 .
Figure 5. RMSE image resulting from the fully constrained LSMA.

Figure 6 .
Figure 6.(a) A 2 m resolution satellite image of the study area on 2 October 2013; (b) one of 100 randomly distributed validation samples (show in red rectangle) on the Google Earth image; (c) interpretation of one validation samples.

Figure 7 .
Figure 7. Scatter plots of LSMA modeled surface component fractions against those extracted from Google Earth image: (a) vegetation fraction; (b) soil fraction.

Figure 7 .
Figure 7. Scatter plots of LSMA modeled surface component fractions against those extracted from Google Earth image: (a) vegetation fraction; (b) soil fraction.

Figure 8 .
Figure 8.The relationship between the Planck blackbody radiation and temperature for ASTER thermal infrared bands.

Figure 8 .
Figure 8.The relationship between the Planck blackbody radiation and temperature for ASTER thermal infrared bands.

Figure 9 .
Figure 9. Emissivity spectral curve for each component.These curves were the average of typical urban ground objects listed in the ASTER Spectral Library (Versions 2.0).

Figure 9 .
Figure 9. Emissivity spectral curve for each component.These curves were the average of typical urban ground objects listed in the ASTER Spectral Library (Versions 2.0).

Figure 10 .
Figure 10.Flowchart of the second step (Step 2) of our model.

Figure 10 .
Figure 10.Flowchart of the second step (Step 2) of our model.

) 3 . 4 .
Retrieval of Land Surface Component Net Radiant Flux The flowchart for the third step of our modified model (Sections 3.4-3.6) is shown in Figure 11.This figure demonstrates the estimation of component net radiant flux, internal heat flux and the final ET result based on the intermediate parameters from Step 1 and Step 2.

_
are component surface net radiant flux of vegetation, soil, high-and low-albedo impervious surfaces;, , _ , and _ are the surface albedo of these components.

Figure 12 .
Figure 12.The spatial variability of daily evapotranspiration of Hefei's urban built-up area on 13 October 2013, inferred from our modified model.

Figure 12 .
Figure 12.The spatial variability of daily evapotranspiration of Hefei's urban built-up area on 13 October 2013, inferred from our modified model.

Figure 13 .
Figure 13.(a) Scatter plot of the modified model derived ET against the China-ASEAN ET product, based on the entire 552 sampled points; (b) scatter plot of the modified model derived ET against the China-ASEAN ET product, for each land cover (124 sampled points for vegetation, 130 for soil, and 102 for impervious surfaces).

Figure 13 .
Figure 13.(a) Scatter plot of the modified model derived ET against the China-ASEAN ET product, based on the entire 552 sampled points; (b) scatter plot of the modified model derived ET against the China-ASEAN ET product, for each land cover (124 sampled points for vegetation, 130 for soil, and 102 for impervious surfaces).

Figure 14 .
Figure 14.Sensitivity of ET with the variation of vegetation fraction (a) and soil fraction (b).

Figure 14 .
Figure 14.Sensitivity of ET with the variation of vegetation fraction (a) and soil fraction (b).

Figure 15 .
Figure 15.(a) The spatial variability of daily evapotranspiration of Hefei's urban built-up area on 13 October 2013, inferred from Zheng's model; (b) the difference between the results of our model and Zheng's model.

Figure 15 .
Figure 15.(a) The spatial variability of daily evapotranspiration of Hefei's urban built-up area on 13 October 2013, inferred from Zheng's model; (b) the difference between the results of our model and Zheng's model.

Figure 16 .
Figure 16.(a) Scatter plot of Zheng's model derived ET against the China-ASEAN ET product, based on the entire 552 sampled points; (b) scatter plot of Zheng's model derived ET against the China-ASEAN ET product, for each land cover (136 sampled points for vegetation, 117 for soil, and 93 for impervious surfaces).

Figure 16 .
Figure 16.(a) Scatter plot of Zheng's model derived ET against the China-ASEAN ET product, based on the entire 552 sampled points; (b) scatter plot of Zheng's model derived ET against the China-ASEAN ET product, for each land cover (136 sampled points for vegetation, 117 for soil, and 93 for impervious surfaces).

Figure 17 .
Figure 17.Performance of the inversion values inferred from the modified model and Zheng's model (the statistics in the Taylor diagram); an ideal model would have a standard deviation ratio ( ) of 1.0 and a correlation coefficient of 1.0 (REF is the reference point).

Figure 17 .
Figure 17.Performance of the inversion values inferred from the modified model and Zheng's model (the statistics in the Taylor diagram); an ideal model would have a standard deviation ratio (σ norm ) of 1.0 and a correlation coefficient of 1.0 (REF is the reference point).

Figure 18 .
Figure 18.(a) The testing area of Wuhan in Hubei province; (b) the spatial variability of daily evapotranspiration of Wuhan's urban built-up area on 8 August 2013, inferred from our modified model.

Figure 19 .
Figure 19.(a) Scatter plot of the modified model derived ET against the China-ASEAN ET product in Wuhan testing area, based on the entire 1446 sampled points; (b) scatter plot of the modified model derived ET against the China-ASEAN ET product, for each land cover (229 sampled points for vegetation, 149 for soil and 160 for impervious surfaces).

Figure 18 .
Figure 18.(a) The testing area of Wuhan in Hubei province; (b) the spatial variability of daily evapotranspiration of Wuhan's urban built-up area on 8 August 2013, inferred from our modified model.

Figure 18 .
Figure 18.(a) The testing area of Wuhan in Hubei province; (b) the spatial variability of daily evapotranspiration of Wuhan's urban built-up area on 8 August 2013, inferred from our modified model.

Figure 19 .
Figure 19.(a) Scatter plot of the modified model derived ET against the China-ASEAN ET product in Wuhan testing area, based on the entire 1446 sampled points; (b) scatter plot of the modified model derived ET against the China-ASEAN ET product, for each land cover (229 sampled points for vegetation, 149 for soil and 160 for impervious surfaces).

Figure 19 .
Figure 19.(a) Scatter plot of the modified model derived ET against the China-ASEAN ET product in Wuhan testing area, based on the entire 1446 sampled points; (b) scatter plot of the modified model derived ET against the China-ASEAN ET product, for each land cover (229 sampled points for vegetation, 149 for soil and 160 for impervious surfaces).

Table 2 .
The linear fitting coefficients of Planck function for ASTER thermal infrared bands 1 .

Table 2 .
The linear fitting coefficients of Planck function for ASTER thermal infrared bands 1 .

Table 3 .
Emissivity for each component at different ASTER thermal infrared bands.

Table 3 .
Emissivity for each component at different ASTER thermal infrared bands.

Table 5 .
Meteorological parameters required for calculating average atmospheric water vapor [58].
•K −1 ); T air is the air temperature (K), which can be obtained from a daily dataset of China's surface climate[58]; and r ahk is the aerodynamic resistance to heat transfer of component k (s•m −1 ).

Table 6 .
Average albedo of the four components for the study area.

Table 7 .
Linear relationship between our modeled ET and the China-ASEAN ET product for each land cover type (vegetation, soil, and impervious surfaces).

Table 7 .
Linear relationship between our modeled ET and the China-ASEAN ET product for each land cover type (vegetation, soil, and impervious surfaces).

Table 8 .
ET error between the initial ET and the corresponding ET with the variation of f v (or f s ).

Table 9 .
Linear relationship between the ET estimated from Zheng's model and the China-ASEAN ET product for each land cover type (vegetation, soil, and impervious surfaces).

Table 9 .
Linear relationship between the ET estimated from Zheng's model and the China-ASEAN ET product for each land cover type (vegetation, soil, and impervious surfaces).

Table 10 .
Indicators of the Taylor diagram.

Table 11 .
Linear relationship between our modeled ET and the China -ASEAN ET product for each land cover type (vegetation, soil, and impervious surfaces) in Wuhan testing area.

Table 10 .
Indicators of the Taylor diagram.

Table 11 .
Linear relationship between our modeled ET and the China -ASEAN ET product for each land cover type (vegetation, soil, and impervious surfaces) in Wuhan testing area.

Table A1 .
Cont.Reflectance for band b in a pixel R k,b Reflectance of endmember k for band b in a pixel L λ Thermal radiation of land surface at wavelength λ ε λ Land surface emissivity at wavelength λ B λ (T sur ) Planck blackbody radiation when the land surface temperature is T sur B λ,sens At-sensor thermal radiation at wavelength λ τ λ Atmospheric transmittance at wavelength λ T air Near-surface average atmospheric temperature ε λk Emissivity of k component at wavelength λ B λ (T k ) Planck blackbody radiation of k component at wavelength λ when the k component temperature is T k f v Temperature scale with T * = T sur − T air u * *

Table A1 .
Cont.Broad band albedo of HJ-1A CCD bands G imp_h Internal heat flux of high-albedo impervious surface G imp_l Internal heat flux of low-albedo impervious surface H v Sensible heat flux of vegetation H s Sensible heat flux of soil H imp_h Sensible heat flux of high-albedo impervious surface H imp_h Sensible heat flux of low-albedo impervious surface L

Table A2 .
Parameters of satellite image data and meteorological data for Wuhan testing area.