Mapping Regional Turbulent Heat Fluxes via Assimilation of MODIS Land Surface Temperature Data into an Ensemble Kalman Smoother Framework

Estimation of turbulent heat fluxes via variational data assimilation (VDA) approaches has been the subject of several studies. The VDA approaches need an adjoint model that is difficult to derive. In this study, remotely sensed land surface temperature (LST) data from the Moderate Resolution Imaging Spectroradiometer (MODIS) are assimilated into the heat diffusion equation within an ensemble Kalman smoother (EnKS) approach to estimate turbulent heat fluxes. The EnKS approach is tested in the Heihe River Basin (HRB) in northwest China. The results show that the EnKS approach can estimate turbulent heat fluxes by assimilating low temporal resolution LST data from MODIS. The findings indicate that the EnKS approach performs fairly well in various hydrological and vegetative conditions. The estimated sensible (H) and latent (LE) heat fluxes are compared with the corresponding observations from large aperture scintillometer systems at three sites (namely, Arou, Daman, and Sidaoqiao) in the HRB. The turbulent heat flux estimates from EnKS agree reasonably well with the observations, and are comparable to those of the VDA approach. The EnKS approach also provides statistical information on the H and LE estimates. It is found that the uncertainties of H and LE estimates are higher over wet and/or densely vegetated areas (grassland and forest) compared to the dry and/or slightly vegetated areas (cropland, shrubland, and barren land).


Introduction
Sensible (H) and latent (LE) heat fluxes play a crucial role in the exchange of energy, water, and carbon between the land surface and overlying atmosphere. These fluxes are required in many disciplines such as meteorology, hydrology, ecology, agronomy, and water resources planning and management (Dickinson, 1987;Mascart et al., 1991;Murray et al., 2012;Wang, & Dickinson, 2012).
Retrieval-based techniques invert a physical or empirical model that relates the measured quantity to the variable of interest (Margulis et al., 2005). The needed inputs and types of retrieval-based models vary significantly and range from simple regression approaches to complex physical models. These techniques often generate turbulent heat fluxes estimates only for instances in which observation inputs are available. Data assimilation approaches merge sequences of remote sensing observations with physically based models to estimate turbulent heat fluxes (Caparrini et al., 2004a(Caparrini et al., , 2004b. In contrast to the retrieval-based techniques that can obtain turbulent heat fluxes only at the temporal resolution of remotely sensed observations, the data assimilation approaches provide H and LE estimates between the observations. The data assimilation approaches can be categorized into two main groupings: 1) variational schemes, and 2) ensemble schemes.
Due to the abovementioned drawbacks of VDA approaches, ensemble-based data assimilation techniques (e.g., ensemble Kalman smoother, EnKS) have been used widely in hydrology (Bateni & Entekhabi, 2012b;Dunne & Entekhabi, 2005Lei et al., 2014;Reichle et al., 2008;Xu, Bateni, et al., 2018). Compared to the VDA approaches, the EnKS methods have some unique characteristics: (1) they do not need an adjoint model, and thus they are relatively easy to implement; (2) they can easily generate background-error covariance and directly provide uncertainty of their estimates, and (3) they are able to account for a wide range of possible model and measurement errors (Bateni & Entekhabi, 2012b;Kalnay et al., 2007;Reichle et al., 2002;Whitaker et al., 2009). Bateni and Entekhabi (2012b) assimilated half-hourly ground-based LST measurements into the EnKS system to estimate turbulent heat fluxes at the First International Satellite Land Surface Climatology Project Field Experiment (FIFE) site. In a follow-up study, Xu, Bateni, et al. (2018) estimated H and LE by assimilating hourly LST data from the Geostationary Operational Environmental Satellite into the EnKS system. LST observations with high temporal resolution were used in both of the abovementioned works because they can capture the diurnal cycle of LST. However, high frequency LST data from Geostationary Operational Environmental Satellite have relatively low spatial resolution and are often unavailable in high latitudes . In contrast, LST data from the Moderate Resolution Imaging Spectroradiometer (MODIS) have high spatial resolution (1 km × 1 km) as well as global coverage.
In this study, MODIS LST data (from Aqua and Terra platforms) are assimilated into the heat diffusion equation via the EnKS system to estimate turbulent heat fluxes. The two key unknown parameters of the EnKS approach are the neutral heat transfer coefficient (C HN ) and evaporative fraction (EF). C HN governs the sum of turbulent heat fluxes (H + LE), and EF represents their partitioning. The EnKS approach is tested in the Heihe River Basin (HRB), which samples various hydrological conditions. The turbulent heat fluxes observations from the large aperture scintillometer (LAS) at three sites (namely, Arou, Daman, and Sidaoqiao) in the HRB are used to validate the EnKS estimates. The uncertainty of estimated turbulent heat fluxes from EnKS is also evaluated under different environmental conditions. Finally, performance of EnKS is compared with those of the ensemble open loop (EnOL) (without assimilating MODIS LST data) and VDA approaches. The readers are referred to Xu, He, et al. (2019) for detailed information on the VDA approach.
The objectives of our study are to (1) estimate H and LE by assimilating low temporal resolution LST data from polar orbiting satellites into an EnKS framework; (2) investigate performance of the EnKS framework in the HRB and validate the turbulent heat fluxes estimates with the large aperture scintillometer (LAS)

10.1029/2019EA000705
Earth and Space Science observations; (3) compare the results of the EnKS approach with those of Xu, He, et al. (2019) VDA approach;and (4) assess the uncertainty of EnKS H and LE estimates in different vegetative and hydrologic conditions. This paper is organized as follows. Sections 2 and 3 explains the methods and models including the heat diffusion equation, SEB equation, and the EnKS algorithm. Section 4 describes the experiment sites and data sets. Section 5 provides the results and the discussion. Finally, conclusions are reported in section 6.

System Model
The SEB equation can be written as where H and LE represent the sensible and latent heat fluxes (W m −2 ), and G is the ground heat flux (W m −2 ). R n is the net radiation (W m −2 ) and can be defined as where R ↓ S is the downward shortwave radiation; α is the surface albedo (-), and R ↓ L and R ↑ L are the downward and upward longwave radiations, respectively (W m −2 ).
Sensible heat flux can be represented in terms of the gradient between the LST (T) and air temperature (T a ), where ρ is the air density (kg m −3 ); C P is the air heat capacity (1,012 J kg −1 K −1 ), and U is the referenceheight wind speed (m s −1 ). The bulk heat transfer coefficient (C H ) is related to the neutral bulk heat transfer coefficient (C HN ) and the atmospheric stability correction function (f) via (Abdolghafoorian et al., 2017;Bateni, Entekhabi, & Castelli, 2013;Xu et al., 2014Xu et al., , 2016) where Ri is the Richardson number. C HN depends on the vegetation phenology and is set to vary monthly (Caparrini et al., 2004a(Caparrini et al., , 2004bCrow & Kustas, 2005). C HN constitutes the first unknown parameter of the EnKS system.
The second unknown parameter of EnKS is EF, which is defined as EF is nearly constant for near-peak radiation hours on days without precipitation (Tang & Li, 2017). In this study, EF is assumed to be constant during the assimilation window (0900-1800 LT) in each day, but it can vary from day to day.
The EnKS approach begins with the a priori EF estimate, and improves it via assimilation of LST. According to Xu, He, et al. (2019), the a priori EF estimate is obtained by where EF max and EF min are the maximum and minimum EF values for a specific land cover type (e.g., barren land, grassland, cropland, forest, etc.); φ is the calibration coefficient, and τ is the environmental index. Following Xu, He, et al. (2019), τ is obtained from a leaf area index (LAI) (for vegetated land) or apparent thermal inertia (ATI) (for barren land). The readers are referred to Xu, He, et al. (2019) for more information on τ. The three unknown parameters in equation (6) (i.e., EF max , EF min , and φ) are estimated by the least-square method using EF measurements from 12 eddy covariance (EC) flux towers in the HRB, and corresponding τ values from LAI or ATI. The estimated EF max , EF min , and φ values for four land cover types (i.e., cropland, grassland, forest, and barren land) in the HRB are reported in Table 1. The cropland is mainly covered by seeded corns. The grassland consists mainly of alpine meadow. The forest is a mixture of deciduous broadleaf forest, evergreen needleleaf forests, mixed forests,

10.1029/2019EA000705
Earth and Space Science and shrub. The barren land consists of barren or sparsely vegetated lands. Detailed information on the 12 EC flux towers in the HRB can be found in Xu, He, et al. (2019).

State-Augmented Ensemble Kalman Smoother
The state-augmented EnKS technique contains a forecast step and an update step. In the forecast step, the heat diffusion equation is used to estimate the dynamics of LST. Then, in the update step, the augmented state [i.e., estimated LST values from the heat diffusion equation, and the a priori EF estimates from equation (6)] is improved by assimilating MODIS LST data.

State and Parameter Propagation Models (Forecast Step)
In the forecast step, EnKS propagates an ensemble of LSTs by using the heat diffusion equation. The dynamics of soil temperature at depth z and time t (T[z, t]) is governed by the heat diffusion equation, which is given by where C and K are the soil volumetric heat capacity (J m −3 K −1 s −1 ) and thermal conductivity (J m −1 K −1 ), respectively, and ω is the model error.
Applying the heat diffusion equation requires specification of boundary conditions at the top and bottom of the soil column. The boundary condition at the top of the soil column is obtained by −KdT(z = 0,t)/dz = G(t).

Earth and Space Science
The bottom boundary condition is dT(z = l, t) /dz = 0 (where l is the depth of bottom boundary). In this study, l is taken equal to 0.5 m because the soil temperature at the depth of 0.3-0.5 m is almost constant (Hu & Islam, 1995). The initial soil temperature (i.e., T(z, t = τ 0 ), where τ 0 is the initial time) in the soil column is obtained from reanalysis product (Shi et al., 2014).
To characterize the model error term, the uncertain input variables are perturbed by adding normally distributed random errors (with zero mean and specified variance) to them and generating ensembles. Random errors are created in a physically reasonable range to disturb uncertain inputs when creating ensembles, and reflect uncertainties in measurements, boundary conditions, and propagated states. It is hypothesized that the uncertain inputs are initial profile of soil temperature, heat diffusion coefficient (D = K/C), air temperature, wind speed, downward shortwave radiation, and augmented state vector variables (LST and EF) (Bateni & Entekhabi, 2012b). Following Bateni and Entekhabi (2012b) and Xu, Bateni, et al. (2018), normally distributed errors with a mean zero and a standard deviation of 2 K, 0.1 m 2 s −1 , 1 K, 0.1 m s −1 , 30 W m −2 , and 0.05 are added to the initial soil temperature in the soil column, heat diffusion coefficient, air temperature, wind speed, downward shortwave radiation, and EF, respectively. These numbers are chosen based on the order of magnitude of uncertain inputs (i.e., initial profile of soil temperature, heat diffusion coefficient (D = K/C), air temperature, wind speed, downward shortwave radiation, and augmented state vector variables (LST and EF). For example, uncertainties in the heat diffusion coefficient (D = 0.63/2.58 × 10 −6 = 0.24 × 10 −6 ) controls the heat diffusivity through the soil slab. To account for uncertainties in D and consequently heat diffusivity through the soil slab, a normally distributed noise with a mean of zero and a standard deviation of 0.1 × 10 −6 m 2 s −1 is added to D (Bateni & Entekhabi, 2012b;Xu, Bateni, et al., 2018).
An ensemble of initial conditions T j (z, τ 0 ) (j = 1, … , N e , where j denotes the j th ensemble member, and N e is the number of ensembles) is generated by adding an ensemble of N e model errors to the initial profile of soil temperature. This ensemble is integrated forward in time by the heat diffusion equation to generate the a priori estimates of LST (the so-called forecast or EnOL LST estimates) at each time step in the assimilation window (i.e., 0900-1800 LT). In this study, N e is set to 100 (Bateni & Entekhabi, 2012b;.
In addition to LST, an ensemble of EF is integrated forward in time. Since EF is assumed to be invariant in the assimilation window, it can be propagated by where ω' represents the uncertainty of EF, which is presumed to have a normal distribution. An ensemble of N e model errors (ω') is added to the a priori EF estimate from equation (6) to generate an ensemble of the a priori EF estimates, EF j (j = 1, … , N e ). This ensemble is propagated forward in time (starting from t = τ 0 ) by equation (8) to generate the a priori EF estimates (the so-called forecast EF estimates) at each time step in the assimilation window.

Joint State-Parameter Update Model (Update Step)
The state (LST) and model unknown parameters (C HN and EF) can be updated by the state-augmented EnKS approach in which EF is added into the system state vector (LST), and artificially treated as an additional state variable. The augmented state vector, T, is then defined as where the jth column of T(t) contains the jth realization of the forecast LST [from equation (7)] and EF [from equation (8)] at time t. The MODIS LST observation at time t, T obs (t), is related to the augmented state vector T(t) by the operator H: where H = [1, 0]. The Gaussian measurement error (ε) with a mean of zero and a covariance of R is created and added to LST observations to account for the contribution of observations error to the posterior covariance and avoid correlation among the ensembles (Burgers et al., 1998). Based on Li et al. (2014) and Lu et al. (2018), the standard deviation of 3 K is chosen for perturbing MODIS LST data.
Each ensemble member of the augmented state vector can be updated by (Evensen, 2009) where the different terms are defined by The EnKS framework [i.e., equation (11)] updates the forecast (or the a priori) LST and EF (denoted by superscript f) and generate the so-called analysis (or the a posteriori) LST and EF (denoted by superscript a). T f j and T a j represent the jth ensemble member of forecast and analysis augmented state vector.
The MODIS LST observation at update time t is used in the EnKS system to update the augmented state vector not only at that update time but also at previous times, t' (Evensen, 2009;Lei et al., 2014;Li, Ryu, et al., 2013). K is the Kalman gain matrix, and R is the measurement error covariance. P f represents the error covariance matrix of the forecast model state (LST) and parameter (EF); the superscript T denotes transpose, andX f t ð Þ represents the mean of forecasted augmented state variables at time t.
The state-augmented EnKS algorithm is run with a number of reasonable C HN values during each monthly period. The C HN value that yields the minimum difference between the estimated and observed LST, and the

Study Domain and Data
The HRB is the second largest basin in northwest China (approximately 37.7 to 42.7°N, and 97.1 to 102.0°E). It has complex biological and environmental characteristics, and various landscapes including ice/frozen

10.1029/2019EA000705
Earth and Space Science soil/snow, forests, grasslands, oases, deserts, and lakes. Air temperature and precipitation show a northsouth gradient over the HRB.
The "Heihe Watershed Allied Telemetry Experimental Research" (HiWATER) was established in 2012 (Li, Cheng, et al., 2013Liu et al., 2011Liu et al., , 2018Xu, Guo, et al., 2018). The half-hourly sensible heat flux was measured by the LAS instrument. The net radiation and ground heat flux were measured every 30 min by the four-component radiometer and ground heat flux plate, respectively. Finally, latent heat flux observations were obtained as the residual of the SEB equation (LE = R N -G − H). The HiWATER experiment data can be downloaded from the Heihe Data Archive (http://card.westgis.ac.cn/).
The hourly micrometeorological data including wind speed, air temperature and humidity, atmospheric pressure, and incoming shortwave and longwave radiation are generated by the Weather Research and Forecasting (WRF) model (5 × 5 km) and used as forcing data in the EnKS approach (Pan et al., 2012). Using the nearest method, the forcing data were resampled from 5 to 1 km to match the size of computational pixels (1 × 1 km) over the study domain. The MODIS LST data are available on a daily basis since the launch of the MODIS sensor in 1999 for Terra platform and in 2002 for Aqua platform. In this study, MODIS LST data from Aqua (MYD11A1) and Terra (MOD11A1) platforms with the spatial resolution of 1 × 1 km and 2-revisits during daytime are assimilated into the EnKS system to estimate turbulent heat fluxes. Because of cloud contamination, 65.33% of MODIS LSTs are available during the modeling period. They can be downloaded from the Level-1 and Atmosphere Archive and Distribution System Distributed Active Archive Center archive (https://ladsweb.nascom.nasa.gov/search/). The land cover data were provided by Zhong et al. (2014). LAI and Albedo are obtained from the Global Land Surface Satellite product (Xiao et al., 2014(Xiao et al., , 2016 (http:// glass-product.bnu.edu.cn/). C and K are determined from the soil type and moisture (Chen, 2008;de Vries, 1963). The soil type is found from the HRB Digital Soil Mapping product . The C and K vary from 1.52 × 10 6 to 3.36 × 10 6 J m −3 K −1 s −1 and from 0.37 to 1.21 J m −1 K −1 in the HRB. Soil moisture data are obtained from the Soil Moisture Active Passive (SMAP) L4 product with the spatial resolution of 9 × 9 km (https://search.earthdata.nasa.gov/).

Earth and Space Science
The HRB is discretized into square computational pixels. The size of each computational pixel is 1 km, resulting in 237,573 pixels. The modeling period covers 5 months from 1 May to 30 September 2015 (Day of Year, DOY 121-273), and the daily assimilation window ranges from 0900 to 1800 LT.

Results and Discussions
The retrieved monthly C HN maps from the EnKS approach and corresponding LAI values are shown in Figure 2. As indicated, C HN varies with the vegetation phenology, and its changes are consistent with those of LAI. The C HN values increase from May to July, and decrease from July to September due to the seasonal variations in the vegetation phenology (i.e., LAI). The C HN estimates are generally higher at the vegetated areas (i.e., grassland, cropland, and forest) in the south of HRB than the barren land in the north. In the south of HRB, C HN increases with the growth of the crop and reaches its maximum in July in which LAI is at its peak. There is a sharp drop in LAI in September due to the crop harvest. Analogously, C HN reduces in September. While there is almost no vegetation in the north of HRB and LAI is invariant, the corresponding retrieved C HN values vary in different monthly modeling periods. This is because C HN depends not only on the vegetation phenology (LAI) but also to a lesser extent on the wind speed, friction velocity, and solar elevation (Duynkerke, 1992;Smedman et al., 2007).
Soil moisture controls latent heat flux and is the key indicator of EF (Bateni, Entekhabi, & Castelli, 2013). Therefore, the variations of EF estimates should be consistent with those of soil moisture. Figure 3 shows the spatial patterns of the a priori [obtained from equation (6)] and the a posteriori [updated by the EnKS approach, equation (11)] EF estimates, standard deviation of EnKS EF retrievals, and SMAP soil moisture data over the HRB for DOYs 121, 153, 204, 257, and 273. As shown, the a posteriori EF estimates can capture variations in the SMAP soil moisture more robustly than the a priori EF estimates. For example, the estimated EF values are highest in the southeast of HRB on DOY 204 that soil moisture is at its peak. The higher EF estimates in the south of HRB are due to the high precipitation and dense vegetation. The EF estimates decrease toward the center and north of HRB because of low precipitation and sparse vegetation.
The standard deviations of EnKS EF estimates are shown in Figure 3 as well. The standard deviations of EF retrievals are generally lower in the north of HRB. In the north, evaporation is in its second phase (i.e., water limited); the drying rate is mainly controlled by land surface state variable (i.e., LST) rather than atmospheric factors. Therefore, the coupling between EF and LST is vigorous and the uncertainty of EF estimates is reduced. In contrast, because of the dense vegetation cover and heavy precipitation in the south of HRB, evaporation is in its first phases (i.e., energy limited) and is mainly controlled by the atmospheric state variables (i.e., air temperature and specific humidity). Therefore, the coupling between EF and LST becomes weak, resulting in a higher uncertainty in the EF estimates. By directly providing the standard deviation of the EF estimates, the EnKS approach shows the relative dependence of evaporation on atmospheric factors and surface properties in different environmental conditions. In addition to the different environmental conditions, the standard deviation of the EF estimates also caused by the noise in the forcing data as well as parametric model errors. Figure 4 shows the hourly LST estimates from EnOL (green points) and EnKS (red points) versus the corresponding measurements at the Arou, Daman, and Sidaoqiao sites. The predicted LSTs from the EnKS model are closer to the observations compared to those of the EnOL model. This is because the EnKS approach shifts the LST ensemble members toward the MODIS LST measurements in each update step, and reduces their difference. The statistical indices of hourly LST estimates from the EnKS and EnOL approaches at the three sites are summarized in Table 2. The root mean square errors (RMSEs) of EnOL LST estimates at the Arou, Daman, and Sidaoqiao sites are 4.12, 3.25, and 2.99 K, respectively. The EnKS approach reduces the

10.1029/2019EA000705
Earth and Space Science abovementioned RMSEs by 10.19%, 19.38%, and 28.43%, respectively. It is worth mentioning that the EnKS LST estimates have the spatial resolution of 1 × 1 km (i.e., the size of computational pixels in the HRB), while LST data at the Arou, Daman, and Sidaoqiao site are in situ point-scale measurements. The abovementioned inconsistency in the spatial resolutions, simplistic assumptions (e.g., constant daily EF, constant monthly C HN , constant soil thermal conductivity and heat capacity, etc.), and measurement errors cause misfits   (Zupanski et al., 2005). The histogram of normalized innovations can be used as a means to evaluate the performance of EnKS, and realize whether the observation and model errors are selected properly (De Lannoy et al., 2010;Kumar et al., 2008). If the EnKS's Gaussian assumption is met, the normalized innovations will have the standard normal distribution with the zero mean and unit variance, N(0, 1) (Crow, 2003;Margulis et al., 2002). The deviation of normalized innovations from the standard normal distribution implies that EnKS is suboptimal (Bateni & Entekhabi, 2012b).
The histogram of normalized innovations over the HRB is compared with N(0, 1) in Figure 5. The distribution of LST normalized innovations has a small negative bias over the HRB. This nonzero mean indicates that the input and output error statistics are not consistent. Hence, it suggests that the defined model and observation errors are not optimal. It might be also related to auto-correlated observation errors. The histogram of normalized innovations in Figure 5 also shows that the forecast error covariance (HP f H T +R) −1/2 is slightly overestimated. The overestimated forecast error covariance leads to a slightly higher peak for normalized innovations compared to the normal distribution (Reichle et al., 2002;Zupanski et al., 2005). This is most likely due to the inclusion of nonoptimal measurements and model error parameters. Overall, the normalized innovations of LST are close to the white noise, N(0, 1), over the HRB, implying that the spread of ensembles is adequate, the model estimates are reliable and almost near optimal (Hol et al., 2008;De Lannoy et al., 2010;Zupanski et al., 2005). At the same time, the small inconsistency between the normalized innovations histogram and the standard normal distribution indicates that the EnKS model is not truly optimal. Figure 6 shows time series of daytime-averaged (0900-1800 LT) observed (open circles) and estimated sensible and latent heat fluxes from the EnKS (blue solid lines) and EnOL (red dashed lines) approaches at the Arou, Daman, and Sidaoqiao sites during DOY 121-273, 2015. The standard deviations of turbulent heat fluxes estimates from the EnKS and EnOL methods are also shown by blue and pink bands, respectively. As indicated, the EnKS H and LE estimates are closer to the observations compared to those of EnOL. This shows that the EnKS approach can extract the implicit information contained in the MODIS LST data to improve the EnOL H and LE retrievals. Also, lower H and LE uncertainties are obtained by EnKS, implying that EnKS can move the ensemble members toward their true values. Figure 6 also indicates that uncertainties of H and LE estimates are larger in the wet period (DOYs 180-240) than the dry period (DOYs

10.1029/2019EA000705
Earth and Space Science fairly well with the observations and mainly fall around the 1:1 line. The misfit between the model estimates and observations might be due to the errors in the forcing data (e.g., air temperature, humidity, incoming solar radiation, and wind speed), misspecification of model parameters (e.g., albedo, emissivity, soil thermal conductivity, soil heat capacity, etc.), and simplistic model assumptions (i.e., daily constant EF and monthly constant C HN ). Figure 7 also shows that LE estimates are more scattered compared to H retrievals. This is because of more sources of errors in the LE estimates as discussed above. LE observations are also obtained as the residual of the SEB equation (i.e., LE = R N − H − G), which can yield erroneous LE data .
The statistical indices of hourly turbulent heat fluxes estimates from the EnOL, EnKS, and VDA approaches at the three sites are summarized in Table 3

Earth and Space Science
The standard deviation of H and LE estimates from the EnOL and EnKS approaches at the three experimental sites are summarized in Table 4. As indicated, the standard deviation of estimated H and LE values from the EnKS approach is lower than that of EnOL, implying that the EnKS approach is able to use information in the MODIS LST data in order to generate less uncertain H and LE values. Table 4 also shows that the uncertainty of H and LE estimates is larger at the Arou site (wet) compared to the Sidaoqiao site (semiarid). Figure 8 shows the monthly mean H and LE estimates from the EnKS method over the HRB during the growing season (May-September). LE increases from May to July and then decreases from July to September over the HRB, which is consistent with the temporal variations in precipitation and LAI. As can be seen, there is a sharp variation in the H and LE estimates from north to south of HRB. The H (LE) values in barren areas in the north of HRB are higher (lower) than those of grasslands, croplands, and forest in the south. LE is low in the north and center of HRB due to rare precipitation and sparse vegetation cover, but it increases significantly toward the south because of high precipitation and dense vegetation cover. High LE values in the center of HRB are associated with crop growth and irrigation in the cropland areas. Figure 9 shows the standard deviations of H and LE estimates from EnKS over the HRB during the growing season. As expected, the standard deviations of H and LE estimates are lower in the north and center of the HRB due to low precipitation and sparse vegetation cover, but they increase toward the south of the HRB because of high precipitation and dense vegetation cover. Figure 10 shows the monthly mean estimated LST from EnKS (first row) and observed LST from MODIS (second row) over the HRB during the growing season (May-September). Monthly mean MODIS LST maps are obtained by averaging MODIS LST observations. MODIS LST data are available twice daily only in clearsky conditions (Bisht et al., 2005). The estimated LST values are consistent with the observations. Because of seasonal variation in solar radiation (not shown here), LST increases from May to July and decreases from July to September. The higher latent heat flux in the south keeps the surface relatively cool whereas the areas to the center and north with lower latent heat flux show higher temperature.
The third row in Figure 10 shows the relative difference of monthly mean estimated and observed LST. In the north and center of the HRB (dry and slightly vegetated conditions), the relative difference of estimated and observed LST is small, but it increases toward the south of the HRB (wet and densely vegetated conditions). The higher relative difference in the south is because of the heavy precipitation and dense vegetation that make the retrieval of EF from MODIS LST data more uncertain.
The HRB contains various vegetative and hydrological conditions, which make it possible to assess the uncertainty of EnKS estimates in different environmental conditions. Figure 11 shows the statistical boxplots of area-averaged standard deviations of evaporative fraction (EF Std.), LST (LST Std.), sensible heat flux (H Std.), and latent heat flux (LE Std.) over five land cover types (cropland, forest, grassland, shrub land, and barren land) in the HRB. The statistical boxplots of LAI and soil moisture data over the same five land over types are also shown in Figure 11. The lower and upper box edges are the 1st and 3rd quartiles in the sample distribution, respectively. The median position is marked within the box. Lines extending from the box ends indicate the 5th and 95th percentiles. As expected, the standard deviations (uncertainties) of EF, LST, H, and LE estimates for grassland (with higher soil moisture) and forest (with higher LAI) are larger than those of cropland, shrubland, and barren land. The lowest standard deviations can be found in the barren land (with the lowest soil moisture and LAI).
The MODIS global ET product provides ET data with the spatial resolution of 1 × 1 km every 8 days (Mu et al., 2007). Figure 12 compares the 8-day ET estimates from EnKS (red points) and MODIS (green Since the MODIS ET product has a 8-day temporal resolution, the daily LAS observations and EnKS estimates are summed to produce 8-day ET values. The results show that the EnKS ET estimates are closer to the in situ measurements. The statistical indices of 8-day ET retrievals from EnKS and MODIS product at the three sites are summarized in Table 5. The three-site-averaged MAPE (RMSE) from the EnKS approach is 16.45% (7.18 mm), which is 65.54% (48.53%) lower than the MAPE (RMSE) of 47.73% (13.95 mm) from the MODIS ET product.

Conclusion
Remotely sensed LST data from Moderate Resolution Imaging Spectroradiometer (MODIS) are assimilated into the heat diffusion equation within the ensemble Kalman smoother (EnKS) approach to estimate turbulent heat fluxes. The EnKS approach is tested over the HRB with contrasting vegetative and hydrological conditions. The key unknown parameters of the EnKS approach are the neutral heat transfer coefficient (C HN ) and evaporative fraction (EF). The results show that the spatiotemporal patterns of C HN estimates agree well with those of vegetation phenology. The patterns in the estimated evaporative fraction (EF) maps resemble those of the SMAP soil moisture product.
The EnKS H and LE estimates are validated with the LAS observations at the Arou (grassland), Daman (cropland), and Sidaoqiao (shrub-forest) sites in the south, center, and north of HRB, respectively. The hourly and daytime-averaged sensible and latent heat flux estimates from the EnKS model agree relatively well with observations at the three experimental sites. For the EnKS approach, the three-site-averaged RMSEs of hourly sensible and latent heat fluxes are 46.27 W m −2 and 103.82 W m −2 , respectively, which are 24.35% and 23.57% lower than those of EnOL. These outcomes imply that the EnKS approach can take advantage of implicit information in the MODIS LST data to improve the EnOL H and LE estimates. The performance of EnKS is compared with the VDA approach. It is observed that the VDA approach performs slightly better than EnKS.
EnKS can directly generate the uncertainties of its estimates. The uncertainties of H and LE estimates over the wet and densely vegetated areas in the south of HRB are higher than the dry and slightly vegetated areas in the center and north. In dry and slightly vegetated areas, the drying rate is mainly controlled by the land surface state variable (i.e., LST). Therefore, the coupling between EF and LST becomes vigorous and the uncertainties of turbulent heat fluxes estimates are reduced. In contrast, in the wet and densely vegetated sites, the drying rate is mainly affected by atmospheric state variables (i.e., air temperature and specific humidity) and coupling between EF and LST becomes weak, resulting in higher H and LE uncertainties.