Abstract

Estimation of turbulent heat fluxes by assimilating sequences of land surface temperature (LST) observations into a variational data assimilation (VDA) framework has been the subject of numerous studies. The VDA approaches are focused on the estimation of two key parameters that regulate the partitioning of available energy between sensible and latent heat fluxes. These parameters are neutral bulk heat transfer coefficient CHN and evaporative fraction (EF). The CHN mainly depends on the roughness of the surface and varies on the time scale of changing vegetation phenology. The existing VDA methods assumed that the variations in vegetation phenology over the period of one month are negligible and took CHN as a monthly constant parameter. However, during the growing season, bare soil may turn into a fully vegetated surface within a few weeks. Thus, assuming a constant CHN may result in a significant error in the estimation of surface fluxes, especially in regions with a high temporal variation in vegetation cover. In this study the VDA approach is advanced by taking CHN as a function of leaf area index (LAI). This allows the characterization of the dynamic effect of vegetation phenology on CHN. The performance of the new VDA model is tested over three sites in the United States and one site in China. The results show that the new model outperforms the previous one and reduces the root-mean-square error (and bias) in sensible and latent heat flux estimates across the four sites on average by 31% (61%) and 21% (37%), respectively.

1. Introduction

Accurate estimation of turbulent heat fluxes (i.e., sensible and latent heat fluxes) is of significant importance in numerous atmospheric, hydrologic, meteorological, climatological, and environmental processes (Alfieri et al. 2009; Stephens et al. 2012). The sum of sensible H and latent (LE) heat fluxes (i.e., H + LE) and their relative partitioning (i.e., H/LE) affects the development of the boundary layer and forces the dynamics of the troposphere. The turbulent fluxes impact weather and climate by modifying the spatiotemporal variability of rainfall (Andrews et al. 2009) and influencing the development of boundary layer clouds (Lewellen and Lewellen 2002).

In situ measurements of heat and moisture fluxes are costly (Baldocchi et al. 2001) and therefore are available only over a limited number of field experiment sites [e.g., Southern Great Plains (SGP), Heihe Watershed Allied Telemetry Experimental Research (HiWATER), First International Satellite Land Surface Climatology Project (ISLSCP) Field Experiment (FIFE) and flux tower stations (e.g., AmeriFlux, Fluxnet, and EuroFlux)].

Among numerous techniques developed to make quantitative estimates of turbulent heat fluxes (e.g., Anderson et al. 1997; Bastiaanssen et al. 1998; Allen et al. 2011; Kalma et al. 2008; Polonio and Soler 2000; Yilmaz et al. 2014), the variational data assimilation (VDA) approaches (Peters-Lidard et al. 2011; Bateni and Entekhabi 2012; Sini et al. 2008; Reichle et al. 2001b) have recently gained a substantial amount of attention. The VDA approaches estimate land surface turbulent heat fluxes by assimilating land surface temperature (LST) observations into the physical land surface models. The uncertainty of the estimated fluxes can be quantified using a Hessian-based approach (Abdolghafoorian and Farhadi 2016). The main characteristics of the VDA schemes that make them suitable for estimation of turbulent heat fluxes are that 1) they do not need any calibration and empirical fitting (Caparrini et al. 2003), 2) they estimate heat fluxes even for instances in which LST observations are not available or there are data gaps (Lewis et al. 2012), and 3) they retrieve ground heat flux without applying the empirical formulations that parameterize ground heat flux as a fraction of net radiation (Campo et al. 2009).

The two main unknown parameters of the VDA approaches for estimating turbulent fluxes are neutral bulk heat transfer coefficient CHN and evaporative fraction (EF). The CHN scales the sum of turbulent heat fluxes and EF scales their partitioning. The existing VDA methods (e.g., Caparrini et al. 2004a, 2003; Bateni et al. 2013b; Sini et al. 2008; Xu et al. 2014) assume that changes in vegetation dynamics over the period of one month are negligible, and thus take CHN as a monthly constant parameter. However, bare soil may turn into a fully vegetated surface in only a few weeks during the growing season, undermining the assumption of constant monthly CHN. In this study the VDA approach for land surface flux estimation is advanced by taking CHN as a function of vegetation phenology [i.e., leaf area index (LAI)]. This allows the characterization of the variations in CHN as LAI changes and hence improves the performance of the VDA approach, especially at sites in which temporal variation of vegetation phenology is significant. Moreover, the period of assimilation of the VDA approach is no longer constrained by the monthly periods in which CHN is assumed to be constant.

Using in situ measurements, Sugita and Brutsaert (1990), Kubota and Sugita (1994), Betts and Ball (1998), and Mahrt and Vickers (2004) found a strong relationship between CHN and LAI. In addition to these experimental studies, a variety of physically based models have been developed that describe the influence of LAI on CHN (Massman 1999; Hasager and Jensen 1999; Su 2002; Timmermans et al. 2013). The results of the previous VDA models also indicate a relationship between CHN estimates and seasonal land-cover phenology (Caparrini et al. 2004a; Bateni et al. 2014). All of these studies show that the evolution of LAI has a significant effect on CHN.

This study augments the Bateni et al. (2013b) VDA model (hereafter called VDAold) by taking CHN as a function of LAI rather than assuming it to be a constant monthly parameter. The CHN–LAI relationship introduced by Farhadi et al. (2014, 2016) and implemented in their proposed methodology of conditional sampling of surface states for estimation of turbulent heat fluxes is used in this study to characterize variations in CHN. The applicability and performance of the new VDA approach (hereafter called VDAnew) is tested at point scale using field site information at three sites in the United States (Audubon, Arizona; Brookings, South Dakota; and Bondville, Illinois) and one site in China (Daman). These sites are chosen because they represent various vegetation types (i.e., grassland, cropland, and shrubland), and LAI at these sites varies significantly in time (i.e., from 0.4 to 4.3). The wide range of LAI values and various vegetation types provide a unique opportunity to assess the effectiveness of taking CHN as a function of LAI in the VDA approach for estimation of turbulent heat fluxes.

This paper is organized as follows. Section 2a shows different components of the surface energy balance (SEB) equation and its key unknown parameters (i.e., CHN and EF). Section 2b describes the heat diffusion equation as the physical constraint of the model. Section 2c introduces the VDA framework for estimation of the key unknown parameters. Section 3 explains the data and sites in which the VDAnew scheme is tested. Results of this study and conclusions are given in sections 4 and 5, respectively.

2. Methodology

a. SEB scheme

The core of the model is based on the combined source SEB model, which treats soil and vegetation as only one effective medium, and is given by

 
formula

where H, LE, and G are the sensible, latent, and ground heat fluxes, respectively, and Rn is the net radiation. Sensible heat flux can be represented in terms of the gradient of temperature between the land surface temperature (composite of soil and vegetation temperature at the ground surface) and air above it Ta:

 
formula

where U is the wind speed, cp is the specific heat of air, ρ is the air density, CHN is the bulk heat transfer coefficient under neutral atmospheric condition, and f is the atmospheric stability correction function, which itself is a function of Richardson number (Ri). The stability correction function proposed by Caparrini et al. (2004b), which has performed satisfactorily in previous studies (e.g., Crow and Kustas 2005; Farhadi et al. 2014; Bateni et al. 2013b; Farhadi et al. 2016), is used in this study. The CHN captures the surface boundary influence on near-surface turbulence and is essential for accurate estimation of sensible heat flux. It can be related to the roughness length scales for heat ZOH and momentum ZOM via (Andreas and Murphy 1986; Ann‐Sofi et al. 2007):

 
formula

where k is the von Kármán constant and Zref is the reference height. Variables ZOH and ZOM are related via and B is the Stanton number (Verhoef et al. 1997). Substituting into (3) leads to

 
formula

Brutsaert (1979) parameterized ZOH and showed its dependency on LAI. Sugita and Brutsaert (1990) obtained an exponential relationship between ZOH and LAI over the FIFE site. A similar exponential relationship was found by Kubota and Sugita (1994) in the energy and water balance experimental field of the Environmental Research Center (ERC) at the University of Tsukuba. Using experimental data, Matsushima (2005) parameterized ZOH and showed that it increases as LAI grows. Variable kB−1 also has been the subject of several theoretical and observational studies. A number of studies (e.g., McNaughton and Hurk 1995; Lhomme et al. 2000; Mahrt and Vickers 2004) showed empirically that kB−1 is mainly a function of LAI. In addition to the empirical studies, several physically based models were developed to characterize the influence of canopy (LAI) on kB−1 (Hasager and Jensen 1999; Su 2002; Timmermans et al. 2013). All of these studies along with (4) indicated that CHN principally depends on LAI. Recently, Farhadi et al. (2014, 2016) parameterized the influence of LAI on CHN via

 
formula

The parameterization is limited to the first-order seasonality effects with constant parameters a and b but seasonally varying LAI. Parameter b is indicative of neutral bulk heat transfer coefficient for bare soil (when LAI is zero), and a characterizes the effect of vegetation density on CHN. It should be noted that the proposed CHN–LAI relationship is consistent with the literature (e.g., Kubota and Sugita 1994; Jensen and Hummelshøj 1995; Qualls and Brutsaert 1996) and performed well in estimation of heat fluxes at several field sites (Farhadi et al. 2014) and in mapping of heat fluxes at regional scale using remote sensing data (Farhadi et al. 2016). This relationship is used herein to characterize the effect of LAI on CHN. The coefficients a and b are the constant unknown parameters of the SEB and will be obtained via the VDAnew approach.

Evaporative fraction is the third unknown parameter of SEB scheme and is defined as the ratio of latent heat flux to the sum of turbulent heat fluxes:

 
formula

Evaporative fraction is almost constant for near-peak radiation hours on days without precipitation (Cragoa and Brutsaert 1996; Gentine et al. 2007). Using (6), LE can be estimated as a function of EF and H.

b. Heat diffusion equation

The heat diffusion equation in a combined source SEB model represents the time evolution of ground temperature at different depths z within the column of soil over time t and is given by

 
formula

where c is the soil volumetric heat capacity and p is the soil thermal conductivity (Campbell 1985). These variables are assumed to be constant during the assimilation period, with values that are determined based on soil texture and soil moisture at each investigated site (Hopmans et al. 2002; Campbell 1985). The boundary condition at the top of the soil column, T (z = 0, t), is obtained by solving the surface forcing equation , where G is the ground heat flux. The ground temperature at the surface T (z = 0, t) represents the composite of soil and vegetation temperature (i.e., LST). The diurnal variation of heat wave is negligible beyond depths 0.3–0.5 m (Hirota et al. 2002); thus, a Neumann boundary condition is implemented at the bottom of the soil column [ is the depth of the lower boundary of the soil column] (Bateni et al. 2013b).

c. VDA system

The key unknown parameters of the VDAnew system (i.e., EF, a, and b) are obtained by assimilating land surface temperature observations Tobs inside the assimilation window [t1 = 0900 local time (LT), t2 = 1600 LT] into the VDA framework. Evaporative fraction varies on a daily basis, while a and b are invariant during the assimilation period of N days. The VDA approach is based on the minimization of the cost function J that represents the aggregated error over the entire assimilation period:

 
formula

The first term in the cost function accounts for the misfit between the LST estimates from the heat diffusion equation and observations over the assimilation period (N days). The second, third, and fourth terms are the squared error of the unknown parameters EF, a, and b with respect to their last best estimates (represented by primed variables), respectively. The last term in the cost function is the physical constraint (i.e., the heat diffusion equation), which is adjoined to the model via the Lagrange multiplier λ. Parameters CT, Ca, Cb, and CEF are constant coefficients that weigh the terms in the cost function and control the rate of convergence (Reichle et al. 2001a). The optimal values of the unknown parameters (i.e., EF, a, and b) are found by minimizing the cost function. The minimum value of J is obtained by setting its first variation equal to zero (δJ = 0). Setting δJ to zero leads to the following set of the so-called Euler–Lagrange equations, which should be solved simultaneously:

 
formula
 
formula
 
formula
 
formula
 
formula
 
formula
 
formula

where σ is the Stefan–Boltzmann constant and εs is the surface emissivity.

The optimal values of the unknown parameters are found via an iterative process, which is initiated by solving the heat diffusion (7) forward in time to estimate the LST during the assimilation window (t1, t2) in which EF is almost constant. It continues by solving the adjoint model (9a) backward in time using the terminal and boundary conditions [(9b)(9d)] to obtain . Once the state variable (i.e., T) and the Lagrange multiplier (i.e., λ) are calculated, the VDAnew scheme updates the unknown parameters via (10)(12). The iterative process is terminated when the unknown parameters converge to their optimal values. Parameters CT, Ca, Cb, and CEF are taken equal to 100 K2, 0.001, 0.001, and 0.001 based on several sensitivity analysis tests and previous relevant studies (Bateni et al. 2013b; Xu et al. 2014). Setting CT = 100 K2 and Ca = Cb = CEF = 0.001 leads to a stable iterative process that converges to an optimal solution after a reasonable iteration number. The initial guess of b can be chosen based on the bare soil roughness from the acceptable range (i.e., −7 < b < −5) reported in the literature (Farhadi et al. 2016; Bateni et al. 2013a). The initial guess of a is chosen from the range 0 < a < 1 depending on vegetation density. The results of sensitivity analysis (not shown) revealed that as long as the initial guesses of a and b are chosen from the physically reasonable range, the VDA approach will generate reliable results.

Although the use of the LAI measurements as input forcing is key to the VDAnew scheme, the variation of LAI (not its magnitude) during the assimilation period is crucial for characterizing the effect of vegetation dynamics on the CHN. The optimized parameters a and b in the VDA scheme compromise for the instrumental errors of the measured LAI.

3. Study locations and data

Performance of the VDAnew model is tested at four sites: Audubon, Brookings, and Bondville in the United States and Daman in China. These sites are chosen because they represent a diverse range of vegetation types (i.e., grassland, cropland, and shrubland), and LAI in these sites varies significantly within a short period of time. Audubon, Brookings, and Bondville are part of the AmeriFlux network (http://ameriflux.ornl.gov/). Audubon is an open shrubland with sandy clay loam soil type. The region has a temperate arid climate. Brookings has humid continental climate with hot summers and cold winters. The vegetation cover is a mixture of C3 and C4 grazed grassland and the soil type is well drained clay loams. Bondville has temperate continental climate and is covered with corn (C4) and soybeans (C3). The soil type of this site is silt loam. The AmeriFlux archive contains half-hourly micrometeorological data (e.g., incoming solar radiation, air temperature, specific humidity, and wind speed) and surface energy fluxes. Sensible and latent heat fluxes were measured with the eddy covariance approach. Surface net radiation and soil heat flux were obtained using net radiometer and soil heat flux plate, respectively (Richardson et al. 2006). Surface temperature is measured by an infrared thermometer.

The LAI data with spatial resolution of 1 km and temporal revisit of 8 days were obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS), available from the Oak Ridge National Laboratory Distributed Active Archive Center (http://daac.ornl.gov/MODIS/modis.html). To reduce the geolocation and pixel-shift errors, the spatial aggregation from 7 × 7 pixels around the tower pixel rather than the mere usage of the tower pixel is used to obtain LAI values. The average spatial variations of LAI for Audubon, Brookings, and Bondville are 0.04, 0.26, and 0.19, respectively. Validation of the VDAnew model is done over the periods in which more than 95% of the MODIS pixels have acceptable quality (see Table 1).

Table 1.

Characteristics of the study sites (for more information on U.S. sites, see http://fluxnet.ornl.gov/ and http://ameriflux.ornl.gov/).

Characteristics of the study sites (for more information on U.S. sites, see http://fluxnet.ornl.gov/ and http://ameriflux.ornl.gov/).
Characteristics of the study sites (for more information on U.S. sites, see http://fluxnet.ornl.gov/ and http://ameriflux.ornl.gov/).

Daman is a wet, densely vegetated site covered by seeded corn. It is located in the middle reach of Heihe River basin, China. In 2012, a multiscale observation experiment was carried out from June to September in the framework of HiWATER (Li et al. 2013). The automatic meteorological stations recorded shortwave and longwave radiation (in both upward and downward directions), net radiation, air temperature, precipitation, wind speed, atmospheric pressure, soil heat flux, soil moisture, and land surface temperature (measured by vertically downward infrared temperature sensors) every 30 min (Wu et al. 2015). The sensible and latent heat fluxes are obtained from the Daman superstation eddy covariance system. At the Daman site, LAI was measured by an LAI-2000 scanner during the field experiment with a spatial resolution of about 1 km (Li et al. 2013). LAI varies rapidly from 0.44 to 4.3 during early growing stages of maize (Julian days 159–188). During the next last two months of the HiWATER experiment, LAI variation becomes negligible. Hence, at the Daman site, the VDAnew system is applied only to the early stages of growing season (Julian days 159–188) when LAI prominently affects CHN.

Many studies (e.g., Behera et al. 2010; Fang et al. 2012; Fensholt et al. 2004) have evaluated the MODIS LAI product and LAI-2000 Plant Canopy Analyzer measurements against direct field measurements, where LAI values are derived by destructive sampling. These studies conclude that both MODIS and LAI-2000 scanner retrieve LAI with good accuracy.

The energy balance ratio [EBR = (H + LE)/(RnG)] for each study site is calculated over the assimilation period. The daytime (0900–1600 LT) EBR for Audubon, Brookings, Bondville, and Daman are 0.758, 0.783, 0.782, and 0.90, respectively. Whenever the EBR is not within the accepted range, the Bowen ratio–energy balance method (Perez et al. 1999) is applied to partition available energy between sensible and latent heat.

The location, vegetation type, and range of LAI at the four field sites are listed in Table 1. The growing season at the Daman site and both the growing and senescence seasons at the Audubon, Brookings, and Bondville sites are examined. Among the four sites, the Daman (Audubon) site experiences the most (least) variation in vegetation.

The observed ln(CHN), which is obtained by substituting the observed sensible heat flux, wind speed, land surface, and air temperature in (3), versus LAI is shown in Fig. 1 for the four investigated sites. Since MODIS LAI products are available every 8 days, daily LAI data are computed via linear interpolation. As expected, CHN increases when vegetation becomes denser. The correlation coefficient between observed ln(CHN) and LAI (herein called robs) is shown in Fig. 1 for each site. A number of studies showed that ZOH and kB−1 and, consequently, CHN depend on LAI and, to a lesser extent, on wind speed, friction velocity, and solar elevation (e.g., Sugita and Brutsaert 1990; Su 2002; Timmermans et al. 2013; Colaizzi et al. 2006). It is hypothesized that these factors explain why the correlation coefficient in Fig. 1 varies among the sites. Overall, the positive relationship between ln(CHN) and LAI is consistent with the findings of other studies (Kubota and Sugita 1994; Farhadi et al. 2014; Massman 1999; Boulet et al. 2012; Farhadi et al. 2016). Figure 1 also indicates that CHN varies significantly at some of the investigated sites only over a short time period. For example, at the Daman site, CHN increases by a factor of about 20 from 0.001 to 0.02 during the growing season. These large variations in CHN occur over only a 30-day period (Julian days 159–188) in which LAI varies from 0.4 to 4.3. This reveals that the assumption of monthly constant CHN may lead to large errors and adversely affect performance of the VDA scheme.

Fig. 1.

(left) Time series of LAI (solid line), observed ln(CHN) (dashed line), and raw LAI measurement (asterisks). (right) Time series of observed ln(CHN) vs LAI at the four study sites.

Fig. 1.

(left) Time series of LAI (solid line), observed ln(CHN) (dashed line), and raw LAI measurement (asterisks). (right) Time series of observed ln(CHN) vs LAI at the four study sites.

4. Results

a. Land surface temperature

Figure 2 shows the root-mean-square error (RMSE) of LST estimates TRMSE from the VDAnew model versus iteration number at the four study sites. As indicated in Fig. 2, the error between observed and estimated LST decreases rapidly with iteration number and TRMSE reaches an asymptotic value at a certain iteration number. Theoretically, the misfit between observed and estimated LST should approach zero after a certain number of iterations. However, there is a deviation from zero that is due to the simplifying assumptions (e.g., constant daily EF, constant soil thermal conductivity, and constant heat capacity as well as errors in input data) and/or measurement errors.

Fig. 2.

Values of TRMSE vs iteration number.

Fig. 2.

Values of TRMSE vs iteration number.

b. Neutral bulk heat transfer coefficient and evaporative fraction

As mentioned above, CHN is one of the key unknown parameters of the SEB equation. In the VDAnew approach, CHN is treated as a function of LAI and thus varies on a daily basis [see (5)], while the VDAold approach assumes CHN is constant over the period of one month. The estimated CHN values by the VDAold and VDAnew approaches for all the investigated sites are shown in Table 2. Unlike the VDAold that provides only an average CHN value during the assimilation period, this study (VDAnew approach) accounts for the evolution of CHN within the assimilation period. As indicated in Table 2, CHN estimates from VDAold fall within the range of variability of CHN retrievals from VDAnew. This implies that VDAnew robustly captures variations in CHN as LAI changes via reliable estimation of the unknown parameters of CHN (i.e., a and b). The estimated values of a and b are also presented in Table 3.

Table 2.

Estimated CHN values from the VDAnew and VDAold approaches for the growing and senescence seasons.

Estimated CHN values from the VDAnew and VDAold approaches for the growing and senescence seasons.
Estimated CHN values from the VDAnew and VDAold approaches for the growing and senescence seasons.
Table 3.

Estimated b and a values from VDAnew approaches.

Estimated b and a values from VDAnew approaches.
Estimated b and a values from VDAnew approaches.

The other unknown parameter of the SEB equation is EF. Figure 3 shows time series of EF estimates from the VDAnew model at all the sites. For comparison, retrieved EF values from the VDAold model and observations are also shown in the same Fig. 3. The retrieved EF values from VDAnew are in good agreement with observations in terms of both magnitude and day-to-day dynamics. The VDAnew model also yields more accurate EF estimates compared to VDAold. This is due to the fact that CHN (which is simultaneously retrieved with EF) is estimated more accurately by the VDAnew method since its variation with vegetation phenology is captured. This in turn results in more accurate estimation of the second unknown parameter of the energy balance equation (i.e., EF). A lower LAI (e.g., at early stages of growing season) corresponds to a lower CHN estimate in the VDAnew model and hence smaller sensible heat flux. This yields higher EF estimates that are closer to observations than those of the VDAold model. The assumption of constant monthly CHN in VDAold results in overestimation (underestimation) of H in sparse (dense) vegetation conditions.

Fig. 3.

Time series of daily EF estimates from the VDAnew and VDAold models at the Daman, Brookings, Bondville, and Audubon sites. Corresponding observations are shown by circles. In VDAold, CHN is assumed constant during the interval between two vertical dotted lines.

Fig. 3.

Time series of daily EF estimates from the VDAnew and VDAold models at the Daman, Brookings, Bondville, and Audubon sites. Corresponding observations are shown by circles. In VDAold, CHN is assumed constant during the interval between two vertical dotted lines.

The RMSE and bias of EF estimates from the VDAnew model are compared to those of the VDAold model in Table 4. At almost all of the sites, RMSE and bias of EF estimates from VDAnew are significantly less than those of VDAold. The four-site average of RMSE (absolute bias) of estimated EF from VDAnew and VDAold are 0.10 (0.01) and 0.14 (0.03), respectively. Reduction in misfit between the measured and estimated daily EF is interpreted as an improvement in the model performance. To indicate the effectiveness of VDAnew in improving EF estimates, the percent improvement in RMSE of EF estimates [RMSE–EFimp (%)] is shown in Table 5 as a function of robs and the average daily variation of LAI . The value of for each site is approximated as follows:

 
formula
Table 4.

Comparing the performance of VDAold and VDAnew at the four study sites. RMSE and bias from VDAold are shown in parentheses.

Comparing the performance of VDAold and VDAnew at the four study sites. RMSE and bias from VDAold are shown in parentheses.
Comparing the performance of VDAold and VDAnew at the four study sites. RMSE and bias from VDAold are shown in parentheses.
Table 5.

Relationship between percent of RMSE reduction, robs, and . Percent of biases reduction are shown in parentheses. Asterisk indicates that bias is very small for both VDAold and VDAnew.

Relationship between percent of RMSE reduction, robs, and . Percent of biases reduction are shown in parentheses. Asterisk indicates that bias is very small for both VDAold and VDAnew.
Relationship between percent of RMSE reduction, robs, and . Percent of biases reduction are shown in parentheses. Asterisk indicates that bias is very small for both VDAold and VDAnew.

The percent improvement in RMSE of estimated variable X for each site is obtained as follows:

 
formula

The value of RMSE–EFimp (%) is obtained from (14) by replacing X with EF. As expected, sites with higher robs and tend to show more improvement in EF estimates when VDAnew is used. For example, at the Daman site, which has the highest robs and values compared to the other field sites, taking CHN as a function of LAI (rather than assuming it to be constant) showed the most improvement in EF retrieval. The performance of the data assimilation system is significantly improved at this site and a reduction of 50% in RMSE of EF estimates is obtained by using VDAnew instead of VDAold. On the other hand, since the vegetation variation at the Audubon site during the assimilation period is not considerable, adding LAI information to VDA only slightly improves the RMSE of EF estimates.

By replacing the variable “RMSE” with “bias” in (14), the improvement in bias for the estimated variable X by the new VDA system is obtained. The four-site average of bias–EFimp (%) is 63% and the largest improvement in terms of bias of EF estimates [bias–EFimp (%)] is found at the Daman site (108%).

c. Sensible and latent heat fluxes

Time series of daytime-averaged estimated sensible and latent fluxes from the VDAold and VDAnew models as well as corresponding measurements are shown in Figs. 4 and 5, respectively. The VDAnew model considers daily variation of surface roughness by taking CHN as a function of LAI, while the VDAold model assumes vegetation phenology is constant over the period of one month. Hence, VDAnew is expected to outperform VDAold in estimation of sensible and consequently latent heat fluxes, especially for sites in which LAI varies significantly over the assimilation period. For example, vegetation shows a rapid growth in the beginning of the assimilation period at the Brooking, Daman, and Bondville sites (see Fig. 1). Correspondingly, the VDAnew model yields much better estimates of CHN and EF and, consequently, sensible and latent heat fluxes compared to the VDAold model, in that period at the abovementioned sites.

Fig. 4.

Time series of daytime-averaged sensible heat fluxes at the Audubon, Daman, Brookings, and Bondville sites. In VDAold, CHN is assumed constant during the interval between two vertical dotted lines.

Fig. 4.

Time series of daytime-averaged sensible heat fluxes at the Audubon, Daman, Brookings, and Bondville sites. In VDAold, CHN is assumed constant during the interval between two vertical dotted lines.

Fig. 5.

Time series of daytime-averaged latent heat fluxes at the Audubon, Daman, Brookings, and Bondville sites. In VDAold, CHN is assumed constant during the interval between two vertical dotted lines.

Fig. 5.

Time series of daytime-averaged latent heat fluxes at the Audubon, Daman, Brookings, and Bondville sites. In VDAold, CHN is assumed constant during the interval between two vertical dotted lines.

The RMSE and bias of H and LE estimates from VDAnew are compared with those of VDAold for each site in Table 4. For sensible heat flux, the four-site averages of RMSE (absolute bias) from VDAnew and VDAold are 39.7 (7.2) and 55.9 (15.5) W m−2, respectively. Corresponding RMSE (absolute bias) values for LE estimates are 58.4 (23.7) and 72.4 (35.2) W m−2. The reduction in misfit between observed and estimated H and LE clearly demonstrates the importance of characterizing the effect of vegetation dynamics on CHN for estimation of heat and moisture fluxes. The larger bias of H and LE estimates at the Daman site compared to those of the other sites is mainly due to the specific environmental conditions of this site. Daman has a denser vegetation cover and wetter land surface (compared to the other field sites) that degrades the performance of the VDA approach at this site [as discussed in previous studies (Abdolghafoorian and Farhadi 2016; Crow and Kustas 2005; Xu et al. 2014)]. In fact, during the energy-limited or first stage of evaporation (when soil is wet and EF is at or close to potential EF), the drying rate is mainly controlled by atmospheric conditions rather than land surface temperature (Shokri et al. 2008). Therefore, under these conditions, the coupling between EF and LST observations is weak, estimating parameters from the LST measurements becomes more uncertain, and consequently VDA performs less robustly.

RMSE−Himp (%) and RMSE−LEimp (%) [obtained from (14) where X = H and X = LE, respectively] represent the percentage of improvement in the performance of the VDAnew model in estimating fluxes of heat and moisture. These improvements are presented in Table 5 as a function of robs and . As expected, sites with higher robs and tend to demonstrate more improvements in H and LE estimates when VDAnew is applied. In fact, in these sites the assumption of constant CHN during the assimilation period (as done in the VDAold scheme) introduces more error to the estimation methodology. Maximum improvements in H and LE estimates occur at the Daman site and are 54% and 30%, respectively. The four-site averages of bias–Himp (%) and bias–LEimp (%) [obtained from (14) where X = H and X = LE and variable RMSE is replaced by bias] are 61% and 37%, respectively.

5. Conclusions

Land surface temperature (LST) observations are assimilated into a new variational data assimilation (VDA) framework to estimate the key unknown parameters of heat and moisture fluxes [i.e., evaporative fraction (EF) and neutral bulk heat transfer coefficient CHN]. Unlike the previous VDA models (e.g., Caparrini et al. 2004b; Sini et al. 2008; Bateni et al. 2013a), which consider CHN as a constant parameter during the assimilation period, our proposed VDA model takes CHN as a function of vegetation phenology [i.e., leaf area index (LAI)] to improve turbulent heat flux estimates. For this purpose, the functional relationship proposed by Farhadi et al. (2014, 2016) is used to characterize the effect of LAI on CHN.

The relationship between CHN and LAI is examined at four sites, namely, Audubon, Daman, Brookings, and Bondville, which sample a fairly wide range of hydrological and vegetative conditions. The correlations between observed CHN and LAI at these field sites are consistent with the CHN–LAI relationship proposed by Farhadi et al. (2014, 2016) and findings of other studies (Kubota and Sugita 1994; Timmermans et al. 2013). Variation of CHN over the assimilation period, especially at the Daman and Bondville sites, reveals that the assumption of constant CHN over the assimilation period can adversely affect the performance of the VDA scheme. The performance of our proposed model (VDAnew) is compared to that of the Bateni et al. (2013b) model (VDAold). Evaporative fraction estimates from the VDAnew model are in agreement with the observations in terms of magnitude and day-to-day dynamics even though no rainfall or soil moisture data are used in the model. Taking CHN as a function of vegetation phenology (i.e., LAI) improves EF estimates (in comparison with the VDAold scheme). The CHN in the VDAnew model captures the daily variation of surface roughness that influences the exchange of turbulent fluxes between the surface and the atmosphere. Thus, sensible heat flux H and consequently latent heat flux (LE) estimates from the VDAnew model are more accurate than those of VDAold. Results show that the new model outperforms the VDAold model and decreases the RMSE (and bias) in sensible and latent heat flux estimates across the four sites on average by 31% (61%) and 21% (37%), respectively. As expected, a larger improvement in estimation of EF, H, and LE is observed at sites with stronger correlation between observed ln(CHN) and LAI (i.e., robs) and during the periods in which temporal variation of vegetation phenology is significant (high ). In fact, in these sites the assumption of constant CHN during the assimilation period (as done in the VDAold scheme) introduces more error to the estimation methodology.

Acknowledgments

Funding for this research project is provided by the George Washington University (CEE 176104). The authors also thank the National Science Foundation of China (41671335) for publication support. The authors would like to thank all the scientists, engineers, and students who participated in the HiWATER field campaigns. The AmeriFlux data used in this paper can be downloaded at http://ameriflux.lbl.gov/. Tilden Meyers is the primary investigator at the Audubon Research Ranch, Bondville, and Brookings. The authors appreciate access to these data.

REFERENCES

REFERENCES
Abdolghafoorian
,
A.
, and
L.
Farhadi
,
2016
:
Uncertainty quantification in land surface hydrologic modeling: Toward an integrated variational data assimilation framework
.
IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.
,
9
,
2628
2637
, doi:.
Alfieri
,
J. G.
,
X.
Xiao
,
D.
Niyogi
,
R. A.
Pielke
,
F.
Chen
, and
M. A.
LeMone
,
2009
:
Satellite-based modeling of transpiration from the grasslands in the Southern Great Plains, USA
.
Global Planet. Change
,
67
,
78
86
, doi:.
Allen
,
R.
,
A.
Irmak
,
R.
Trezza
,
J. M. H.
Hendrickx
,
W.
Bastiaanssen
, and
J.
Kjaersgaard
,
2011
:
Satellite-based ET estimation in agriculture using SEBAL and METRIC
.
Hydrol. Processes
,
25
,
4011
4027
, doi:.
Anderson
,
M. C.
,
J. M.
Norman
,
J. R.
Diak
,
W. P.
Kustas
, and
J. R.
Mecikalski
,
1997
:
A two-source time-integrated model for estimating surface fluxes using thermal infrared remote sensing
.
Remote Sens. Environ.
,
60
,
195
216
, doi:.
Andreas
,
E. L
, and
B.
Murphy
,
1986
:
Bulk transfer coefficients for heat and momentum over leads and polynyas
.
J. Phys. Oceanogr.
,
16
,
1875
1883
, doi:.
Andrews
,
T.
,
P. M.
Forster
, and
J. M.
Gregory
,
2009
:
A surface energy perspective on climate change
.
J. Climate
,
22
,
2557
2570
, doi:.
Ann-Sofi
,
S.
,
S.
Ulf
,
H.
Erik
, and
J.
Cecilia
,
2007
:
Critical re-evaluation of the bulk transfer coefficient for sensible heat over the ocean during unstable and neutral conditions
.
Quart. J. Roy. Meteor. Soc.
,
133
,
227
250
, doi:.
Baldocchi
,
D.
, and Coauthors
,
2001
:
FLUXNET: A new tool to study the temporal and spatial variability of ecosystem-scale carbon dioxide, water vapor, and energy flux densities
.
Bull. Amer. Meteor. Soc.
,
82
,
2415
2434
, doi:.
Bastiaanssen
,
W. G. M.
,
M.
Menenti
,
R. A.
Feddes
, and
A. A. M.
Holtslag
,
1998
:
A remote sensing Surface Energy Balance Algorithm for Land (SEBAL). 1. Formulation
J. Hydrol.
,
212–213
,
198
212
, doi:.
Bateni
,
S. M.
, and
D.
Entekhabi
,
2012
:
Surface heat flux estimation with the ensemble Kalman smoother: Joint estimation of state and parameters
.
Water Resour. Res.
,
48
,
W08521
, doi:.
Bateni
,
S. M.
,
D.
Entekhabi
, and
F.
Castelli
,
2013a
:
Mapping evaporation and estimation of surface control of evaporation using remotely sensed land surface temperature from a constellation of satellites
.
Water Resour. Res.
,
49
,
950
968
, doi:.
Bateni
,
S. M.
,
D.
Entekhabi
, and
D.-S.
Jeng
,
2013b
:
Variational assimilation of land surface temperature and the estimation of surface energy balance components
.
J. Hydrol.
,
481
,
143
156
, doi:.
Bateni
,
S. M.
,
D.
Entekhabi
,
S. A.
Margulis
,
F.
Castelli
, and
L.
Kergoat
,
2014
:
Coupled estimation of surface heat fluxes and vegetation dynamics from remotely sensed land surface temperature and fraction of photosynthetically active radiation
.
Water Resour. Res.
,
50
,
8420
8440
, doi:.
Behera
,
S. K.
,
P.
Srivastava
,
U. V.
Pathre
, and
R.
Tuli
,
2010
:
An indirect method of estimating leaf area index in Jatropha curcas L. using LAI-2000 Plant Canopy Analyzer
.
Agric. For. Meteor.
,
150
,
307
311
, doi:.
Betts
,
A. K.
, and
J. H.
Ball
,
1998
:
FIFE surface climate and site-average dataset 1987–89
.
J. Atmos. Sci.
,
55
,
1091
1108
, doi:.
Boulet
,
G.
,
A.
Olioso
,
E.
Ceschia
,
O.
Marloie
,
B.
Coudert
,
V.
Rivalland
,
J.
Chirouze
, and
G.
Chehbouni
,
2012
:
An empirical expression to relate aerodynamic and surface temperatures for use within single-source energy balance models
.
Agric. For. Meteor.
,
161
,
148
155
, doi:.
Brutsaert
,
W.
,
1979
:
Heat and mass transfer to and from surfaces with dense vegetation or similar permeable roughness
.
Bound.-Layer Meteor.
,
16
,
365
388
, doi:.
Campbell
,
G. S.
,
1985
: Soil Physics with BASIC: Transport Models for Soil-Plant Systems. Elsevier, 149 pp.
Campo
,
L.
,
F.
Castelli
,
D.
Entekhabi
, and
F.
Caparrini
,
2009
:
Land–atmosphere interactions in an high resolution atmospheric simulation coupled with a surface data assimilation scheme
.
Nat. Hazards Earth Syst. Sci.
,
9
,
1613
1624
, doi:.
Caparrini
,
F.
,
F.
Castelli
, and
D.
Entekhabi
,
2003
:
Mapping of land–atmosphere heat fluxes and surface parameters with remote sensing data
.
Bound.-Layer Meteor.
,
107
,
605
633
, doi:.
Caparrini
,
F.
,
F.
Castelli
, and
D.
Entekhabi
,
2004a
:
Variational estimation of soil and vegetation turbulent transfer and heat flux parameters from sequences of multisensor imagery
.
Water Resour. Res.
,
40
,
W12515
, doi:.
Caparrini
,
F.
,
F.
Castelli
, and
D.
Entekhabi
,
2004b
:
Estimation of surface turbulent fluxes through assimilation of radiometric surface temperature sequences
.
J. Hydrometeor.
,
5
,
145
159
, doi:.
Colaizzi
,
P. D.
,
S. R.
Evett
,
T. A.
Howell
, and
J. A.
Tolk
,
2006
:
Comparison of five models to scale daily evapotranspiration from one-time-of-day measurements
.
Trans. ASABE
,
49
,
1409
1417
, doi:.
Cragoa
,
R.
, and
W.
Brutsaert
,
1996
:
Daytime evaporation and the self-preservation of the evaporative fraction and the Bowen ratio
.
J. Hydrol.
,
178
,
241
255
, doi:.
Crow
,
W. T.
, and
W. P.
Kustas
,
2005
:
Utility of assimilating surface radiometric temperature observations for evaporative fraction and heat transfer coefficient retrieval
.
Bound.-Layer Meteor.
,
115
,
105
130
, doi:.
Fang
,
H.
,
S.
Wei
, and
S.
Liang
,
2012
:
Validation of MODIS and CYCLOPES LAI products using global field measurement data
.
Remote Sens. Environ.
,
119
,
43
54
, doi:.
Farhadi
,
L.
,
D.
Entekhabi
,
G.
Salvucci
, and
J.
Sun
,
2014
:
Estimation of land surface water and energy balance parameters using conditional sampling of surface states
.
Water Resour. Res.
,
50
,
1805
1822
, doi:.
Farhadi
,
L.
,
D.
Entekhabi
, and
G.
Salvucci
,
2016
:
Mapping land water and energy balance relations through conditional sampling of remote sensing estimates of atmospheric forcing and surface states
.
Water Resour. Res.
,
52
,
2737
2752
, doi:.
Fensholt
,
R.
,
I.
Sandholt
, and
M. S.
Rasmussen
,
2004
:
Evaluation of MODIS LAI, fAPAR and the relation between fAPAR and NDVI in a semi-arid environment using in situ measurements
.
Remote Sens. Environ.
,
91
,
490
507
, doi:.
Gentine
,
P.
,
D.
Entekhabi
,
A.
Chehbouni
,
G.
Boulet
, and
B.
Duchemin
,
2007
:
Analysis of evaporative fraction diurnal behaviour
.
Agric. For. Meteor.
,
143
,
13
29
, doi:.
Hasager
,
C. B.
, and
N. O.
Jensen
,
1999
:
Surface-flux aggregation in heterogeneous terrain
.
Quart. J. Roy. Meteor. Soc.
,
125
,
2075
2102
, doi:.
Hirota
,
T.
,
J. W.
Pomeroy
,
R. J.
Granger
, and
C. P.
Maule
,
2002
:
An extension of the force-restore method to estimating soil temperature at depth and evaluation for frozen soils under snow
.
J. Geophys. Res.
,
107
,
4767
, doi:.
Hopmans
,
J. W.
,
J.
Šimunek
, and
K. L.
Bristow
,
2002
:
Indirect estimation of soil thermal properties and water flux using heat pulse probe measurements: Geometry and dispersion effects
.
Water Resour. Res.
,
38
,
1006
, doi:.
Jensen
,
N. O.
, and
P.
Hummelshøj
,
1995
:
Derivation of canopy resistance for water vapour fluxes over a spruce forest, using a new technique for the viscous sublayer resistance
.
Agric. For. Meteor.
,
73
,
339
352
, doi:.
Kalma
,
J. D.
,
T. R.
McVicar
, and
M. F.
McCabe
,
2008
:
Estimating land surface evaporation: A review of methods using remotely sensed surface temperature data
.
Surv. Geophys.
,
29
,
421
469
, doi:.
Kubota
,
A.
, and
M.
Sugita
,
1994
:
Radiometrically determined skin temperature and scalar roughness to estimate surface heat flux. Part I: Parameterization of radiometric scalar roughness
.
Bound.-Layer Meteor.
,
69
,
397
416
, doi:.
Lewellen
,
D. C.
, and
W. S.
Lewellen
,
2002
:
Entrainment and decoupling relations for cloudy boundary layers
.
J. Atmos. Sci.
,
59
,
2966
2986
, doi:.
Lewis
,
P.
,
J.
Gómez-Dans
,
T.
Kaminski
,
J.
Settle
,
T.
Quaife
,
N.
Gobron
,
J.
Styles
, and
M.
Berger
,
2012
:
An Earth Observation Land Data Assimilation System (EO-LDAS)
.
Remote Sens. Environ.
,
120
,
219
235
, doi:.
Lhomme
,
J. P.
,
A.
Chehbouni
, and
B.
Monteny
,
2000
:
Sensible heat flux–radiometric surface temperature relationship over sparse vegetation: Parameterizing B-1
.
Bound.-Layer Meteor.
,
97
,
431
457
, doi:.
Li
,
X.
, and Coauthors
,
2013
:
Heihe watershed allied telemetry experimental research (HiWater) scientific objectives and experimental design
.
Bull. Amer. Meteor. Soc.
,
94
,
1145
1160
, doi:.
Mahrt
,
L.
, and
D.
Vickers
,
2004
:
Bulk formulation of the surface heat flux
.
Bound.-Layer Meteor.
,
110
,
357
379
, doi:.
Massman
,
W. J.
,
1999
:
A model study of kBH−1 for vegetated surfaces using “localized near-field” Lagrangian theory
.
J. Hydrol.
,
223
,
27
43
, doi:.
Matsushima
,
D.
,
2005
:
Relations between aerodynamic parameters of heat transfer and thermal-infrared thermometry in the bulk surface formulation
.
J. Meteor. Soc. Japan
,
83
,
373
389
, doi:.
McNaughton
,
K. G.
, and
B. J. J. M.
Hurk
,
1995
:
A “Lagrangian” revision of the resistors in the two-layer model for calculating the energy budget of a plant canopy
.
Bound.-Layer Meteor.
,
74
,
261
288
, doi:.
Perez
,
P. J.
,
F.
Castellvi
,
M.
Ibañez
, and
J. I.
Rosell
,
1999
:
Assessment of reliability of Bowen ratio method for partitioning fluxes
.
Agric. For. Meteor.
,
97
,
141
150
, doi:.
Peters-Lidard
,
C. D.
,
S. V.
Kumar
,
D. M.
Mocko
, and
Y.
Tian
,
2011
:
Estimating evapotranspiration with land data assimilation systems
.
Hydrol. Processes
,
25
,
3979
3992
, doi:.
Polonio
,
D.
, and
M. R.
Soler
,
2000
:
Surface fluxes estimation over agricultural areas. Comparison of methods and the effects of land surface inhomogeneity
.
Theor. Appl. Climatol.
,
67
,
65
79
, doi:.
Qualls
,
R.
, and
W.
Brutsaert
,
1996
:
Effect of vegetation density on the parameterization of scalar roughness to estimate spatially distributed sensible heat fluxes
.
Water Resour. Res.
,
32
,
645
652
, doi:.
Reichle
,
R. H.
,
D.
Entekhabi
, and
D. B.
McLaughlin
,
2001a
:
Downscaling of radio brightness measurements for soil moisture estimation: A four-dimensional variational data assimilation approach
.
Water Resour. Res.
,
37
,
2353
2364
, doi:.
Reichle
,
R. H.
,
D. B.
McLaughlin
, and
D.
Entekhabi
,
2001b
:
Variational data assimilation of microwave radio brightness observations for land surface hydrology applications
.
IEEE Trans. Geosci. Remote Sens.
,
39
,
1708
1718
, doi:.
Richardson
,
A. D.
, and Coauthors
,
2006
:
A multi-site analysis of random error in tower-based measurements of carbon and energy fluxes
.
Agric. For. Meteor.
,
136
,
1
18
, doi:.
Shokri
,
N.
,
P.
Lehmann
,
P.
Vontobel
, and
D.
Or
,
2008
:
Drying front and water content dynamics during evaporation from sand delineated by neutron radiography
.
Water Resour. Res.
,
44
,
W06418
, doi:.
Sini
,
F.
,
G.
Boni
,
F.
Caparrini
, and
D.
Entekhabi
,
2008
:
Estimation of large-scale evaporation fields based on assimilation of remotely sensed land temperature
.
Water Resour. Res.
,
44
,
W06410
, doi:.
Stephens
,
G. L.
, and Coauthors
,
2012
:
An update on Earth’s energy balance in light of the latest global observations
.
Nat. Geosci.
,
5
,
691
696
, doi:.
Su
,
Z.
,
2002
:
The Surface Energy Balance System (SEBS) for estimation of turbulent heat fluxes
.
Hydrol. Earth Syst. Sci.
,
6
,
85
100
, doi:.
Sugita
,
M.
, and
W.
Brutsaert
,
1990
:
Regional surface fluxes from remotely sensed skin temperature and lower boundary layer measurements
.
Water Resour. Res.
,
26
,
2937
2944
, doi:.
Timmermans
,
J.
,
Z.
Su
,
C.
van der Tol
,
A.
Verhoef
, and
W.
Verhoef
,
2013
:
Quantifying the uncertainty in estimates of surface–atmosphere fluxes through joint evaluation of the SEBS and SCOPE models
.
Hydrol. Earth Syst. Sci.
,
17
,
1561
1573
, doi:.
Verhoef
,
A.
,
H. R.
De Bruin
, and
B. J. J. M.
Van Den Hurk
,
1997
:
Some practical notes on the parameter kB−1 for sparse vegetation
.
J. Appl. Meteor.
,
36
,
560
572
, doi:.
Wu
,
X.
,
J.
Zhou
,
H.
Wang
,
Y.
Li
, and
B.
Zhong
,
2015
:
Evaluation of irrigation water use efficiency using remote sensing in the middle reach of the Heihe River, in the semi-arid northwestern China
.
Hydrol. Processes
,
29
,
2243
2257
, doi:.
Xu
,
T. R.
,
S. M.
Bateni
,
S.
Liang
,
D.
Entekhabi
, and
K. B.
Mao
,
2014
:
Estimation of surface turbulent heat fluxes via variational assimilation of sequences of land surface temperatures from Geostationary Operational Environmental Satellites
.
J. Geophys. Res. Atmos.
,
119
,
10 780
10 798
, doi:.
Yilmaz
,
M. T.
,
M. C.
Anderson
,
B.
Zaitchik
,
C. R.
Hain
,
W. T.
Crow
,
M.
Ozdogan
,
J. A.
Chun
, and
J.
Evans
,
2014
:
Comparison of prognostic and diagnostic surface flux modeling approaches over the Nile River basin
.
Water Resour. Res.
,
50
,
386
408
, doi:.