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

Abstract

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.

1. Introduction

Accurate estimation of turbulent heat fluxes (sensible and latent heat fluxes) is required in agricultural, hydrological, and meteorological applications [1,2,3]. Flux towers have been installed across the world to measure turbulent heat fluxes. However, they are expensive and have a sparse distribution (i.e., limited spatial representation). Consequently, a number of models have been developed to estimate turbulent heat fluxes based on in-situ and remotely sensed observations [4,5,6,7,8,9,10,11,12,13].

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.

2. Methodology

2.1. Heat Diffusion Equation

The heat diffusion equation describes the variation in soil temperature (T) with depth z and time t, and is given by,

$$C\frac{\partial T\left(z,t\right)}{\partial t}=K\frac{{\partial }^{2}T\left(z,t\right)}{\partial z^2},$$

(1)

where C and K are the soil volumetric heat capacity (J·m⁻³·K⁻¹) and thermal conductivity (W·m⁻¹·K⁻¹), 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,

$$\frac{-KdT\left(z=0,t\right)}{dz}=G\left(t\right),$$

(2)

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

$$\frac{\partial T\left(l,t\right)}{\partial z}=0,$$

(3)

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 (τ₀),

$$f\left(z,{\tau }_0\right)=f_1\left(z\right),$$

(4)

where f₁(z) represents the soil temperature profile at the initial time (τ₀). The readers are referred to Bateni et al. [35] for a detailed description of obtaining f₁(z).

2.2. Surface Energy Balance (SEB) equation

The model error term (ω) was added to the SEB equation as follows,

$$Rn=H+\mathrm{LE}+ G+\omega \left(t\right),$$

(5)

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⁻²) that can be defined as,

$$Rn=\left(1-\alpha \right)R_S^{\downarrow }+R_L^{\downarrow }-\epsilon \sigma T^4,$$

(6)

where α is the surface albedo (–), $R_S^{\downarrow }$ and $R_L^{\downarrow }$ are the incoming shortwave and longwave radiations (W·m⁻²), respectively, ε is the surface emissivity (–), σ is the Stefan-Boltzmann constant (5.67 × 10⁻⁸ W· m⁻²·K⁻⁴).

The sensible heat flux can be calculated via

$$H=\rho c_P{C}_{H}U\left(T-T_a\right),$$

(7)

where ρ is the air density (kg·m⁻³), C_P is the air heat capacity (1012 J·kg⁻¹·K⁻¹), and U and T_a are the reference-height wind speed (m·s⁻¹) 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],

$$C_H=C_{HN}f\left(Ri\right)=C_{HN}\left[1+2\left(1-e^{10Ri}\right)\right],$$

(8)

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),

$$\mathrm{EF}=\mathrm{LE}/\left(H+\mathrm{LE}\right),$$

(9)

2.3. 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,

$$\begin{array}{ll}J\left(T,R,\mathrm{EF},\lambda \right)=& {{\displaystyle \sum }}_{i=1}^N{{\displaystyle \int }}_{t_0}^{t_1}{\left[T_{OBS,i}\left(t\right)-T_i\left(t\right)\right]}^T{R}_T^{-1}\left[T_{OBS,i}\left(t\right)-T_i\left(t\right)\right]dt\\ & +{\left(R-R^{\prime }\right)}^T{B}_R^{-1}\left(R-R^{\prime }\right)+{{\displaystyle \sum }}_{i=1}^N{\left({\mathrm{EF}}_i-{\mathrm{EF}}_i^{\prime }\right)}^T{B}_{\mathrm{EF}}^{-1}\left({\mathrm{EF}}_i-{\mathrm{EF}}_i^{\prime }\right)\\ & +2{{\displaystyle \sum }}_{i=1}^N{{\displaystyle \int }}_{t_0}^{t_1}{{\displaystyle \int }}_0^l{\lambda }_i\left(z,t\right)\left[\frac{\partial T_i\left(z,t\right)}{\partial t}-\left(\frac{K}{C}\right)\frac{\partial T_i{\left(z,t\right)}^2}{\partial z^2}\right]dzdt\\ & +{{\displaystyle \sum }}_{i=1}^N{{\displaystyle \int }}_{t_0}^{t_1}{{\displaystyle \int }}_{t_0}^{t_1}{\omega }_i\left(t\prime \right)Q_{\omega }^{-1}\left(t^{\prime },t\prime \prime \right){\omega }_i\left(t\prime \prime \right)d{t}^{\prime }d{t}^{″},\end{array}$$

(10)

where N is the total number of LST observations during the assimilation period, and $R_T^{-1}$ is the inverse covariance matrix of LST observation errors, $B_R^{-1}$ and $B_{\mathrm{EF}}^{-1}$ are the inverse background error covariance matrices of R and EF, respectively, and $Q_{\omega }^{-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^{\prime }$ and ${\mathrm{EF}}_i^{\prime }$). 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_T^{-1}$, $B_R^{-1}$, and $B_{\mathrm{EF}}^{-1}$ 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_T^{-1}$, $B_R^{-1}$, and $B_{\mathrm{EF}}^{-1}$ are set to 0.01 K⁻², 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,

$$Q_{\omega }^{-1}\left(t\prime ,t\prime \prime \right)={\sigma }_{\omega }^2\exp (\left|t\prime -t\prime \prime \right|/\tau ).$$

(11)

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⁻² 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,

$$BIAS=\frac{{{\displaystyle \sum }}_{i=1}^N\left(P_i-O_i\right)}{N},$$

(12)

$$RMSE=\sqrt{\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N{(P_i-O_i)}^2},$$

(13)

where $P_i$ and $O_i$ are the predicted and observed values at time step i, respectively.

3. Study Domain and Data

The WC-VDA approach was tested at the Mead, Audubon, Brookings, and Willow Creek sites (Day of Year (DOY) 121–273, 2005) in the US (http://ameriflux.lbl.gov/) [48], as well as the Arou and Desert sites (DOY 121–273, 2015) in China [3,49,50] (http://card.westgis.ac.cn/). Characteristics of these sites are summarized in

Table 1. Characteristics of the six study sites.

Site	Location	Land Cover	LAI	fc	SM	Elevation (m)
Mead	41.18N, 96.44W	Cropland	1.75	0.53	0.25	363
Arou	38.04N, 100.46E	Grassland	2.53	0.72	0.34	3033
Audubon	31.59N, 110.51W	Grassland	0.39	0.17	0.07	1469
Brookings	44.34N, 96.83W	Grassland	2.24	0.67	0.45	510
Willow Creek	45.80N, 90.08W	Forest	4.8	0.91	0.30	520
Desert	42.11N, 100.99E	Barren land	0	0	0.01	1000

LAI, f_c, and SM represent the mean leaf area index (m²·m⁻²), vegetation cover fraction (-), and soil moisture (m³·m⁻³) over the growing period (DOY 121–273), respectively.

. The experimental sites have four vegetation cover types: grassland, cropland, forest, and barren land. Mead is an irrigated cropland. Willow Creek is a deciduous broadleaf forest with dense vegetation cover (high leaf area index, LAI). Arou, Audubon, and Brookings are grasslands with a large range of soil moisture (0.07–0.45 m³·m⁻³) and LAI (0.39–2.53 m²·m⁻²). Desert is barren land with sparse vegetation and dry soil. The wide range of soil moisture and LAI values allows us to assess performance of the WC-VDA approach and variations of the model error term under different environmental conditions.

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].

4. Results and Discussions

4.1. Neutral Bulk Heat Transfer Coefficient

Table 2. Neutral bulk heat transfer coefficient (C_HN) estimates from the weak constraint-variational data assimilation (WC-VDA) model, and corresponding leaf area index (LAI) values.

Site		121–151	152–181	182–212	213–243	244–273
Mead	CHN	0.22 × 10−2	0.40 × 10−2	1.22 × 10−2	1.48 × 10−2	0.78 × 10−2
	LAI	0.77	1.28	2.42	2.70	1.54
Arou	CHN	0.12 × 10−2	0.46 × 10−2	3.35 × 10−2	5.12 × 10−2	1.98 × 10−2
	LAI	0.46	1.77	3.79	4.22	2.42
Audubon	CHN	0.31 × 10−2	0.35 × 10−2	0.42 × 10−2	0.31 × 10−2	0.25 × 10−2
	LAI	0.23	0.30	0.35	0.54	0.53
Brookings	CHN	0.72 × 10−2	1.54 × 10−2	3.12 × 10−2	1.55 × 10−2	0.62 × 10−2
	LAI	1.30	2.37	2.85	2.74	1.94
Willow Creek	CHN	6.33 × 10−2	7.45 × 10−2	8.61 × 10−2	14.35 × 10−2	8.13 × 10−2
	LAI	2.64	5.35	5.88	5.73	4.40
Desert	CHN	0.13 × 10−3	0.14 × 10−3	0.18 × 10−3	0.19 × 10−3	0.12 × 10−3
	LAI	0	0	0	0	0

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, 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].

4.2. 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.

4.3. Land Surface Temperature

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. Bias and Root Mean Square Error (RMSE) of half-hourly land surface temperature (LST) estimates from the weak constraint-variational data assimilation (WC-VDA) and strong constraint-variational data assimilation (SC-VDA) approaches at the six experimental sites.

	Bias (K)		RMSE (K)
Study Sites	WC-VDA	SC-VDA	WC-VDA	SC-VDA
Mead	0.82	−1.23	1.55	1.71
Arou	0.30	0.5	1.97	2.07
Audubon	0.06	0.23	2.37	2.75
Brookings	−0.76	−0.86	1.71	1.72
Willow Creek	0.21	0.56	0.99	1.06
Desert	−0.41	−0.92	2.05	2.53
Six-sites-average	0.04	−0.29	1.77	1.97

. 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.

4.4. Sensible and Latent Heat Fluxes

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⁻². At the Mead, Audubon, and Brookings sites, the WC-VDA model tends to overestimate LE when it is larger than 200 W·m⁻². 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) as follows: LE = EF/(1 − EF) × H. Hence, the uncertainty of LE values is due to the errors in the EF, C_HN, and LST retrievals. More sources of errors increase scattering of LE estimates.

The Bias and RMSE of half-hourly turbulent heat flux estimates from the WC-VDA approach are summarized in

Table 4. a. Statistical indices of half-hourly sensible (H) and latent (LE) estimates from the weak constraint-variational data assimilation (WC-VDA) and strong constraint-variational data assimilation (SC-VDA) approaches at the six experimental sites. b. Statistical indices of daily H and LE estimates from the WC-VDA and SC-VDA approaches at the six experimental sites.

a
Study Site	Method	H (W·m−2)		LE (W·m−2)
		Bias	RMSE	Bias	RMSE
Mead	WC-VDA (SC-VDA)	8.73 (10.21)	64.60 (69.25)	−3.81 (−8.65)	63.07 (75.55)
Arou	WC-VDA (SC-VDA)	−11.64 (−15.75)	52.28 (55.63)	−0.49 (5.32)	70.54 (79.88)
Audubon	WC-VDA (SC-VDA)	16.25 (20.21)	56.44 (67.82)	7.95 (12.19)	54.37 (57.55)
Brookings	WC-VDA (SC-VDA)	4.77 (10.88)	42.72 (45.35)	−7.31 (−11.60)	81.36 (83.86)
Willow Creek	WC-VDA (SC-VDA)	7.67 (22.17)	86.74 (91.18)	19.88 (32.49)	110.82 (125.60)
Desert	WC-VDA (SC-VDA)	7.08 (12.68)	54.30 (68.21)	0.95 (1.25)	24.45 (27.91)
b
Study Site	Method	H (W·m−2)		LE (W·m−2)
		Bias	RMSE	Bias	RMSE
Mead	WC-VDA (SC-VDA)	5.32 (9.45)	40.89 (43.28)	−3.05 (−7.51)	40.75 (47.62)
Arou	WC-VDA (SC-VDA)	−0.37 (−3.43)	31.56 (38.69)	3.56 (17.74)	32.71 (50.58)
Audubon	WC-VDA (SC-VDA)	12.04 (16.15)	31.95 (42.09)	0.51 (−7.32)	28.86 (22.98)
Brookings	WC-VDA (SC-VDA)	2.50 (8.82)	33.91 (36.54)	−13.13 (−17.61)	60.52 (64.87)
Willow Creek	WC-VDA (SC-VDA)	5.38 (19.88)	52.46 (57.23)	19.12 (25.76)	56.59 (73.18)
Desert	WC-VDA (SC-VDA)	7.51 (12.69)	33.30 (49.64)	0.43 (0.97)	10.04 (12.71)

a. 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⁻², which is 10.16% lower than the RMSE of 66.24 W·m⁻² from SC-VDA. For latent heat flux, the six-site-averaged RMSE is 75.06 W·m⁻² for the SC-VDA approach, and 67.44 W·m⁻² 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⁻² and 45.32 W·m⁻², 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.

4.5. Model Error Analysis

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).

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⁻² to 30 W·m⁻² 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⁻². 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).

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⁻² (−30 W·m⁻²), and 1 m·s⁻¹ (−1 m·s⁻¹) and standard deviation of 2 K, 1 K, 0.02, 10 W·m⁻², and 0.2 m s⁻¹ 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.

5. 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⁻² and 67.44 W·m⁻², 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.