Evaluation of the Weak Constraint Data Assimilation Approach for Estimating Turbulent Heat Fluxes at Six Sites

A number of studies have estimated turbulent heat fluxes by assimilating sequences of land surface temperature (LST) observations into the strong constraint-variational data assimilation (SC-VDA) approaches. The SC-VDA approaches do not account for the structural model errors and uncertainties in the micrometeorological variables. In contrast to the SC-VDA approaches, the WC-VDA approach (the so-called weak constraint-VDA) accounts for the effects of structural and model errors by adding a model error term. In this study, the WC-VDA approach is tested at six study sites with different climatic and vegetative conditions. Its performance is also compared with that of SC-VDA at the six study sites. The results show that the WC-VDA produces 10.16% and 10.15% lower root mean square errors (RMSEs) for sensible and latent heat flux estimates compared with the SC-VDA approach. The model error term can capture errors in the turbulent heat flux estimates due to errors in LST and micrometeorological measurements, as well as structural model errors, and does not allow those errors to adversely affect the turbulent heat flux estimates. The findings also indicate that the estimated model error term varies reasonably well, so as to capture the misfit between predicted and observed net radiation in different hydrological and vegetative conditions. Finally, synthetically generated positive (negative) noises are added to the hydrological input variables (e.g., LST, air temperature, air humidity, incoming solar radiation, and wind speed) to examine whether the WC-VDA approach can capture those errors. It was found that the WC-VDA approach accounts for the errors in the input data and reduces their effect on the turbulent heat flux estimates.

Land surface temperature (LST) lies at the heart of the surface energy balance (SEB) equation, and has been used in many studies to estimate turbulent heat fluxes [14]. The existing methods for estimating turbulent heat fluxes based on LST observations are mainly divided into five groups. The first group (known as the triangle method) uses the empirical relationship between LST and vegetation indices to predict sensible (H) and latent (LE) heat fluxes [7,[15][16][17][18][19][20]. The second group (called the diagnostic approach) obtains sensible and latent heat fluxes by solving the SEB equation [5,6,8,11,13,21,22]. The third group (known as the surface temperature initiated closure method) estimates turbulent heat fluxes by integrating LST observations into the Penman-Monteith equation [23,24]. The fourth group combines hydrological variables (e.g., soil moisture and temperature, vegetation properties, topography, and meteorological forcing data) with a land surface model (e.g., Noah) within an Ensemble Kalman Filter (EnKF) method to predict sensible and latent heat fluxes [25][26][27][28][29][30][31].
The fifth group assimilates sequences of LST observations into the four dimensional variational data assimilation (4D-VDA) frameworks to retrieve H and LE [32][33][34][35][36][37][38][39][40]. This group estimates sensible and latent heat fluxes even for instances in which there are no LST data. The VDA approaches can obtain the key parameters of the land surface control on the partitioning of the available energy between H and LE [i.e., neutral bulk heat transfer coefficient (C HN ) and evaporative fraction (EF)]. C HN scales the sum of the turbulent heat fluxes (i.e., H + LE), is mainly a function of vegetation phenology, and is assumed to vary on a monthly time scale. EF scales partitioning between the turbulent heat fluxes (i.e., EF = LE/(H + LE)), and is assumed to be invariant during the assimilation window (i.e., 0900-1800 LT) in each day.
The studies in the fifth group often have used the strong constraint-VDA (SC-VDA) approach that presumes the model structure is perfect and the input data have no errors. However, the simplistic assumptions (e.g., constant daily EF and constant monthly C HN ) cause structural model errors in the VDA approaches. Moreover, there are errors in the micrometeorological input variables, and uncertainties in the specifications of the input parameters (e.g., albedo, soil emissivity, and soil thermal properties). Bateni et al. [35] and Xu et al. [38] added a model error term (ω) to the SEB equation within the VDA system [the so-called weak constraint-VDA (WC-VDA) methodology] to account for the errors resulting from uncertainties in the micrometeorological input variables, inaccurate model parameters, and simplistic model parameterizations. They tested the WC-VDA approach at the First International Satellite Land Surface Climatology Project Field Experiment (FIFE) and Huazhaizi sites. Both of these sites are subhumid grasslands.
The objectives of this study were to (1) evaluate the performance of the WC-VDA approach at six sites with different climatic and vegetative conditions, (2) compare the results of WC-VDA to those of SC-VDA, and (3) assess behavior of the model error term under different hydrological and vegetative conditions. This paper is organized as follows: Section 2 introduces the methodology including the heat diffusion equation, SEB equation, and adjoint state formulation; Section 3 explains the study area and data; results and discussions are given in Section 4; and finally, conclusions are reported in Section 5.

Heat Diffusion Equation
The heat diffusion equation describes the variation in soil temperature (T) with depth z and time t, and is given by, where C and K are the soil volumetric heat capacity (J·m −3 ·K −1 ) and thermal conductivity (W·m −1 ·K −1 ), respectively. To solve the heat diffusion equation, the boundary conditions at the top and bottom of the soil column need to be specified. The boundary condition at the top of the soil column is obtained by the land surface forcing, where G(t) is the ground heat flux at the land surface at time t. At the bottom boundary of the soil column, a Neumann boundary condition, is implemented as where l is the depth of the bottom boundary. According to Hu and Islam [41], the soil temperature at the depth of 0.3-0.5 m is almost invariant in a daily time-scale. Hence, l = 0.5 m was used in this study. The soil temperature profile is estimated by integrating the heat diffusion equation from a starting time (τ 0 ), where f 1 (z) represents the soil temperature profile at the initial time (τ 0 ). The readers are referred to Bateni et al. [35] for a detailed description of obtaining f 1 (z).

Surface Energy Balance (SEB) equation
The model error term (ω) was added to the SEB equation as follows, where H and LE are the sensible and latent heat fluxes, respectively. ω(t) is the model error term, which accounts for the model and observation errors, and Rn is the net radiation (W·m −2 ) that can be defined as, where α is the surface albedo (-), R ↓ S and R ↓ L are the incoming shortwave and longwave radiations (W·m −2 ), respectively, ε is the surface emissivity (-), σ is the Stefan-Boltzmann constant (5.67 The sensible heat flux can be calculated via where ρ is the air density (kg·m −3 ), C P is the air heat capacity (1012 J·kg −1 ·K −1 ), and U and T a are the reference-height wind speed (m·s −1 ) and air temperature (K), respectively. Bulk heat transfer coefficient (C H ) is related to the heat transfer coefficient under neutral atmospheric condition (C HN ) and the atmospheric stability correction function (f (Ri)) via [36,38,40,42], Remote Sens. 2018, 10, 1994 4 of 23 where Ri is the Richardson number. C HN depends on the geometry of the land surface and constitutes the first unknown parameter of the WC-VDA scheme. C HN is assumed to vary on a monthly (30-day) time scale. Evaporative fraction (EF) is the second unknown parameter of the WC-VDA approach. It is almost constant for near-peak radiation hours on days without precipitation [43]. In this study, EF was assumed to be constant during the assimilation window (0900-1800 LT),

Adjoint State Formulation
The key unknown parameters (C HN and EF) and the model error term (ω) in the WC-VDA approach were obtained by minimizing the cost function J, which is defined as follows, where N is the total number of LST observations during the assimilation period, and R −1 T is the inverse covariance matrix of LST observation errors, B −1 R and B −1 EF are the inverse background error covariance matrices of R and EF, respectively, and Q −1 ω is the inverse model error covariance matrix. The first term on the right-hand side of the cost function represents the misfit between the LST observations and predictions. The second and third terms represent the misfit between the parameters estimates (R and EF) and their prior values (R and EF i ). C HN is related to R via C HN = exp(R) to make it strictly positive and meaningful. The key unknown parameters (C HN and EF) and the model error term (ω) were obtained by minimizing the difference between the LST observations and predictions from the heat diffusion equation. Following Bateni et al. [35], R −1 T , B −1 R , and B −1 EF are assumed to be diagonal matrices of numerically constant values whose relative magnitude controls the rate of convergence of the iterative minimization procedure explained in Appendix A. The magnitudes of the diagonal elements in R −1 T , B −1 R , and B −1 EF are set to 0.01 K −2 , 1000, and 1000, respectively [34,35]. The fourth term is the heat diffusion equation (the physical constraint), which is adjoined to the cost function via the Lagrange multiplier Λ. The last term adds the model error (ω(t)) to the objective function to account for the uncertainties in the forcing as well as the parametric and structural model errors. Following Bateni et al. [35], Xu et al. [37], and Reichle [44], the inverse model error covariance can be expressed as, This type of error covariance model is often used when the covariance structure of errors is poorly known [44,45]. More information about the weak constraint VDA approach and the model error term based on the Bayes' theorem can be found in Rodgers [46] and Tremolet [47]. The σ ω and τ are the standard deviation and decorrelation timescale, which are set to 100 W·m −2 and 6 h, respectively [35]. The strong constraint-VDA (SC-VDA) approach does not include the model error term and treats the physical model as perfect. While, the weak constraint-VDA incorporates the model error term (ω) into the cost function J to account for model uncertainties.
The optimal values of C HN , EF, and ω were obtained by setting the first variation of J to zero (i.e., δJ = 0). Imposing δJ = 0 leads to a number of Euler-Lagrangian equations (shown in Appendix B), which should be solved simultaneously through an iterative procedure to obtain optimal values for C HN , EF, and ω. The Bias and Root Mean Square Error (RMSE) metrics were used to access the performances of the VDA approach, where P i and O i are the predicted and observed values at time step i, respectively.

Study Domain and Data
The WC-VDA approach was tested at the Mead, Audubon, Brookings, and Willow Creek sites  The half-hourly micrometeorological data (wind speed, air temperature and humidity, incoming shortwave and longwave radiations) were measured at the six sites by Automatic Weather Stations (AWS). At all the sites, the sensible and latent heat fluxes were measured by the eddy covariance (EC) instrument. Albedo and LAI data were obtained from the Global Land Surface Satellite (GLASS) product [51,52] (http://glass-product.bnu.edu.cn/). The soil heat capacity (C) and conductivity (K) depend on the soil moisture and texture [53][54][55]. In this study, soil texture and moisture at each site were obtained from the Harmonized World Soil Database (HWSD) and Heihe River Basin (HRB) Digital Soil Mapping product, respectively [56]. Table 2 represents the neutral bulk heat transfer coefficient (C HN ) estimates and LAI values for the six experimental sites. As shown, changes in the C HN estimates are consistent with those of LAI. The C HN values are generally higher at the densely vegetated sites (Arou, Brookings, and Willow Creek), and lower at the sparsely vegetated sites (Audubon and Desert). At the Mead site, C HN increases with the growth of the crop and reaches its maximum in the fourth period (DOY 213-243) where LAI is at its peak value. The higher C HN value at the Willow Creek site is attributed to its dense vegetation cover (LAI = 4.8). There is no vegetation coverage at the Desert site (LAI is zero), and hence the C HN estimates at Desert are lower than those of the other sites. Although LAI is invariant at the Desert site, Remote Sens. 2018, 10, 1994 6 of 23 the C HN estimates vary in different monthly modeling periods. This is because C HN depends on not only vegetation phenology (LAI) but also wind speed, friction velocity, and solar elevation [57,58]. Desert

Evaporative Fraction
The estimated evaporative fraction (EF) values from the WC-VDA approach are compared with the observations in Figure 1. As shown, the estimated EF values are in good agreement with the observations, and show a characteristic response to the wetting and dry down events although rainfall measurements are not used in the WC-VDA approach. For example, at the Audubon site, the estimated EF values increase sharply when precipitation occurs (e.g., DOYs 145, 201, 206, 226, 251, and 236) and decrease in the following days (e.g., DOYs 148, 203, 207, 228, 253, and 237) when the land surface becomes dry. When the soil is very wet, the sensible heat flux can become negative, leading to EF observation values of larger than 1.0 [43]. This can be seen in the Mead, Brookings, and Willow Creek sites. To avoid the numerical instability, EF in the WC-VDA approach is set to be less than 0.99, which results in the EF estimations being smaller than the observations.  Figure 2 shows the estimated LST from the WC-VDA system versus measurements at the six sites. As indicated, the LST estimations are in good agreement with the observations and mainly fall around the 1:1 line. The Bias and RMSE of estimated LST by the WC-VDA approach are shown in Table 3. For comparison, the corresponding values from the SC-VDA approach are also shown in the same Table. The Bias and RMSE of LST estimates from WC-VDA are smaller than those of SC-VDA in all the sites. The six-site-averaged RMSE of LST predictions decreases from 1.97 K to 1.77 K by the addition of the model error term within the VDA system. The results show that the WC-VDA approach can update CHN and EF more effectively by absorbing spike errors in LST and forcing data and not allowing measurement errors to affect the VDA performance.  Figure 2 shows the estimated LST from the WC-VDA system versus measurements at the six sites. As indicated, the LST estimations are in good agreement with the observations and mainly fall around the 1:1 line. The Bias and RMSE of estimated LST by the WC-VDA approach are shown in Table 3. For comparison, the corresponding values from the SC-VDA approach are also shown in the same Table. The Bias and RMSE of LST estimates from WC-VDA are smaller than those of SC-VDA in all the sites. The six-site-averaged RMSE of LST predictions decreases from 1.97 K to 1.77 K by the addition of the model error term within the VDA system. The results show that the WC-VDA approach can update C HN and EF more effectively by absorbing spike errors in LST and forcing data and not allowing measurement errors to affect the VDA performance.     Figure 3a,b shows the daytime-averaged (0900-1800 LT) sensible and latent heat flux estimates at the six sites during DOY 121-273. The estimated turbulent heat fluxes agree well with the ground measurements, implying that the WC-VDA approach can yield accurate turbulent heat flux estimates. However, the turbulent heat flux estimates degrade in wet periods (e.g., DOY 180-200 at the Brookings site and DOY 180-220 at the Willow Creek site). The Bias and RMSE of turbulent heat flux estimates from the WC-VDA approach at the Willow Creek site (higher precipitation and denser vegetation cover) are larger than those of Audubon and Desert sites (lower precipitation and sparser canopy cover) (Figure 4a,b). This is because drying rate is mainly controlled by the land surface state variable (i.e., LST) in dry and slightly vegetated sites. But, in wet and densely vegetated sites, evaporation is mainly influenced by the atmospheric state variables (i.e., air temperature and specific humidity) [59,60]. Figure 4a,b compares the half-hourly sensible and latent heat flux estimates from the WC-VDA method with the EC measurements at the six study sites. The estimated sensible and latent heat fluxes from the WC-VDA approach are in good consistency with the observations, and the scatterplots mainly fall around the 1:1 line. The WC-VDA model tends to overestimate (underestimate) H when it is smaller (larger) than 100 W·m −2 . At the Mead, Audubon, and Brookings sites, the WC-VDA model tends to overestimate LE when it is larger than 200 W·m −2 . Also, LE estimates are more sparsely distributed around the 1:1 line at the Willow Creek site (forest). This happens because performance of the WC-VDA degrades in heavily vegetated sites as LE is controlled dominantly by the atmospheric variables rather than the land surface processes. Similar results were also found in Xu et al. [36,39]. The misfits between the model estimates and observations are mainly due to the simplistic model assumptions (i.e., daily constant EF and monthly constant C HN ), and measurement errors. Figure 4a,b also shows that LE estimates are more scattered around the 45-degree line than H. This is because the uncertainty of the H estimates comes from the errors in the C HN and LST estimates (see Equation (7)). LE is obtained by re-writing Equation (9)      The Bias and RMSE of half-hourly turbulent heat flux estimates from the WC-VDA approach are summarized in Table 4a. For comparison, corresponding values from SC-VDA are also shown in the same table. For half-hourly sensible heat flux, the six-site-averaged RMSE from the WC-VDA approach is 59.51 W·m −2 , which is 10.16% lower than the RMSE of 66.24 W·m −2 from SC-VDA. For latent heat flux, the six-site-averaged RMSE is 75.06 W·m −2 for the SC-VDA approach, and 67.44 W·m −2 for the WC-VDA approach. In a similar effort, the statistical indices for the daily H and LE estimates from the WC-VDA and SC-VDA approaches are shown in Table 4b. The six-site-averaged RMSE of estimated daily H and LE by SC-VDA are 44.58 W·m −2 and 45.32 W·m −2 , respectively. The WC-VDA approach reduces the aforementioned RMSEs by 16.22% and 15.60%, implying that the WC-VDA method can absorb noise in the LST and micrometeorological forcing data, and does not allow measurement errors to propagate into EF and C HN , and consequently turbulent heat fluxes. Overall, the results in Table 4a,b show that WC-VDA outperforms SC-VDA.   Figure 5 shows time series of daily-averaged model error term (ω) and misfit between the predicted and observed net radiation (Rn pre − Rn obs ) from the WC-VDA approach at the six study sites. The misfit between the predicted and observed R n is due to the difference between T obs and T pre , which itself is mainly because of the simplistic assumptions (i.e., constant soil thermal conductivity (P) and heat capacity (C) in the heat diffusion equation, monthly constant C HN , and daily constant EF). As shown in Figure 5, the estimated ω values show good consistency with the (Rn pre − Rn obs ). The model error term increases or decreases to capture the misfit between Rn pre − Rn obs . For example, the model error term estimates at the Desert site with larger Rn pre − Rn obs are higher than those at the Brookings and Willow Creek sites with smaller Rn pre − Rn obs . The positive model error term at the Audubon (e.g., DOY 159-163) and Desert (e.g., DOY 196-200) sites is consistent with the positive (Rn pre − Rn obs ). For positive (Rn pre − Rn obs ), the model error term becomes positive to capture the overestimated Rn, thereby not allowing the errors in Rn pre to propagate into the sensible and latent heat fluxes. In contrast, when (Rn pre − Rn obs ) is negative (i.e., Rn is underestimated), the model error term becomes negative to balance the SEB equation. At some of the sites (e.g., Mead and Arou), there is a weaker consistency between ω and (Rn pre − Rn obs ). This is because the noise in the forcing data can affect the patterns of ω and weaken the consistency between ω and (Rn pre − Rn obs ).

Model Error Analysis
contrast, when (Rnpre − Rnobs) is negative (i.e., Rn is underestimated), the model error term becomes negative to balance the SEB equation. At some of the sites (e.g., Mead and Arou), there is a weaker consistency between ω and (Rnpre − Rnobs). This is because the noise in the forcing data can affect the patterns of ω and weaken the consistency between ω and (Rnpre − Rnobs). To further evaluate the model robustness, the frequency histograms of the model error term (ω) and misfit between the predicted and observed net radiation (Rnpre − Rnobs) are shown in Figure 6a,b, respectively. As illustrated, the frequency histograms of ω and (Rnpre − Rnobs) are comparable, and they range from −30 W·m −2 to 30 W·m −2 at most of the sites. At the Mead, Brookings, and Desert sites, the distribution of (Rnpre − Rnobs) is negatively-biased. Most likely this happens due to the inexact soil To further evaluate the model robustness, the frequency histograms of the model error term (ω) and misfit between the predicted and observed net radiation (Rn pre − Rn obs ) are shown in Figure 6a,b, respectively. As illustrated, the frequency histograms of ω and (Rn pre − Rn obs ) are comparable, and they range from −30 W·m −2 to 30 W·m −2 at most of the sites. At the Mead, Brookings, and Desert sites, the distribution of (Rn pre − Rn obs ) is negatively-biased. Most likely this happens due to the inexact soil thermal properties (C and K), and the simplistic assumptions (e.g., constant daily EF and constant monthly C HN ). Remarkably, the ω estimates at these sites have a negative bias to capture the negative misfit in (Rn pre − Rn obs ). At the Arou and Willow Creek sites, the frequency histograms of ω and Rn pre − Rn obs are unbiased, and they both distribute almost symmetrically around 0 W·m −2 . At the Audubon site, there is a weaker consistency between the frequency histograms of ω and Rn pre − Rn obs . A possible explanation for this is because the model error term accounts for not only the structural model deficiencies (e.g., constant daily EF, constant monthly C HN , etc.), which leads to the underestimation or overestimation of Rn, but also the noise in the forcing variables, as well as parametric model errors. It is apparent that the noisy forcing and parametric model errors (pertains to having the inaccurate values of input parameters) can affect the distribution of ω and weaken the consistency between ω and (Rn pre − Rn obs ).
Remote Sens. 2018, 10, x FOR PEER REVIEW 16 of 23 thermal properties (C and K), and the simplistic assumptions (e.g., constant daily EF and constant monthly CHN). Remarkably, the ω estimates at these sites have a negative bias to capture the negative misfit in (Rnpre − Rnobs). At the Arou and Willow Creek sites, the frequency histograms of ω and Rnpre − Rnobs are unbiased, and they both distribute almost symmetrically around 0 W·m −2 . At the Audubon site, there is a weaker consistency between the frequency histograms of ω and Rnpre − Rnobs. A possible explanation for this is because the model error term accounts for not only the structural model deficiencies (e.g., constant daily EF, constant monthly CHN, etc.), which leads to the underestimation or overestimation of Rn, but also the noise in the forcing variables, as well as parametric model errors.
It is apparent that the noisy forcing and parametric model errors (pertains to having the inaccurate values of input parameters) can affect the distribution of ω and weaken the consistency between ω and (Rnpre − Rnobs).  Two numerical experiments were conducted to examine whether the VDA approach with model uncertainty can capture errors in the meteorological input data. In the first (second) numerical experiment, positive (negative) noises with mean of 4 K (−4 K), 2 K (−2 K), 0.1 (−0.1), 30 W·m −2 (−30 W·m −2 ), and 1 m·s −1 (−1 m·s −1 ) and standard deviation of 2 K, 1 K, 0.02, 10 W·m −2 , and 0.2 m s −1 were added respectively to the measured LST, air temperature, air humidity, incoming solar radiation, and wind speed at the Audubon site. The noises generated by the random Gaussian number generator were used to test the performance of the WC-VDA and SC-VDA approaches. Figure 7a shows the errors of H and LE estimates when the positively noised input variables were used in the WC-VDA and SC-VDA approaches. Analogously, Figure 7b indicates the errors of estimated H and LE values for negatively noised input variables. As shown, the WC-VDA approach can capture both the positive and negative noises in the input variables and reduces their effect on the H and LE estimates. In contrast, the SC-VDA approach is more susceptible to errors in the input data, which can be seen by higher errors in its H and LE estimates. As shown in Figure 7a,b, the time series of errors in the LE estimates from the WC-VDA and SC-VDA approaches indicates a significant increase from early August until mid-September (i.e., around DOYs 212-255). This happens because the Audubon site experiences high rainfall in this period. Performance of the VDA degrades in wet periods because LE is in Stage-I (energy-limited), is controlled mainly by the atmospheric variables, and the coupling between LST and EF weakens [59,60]. In this case, the retrieval of EF from LST becomes more uncertain, which ultimately causes large errors in the LE estimates. From mid-September towards the end of the modeling period, there is no rainfall and consequently the uncertainty of LE retrievals decreases.
Two numerical experiments were conducted to examine whether the VDA approach with model uncertainty can capture errors in the meteorological input data. In the first (second) numerical experiment, positive (negative) noises with mean of 4 K (−4 K), 2 K (−2 K), 0.1 (−0.1), 30 W·m −2 (−30 W·m −2 ), and 1 m·s −1 (−1 m·s −1 ) and standard deviation of 2 K, 1 K, 0.02, 10 W·m −2 , and 0.2 m s −1 were added respectively to the measured LST, air temperature, air humidity, incoming solar radiation, and wind speed at the Audubon site. The noises generated by the random Gaussian number generator were used to test the performance of the WC-VDA and SC-VDA approaches. Figure 7a shows the errors of H and LE estimates when the positively noised input variables were used in the WC-VDA and SC-VDA approaches. Analogously, Figure 7b indicates the errors of estimated H and LE values for negatively noised input variables. As shown, the WC-VDA approach can capture both the positive and negative noises in the input variables and reduces their effect on the H and LE estimates. In contrast, the SC-VDA approach is more susceptible to errors in the input data, which can be seen by higher errors in its H and LE estimates. As shown in Figure 7a,b, the time series of errors in the LE estimates from the WC-VDA and SC-VDA approaches indicates a significant increase from early August until mid-September (i.e., around DOYs 212-255). This happens because the Audubon site experiences high rainfall in this period. Performance of the VDA degrades in wet periods because LE is in Stage-I (energy-limited), is controlled mainly by the atmospheric variables, and the coupling between LST and EF weakens [59,60]. In this case, the retrieval of EF from LST becomes more uncertain, which ultimately causes large errors in the LE estimates. From mid-September towards the end of the modeling period, there is no rainfall and consequently the uncertainty of LE retrievals decreases. (a)

Conclusions
In this study, the weak constraint-variational data assimilation (WC-VDA) method was tested at six sites with contrasting hydrological conditions. The key unknown parameters of the WC-VDA approach are neutral bulk transfer coefficient (CHN) and evaporative fraction (EF). The results showed that the neutral heat transfer coefficient (CHN) increases with the increase of leaf area index (LAI) although no information on vegetation phenology was used as input. The variations in evaporative fraction (EF) estimates are consistent with the drying and wetting events, while no precipitation or soil moisture observations were used in the WC-VDA approach. The estimated sensible and latent heat fluxes from the WC-VDA approach were validated with the eddy covariance (EC) observations at the six study sites, and compared to those of strong constraint-variational data assimilation (SC-VDA). For the WC-VDA scheme, the six-site-averaged root mean square errors (RMSEs) of sensible and latent heat fluxes are 59.51 W·m −2 and 67.44 W·m −2 , respectively, which are 10.16% and 10.15% lower than those of SC-VDA. The findings indicated that the WC-VDA approach outperforms SC-VDA because the model error term can account for errors resulting from uncertainties in the micrometeorological input variables and inaccurate model parameters, as well as simplistic model parameterizations. This study evaluated the performance of the model error term in different hydrological and vegetative conditions. It was found that the model error term shows good consistency with the misfit between the predicted and observed net radiation (Rnpre − Rnobs), implying

Conclusions
In this study, the weak constraint-variational data assimilation (WC-VDA) method was tested at six sites with contrasting hydrological conditions. The key unknown parameters of the WC-VDA approach are neutral bulk transfer coefficient (C HN ) and evaporative fraction (EF). The results showed that the neutral heat transfer coefficient (C HN ) increases with the increase of leaf area index (LAI) although no information on vegetation phenology was used as input. The variations in evaporative fraction (EF) estimates are consistent with the drying and wetting events, while no precipitation or soil moisture observations were used in the WC-VDA approach. The estimated sensible and latent heat fluxes from the WC-VDA approach were validated with the eddy covariance (EC) observations at the six study sites, and compared to those of strong constraint-variational data assimilation (SC-VDA). For the WC-VDA scheme, the six-site-averaged root mean square errors (RMSEs) of sensible and latent heat fluxes are 59.51 W·m −2 and 67.44 W·m −2 , respectively, which are 10.16% and 10.15% lower than those of SC-VDA. The findings indicated that the WC-VDA approach outperforms SC-VDA because the model error term can account for errors resulting from uncertainties in the micrometeorological input variables and inaccurate model parameters, as well as simplistic model parameterizations. This study evaluated the performance of the model error term in different hydrological and vegetative conditions. It was found that the model error term shows good consistency with the misfit between the predicted and observed net radiation (Rn pre − Rn obs ), implying that the model error term can capture the difference between the predicted and observed net radiation in various hydrological environments.
This hinders the errors in the estimated net radiation values to transfer into the turbulent heat flux estimates. Finally, numerical experiments indicated that the WC-VDA approach can capture the synthetically added positive and negative noises to the nominal input variables and reduce their effect on the H and LE estimates.

Acknowledgments:
The authors would like to thank all the scientists, engineers, and students who participated in HiWATER field campaigns. We thank two reviewers for their valuable comments that greatly improved the presentation of this paper.

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

Appendix A. Difference between the 3D-and 4D-VDA Approaches
The principle of three dimensional variational data assimilation (3D-VDA) is to obtain an approximate solution by minimizing the cost function J, where X and X b are the analysis model state and background model state, respectively, y is the vector of observations, H is the observation operator, and B and R are the covariance matrix of the background and observation errors, respectively. Four dimensional variational data assimilation (4D-VDA) is an extension of 3D-VDA, which is used for time series of observations. The 4D-VDA method allows comparison of the model state with observations at the appropriate times. Over a specified assimilation period, the cost function J for the 4D-VDA can be shown as, where y i and X i are the observation and model state at time step i, respectively. R i and H i are the observation error covariance matrix and observation operator at time step i, respectively, and N is the total number of time steps in the assimilation period. The 3D-VDA method collects observations during analysis time and does not require model integration. Therefore, the analysis increment does not evolve over time and requires less computing resources [61]. While, in the 4D-VDA method, all the observations in the assimilation period are incorporated into the model. The 4D-VDA method uses tangent linear and adjoint models to generate the propagation of analytical increments during the assimilation period, and requires more computational resources than 3D-VDA [62][63][64].