Estimation of land surface heat fluxes based on visible infrared imaging radiometer suite data: case study in northern China

Abstract. Evapotranspiration (ET) plays an important role in surface–atmosphere interactions and can be monitored using remote sensing data. The visible infrared imaging radiometer suite (VIIRS) sensor is a generation of optical satellite sensors that provide daily global coverage at 375- to 750-m spatial resolutions with 22 spectral channels (0.412 to 12.05  μm) and capable of monitoring ET from regional to global scales. However, few studies have focused on methods of acquiring ET from VIIRS images. The objective of this study is to introduce an algorithm that uses the VIIRS data and meteorological variables to estimate the energy budgets of land surfaces, including the net radiation, soil heat flux, sensible heat flux, and latent heat fluxes. A single-source model that based on surface energy balance equation is used to obtain surface heat fluxes within the Zhangye oasis in China. The results were validated using observations collected during the HiWATER (Heihe Watershed Allied Telemetry Experimental Research) project. To facilitate comparison, we also use moderate resolution imaging spectrometer (MODIS) data to retrieve the regional surface heat fluxes. The validation results show that it is feasible to estimate the turbulent heat flux based on the VIIRS sensor and that these data have certain advantages (i.e., the mean bias error of sensible heat flux is 15.23  W m−2) compared with MODIS data (i.e., the mean bias error of sensible heat flux is −29.36  W m−2). Error analysis indicates that, in our model, the accuracies of the estimated sensible heat fluxes rely on the errors in the retrieved surface temperatures and the canopy heights.


Introduction
Evapotranspiration (ET) is an important parameter in the interactions among soil, vegetation, and the atmosphere, and it is central part of surface energy and water budgets. 1,2 Approximately 46% to 53% of the incoming solar radiation absorbed by the surface is globally consumed by the ET process as latent heat. 3,4 ET can provide useful information about the heat and vapor exchange between the land surface and the atmosphere. A thorough understanding of ET can improve water resources management at different scales; in particular, detailed ET estimations in agricultural areas can improve the detection of water stresses and aid in irrigation scheduling. 5,6 2 Study Area and Datasets

Study Area
Our study was conducted in the middle stream of the Heihe River Basin (HRB) at the Zhangye oasis, which is located near the city of Zhangye in the arid region of Gansu Province in northwestern China (100.10°E to 100.66°E, 38.68°N to 39.15°N). The Zhangye oasis has a semiarid climate, with an annual precipitation of <200 mm and an annual potential ET of ∼1200 to 1800 mm; the annual average temperature is ∼7°C. 27,28 Areas of the Gobi desert and alpine vegetation are located near the study area in the Qilian Mountains (see Fig. 1). Around the oasis, the land surface is mainly covered by maize, which has a growing season that lasts from May to September. The artificial oasis is highly heterogeneous, which impacts its thermodynamic and hydraulic features. Consequently, its water use efficiency and ET are variable.
The HRB has long served as a test area for integrated watershed studies and land surface and hydrological experiments. 28 One major objective of the Heihe Watershed Allied Telemetry Experimental Research (HiWATER), namely the Multi-Scale Observation Experiment on ET (HiWATER-MUSOEXE), is to capture strong land surface heterogeneities and the associated uncertainties within a watershed. [29][30][31]

Datasets
The data in this study were mainly derived from VIIRS and MODIS. We combined these two sensors datasets with ancillary data and in situ HiWATER-MUSOEXE data to estimate and validate the VIIRS and MODIS heat fluxes.

Remote sensing data
Two kinds of remote sensing data were used in this study, specifically MODIS and VIIRS data. The MODIS sensor launched on December 18, 1999 aboard the Terra satellite has been in operation for nearly 18 years and provides global data that are used for monitoring long-term regional land surface heat and water vapor fluxes as well as other Earth-related changes. 32 The MODIS data have high temporal (i.e., at least once every day) and spectral resolutions (i.e., 36 bands in total). The NPP/VIIRS was launched in late 2011 and has imaging capabilities in 22 bands with a spectral coverage from 0.412 to 12 μm. Among the 22 bands, there are five high-resolution imaging channels (I-bands), 16 moderate resolution channels (M-bands), and 1 day/night band (DNB). The I-bands have spatial resolutions of 375 m at nadir, and M-bands (and DNB) have a resolution of 750 m at nadir. In the cross-track direction, the scan swath width is ∼3040 km. The data for one granule covers ∼85 s of data and comprises 48 scans. 33 The images produced by VIIRS are stable with balanced band performances.
In this paper, we mainly use two kinds of VIIRS operational products [i.e., VIIRS/NPP surface reflectance daily L2G global 1 km and 500 m SIN grid (VNP09GA) and VIIRS/NPP daily gridded land surface temperature 1-km SIN grid day (NPP_DLSTD_L3D)], which can be freely downloaded. 34 The parameters from the MODIS sensor used for estimating heat fluxes, such as normalized difference vegetation index (NDVI), albedo, LSTs, and emissivity, are freely available from the NASA website. We used the MODIS and VIIRS image data that cover the HRB region from 2012. Many algorithms for estimating ET using surface parameters require clear-sky conditions; therefore, we combined the available information on data quality with visual interpretations to select satellite images without clouds. We obtained data for 9 days during the period of ground observations: June 19 and 30, July 8 and 27, August 3, 22, and 29, and September 2 and 14. These two kinds of satellite products used in this study are shown in Table 1.
All of the above data are required for the corresponding preprocessed. The MODIS products that cover the study area were downloaded and corrected for geometric and projection effects. VNP09GA is a level-2G surface reflectance product obtained by adjusting the top-of-atmosphere reflectance to compensate for atmospheric effects. The effects of molecular gases, including ozone and water vapor, and the effects of atmospheric aerosols are corrected. Previous studies showed that VIIRS LST products have reasonable accuracies. 25 The VIIRS data covering the HRB were preprocessed using geometric correction, radiometric calibration, and atmospheric correction. 14 To estimate the regional ET, various surface variables, including the NDVI, the surface albedo, 26 the emissivity, 35 the leaf area index (LAI), and the fractional vegetation coverage (FVC), were retrieved by applying the lumped method to the VIIRS data. 14 These two kinds of satellite data were organized and paired with each other on the Universal Transverse Mercator (UTM) projection (i.e., WGS 84 UTM ZONE 47N).

Data from the HiWATER experiment
The in situ HRB observational data were provided by HiWATER. From June to September 2012, HiWATER designed nested observation matrices of 30 km × 30 km and 5.5 km × 5.5 km within the middle stream oasis.
For the larger observation matrix, four EC systems and one superstation (EC15) were installed in the oasis-desert ecosystem. Each station was supplemented with an automatic meteorological station (AMS) to record the meteorological and soil variables and to monitor the spatiotemporal variations of ET and its associated factors, 28 including air pressure, precipitation, wind speed/direction, four-component radiation, and soil moisture. Each AMS has a different observation date, and its specific observation time and sensor characteristics can be found in Ref. 30. Within the oasis, 17 AMSs and EC towers were used, including one in a vegetable field (EC01), one in a residential area (EC04), one in an orchard (EC17), and 14 in maize fields. Because the HHZ station (shown in Fig. 1) and EC16 station lacked net radiations and soil heat flux observation data, they were excluded from this study. Additionally, four groups of optical LASs (eight sets, with two sets per group) were installed in the 3 × 3 and 2 × 1 MODIS pixels within the kernel experimental area (three groups in three 3 × 1 MODIS pixels, named LAS1 to LAS3, from west to east, respectively, and one group in a 2 × 1 MODIS pixels, LAS4). Considering the observational scale of the LAS systems and the resolution of satellite data used in this paper, we mainly use the EC data to validate the surface variables (i.e., land surface temperature, net radiation, and soil heat flux) and the LAS data to validate the sensible heat flux. The details of these LAS systems are shown in Table 2.
An LAS is a device that derives the turbulent intensity through measuring the refractive index of air, C 2 n ðm −2∕3 Þ. The sensible heat flux (LAS-H) was obtained using Monin-Obukhov similarity theory (MOST) via an iterative method that combines various meteorological data (e.g., wind speed, air temperature, and air pressure). 7,30 Finally, the LAS-H data were averaged across 30-min intervals. The data were carefully screened to ensure the LAS data quality: (i) the data were rejected when C 2 n exceeded the saturated criterion, which was determined according to Ref. 36, (ii) the data were rejected when the demodulated signal was less than the threshold, 30 (iii) the data were rejected when collected on nights with weak turbulence, and (iv) the data were rejected when precipitation occurred. Additionally, to avoid advection conditions, we selected data at 30-min intervals during periods when the H values of the EC and LAS measurements were >10 Wm −2 . The soil heat flux was measured using three soil heat plates at a depth of 6 cm at each site and the surface soil heat flux was calculated based on the soil temperature and moisture above the plates. 14 The surface meteorological variables, such as wind direction, wind speed, air pressure, and relative humidity, were used to interpolate images using inverse distance weighting. The data obtained during the HiWATER experiment can be freely obtained from the Heihe Plan Data Management Center. 37

Methodology
We have developed a single-source ET estimation algorithm that combines sensor-derived surface variables with meteorological data over the northwestern China HRB area. For comparison, we use the MODIS-derived surface variables to calculate the regional ET values and the two ETs are compared with ground-based observations; the errors in the model are then analyzed and discussed.

Pixel Evapotranspiration Algorithm
The SEB describes the transfer of energy between the land surface and the atmosphere. The energy budget is commonly expressed as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 2 2 4 , and LE is the latent heat flux (W m −2 ). Sensible heat flux is the flux caused by the temperature gradient between the land and the air, whereas LE is the change in energy due to ET or condensation. Atmospheric turbulence at the boundary layer causes the heat transfer at the land surface and atmospheric boundary. LE can be derived from the residual of the energy balance equation when R n , G, and H are known. Figure 2 shows the flowchart for ET retrieval using VIIRS data. The models that permit the estimations of R n , G, and H using remote sensing data are introduced. R n , also known as the radiation balance or radiation budget, is the main energy source of land surface heat and mass transfers and exchanges, such as evaporation, and is calculated as the difference between the incoming and outgoing radiation where R s ↓ is the downward shortwave radiation (W m −2 ), R s ↑ is the shortwave radiation reflected by the surface (W m −2 ), R L ↓ is the downward atmospheric longwave radiation (W m −2 ), and R L ↑ is the upward longwave radiation emitted by the surface (W m −2 ). α is the surface albedo, ε s is the emissivity of the land surface, R L ↓¼ ε a σT 4 a , ε a is the effective emissivity of the atmosphere, σ ¼ 5.67 × 10 −8 W m −2 K −4 is the Stefan-Boltzmann constant, T a is the air temperature, and T s is the surface radiation temperature.
G is the heat exchange between the surface and the deeper soils and is related to the exchanges of water and heat within the soil. Generally, heat is transferred from deeper depths to the surface at night. Studies have found that soil heat flux is correlated with R n and can be derived using land surface parameters. 15 Because canopies exert significant influences on G, the fractional canopy coverage f c is used to determine the ratio of G to R n 14 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 3 5 3 where Γ c ¼ 0.05 and Γ s ¼ 0.315 represent the ratio of G and R n with full canopy coverage and bare soil.
H is the turbulent heat transfer between the surface and the atmosphere that is derived using the temperature gradient between them. According to gradient diffusion theory E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 2 7 4 where ρ is the air density (kg m −3 ), C p is the specific heat of air at a constant pressure, and T aero is the radiative surface temperature (K). T a is the air temperature at a reference height (K), and r ac is the aerodynamic resistance (s −1 m), which is determined using the zero-plane displacement height and the roughness length. r ac was calculated based on MOST using a stability correction function. [38][39][40] Because remote sensing-based methods cannot obtain T aero , the value of T aero is usually replaced by the radiative surface temperature (T s ). It should be noted that T s is not strictly equal to the T aero in Eq. (4). The difference between these terms for homogeneous and full-coverage vegetation is ∼1 to 2 K 41 and it can reach 10 K in sparsely vegetative areas. 42 The method used by Chen et al. 43 corrects for this discrepancy by adding an "excess resistance" r bh and thus, was used in this study. The modified H can be specified as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 1 1 1

Accuracy Evaluation Index
To control model inputs and analyze error sources, the land surface temperature, net radiation, and soil heat flux were evaluated using in situ data. The analyses were performed using the root mean square error (RMSE), the mean bias error (MBE), and the coefficient of determination (R 2 ) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 6 8 0 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 6 2 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 5 7 8 where m i represents the observed value, e i represents the estimated value, and n represents the number of observations. The ground-based land surface temperature (T s ) was calculated using the Stefan-Boltzmann law from the AMS measurements of the longwave radiation fluxes 25 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 4 8 3 where R L ↑and R L ↓are the in situ surface upwelling and atmospheric downwelling longwave radiations, respectively, and ε S is the surface broadband emissivity.

Evaluation of Surface Variables
The station validation results of two LSTs are shown in Tables 3 and 4. The R 2 , MBE, and RMSE values of the VIIRS LST data are 0.77, 0.55 K, and 2.40 K, respectively. As shown in Table 3, the VIIRS MBE varies from −6.67 to 1.95 K; the RMSE varies from 1.16 to 7.10 K. The VIIRS LSTs are warmer than the ground LSTs for most stations. From Table 4, the MBE values for There are two possible reasons for the large Terra MODIS LST errors over barren surface sites. First, the LST product algorithm does not cover the wide range of LSTs, i.e., the daytime LSTs are greater than the air temperature at the surface level by >16 K, which is beyond the upper limit. Second, the large errors are present in the surface emissivity values in MODIS bands 31 and 32 estimated from land cover types. 44 As shown in Tables 3 and 4, the EC4 accuracy of both LST products is low. The main causes of the large errors are as follows: (1) buildings and soil/vegetation are distinct materials, and the emissivity algorithm may not be suitable for buildings and (2) the EC4 foundation is nonuniform and is not suitable for validation. After removing the EC4 data, the R 2 , MBE, and RMSE values of the VIIRS LST data are 0.86, 0.75 K, and 1.97 K, respectively, and the R 2 , MBE, and RMSE values of the MODIS LST data are 0.82, −0.87 K, and 2.77 K, respectively. The SD errors for both were large due to large errors on particular days. For example, although it was briefly cloudy above station SD on July 27, this area was not identified as cloudy in the cloud detection process.
The validation of net radiation and soil heat flux based on scatterplots is shown in Fig. 3. As shown in this figure, the R 2 , MBE, and RMSE values of the R n calculated using the VIIRS data were 0.68, −8.44 W m −2 , and 40.80 W m −2 , respectively, whereas the errors in the R n estimates obtained using the MODIS data were 0.64, −33.87 W m −2 , and 55.25 W m −2 , respectively. According to the sensitivity analysis of Eq. (2), R s ↓ and R L ↓ are highly sensitive variables when calculating R n , whereas the albedo, LST are not as sensitive. However, the R s ↓ and R L ↓ values used in the ET model are derived from ground observations; thus, the R s ↓ and R L ↓ are relatively accurate. Therefore, LST and albedo are the main sources of error and are discussed in Sec. 4.3. A 1 K bias in LST would result in an R n error of −0.96%, but an albedo bias of 0.03 can lead to an R n error of ∼20 W m −2 when the incoming solar radiation is large. Studies have shown that land surface albedos retrieved using VIIRS data are better than those retrieved using MODIS; 26 thus, the R n values calculated using the VIIRS data show better results.
The R 2 , MBE, and RMSE values of G calculated using the VIIRS data were 0.58, −11.80 W m −2 , and 30.84 W m −2 , respectively, whereas the errors in the G values calculated using the MODIS data were 0.52, −20.47 W m −2 , and 40.98 W m −2 , respectively. Since G is calculated using R n , G and R n show the same sensitivity results. As shown in the validation of G in Fig. 3, the sensor-retrieved G is overestimated at many sites, which may be related to the soil temperature and moisture above the soil heat plates at each site. For example, at EC5, the soil temperature and moisture were the same at different depths after July 19, which resulted in the surface G being equal to the G at a depth of 6 cm. However, the G values below the surface are usually less than those at the soil surface; thus, the validation results of G at EC5 indicate that G was overestimated. However, G is still underestimated as a whole (MBE < 0) because R n is also underestimated. Because the barren areas in the southwest border of the Qilian Mountains, the groundwater levels are high due to snow melt and the downward movement of water, leading to a soil moisture content of ∼30% according to in situ measurements at a depth of 2 cm. 14 As shown in Figs. 4(e) and 4(f), some areas in the H distribution map have values less than zero due to inversion from the oasis effect or irrigation. The HiWATER soil moisture data show that irrigation occurred on June 19, 2012. Irrigation is the main source of water within the oasis and cools the land surface to temperatures below that of the air. To further analyze the heat fluxes retrieved from the two sensors, the H values were evaluated using in situ data. In recent years, many researchers have employed a footprint model to validate remote sensing-based ET estimates (specifically only H and LE). When combined with the flux contribution source areas, this model solves the problems associated with mismatches between the observed values and the remote sensing-based estimates and leads to better results. 45,46 In this study, we employ the path-weighting function of the LAS and the footprint model for point fluxes to calculate the LAS source area. 30,46 The station-validated results for the H values are shown in Fig. 5 Fig. 3 Comparison of the R n and G data derived from the two sensors with measured data.

Error and Sensitivity Analysis
Many sources of error are unavoidable, including the calibration coefficients of the onboard satellite sensors; the drifts in the sensors caused by the aging of the instruments; mismatches in the timescales of the observational data, which are averaged over time, and the remotely sensed turbulent heat fluxes, which are instantaneous; and the geometric corrections cause a half-pixel bias less than or equal to the deviation of the artificially subjective interpretation. The estimation of H is the most important and difficult step in the remote sensing-based assessment of the SEB; thus, the sensitivity of H was analyzed first. Land surface variables, including LST, canopy height, FVC, and LAI, and meteorological variables, including wind speed, air temperature, and relative humidity, are the major factors in the sensitivity analysis of H. Figure 6 presents a sensitivity analysis of H. In this case, the average LST is 304 K, with a range from 297.5 to 311 K, and has a step size of 0.5 K; meanwhile, LAI ranges from 0.1 to 2.5 and has a step size of 0.1. The canopy height is 1 m, with a range from 0.1 to 1.9 m. The wind speed ranges from 0.50 to 4.95 m s −1 , the air temperature is T a ¼ 297.90 K, FVC ¼ 0.54, and the relative humidity is RH ¼ 40.29%. In addition, the land cover is maize and the reference H is 250.05 W m −2 .
The air pressure is stable over short periods and has little effect on the ET results. Although, H is sensitive to meteorological variables, such as wind speed and air temperature (see Fig. 6), the meteorological data used in the ET model are derived from ground observations; thus, the meteorological factors are relatively accurate. As shown in Fig. 6, LAI, LST, and the canopy height are all influential variables. H is sensitive to LAI when LAI is <1. The momentum roughness length increases as LAI increases and the turbulent exchange is enhanced. However, when LAI is >1, the plant canopy is regarded as a continuum that is not a sensitive variable. Because our study area is dominated by agriculture and the study period extended from June to September, the crops in the middle reaches of the HRB grew quickly; thus, LAI was usually >1. Hence, LST and canopy height are the main sources of error. 1. Errors in LST. As shown in Fig. 6, a 1 K bias in LST would result in 20.5% error in H when H is 250.05 W m −2 , the corresponding errors in R n and G are both only −0.96%. Therefore, in this paper, LST is mainly affected by the impact of H on the role of heat fluxes. However, as we mentioned in Sec. 3, H depends on the temperature gradient between the surface and the atmosphere. Hence, the LST sensitivity is unstable and depends on the strength of the local turbulence. The strength of the turbulence determines the mass and energy transports and the resistance to heat transfer, which influences the sensitivity of H to the LST. A weaker turbulence corresponds to a weaker LST sensitivity and vice versa. 2. Errors in the canopy height. The canopy height was obtained from a phenophase and classification map. Thus, the accuracy of the canopy height depends mainly on the plant growth state and classification accuracies. Even within the same region, the canopy height of a crop can differ due to differences in seeding times and soil attributes, such as soil moisture and fertilization.
As discussed above, the accuracy of the LST and canopy height parameters greatly influences the estimated turbulent heat flux. Since the VIIRS sensor is still running, the algorithm to calculate its LST product is still being explored. 25 Canopy heights are known a priori and are derived from phenophase classifications; these classifications influence the accuracies of the surface roughness calculations. Multisource remote sensing data, such as active microwave and LiDAR data, can be used to obtain canopy heights for future studies, 14 which would lead to more robust canopy height values. Since a footprint model was used in the validation, discrepancies between the in situ measurements and remote sensing pixels still exist in the footprints. 3 Additionally, to correct the discrepancies between the remotely sensed radiative surface temperatures and the aerodynamic temperatures at the source of heat transport, a brief and well-constrained parameterization scheme (under a uniform and flat plant surface) that included excess resistance was used to calculate the aerodynamic resistance of the heat transfer. 14 Since the surface cover of the HRB is very heterogeneous, 4,30 multiple parameterization methods should be compared to select an optimal method. 43 Although these preliminary results show that VIIRS data have certain advantages over MODIS data when calculating surface heat fluxes, the limited ground observation data in this paper use only 9 days in calculating and validating the results. Future work will involve testing this single-source ET estimation algorithm and the practicality of VIIRS data for calculating the heat fluxes in the range of actual environmental conditions.

Conclusion
We use a single-source model and VIIRS data to calculate the components of SEB. To determine the usefulness VIIRS sensors in estimating regional heat fluxes, we also used MODIS data to calculate the regional ET. After data preparation, SEB component inversion, verification (R n , G, and H), and error analysis, the following conclusions can be drawn: 1. We have developed a single-source ET estimation algorithm that combines sensorderived surface variables with meteorological data over the HRB area of northwestern China. This method can be used to estimate the turbulent heat flux based on VIIRS sensor data, which has certain advantages (i.e., the MBE of the sensible heat flux is 15.23 W m −2 ) compared with using MODIS data (i.e., the MBE of the sensible heat flux is −29.36 W m −2 ). 2. To control model inputs and analyze error sources, the model parameters (including the land surface temperature, net radiation, and soil heat flux) were evaluated using in situ data. The validation results show that the VIIRS LST product exhibits a more consistent agreement with in situ measurements on barren surfaces compared with the Terra MODIS LST product. 3. To estimate ET using a single-source model and remote sensing data, the sensible heat flux is the most important component, but the calculation of this component is also complex. In our model, the accuracy of the sensible heat flux mainly depends on the inversion accuracy of the surface temperature and the canopy height.
The VIIRS data have considerable advantages due to their fine temporal resolutions and free access. Moreover, the expansion and improvement of MODIS data and the next-generation of JPSS platforms will provide a broader application space. Because the Suomi-NPP/VIIRS has been in orbit for several years, the long-term data are promising for applications in monitoring energy budgets.