Feasibility of Estimating Turbulent Heat Fluxes via Variational Assimilation of Reference-Level Air Temperature and Speciﬁc Humidity Observations

: This study investigated the feasibility of partitioning the available energy between sensible ( H ) and latent ( LE ) heat ﬂuxes via variational assimilation of reference-level air temperature and speciﬁc humidity. For this purpose, sequences of reference-level air temperature and speciﬁc humidity were assimilated into an atmospheric boundary layer model (ABL) within a variational data assimilation (VDA) framework to estimate H and LE . The VDA approach was tested at six sites (namely, Arou, Audubon, Bondville, Brookings, Desert, and Willow Creek) with contrasting climatic and vegetative conditions. The unknowns of the VDA system were the neutral bulk heat transfer coe ﬃ cient ( C HN ) and evaporative fraction ( EF ). EF estimates were found to agree well with observations in terms of magnitude and day-to-day ﬂuctuations in wet / densely vegetated sites but degraded in dry / sparsely vegetated sites. Similarly, in wet / densely vegetated sites, the variations in the C HN estimates were found to be consistent with those of the leaf area index (LAI) while this consistency deteriorated in dry / sparely vegetated sites. The root mean square errors (RMSEs) of daily H and LE estimates at the Arou site (wet) were 25.43 (Wm − 2 ) and 55.81 (Wm − 2 ), which are respectively 57.6% and 45.4% smaller than those of 60.00 (Wm − 2 ) and 102.21 (Wm − 2 ) at the Desert site (dry). Overall, the results show that the VDA system performs well at wet / densely vegetated sites (e.g., Arou and Willow Creek), but its performance degrades at dry / slightly vegetated sites (e.g., Desert and Audubon). These outcomes show that the sequences of reference-level air temperature and speciﬁc humidity have more information on the partitioning of available energy between the sensible and latent heat ﬂuxes in wet / densely vegetated sites than dry / slightly vegetated sites.


Sensible and Latent Heat Fluxes
The sensible heat flux is given by: where ρ is the air density, c p is the specific heat capacity of air, U is the wind speed at the reference level (z re f = 2 m), T is the land surface temperature, T a is the air temperature at z re f , and f is the atmospheric stability correction function, which is a function of the Richardson number (Ri). Ri is a measure of atmospheric stability. The stability correction function (f ) proposed by Caparrini et al. [45] performed well in previous studies [49][50][51]53,57,[81][82][83], and thus was used in this study. It is given by f (Ri) = 1 + 2(1 -e 10Ri ). C HN is the first unknown of the VDA approach, which depends on the characteristics of the landscape, and is assumed to be constant during each month [49,50,53,57]. C HN is mainly a function of LAI and to a lesser extent wind speed, friction velocity, solar elevation, and the structure and shape of vegetation (i.e., crown density and vertical distribution of foliage elements) [50,[84][85][86][87][88].
EF is the second unknown of the VDA approach that represents partitioning between the turbulent heat fluxes, and is defined as the ratio of latent heat flux to the sum of turbulent heat fluxes: Crago [89], Caparrini et al. [46,47], and Bateni et al. [49][50][51] showed that EF is almost constant for near-peak radiation hours [09:00-16:00 LT] on days without precipitation. To estimate LE, Equation (2) can be re-written as follows: It is worth mentioning that the assumption of daily constant EF and monthly constant C HN are widely used in VDA studies [42,44,46,47,49,50,53,55,56,[81][82][83] to estimate H and LE.

Atmospheric Boundary Layer (ABL) Model
A mixed-layer model, which performs well compared to the complicated large eddy simulation models, is used to simulate ABL processes [81][82][83]90,91]. The profiles of potential temperature (θ) and specific humidity (q) are assumed to be constant with height throughout the mixed layer. As shown in Figure 1, the thin layer between the ground and the mixed layer is called the surface layer (SL). The surface layer is assumed to be 10% of the mixed-layer top, h (i.e., z SL = 0.1h), and is convectively unstable during the growth phase of the boundary layer [73,79,81,83].

Energy and Moisture Budget Equations
The potential temperature and specific humidity in the mixed layer are given by: where H top and LE top are the entrainment sensible and latent heat fluxes from above the mixed layer, respectively; ε m is the mixed-layer bulk emissivity; t is the time; and L v is the latent heat of vaporization. R ad and R gu are the downward longwave radiation from the free atmosphere and upward longwave radiation from the ground into the mixed layer, respectively. R Ad and R Au are the downward and upward longwave radiative fluxes from within the mixed layer, respectively. The potential temperature at the level of 1000 mb (i.e., θ a ) is obtained from the reference-level air temperature (T a ) via θ a = T a (P 0 /P s ) R d /c p (where R d is the gas constant of dry air, P s is the surface pressure, and P 0 is 1000 mb) [92]. The variations of the potential temperature and specific humidity profiles within the surface layer can be explained by the Monin-Obukhov similarity theory (MOST) [78,93]. Hence, θ SL and q SL are found from θ a and q a using the MOST (see Appendix B). As q and θ are uniform throughout the mixed layer, the initial conditions for Equations (4) and (5) (i.e., θ(t = t o ) and q(t = t o )) are set to be equal to θ SL (t = t o ) and q SL (t = t o ), respectively.

Radiative Fluxes
The downwelling longwave radiation from the overlying atmosphere (R ad ) and the upwelling longwave radiative flux from the land surface beneath (R gu ) are given by: where σ is the Stefan-Boltzman constant, ε ad is the effective emissivity above the mixed layer [73,81,83,[94][95][96], T h+ is the air temperature exactly above the mixed layer [95,96], and ε s is the surface emissivity. The downward and upward longwave radiative fluxes from within the mixed layer (R Ad and R Au ) are presented by: where ε d and ε u are the mixed-layer downward and upward emissivity, respectively [95,96].

Mixed-Layer Height
The daytime ABL height growth is calculated by [81,83,[95][96][97][98][99]: where the various terms in Equation (10) are given by: where u * is the friction velocity; u SL is the wind speed at the top of the surface layer, which is obtained from the reference-level wind speed (U) via the MOST (the readers are referred to Appendix B); G * is the production of mechanical turbulent energy; g is gravitational acceleration; ϕ is the mechanical turbulence dissipation parameter, which is set to 0.01 [73,81,83,96,99]; H v is the virtual heat flux at the surface; E is the rate of evaporation from ground; and δ θ is the potential temperature inversion strength at the top of the mixed layer.

Inversion Strengths of θ and q
There are instantaneous jumps (inversion strengths) in temperature and humidity (δ θ and δ q ) at the top of the mixed layer ( Figure 1). As shown in Equation (15), the inversion strength of θ increases Remote Sens. 2020, 12, 1065 5 of 29 as the boundary layer grows, but it decreases as the boundary layer warms up. Similarly, the inversion strength of q increases as the boundary layer develops, and it increases (becomes more negative) as the boundary layer becomes moister (Equation (16)). The magnitudes of these jumps can be calculated by the following prognostic equations [81,83,100]: where γ θ and γ q are the lapse rates in the potential temperature and specific humidity above the mixed layer.

Entrainment Fluxes
As the boundary layer grows, warm dry air from the free atmosphere enters the mixed layer, and results in the presence of sensible and latent heat fluxes entrainment (H top and LE top ). H top heats up the ABL and LE top reduces its humidity [101]: A typical value of 0.2 is used for A [81,83,90,91].
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 28 As the boundary layer grows, warm dry air from the free atmosphere enters the mixed layer, and results in the presence of sensible and latent heat fluxes entrainment ( and ). heats up the ABL and reduces its humidity [101]: = .

Variational Data Assimilation (VDA) Approach
The unknowns of the VDA system (i.e., CHN and EF) are obtained by minimizing the objective function J. Two different integral time scales are used in the objective function J. The first one covers the entire assimilation period in which CHN is assumed to be constant (N = 30 days). The second time scale constitutes the assimilation window in which EF is presumed to be constant [t0,t1] = [09:00-16:00 LT]. The objective function J is expressed as: , , , , , =

Variational Data Assimilation (VDA) Approach
The unknowns of the VDA system (i.e., C HN and EF) are obtained by minimizing the objective function J. Two different integral time scales are used in the objective function J. The first one covers the entire assimilation period in which C HN is assumed to be constant (N = 30 days). The second time scale constitutes the assimilation window in which EF is presumed to be constant [t 0 ,t 1 ] = [09:00-16:00 LT]. The objective function J is expressed as: The first term in the right-hand side of Equation (19) is the quadratic of difference between the top of the surface layer potential temperature (θ SL ) and the ABL potential temperature (θ) estimates from Equation (4). Similarly, the second term is the quadratic of the misfit between the top of the surface layer specific humidity (q SL ) and the ABL specific humidity (q) estimates from Equation (5). Using the MOST (see Appendix B), θ SL and q SL are retrieved from the reference-level air temperature (T a ) and specific humidity (q a ), respectively. To make C HN strictly positive, it is transformed to R via C HN = e R . The third and fourth terms are the quadratic errors of the unknown parameters (i.e., R and EF) with respect to their prior values (R and EF ). The last two terms are the physical constraints adjoined to the model through the Lagrange multipliers λ 1 and λ 2 . R −1 θ and R −1 q are the inverse error covariance matrices of θ and q, respectively. B −1 R and B −1 EF are the inverse background error covariance matrices of R and EF, respectively. Following Tajfar et al. [81,83], the diagonal elements of R −1 θ , R −1 q , B −1 R , and B −1 EF are set to 10 −1 K −2 , 10 5 (kg/kg) −2 , 10 8 , and 10 9 , respectively. To minimize the objective function, its first variation (δJ) with respect to the independent variables θ, q, λ 1 , λ 2 , R, and EF is set to zero. This leads to a set of the so-called Euler-Lagrange equations that should be solved simultaneously. The Euler-Lagrange equations are presented in Appendix C.

Study Sites
The ability of the VDA approach to partition the available energy between the turbulent heat fluxes was tested at six sites (namely, Arou, Audubon, Bondville, Brookings, Desert, and Willow Creek) with contrasting climatic and vegetative conditions. Audubon, Bondville, Brookings, and Willow Creek are located in the United States ( Figure 2). Arou and Desert are situated in the middle reach of Heihe River basin (HRB) in the Gobi desert in Northern China (http://card.westgis.ac.cn/) ( Figure 3). Audubon (in Arizona) is a grassland water-limited monsoonal site with a temperate arid climate (https://ameriflux.lbl.gov/sites/siteinfo/US-Aud). Bondville (in Illinois) is a cropland that has a humid continental climate with hot summers (https://ameriflux.lbl.gov/sites/siteinfo/US-Bo1). Brookings (in South Dakota) is a grassland with a temperate continental climate and no dry season (https://ameriflux. lbl.gov/sites/siteinfo/US-Bkg). Willow Creek (in Wisconsin) is a deciduous broad-leaf forest with dense vegetation cover and significant precipitation (https://ameriflux.lbl.gov/sites/siteinfo/US-WCr). Arou (in Qinghai) is a dense grassland with high soil moisture (http://data.tpdc.ac.cn/en/data/9074b105-549c-4f15-b8a5-87602a45134d/). Desert (in Inner Mongolia) is barren soil and has low precipitation and sparse canopy cover (http://data.tpdc.ac.cn/en/data/77c7bfd5-38d3-41e2-93ce-cb4d8c0d475a/). Leaf area index (LAI) ranges from 0.00 in Desert (barren land) to 5.67 in Willow Creek (dense forest). Soil moisture (SM) varies from 0.03 in Desert (dry) to 0.36 in Arou (wet). These sites were chosen to sample different land cover types (forest, grassland, cropland, and barren land), and cover a wide range of vegetation density and soil moisture conditions (Table 1).

Results and Discussions
As mentioned previously, CHN and EF are the two key unknowns of the VDA system. Table 3 shows the estimated CHN values from the VDA approach and corresponding LAI values at the six study sites. CHN is mainly affected by vegetation phenology [44,46,47,50,52,55]. As shown in Table 3, the variations in CHN estimates are consistent with those of LAI in the wet and/or densely vegetated sites. For instance, at the Willow Creek site, CHN and LAI decrease slightly during the modeling periods. At the Arou and Brookings sites, CHN reaches its highest value in the second monthly period

Results and Discussions
As mentioned previously, CHN and EF are the two key unknowns of the VDA system. Table 3 shows the estimated CHN values from the VDA approach and corresponding LAI values at the six study sites. CHN is mainly affected by vegetation phenology [44,46,47,50,52,55]. As shown in Table 3, the variations in CHN estimates are consistent with those of LAI in the wet and/or densely vegetated sites. For instance, at the Willow Creek site, CHN and LAI decrease slightly during the modeling periods. At the Arou and Brookings sites, CHN reaches its highest value in the second monthly period  LAI and SM are the key factors affecting the partitioning of available energy between the sensible and latent heat fluxes [50,58,102]. Hence, testing the VDA approach at the abovementioned six sites provides insights into how much information is contained in the sequences of reference-level air temperature and specific humidity for partitioning the available energy between H and LE in various climatic and vegetative conditions. Soil moisture and half-hourly meteorological data (air temperature and specific humidity, wind speed, and atmospheric pressure) at the Audubon, Bondville, Brookings, and Willow Creek sites were retrieved from the AmeriFlux archive (http://www.ameriflux.lbl.gov). The sensible and latent heat fluxes at these sites were measured by eddy covariance (EC) stations, and are available on the Ameriflux database. LAI data were obtained from the Global LAnd Surface Satellites (GLASS) product [103,104]. This product is accessible on the University of Maryland global land cover facility archive (http://www.un-spider.org/links-and-resources/data-sources/global-land-cover-facilityuniversity-maryland-nasagofc-gold).
The Multiscale Observation Experiment on Evapotranspiration over the Heihe Watershed Allied Telemetry Experimental Research (HiWATER-MUSOEXE) provides half-hourly measurements of wind speed, air temperature, soil moisture, relative humidity, and pressure at the Arou and Desert sites (http://card.westgis.ac.cn/hiwater) [9,11,105,106]. The sensible and latent heat fluxes were recorded at these sites by an EC station every 30 min (http://card.westgis.ac.cn/data/e9e38ff4-5a10-4977-9de7-61e9c25bd333). LAI data at the Arou and Desert sites were also obtained from the GLASS product. Vegetation height was obtained by in situ measurements [105].
The initial condition for h is required to solve Equation (10). Similarly, the initial conditions for δ θ and δ q , and the magnitudes of γ θ and γ q are needed to integrate Equations (15) and (16) forward in time. Following Tajfar et al. [81,83], the initial conditions for h, δ θ , and δ q , and the magnitudes of γ θ and γ q were varied from 100 to 500 m, 2 to 6 K, −4.8 × 10 −3 to −0.5 × 10 −3 kg kg −1 , 2 to 8 K km −1 , and −7 × 10 −3 to −0.5 × 10 −3 kg kg −1 km −1 with the increment of 100 m, 0.4 K, 0.4 × 10 −3 kg kg −1 , 0.5 K km −1 , and 0.5 × 10 −3 kg kg −1 km −1 , respectively. For each study site, the magnitudes of h (t = t o ), δ θ (t = t o ), and δ q (t = t o ), γ θ , and γ q that lead to a minimum cost function (J) are reported in Table 2. Table 2. The magnitudes of the initial conditions for h, δ θ , and δ q , as well as γ θ and γ q for the study sites. Site The VDA was applied to the Desert ( and Arou (DOYs 170-259, 2015). These modeling periods were chosen because they have a minimum data gap in H and LE measurements. Linear interpolation was used to fill the data gap [107,108].

Results and Discussions
As mentioned previously, C HN and EF are the two key unknowns of the VDA system. Table 3 shows the estimated C HN values from the VDA approach and corresponding LAI values at the six study sites. C HN is mainly affected by vegetation phenology [44,46,47,50,52,55]. As shown in Table 3, the variations in C HN estimates are consistent with those of LAI in the wet and/or densely vegetated sites. For instance, at the Willow Creek site, C HN and LAI decrease slightly during the modeling periods. At the Arou and Brookings sites, C HN reaches its highest value in the second monthly period and decreases in the third month. A similar trend is observed in LAI. At the Bondville site, the LAI and C HN estimates remain almost constant during the modeling periods (DOYs 182-271).
The consistency between the C HN retrievals and LAI weakens in dry/sparsely vegetated sites. At Audubon, LAI values increase continuously during the modeling periods. However, C HN estimates decrease slightly in the second monthly period. At the Desert site, LAI is invariant, while C HN shows slight variations in different assimilation periods. Among the study sites, Desert has the lowest C HN estimates because it has no vegetation (LAI = 0), while Willow Creek has the highest C HN values due to its dense vegetation cover (LAI = 5.67). In general, the study sites with higher LAI values (e.g., Arou and Willow Creek) have higher C HN estimates compared to the sites with lower LAI (e.g., Desert and Audubon).
Overall, the results show that the C HN estimates agree well with the vegetation phenology at wet and/or densely vegetated sites, although no information on vegetation density is used in the VDA approach. These findings indicate that the VDA approach can extract the implicit information contained in the air temperature and specific humidity measurements to estimate C HN at wet and/or densely vegetated sites, but its performance degrades in dry and/or slightly vegetated sites.  Figure 4 shows the time series of EF estimates at the six study sites. EF values from the measured H and LE are also shown by open circles in Figure 4. The EF retrievals agree well with the observations in terms of both the magnitude and day-to-day fluctuations in wet and/or densely vegetated sites (e.g., Arou and Willow Creek). At these sites, the EF values are consistent with the wetting and dry down events, although no information on rainfall is used in the VDA approach. For example, sharp jumps are seen in the EF estimates after rainfall events at the Willow Creek site (e.g., DOY 181, 201, 206, 209, 221, 231, and 255), followed by reductions in the EF estimates during dry down events (e.g., DOY 171, 234, 245, 249, and 252). This shows that the sequences of air temperature and humidity have a significant amount of information on the partitioning of available energy between the turbulent heat fluxes in wet and/or densely vegetated sites.
The drying rate of land surface is controlled by both land surface properties and atmospheric factors [109]. At dry/sparsely vegetated sites (e.g., Desert and Audubon), evapotranspiration is at its second stage (water limited), and thus is mainly controlled by the land state variable (i.e., LST) rather than the atmospheric state variables (i.e., air temperature and specific humidity). Consequently, the coupling between EF and the atmospheric state variables is weak, and the retrieval of EF from air temperature and specific humidity becomes more uncertain. The Desert site is dry on DOYs 170-245. In this period, the VDA approach performs poorly and the EF estimate show unreasonable spikes. Towards the end of the modeling period (i.e., DOYs 246-259), the Desert site becomes wet due to rainfall events, and the EF estimates can capture the variations in the observations. Similarly, the Audubon site is mostly dry on DOYs 170-210, and thus the VDA cannot robustly estimate EF. However, it becomes moister on DOYs 210-259, and the EF retrievals follow the observations more closely. These results indicate that the sequences of air temperature and humidity in dry/sparsely vegetated sites do not have enough information to constrain EF.
The high soil moisture can result in negative sensible heat flux, which makes EF observations larger than one [44,110]. This occurs on DOYs 184, 201, 206, 209, and 214 in the Willow Creek; DOY 201 in the Brookings; and DOYs 206 and 242 in the Bondville. Nonetheless, in order to prevent numerical instabilities, EF estimates are set to be less than 0.99 in the VDA approach, causing the estimations to deviate from observations when they are above unity. Table 4 shows the MAE and RMSE of the EF estimates at the six study sites. As anticipated, the MAE and RMSE of the EF estimates are lower in wet and/or densely vegetated sites compared to dry and/or slightly vegetated sites. The EF estimates have the highest MAE and RMSE of 0.152 and 0.198 at Desert (dry and barren land). The lowest MAE and RMSE values of 0.039 and 0.053 are observed at Arou (wet and dense vegetation cover).
To evaluate the performance of the VDA approach in various vegetative and climatic conditions, the scatterplot of half-hourly H and LE estimates was compared with the observations in Figures 5 and 6. As shown, the VDA approach performs better at wet/densely vegetated sites. The most and least accurate turbulent heat fluxes estimates are obtained at Arou and Desert, respectively. H and LE estimates mostly fall around the 45-degree line at Arou and Willow Creek. While at Desert, H is overestimated for H values of higher than~200 Wm −2 and LE retrievals are far from the 45-degree line These outcomes show that the information content of the atmospheric state variables (e.g., air temperature and specific humidity) for estimating the turbulent heat fluxes significantly reduces in dry/sparsely vegetated sites (e.g., Audubon and Desert). H and LE estimates at Bondville and Brookings are comparable, and agree fairly well with the observations. Overall, the results indicate that the sequences of air temperature and specific humidity have a significant amount of information for partitioning the available energy between turbulent heat fluxes at sites with high SM and/or LAI values, but their information content significantly reduces at sites with low SM and/or LAI. temperature and specific humidity becomes more uncertain. The Desert site is dry on DOYs 170-245. In this period, the VDA approach performs poorly and the EF estimate show unreasonable spikes. Towards the end of the modeling period (i.e., DOYs 246-259), the Desert site becomes wet due to rainfall events, and the EF estimates can capture the variations in the observations. Similarly, the Audubon site is mostly dry on DOYs 170-210, and thus the VDA cannot robustly estimate EF. However, it becomes moister on DOYs 210-259, and the EF retrievals follow the observations more closely. These results indicate that the sequences of air temperature and humidity in dry/sparsely vegetated sites do not have enough information to constrain EF. The high soil moisture can result in negative sensible heat flux, which makes EF observations larger than one [44,110]. This occurs on DOYs 184, 201, 206, 209, and 214 in the Willow Creek; DOY 201 in the Brookings; and DOYs 206 and 242 in the Bondville. Nonetheless, in order to prevent numerical instabilities, EF estimates are set to be less than 0.99 in the VDA approach, causing the estimations to deviate from observations when they are above unity. Table 4 shows the MAE and RMSE of the EF estimates at the six study sites. As anticipated, the MAE and RMSE of the EF estimates are lower in wet and/or densely vegetated sites compared to dry and/or slightly vegetated sites. The     Creek. While at Desert, H is overestimated for H values of higher than ~200 Wm −2 and LE retrievals are far from the 45-degree line These outcomes show that the information content of the atmospheric state variables (e.g., air temperature and specific humidity) for estimating the turbulent heat fluxes significantly reduces in dry/sparsely vegetated sites (e.g., Audubon and Desert). H and LE estimates at Bondville and Brookings are comparable, and agree fairly well with the observations. Overall, the results indicate that the sequences of air temperature and specific humidity have a significant amount of information for partitioning the available energy between turbulent heat fluxes at sites with high SM and/or LAI values, but their information content significantly reduces at sites with low SM and/or LAI.
At all the study sites, LE estimates are more sparsely scattered around the 1:1 line compared to H estimates. This is due to the fact that errors in H estimates are due to uncertainties in CHN estimates (see Equation (1)) but errors in LE retrievals stem from uncertainties in both H and EF estimates (see Equation (3)). More sources of errors increase the scattering of LE estimates [44,81,83].
The discrepancy between the EF, H, and LE estimates and corresponding observations might be due to the measurement errors and simplistic assumptions, such as the constant monthly CHN,  At all the study sites, LE estimates are more sparsely scattered around the 1:1 line compared to H estimates. This is due to the fact that errors in H estimates are due to uncertainties in C HN estimates (see Equation (1)) but errors in LE retrievals stem from uncertainties in both H and EF estimates (see Equation (3)). More sources of errors increase the scattering of LE estimates [44,81,83].
The discrepancy between the EF, H, and LE estimates and corresponding observations might be due to the measurement errors and simplistic assumptions, such as the constant monthly C HN , constant daily EF, insignificant advection, and convectively well mixed boundary layer, which results in constant profiles of the potential temperature and specific humidity with height. Figure 7 compares the time series of daily H estimates from the VDA and open loop approaches with the measurements at the Arou, Willow Creek, Brookings, Bondwille, Audubon, and Desert, sites. Similarly, Figure 8 compares the time series of daily LE retrievals with the observations. As expected, the VDA approach generates more accurate H and LE values at the energy-limited sites (e.g., Arou, Willow Creek, Bondville, and Brookings). At these sites, H and LE estimates can capture the day-to-day fluctuations in the observations. At the water-limited sites (e.g., Audubon and specifically Desert), VDA cannot capture the daily variations in LE. The VDA approach overestimates (underestimates) H at Desert (Audubon). In the open loop approach, air temperature and specific humidity measurements are not assimilated into the VDA approach, and thus the initial guesses for C HN and EF are not updated. In this approach, the ABL potential temperature (θ) and specific humidity (q) are estimated by integrating Equations (4) and (5) forward in time using the initial guesses of C HN and EF. The large discrepancy between the turbulent heat fluxes estimates from VDA and the open loop at the wet and/or highly vegetated sites shows that the performance of VDA significantly improves at these sites by using the information content of atmospheric state variables.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 28 constant daily EF, insignificant advection, and convectively well mixed boundary layer, which results in constant profiles of the potential temperature and specific humidity with height.    Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 28 constant daily EF, insignificant advection, and convectively well mixed boundary layer, which results in constant profiles of the potential temperature and specific humidity with height.    humidity measurements are not assimilated into the VDA approach, and thus the initial guesses for CHN and EF are not updated. In this approach, the ABL potential temperature ( ) and specific humidity ( ) are estimated by integrating Equations (4) and (5) forward in time using the initial guesses of CHN and EF. The large discrepancy between the turbulent heat fluxes estimates from VDA and the open loop at the wet and/or highly vegetated sites shows that the performance of VDA significantly improves at these sites by using the information content of atmospheric state variables.    Table 5. Similarly, Table 6   The MAE and RMSE of half-hourly sensible and latent heat fluxes estimates from the VDA and open loop approaches at the six study sites are shown in Table 5. Similarly, Table 6 indicates the MAE and RMSE of daily H and LE retrievals. The high MAE and RMSE of turbulent heat fluxes estimates from VDA at the Desert and Audubon sites are attributed to their dry land and/or sparse vegetation. Compared to Desert and Audubon, Arou, Brookings, Bonville, and Willow Creek are wetter and have denser vegetation cover. Hence, the errors of VDA H and LE retrievals at Arou, Brookings, Bonville, and Willow Creek are lower than those of Desert and Audubon. Arou has the lowest MAE and RMSE values because it has a dense vegetation cover and the highest soil moisture.   4%), respectively. These results show that the assimilation of T a and q a leads to a higher improvement in the H and LE estimates at wet/densely vegetated sites.
The VDA approach iteratively improves the C HN and EF estimates by minimizing the difference between the ABL potential temperature and specific humidity estimates from Equations (4) and (5) (i.e., θ and q), and the corresponding values obtained from the reference-level air temperature and specific humidity via the MOST (i.e., θ SL and q SL ) (Appendix B). Thus, a close agreement between the θ and θ SL , and q and q SL (e.g., Arou) shows that the VDA approach can successfully update the initial guesses of C HN and EF and converges to their optimal values. On the other hand, a significant misfit between the θ SL and θ, and q SL and q (e.g., Desert) implies that the VDA cannot effectively improve the initial guesses of C HN and EF and converges to the inaccurate C HN and EF values. Figure 9 shows half-hourly θ estimates from VDA versus θ SL values. Similarly, Figure 10 indicates half-hourly q estimates versus q SL values. As shown, at Arou and Willow Creek, the θ and q estimates from VDA agree well with θ SL and q SL , respectively. This shows that in wet and/or densely vegetated sites, the VDA approach can take advantage of the significant amount of information in the sequences of air temperature and humidity to optimize C HN and EF, and consequently minimize the difference between θ and θ SL , and q and q SL . At the Brookings and Bondville sites, the θ and q estimates are in fairly good agreement with θ SL and q SL , respectively, implying that the time series of the atmospheric state variable have some information to constrain C HN and EF, and retrieve turbulent heat fluxes. At the Desert and Audubon sites, θ and q estimates are more scattered around the 1:1 line, showing the lack of sufficient information in the sequence of air temperature and humidity to accurately tune C HN and EF.
Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 28 Figure 9. Scatterplot of half-hourly potential temperature (θ) estimates from VDA versus at the six experimental sites. The RMSE of the θ (q) estimates at the Desert, Audubon, Bondville, Brookings, Willow Creek, and Arou sites are 2.56 K (0.0021 kg kg −1 ), 2.52 K (0.0015 kg kg −1 ), 1.36 K (0.0012 kg kg −1 ), 1.99 K (0.0012 kg kg −1 ), 1.04 K (0.0011 kg kg −1 ), and 0.95 K (0.0009 kg kg −1 ), respectively. As anticipated, the RMSEs of the θ and q estimates decrease as LAI and/or SM increase. Arou (with the highest SM) has the lowest RMSEs for θ and q. In contrast, Desert (with the lowest LAI and SM) has the highest RMSEs. As mentioned earlier, in dry sites (e.g., Desert and Audubon), evaporation is mainly controlled by the land surface state variable (i.e., LST). Hence, assimilating sequences of air temperature and specific humidity cannot robustly constrain C HN and EF, leading to larger errors in the θ and q estimates. Bondville and Brookings are neither as sparsely vegetated as Desert and Audubon nor as densely vegetated as Willow Creek and Arou. The RMSEs of the θ and q estimates at these sites are smaller than those of Desert and Audubon but larger than those of Willow Creek and Arou. Figure 11 indicates the mean diurnal cycles of the measured and estimated H and LE over the whole modeling period for the six study sites. As anticipated, the magnitude and phase of the diurnal cycles of retrieved sensible and latent heat fluxes from VDA agree well with those of the observations at the Arou and Willow Creek sites. For Audubon and especially Desert, there is a large discrepancy between the diurnal cycles of the estimated and observed H and LE, showing that the performance of the VDA approach significantly degrades in dry and/or sparsely vegetated sites. The results for the Bondville and Brookings sites show that assimilating the state variables of the atmosphere (i.e., θ and q) can capture the diurnal cycles of H and LE relatively well.

Conclusions
In this study, sequences of reference-level air temperature and humidity were assimilated into a VDA approach to partition the available energy between the sensible (H) and latent (LE) heat fluxes at six sites (namely, Arou, Audubon, Bondville, Brookings, Desert, and Willow Creek) with different vegetative and climatic conditions. The VDA approach takes advantage of the implicit information in the time series of reference-level air temperature and specific humidity (as the state variables of

Conclusions
In this study, sequences of reference-level air temperature and humidity were assimilated into a VDA approach to partition the available energy between the sensible (H) and latent (LE) heat fluxes at six sites (namely, Arou, Audubon, Bondville, Brookings, Desert, and Willow Creek) with different vegetative and climatic conditions. The VDA approach takes advantage of the implicit information in the time series of reference-level air temperature and specific humidity (as the state variables of the atmosphere) to estimate the neutral bulk heat transfer coefficient, C HN (that scales the sum of H and LE), and evaporative fraction, EF (that represents the partitioning between H and LE).
The sites with higher leaf area index (LAI) values showed larger C HN estimates. Additionally, at wet and/or densely vegetated sites, the variations in C HN estimates were consistent with those of LAI. This consistency weakened in dry and/or sparsely vegetated sites. Similarly, the EF estimates agreed well with the observations in terms of the magnitude and day-to-day dynamics in wet and/or densely vegetated sites (e.g., Arou and Willow Creek), but this agreement degraded in dry and/or sparsely vegetated sites (e.g., Desert and Audubon). The RMSE of EF estimates at Desert (dry barren land) was 0.198, which is 73.2% higher than that of 0.053 at Arou (wet grassland). These results show that the sequences of air temperature and specific humidity have a significant amount of information on the partitioning of available energy between H and LE in wet and/or densely vegetated sites. This information content decreases by the reduction of soil moisture and/or LAI. ), respectively. The RMSEs of the daily H and LE estimates increased as the site became drier and/or sparser in vegetation density. This is due to the fact that in dry and/or sparsely vegetated sites (e.g., Desert and Audubon), evapotranspiration is mainly controlled by the land surface state variable (i.e., land surface temperature) rather than the atmospheric state variables (i.e., reference-level air temperature and specific humidity). Hence, the coupling between EF and the atmospheric state variables weakens and the estimation of EF from the sequences of air temperature and specific humidity becomes uncertain. In contrast, at wet and/or densely vegetated sites (e.g., Arou and Willow Creek), the evaporative demand is primarily controlled by atmospheric state variables and the VDA system can extract that information to estimate H and LE.
The RMSEs of the atmospheric boundary layer (ABL) potential temperature (θ) estimates at the Desert, Audubon, Bondville, Brookings, Willow Creek and Arou sites were 2.56, 2.52, 1.36, 1.99, 1.04, and 0.95 K, respectively. The corresponding RMSEs for the ABL humidity (q) estimates were 0.0021, 0.0015, 0.0012, 0.0012, 0.0011, and 0.0009 kg kg −1 . The low RMSEs of the θ and q estimates at wet and/or densely vegetated sites (e.g., Arou and Willow Creek) indicated that the VDA approach can effectively update the two key unknowns (i.e., C HN and EF) and obtain their optimal values. In contrast, at dry and/or sparsely vegetated sites (e.g., Desert and Audubon), the VDA system cannot effectively update C HN and EF, leading to high RMSEs for θ and q. The highest and lowest RMSEs of θ and q estimates occurred at Desert and Arou, respectively. The RMSEs of θ and q estimates at the Bondville and Brookings sites fell within those of dry/sparsely vegetated and wet/densely vegetated sites.
The magnitude and phase of the diurnal cycles of the retrieved sensible and latent heat fluxes agreed well with the observations at wet/densely vegetated sites (e.g., Arou and Willow Creek). In contrast, for the dry/lightly vegetated sites (Audubon and especially Desert), a significant difference was found between the diurnal cycles of the retrieved and observed H and LE.
Future studies should focus on the synergistic assimilation of the LST (as the state variable of land surface) and the reference-level air temperature and specific humidity (as the state variables of atmosphere) to improve the turbulent heat flux estimates at the dry and/or sparsely vegetated sites.
In addition, future studies should be advanced by taking C HN as a function of LAI.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix B. Monin-Obukhov Similarity Theory (MOST)
The potential temperature, specific humidity, and wind speed at the top of the surface layer (i.e., θ SL , q SL , and u SL ) are found by expanding θ a , q a , and U from the reference height z re f to the bottom of the mixed layer (z SL ) via the Monin-Obukhov similarity theory (MOST) [77]: where k is the von Karman's constant; d is zero-plane displacement height (d is 2/3 of the vegetation height, z veg ); L v is the latent heat of vaporization; and Ψ h , Ψ m , and Ψ q are the stability functions for heat, momentum, and water vapor and are given by [77,111,112]: where ξ is the stability parameter, which is defined as the ratio of the measurement height to the Monin-Obhukov length z re f /L [77,93]. The Monin-Obukhov length is defined as: where R d is the gas constant for dry air (287 J kg −1 • C −1 ), R v is the gas constant for water vapor (461 J kg −1 o C −1 ), g is the gravitational acceleration, ρ is the air density, and c p is the specific heat capacity of air. Sensible (H) and latent heat fluxes (LE) are obtained respectively from Equations (1) and (3). u * is the friction velocity that is related to the wind speed measurements at the reference level (U) via [77]: where z 0m is the momentum roughness height. The neutral bulk heat transfer coefficient for heat (C HN ) can be related to the roughness length scales for heat (z 0h ) and momentum (z 0m ) via [50]: where z 0h is the roughness length for heat. The roughness length scales for heat and momentum are related through [50,111,112]: where B is the Stanton number. Rigden and Salvucci [77] related kB −1 to the vegetation height (z veg ) through: The following steps explain how θ SL and q SL are obtained from θ a and q a :