A Critical Evaluation on the Role of Aerodynamic and Canopy–Surface Conductance Parameterization in SEB and SVAT Models for Simulating Evapotranspiration: A Case Study in the Upper Biebrza National Park Wetland in Poland

Evapotranspiration (ET) estimation through the surface energy balance (SEB) and soil-vegetation-atmosphere-transfer (SVAT) models are uncertain due to the empirical parameterizations of the aerodynamic and canopy-substrate conductances (gA and gS) for heat and water vapor transfers. This study critically assessed the impact of conductance parameterizations on ET simulation using three structurally different SEB and SVAT models for an ecologically important North-Eastern European wetland, Upper Biebrza National Park (UBNP) in two consecutive years 2015 and 2016. A pronounced ET underestimation (mean bias −0.48 to −0.68 mm day−1) in SEBS (Surface Energy Balance System) was associated with an overestimation of gA due to uncertain parameterization of momentum roughness length and bare soil’s excess resistance to heat transfer (kB−1) under low vegetation cover. The systematic ET overestimation (0.65–0.80 mm day−1) in SCOPE (Soil Canopy Observation, Photochemistry and Energy fluxes) was attributed to the overestimation of both the conductances. Conductance parameterizations in SEBS and SCOPE appeared to be very sensitive to the general ecohydrological conditions, with a tendency of overestimating gA (gS) under humid (arid) conditions. Low ET bias in the analytical STIC (Surface Temperature Initiated Closure) model as compared to SEBS/SCOPE indicated the critical need for calibration-free conductance parameterizations for improved ET estimation.


Introduction
Evapotranspiration (ET) (or latent heat flux, λE) monitoring in the wetlands are important for evaluating the hydrological budget, estimating groundwater discharge, quantifying seasonal water availability, and tracking as well as anticipating drought conditions [1]. Among the hydrological 1.
What is the performance of the three structurally different SEB/SVAT models in the UBNP wetland when simulated with high temporal frequency measurements? 2.
What are the effects of aerodynamic and surface conductances parameterizations and associated state variables in determining the model errors with respect to ET? 3. To what extent the ET modelling errors and conductance parameterizations are impacted by a range of environmental and ecohydrological conditions?
A description of methods including models, study sites, dataset, and data processing is given in Section 2, followed by the results in Section 3. An extended discussion of the results and potential of the SEB and SVAT models for evapotranspiration monitoring in the wetlands is elaborated in Sections 4 and 5, respectively. Symbols used for variables in this study are listed in the in Table S1 (Supplementary Materials).

Model Description
In this study, we used two structurally different SEB models and one SVAT model. The two SEB models are the Surface Energy Balance System (SEBS) [19] and Surface Temperature Initiated Closure (STIC) (version STIC1.2) [56,57], and the SVAT model is the Soil-Canopy-Observation of Photosynthesis and the Energy balance (SCOPE) (version SCOPE1.7) [58]. Among the three models, two are based on bottom-up scaling approach to resolve the surface energy balance components (SEBS and SCOPE1.7) and the other model applies top-down approach (STIC1.2) to estimate the SEB components. The fundamental difference of the three models and their input requirements are detailed  Table 1 and a detailed description of the governing equations are given in Supplementary Materials. A list of symbols and units are presented in Supplementary Materials Table S1. SCOPE1.7 being a SVAT model, tuning of certain input variables was necessary in order to reflect the seasonal variation of canopy biophysical and biochemical attributes. We implemented the Radiative Transfer Model optical (RTMo) model that is incorporated into the SCOPE model framework. The RTMo is composed of an extended versions of the SAIL [59] model and PROSPECT45 [60] model, now referred to as Fluspect. First, an ensemble of the canopy spectra was simulated which was followed by inversion of the RTMo to obtain the best match between the simulated and measured (Landsat 8) spectral reflectance through an iterative optimization approach [61]. An example of selected simulation and corresponding solutions is shown in Figure 1. A detailed overview of the SCOPE1.7 input parameters and variables used in the current study are found in the Table S2 of Supplementary Materials. In the SEBS model, estimation of z 0M is critical. In this study we used NDVI (Normalized Difference Vegetation Index) to compute z 0M [62]. Using five (March, April, June, July, and August) Landsat-8 Near Infra-red (NIR) and Red reflectance spectral bands, we computed NDVI as NDVI = (NIR − Red)/(NIR + Red). We proceeded to interpolate the NDVI into 8-day composites and assigned each day to the nearest distance based on the observed trend. and measured (Landsat 8) spectral reflectance through an iterative optimization approach [61]. An example of selected simulation and corresponding solutions is shown in Figure 1. A detailed overview of the SCOPE1.7 input parameters and variables used in the current study are found in the Table S2 of Supplementary Materials. In the SEBS model, estimation of z0M is critical. In this study we used NDVI (Normalized Difference Vegetation Index) to compute z0M [62]. Using five (March, April, June, July, and August) Landsat-8 Near Infra-red (NIR) and Red reflectance spectral bands, we computed NDVI as NDVI = (NIR − Red)/(NIR + Red). We proceeded to interpolate the NDVI into 8-day composites and assigned each day to the nearest distance based on the observed trend.

Uncertainties in gA and gS Parameterizations in SEB and SVAT
Surface energy balance simulation models use TS to constrain the energy-water fluxes [16,28,63]. The central aspect of contemporary SEB simulation is based on estimating gA and sensible heat flux (H), while solving λE as a residual SEB component. Estimation of gA relies on the Monin-Obukhov Similarity Theory (MOST) or Richardson Number (Ri) criteria, and gA estimates are subjected to uncertainties due to (a) empirical approximation of the displacement height (d0) either as a simple function of vegetation height or leaf area index [64,65], (b) uncertainties in estimating the roughness length of momentum transfer (z0M) and associated surface geometry [66], (c) challenges in accommodating the aerodynamic versus surface temperature inequalities, and (d) complexities in kB −1 parameterization [67,68] to compensate for the differences in the scalar roughness lengths of heat (z0H) and momentum (z0M) transfer. The parameters d0,

Uncertainties in g A and g S Parameterizations in SEB and SVAT
Surface energy balance simulation models use T S to constrain the energy-water fluxes [16,28,63]. The central aspect of contemporary SEB simulation is based on estimating g A and sensible heat flux (H), while solving λE as a residual SEB component. Estimation of g A relies on the Monin-Obukhov Similarity Theory (MOST) or Richardson Number (Ri) criteria, and g A estimates are subjected to uncertainties due to (a) empirical approximation of the displacement height (d 0 ) either as a simple function of vegetation height or leaf area index [64,65], (b) uncertainties in estimating the roughness length of momentum transfer (z 0M ) and associated surface geometry [66], (c) challenges in accommodating the aerodynamic versus surface temperature inequalities, and (d) complexities in kB −1 parameterization [67,68] to compensate for the differences in the scalar roughness lengths of heat (z 0H ) and momentum (z 0M ) transfer. The parameters d 0 , z 0M , z 0H , and kB −1 are the highly variable components in SEB models, and the commonly used empirical response functions of these components to characterize g A have an uncertain transferability in space and time [69,70]. An extended description of the impacts of ambiguous parameterization of d 0 , z 0M , kB −1 on ET estimates is detailed by [67,71,72].
SVAT models represent electrical resistance network analogy through mathematical representation of the detailed physical mechanisms that controls the energy and mass transfers in the soil/vegetation/atmosphere continuum [73]. Surface energy balance is an integral part of the SVAT models and are mainly used to simulate the coupled exchanges of energy-water-carbon to the ecosystem level. In addition to the aerodynamic conductance, SVAT models use different formulations of stomatal (or canopy-surface) conductance and plant hydraulics. However, there is no common consensus on which g S model matches best with the observed surface energy balance fluxes [74]. Therefore, there is a critical need to assess the impact of aerodynamic and canopy-surface conductance representation on the performance of SEB and SVAT models.

Eddy Covariance Estimation of g A and g S for Model Evaluation
Due to the lack of direct g A measurements, a rigorous evaluation of g A from the models is challenging. However, an indirect evaluation of g A retrievals from the three models was performed using the micrometeorological measurements from the EC towers. By using the measured friction velocity (u*) and wind speed (u) at the EC towers and using the equation of [75] we estimated g A as the sum of turbulent conductance and canopy (quasi-laminar) boundary-layer conductance as follows. g A (EC tower) = [(u/u* 2 ) + (2/ku* 2 )(Sc/Pr) 0.67 ] −1 where k is von Karman's constant, 0.4; Sc is the Schmidt Number; Pr is the Prandtl Number and their ratio is generally considered to be unity. Here the conductances of momentum, sensible and latent heat fluxes are assumed to be identical [57].
Similar to g A , a direct evaluation modeled g S could not be performed, as independent ecosystem-scale g S observations are not possible with current measurement techniques. Assuming u*-based g A as the baseline aerodynamic conductance, we estimated g S by inverting the Penman-Monteith equation [76] to evaluate the modeled g S . Here u*-based g A was used in conjunction with the net available energy, λE, air temperature, and humidity measurements from the EC towers [57] as follows.
Here γ is the psychrometric constant, s is the slope of saturation vapor pressure versus air temperature relationship, R N − G is the net available energy, ρ and c p are the density and specific heat of air, respectively.
Among the three models, SEBS does not estimate g S . However, from the estimated g A and λE (of SEBS), an inverse estimation of SEBS-based g S was also performed for comparing with the EC g S estimates in the framework of Penman-Monteith equation.

Study Site: Ecological and Hydrological Significance of UBNP Wetland
The Biebrza is the largest National Park in Poland, and it covers an area of over 59,000 hectares [45]. The area experiences sub-continental/sub-boreal climate with an annual average temperature of 6.8 • C, precipitation ranging between 550 to 700 mm year −1 and evapotranspiration between 460 and 480 mm year −1 [40]. It is located 230 km North East of Warsaw being serviced by the Biebrza River which is a right sided tributary of the Narew River. The UBNP wetland is situated within the upper Biebrza river catchment (22 • 30 -23 • 60 E, 53 • 30 -53 • 75 N). While the lower part of Biebrza experiences a sequence of flood events often after snowmelt in early spring, the upper part of the Biebrza river basin rarely experiences flooding. Most part of the surface runoff takes place within the drainage network [46]. During the dry periods, the wetland is groundwater fed. Most part of the surface runoff in UBNP takes place within the drainage network [46].
Upper Biebrza is the region of large variety of wetland habitats and typically driven by local hydrology [44,79]. The substrate is mainly composed of natural peat supporting highly productive tall sedge, reeds and grass [44]. The ecosystem offers habitat and breeding ground for numerous species of which up to 21 species have been listed as threatened in Europe including but not limited to snipe, aquatic warbler, elk, otters and European beaver [44]. Biebrza valley can be considered as a reference for other similar ecosystems in Western and Central Europe since 96% of the present plant species also belong to the Dutch flora [79,80]. Currently, the UBNP wetland is listed under the RAMSAR convention as one of the most important wetlands worldwide designated under protected river valley [45].

Datasets
In this analysis, we have used eddy-covariance (EC) observations from the Rogozynk EC tower site for the two consecutive years of 2015 and 2016 from spring (March) to mid-summer (July) period. These periods were selected because 2015 and 2016 experienced contrasting precipitation amounts according to the meteorological bureau report and also due to the data availability. Whereas summer of 2015 was actually characterized by very low precipitation (amounting to 50% of the seasonal average computed over the period 1971-2000); the precipitation amounts in spring/summer 2016 were close to the normal trend [7] (Figure 2).
The EC data were available at half-hourly temporal resolution from 5th March to 22nd July for both the years. Data used for this analysis included time series of surface energy balance fluxes (R N , λE, H, and G), shortwave and longwave radiation components (RS in , RS out , RL in , RL out ), and hydrometeorological variables (e.g., T A , R H , u, u*, θ, and P). Assuming nighttime ET to be negligible in the wetland [10,81], daily ET (in mm) was computed by integrating half-hourly λE from sunrise to sunset (corresponds to positive magnitude of RS in ). However, it is also important to mention that wetlands in the arid and semi-arid regions could show substantial amount of ET during night and careful daily ET averaging is needed in such cases [82]. Weekly ET (in mm) was computed by summing daily ET. We did not perform any gap filling, which implies that missing observed or estimated sub-daily or daily ET values were not included in the computation. The energy balance closure of the EC site was 92% and the surface energy balance was closed using the Bowen ratio energy balance method as described in [27,56,83].
Surface emissivity (ε s ) was assumed to be 0.98 and Ts was computed by inverting the longwave radiation measurements as follows.
Albedo computation was based on the ratio between the downwelling and upwelling shortwave radiation α 0 = RS out RS in . The distribution of precipitation (P) and soil moisture (θ) for the year 2015 and 2016 is also shown in the primary y-axis (P) and secondary y-axis (θ), respectively.

Model Evaluation and Comparison
ET estimates from the respective models were compared with measured ET. The relative performance of each model was evaluated using the following statistical indices namely the coefficient of determination (R 2 ), the root mean square error (RMSE), the mean absolute percentage deviation (MAPD), the regression slope, and the regression intercept, respectively.
where, n is the number of data points; oi and pi are observed and estimated λE and ET; � is the mean of observations.

Relationship between ET Modeling Errors, Conductances, and Ecohydrological Factors
To understand the effects of the aerodynamic and canopy-surface conductances and associated variables on λE simulation, we analyzed the relationship between the residual λE errors of individual models with the simulated conductances and associated variables. The degree of relationship was assessed through the correlation (r) values and the significance of correlation was tested by the pvalue. The p-value signifies the probability that we would have found the current results in residual The distribution of precipitation (P) and soil moisture (θ) for the year 2015 and 2016 is also shown in the primary y-axis (P) and secondary y-axis (θ), respectively.

Model Evaluation and Comparison
ET estimates from the respective models were compared with measured ET. The relative performance of each model was evaluated using the following statistical indices namely the coefficient of determination (R 2 ), the root mean square error (RMSE), the mean absolute percentage deviation (MAPD), the regression slope, and the regression intercept, respectively.
where, n is the number of data points; o i and p i are observed and estimated λE and ET; o i is the mean of observations.

Relationship between ET Modeling Errors, Conductances, and Ecohydrological Factors
To understand the effects of the aerodynamic and canopy-surface conductances and associated variables on λE simulation, we analyzed the relationship between the residual λE errors of individual models with the simulated conductances and associated variables. The degree of relationship was assessed through the correlation (r) values and the significance of correlation was tested by the p-value. The p-value signifies the probability that we would have found the current results in residual error analysis if the correlation coefficient is zero (i.e., null hypothesis). If the p-value is less than 5% (i.e., p-value < 0.05), the correlation coefficient is considered as statistically significant.
To examine the influence of ecohydrological conditions on ET prediction errors and biophysical conductances, we further investigated the patterns of mean bias (MB) in weekly ET in comparison to the weekly soil moisture (θ) and climatic aridity or dryness (i.e., annual E P /P) [84] which are considered to represent the general ecohydrological characteristics of ecosystems. UBNP wetland ecosystem generally has high magnitude of soil moisture and low E P /P (low evaporative demand and high precipitation). Therefore, an independent assessment of the effects of E P /P and θ on the predictive capacity of the SEB and SVAT models is crucial. Results of the correlation analysis between MB in weekly ET with weekly θ and weekly E P /P are presented in Section 3.3, and discussions on the effects of θ and weekly E P /P on the biophysical conductance simulations are elaborated in Section 5.

Evaluation of g A and g S Estimates from Models
An evaluation of ecosystem-scale g A is illustrated in Figure 3a combining data from both the years. While the estimated g A values ranged between 0-0.06 m s −1 for STIC1.2 and SEBS, it was 0-0.6 m s −1 for SCOPE1.7, with R 2 ranged between 0.13 to 0.80 between the tower-observed g A and modelled g A . Statistical comparisons between the models revealed that although SCOPE1.7 had a good R 2 , but the magnitude of g A from SCOPE1.7 is ten times higher with respect to the tower observations. As seen in Figure 3a, there appears to be a strong systematic bias in SCOPE1.7 g A which led to very high BIAS and MAPD. Overall, the mean bias and MAPD in modeled g A varied between 0.003-0.12 m s −1 and 81 to 93%, respectively.
Canopy-scale evaluation of half-hourly g S is presented in Figure 3b and the estimated values ranged between 0-0.06 m s −1 for STIC1.2, 0-0.2 m s −1 for SCOPE1.7, and 0-0.01 m s −1 for SEBS. The magnitude of R 2 varied from 0.15-0.53, with mean bias and MAPD of −0.004-0.1 m s −1 and 73-86%, respectively. Like g A , a systematic overestimation of g S was also found in SCOPE1.7 which is consequently revealed in very high mean bias and MAPD.
Water 2018, 10, x FOR PEER REVIEW 9 of 27 error analysis if the correlation coefficient is zero (i.e., null hypothesis). If the p-value is less than 5% (i.e., p-value < 0.05), the correlation coefficient is considered as statistically significant.
To examine the influence of ecohydrological conditions on ET prediction errors and biophysical conductances, we further investigated the patterns of mean bias (MB) in weekly ET in comparison to the weekly soil moisture (θ) and climatic aridity or dryness (i.e., annual EP/P) [84] which are considered to represent the general ecohydrological characteristics of ecosystems. UBNP wetland ecosystem generally has high magnitude of soil moisture and low EP/P (low evaporative demand and high precipitation). Therefore, an independent assessment of the effects of EP/P and θ on the predictive capacity of the SEB and SVAT models is crucial. Results of the correlation analysis between MB in weekly ET with weekly θ and weekly EP/P are presented in Section 3.3, and discussions on the effects of θ and weekly EP/P on the biophysical conductance simulations are elaborated in Section 5.

Evaluation of gA and gS Estimates from Models
An evaluation of ecosystem-scale gA is illustrated in Figure 3a combining data from both the years. While the estimated gA values ranged between 0-0.06 m s −1 for STIC1.2 and SEBS, it was 0-0.6 m s −1 for SCOPE1.7, with R 2 ranged between 0.13 to 0.80 between the tower-observed gA and modelled gA. Statistical comparisons between the models revealed that although SCOPE1.7 had a good R 2 , but the magnitude of gA from SCOPE1.7 is ten times higher with respect to the tower observations. As seen in Figure 3a, there appears to be a strong systematic bias in SCOPE1.7 gA which led to very high BIAS and MAPD. Overall, the mean bias and MAPD in modeled gA varied between 0.003-0.12 m s −1 and 81 to 93%, respectively.
Canopy-scale evaluation of half-hourly gS is presented in Figure 3b

Evaluation of λE (ET) Estimates from Models
Given the structural differences between the three models in computing the SEB components, this section describes the performance of the models through evaluating the predictive capacity of half-hourly λE for spring, summer, and entire growth period for the years 2015 and 2016. The statistical evaluation of daily ET (integration of half-hourly λE for the daytime) is also presented afterwards.
Scatterplots of predicted versus observed λE revealed SCOPE1.7 and STIC1.  (Table 2) revealed SEBS to produce the highest RMSE (62-75 W m −2 ) and MAPD (24-46%) (for the entire growing season as well as summer and spring), followed by SCOPE1.7 (37-52 W m −2 , 22-32%) and STIC1.2 (29-31 W m −2 , 18-19%). As evident from the slope of the linear regressions and mean bias (Table 2), while SEBS and STIC1.2 underestimated λE; SCOPE1.7 had overestimated λE in both years, regardless of the season. (a) Comparison between model derived g A estimates with an estimated aerodynamic conductance based on friction velocity (u*) and wind speed (u) according to [75], (b) Comparison between model derived g S estimates with respect to g S computed by inverting the penman-Monteith model where g A estimates from EC tower were used as aerodynamic input in conjunction with tower measurements of λE, radiation and meteorological variables.

Evaluation of λE (ET) Estimates from Models
Given the structural differences between the three models in computing the SEB components, this section describes the performance of the models through evaluating the predictive capacity of half-hourly λE for spring, summer, and entire growth period for the years 2015 and 2016. The statistical evaluation of daily ET (integration of half-hourly λE for the daytime) is also presented afterwards.
Scatterplots of predicted versus observed λE revealed SCOPE1.7 and STIC1.  (Table 2) revealed SEBS to produce the highest RMSE (62-75 W m −2 ) and MAPD (24-46%) (for the entire growing season as well as summer and spring), followed by SCOPE1.7 (37-52 W m −2 , 22-32%) and STIC1.2 (29-31 W m −2 , 18-19%). As evident from the slope of the linear regressions and mean bias (    Time-series progression of daily ET (mm day −1 ) ( Figure 5) showed large variability in day-to-day ET during the summer months (June-July) (1-7 mm) which is reasonably captured by all the models. Maximum ET was found to occur during June-July (7 mm) and minimum ET was observed during spring (1-3 mm). Among the three models, SEBS revealed consistent underestimation which was more pronounced during the spring. However, SEBS captured the daily ET trend slightly better than the other two models during the dry-down phase of late summer in the year 2016, whereas SCOPE1.7 (STIC1.2) revealed substantial (little) overestimation tendency during this period. Time-series progression of daily ET (mm day −1 ) ( Figure 5) showed large variability in day-today ET during the summer months (June-July) (1-7 mm) which is reasonably captured by all the models. Maximum ET was found to occur during June-July (7 mm) and minimum ET was observed during spring (1-3 mm). Among the three models, SEBS revealed consistent underestimation which was more pronounced during the spring. However, SEBS captured the daily ET trend slightly better than the other two models during the dry-down phase of late summer in the year 2016, whereas SCOPE1.7 (STIC1.2) revealed substantial (little) overestimation tendency during this period.  (Table 3), with a consistent overestimation (underestimation) of daily ET from SCOPE1.7 (STIC1.2). While the MAPD of STIC1.2 varied between 16-21%, it was 38-44% for SCOPE1.7, and both models showed relatively high MAPD in 2016. For SEBS, the RMSE and MAPD in daily ET was 0.74-0.95 mm and 33-40%, respectively (Table 3).   (Table 3), with a consistent overestimation (underestimation) of daily ET from SCOPE1.7 (STIC1.2). While the MAPD of STIC1.2 varied between 16-21%, it was 38-44% for SCOPE1.7, and both models showed relatively high MAPD in 2016. For SEBS, the RMSE and MAPD in daily ET was 0.74-0.95 mm and 33-40%, respectively (Table 3).

Effects of Biophysical Conductance Parameterization on Residual Error of the Models
The effects of biophysical parameterization in the propagation of model errors are presented in the box and violin plots ( Figure 6) which showed the residual λE error (∆ λE ) (= predicted − observed) for different classes of simulated conductances, T 0 and T c for STIC1.2 and SCOPE1.7. Since SEBS does not simulate the surface conductance, the model errors for SEBS are assessed by relating ∆ λE with the aerodynamic conductance, kB −1 and z 0M .
for different classes of simulated conductances, T0 and Tc for STIC1.2 and SCOPE1.7. Since SEBS does not simulate the surface conductance, the model errors for SEBS are assessed by relating ΔλE with the aerodynamic conductance, kB −1 and z0M. Figure 6 revealed that ΔλE from SCOPE1.7 (ΔλESCOPE1.7 hereafter) had a significantly positive relationship with the simulated canopy temperature (Tc) (r = 0.46, p-value < 0.05), and a systematic overestimation of λE (ΔλESCOPE1.7 increased exponentially) with Tc was evident when Tc increased from 10 to 30 °C (Figure 6a). However, ΔλESCOPE1.7 was found to be moderately correlated (r = 0.22 and r = 0.32 for gA and gS respectively) with the two biophysical conductances (Figure 6a). The residual λE error from SEBS (ΔλESEBS hereafter) was inversely related to gA and kB −1 (Figure 6b). Here, a systematic underestimation of λE was evident with increasing gA and kB −1 (Figure 6b) with a correlation of −0.55 and −0.27 (p-value < 0.05, significant), respectively. The mean residual λE error from STIC1.2 (ΔλESTIC1.2 hereafter) showed an increasing pattern with an increase in gA and gS (Figure 6c) having correlation of 0.35 and 0.27 (p-value < 0.05), respectively. However, ΔλESTIC1.2 appeared to be heteroscedastic with an increase in T0, which signifies unequal variability of ΔλESTIC1.2 as the value of T0 increases (Figure 6c).  Figure 6 revealed that ∆ λE from SCOPE1.7 (∆ λE SCOPE1.7 hereafter) had a significantly positive relationship with the simulated canopy temperature (T c ) (r = 0.46, p-value < 0.05), and a systematic overestimation of λE (∆ λE SCOPE1.7 increased exponentially) with T c was evident when T c increased from 10 to 30 • C (Figure 6a). However, ∆ λE SCOPE1.7 was found to be moderately correlated (r = 0.22 and r = 0.32 for g A and g S respectively) with the two biophysical conductances (Figure 6a). The residual λE error from SEBS (∆ λE SEBS hereafter) was inversely related to g A and kB −1 (Figure 6b). Here, a systematic underestimation of λE was evident with increasing g A and kB −1 (Figure 6b) with a correlation of −0.55 and −0.27 (p-value < 0.05, significant), respectively. The mean residual λE error from STIC1.2 (∆ λE STIC1.2 hereafter) showed an increasing pattern with an increase in g A and g S (Figure 6c) having correlation of 0.35 and 0.27 (p-value < 0.05), respectively. However, ∆ λE STIC1.2 appeared to be heteroscedastic with an increase in T 0 , which signifies unequal variability of ∆ λE STIC1.2 as the value of T 0 increases (Figure 6c).

Effects of Environmental and Ecohydrological Factors on the Model Performances
Given the variable performance of the models to homogeneous inputs, this section described the role of environmental (meteorological and radiation) and surface (T S and vegetation index) variables in determining the residual λE error (∆ λE = predicted − observed) of the individual models through principal component regression (PCR) analysis. This will be followed by assessing the role of general ecohydrological conditions on the weekly ET bias from the individual models.

Effects of Environmental and Ecohydrological Factors on the Model Performances
Given the variable performance of the models to homogeneous inputs, this section described the role of environmental (meteorological and radiation) and surface (TS and vegetation index) variables in determining the residual λE error (ΔλE = predicted − observed) of the individual models through principal component regression (PCR) analysis. This will be followed by assessing the role of general ecohydrological conditions on the weekly ET bias from the individual models.   (Figure 7). For all the models, the first principal components (PC1), which is on the horizontal axis, had positive coefficients for T S , T A , φ, and u (exception SCOPE1.7 for year 2016 where u had negative coefficient). Therefore, vectors are directed into the right half of the plot. It is also surprising to see that despite STIC1.2 was not driven by wind speed, there is an apparent relationship between ∆ λE and u in STIC1.2 due to a strong relationship between u and T S . For all the models, maximum PC1 loading was found for T S and T A followed by φ (Figure 7 (Figure 8b); whereas ET underestimation tendency in SEBS was reduced with increasing E p /P ratio (Figure 8b). Contrarily, STIC1.2 showed an overestimation tendency in weekly ET (positive bias) (2-6 mm/week) for low E p /P which progressively declined with increasing E p /P ratio, with a correlation of 0.21 (p-value < 0.05).  (Figure 7). For all the models, the first principal components (PC1), which is on the horizontal axis, had positive coefficients for TS, TA, ϕ, and u (exception SCOPE1.7 for year 2016 where u had negative coefficient). Therefore, vectors are directed into the right half of the plot. It is also surprising to see that despite STIC1.2 was not driven by wind speed, there is an apparent relationship between ΔλE and u in STIC1.2 due to a strong relationship between u and TS. For all the models, maximum PC1 loading was found for TS and TA followed by ϕ (Figure 7

Effects of Model Structure and Biophysical Parameterizations on Residual λE or ET Errors
Overall, STIC1.2 revealed marginal underestimation (negative mean bias) of half-hourly and daily λE and ET during both years (2015 and 2016), whereas SCOPE1.7 showed a tendency of overestimation and SEBS showed a tendency of underestimation especially in the beginning of the spring season. As evident from the residual error analysis, the underestimation tendency of STIC1.2 is due to the underestimation of λE for the low values g A and g S (Figure 6). Uncertainty in the relationship between thermal infrared temperature (T S ) and aggregated moisture availability (M) in STIC1.2 could be a considerable source of underestimation of the conductances which is reflected in the simulation of λE (and ET) through STIC1.2. In STIC1.2, M is modeled as a fraction of the dewpoint temperature difference between the evaporating front and atmosphere (T 0D -T D ) and of infrared temperature-dewpoint differences between surface to atmosphere (T S -T D ), weighted by two different slopes of saturation vapor pressure temperature relationships (s 1 and s 2 ; equation S26 in [56], [M = s 1 (T 0D − T D )/s 2 (T S − T D )]. The estimation of T 0D plays a critical role in the water unlimited wetland ecosystem because estimation of T 0D further requires sound estimation of s 1 . From the definition, s 1 is the slope of the saturation vapor pressure versus temperature curve between the evaporating front and air [s 1 = (e 0 − e A )/(T 0D − T D )]. However, in the wetland, e 0 tends to approach saturation vapor pressure of wet surface and s 1 tend to be the slope of the wet surface to air saturation-vapor pressure versus temperature. In the present case, the estimates of s 1 as a function of air dewpoint temperature (T D ) tend to be lower than the possible s 1 -limits for the water-unlimited surfaces, which is likely to introduce underestimation errors in T 0D [56,57]. Underestimation of s 1 and T 0D would also lead to an underestimation of M (through the denominator in equation S26 in [56], thus leading to an underestimation of the conductances and λE. As demonstrated by [56], the ratio of g S /g A increases with increasing M and the sensitivity of g S to M is substantially higher as compared to the sensitivity of g A to M. An underestimation M would lead to an underestimation of g S , and consequently an overestimation of the denominator [i.e., s + γ(1 + g A /g S )] in the Penman-Monteith equation (because g A /g S in the denominator is overestimated). Thus underestimation errors could be introduced in λE simulation in STIC1.2 for high soil moisture conditions.
In SEBS, the principal source of errors originated from consistent overestimation of g A (as seen in Figure 3a). The impact of g A and associated roughness length parameterization (z 0M and kB −1 ) is therefore evident in the residual λE error in SEBS ( Figure 6). In SEBS, z 0M is estimated as a function of the leaf area index and vegetation index [77], and low values of both the variables in the start of the growing season would tend to an underestimation of z 0M , which would lead to an overestimation of kB −1 and g A , and an underestimation of λE (as see in Figure 6). Despite the strong dependence of kB −1 on T S , radiation, and meteorological variables [64]; no common consensus on a physically-based model for both z 0M and kB −1 is available [67,68]. Empirical parameterization of z 0M and ±50% uncertainties in z 0M can also lead to 25% errors in g A estimation [18,64,77], which would lead to more than 30% uncertainty in ET estimates. This is also evident from the semi-exponential pattern between z 0M and mean ∆ λE SEBS (Figure 6b) that showed a positive correlation (r = 0.26, p-value < 0.05). Major λE differences for kB −1 range greater than 6 (typical range for the wet/saturated surfaces) (Figure 6b) indicates uncertainties in kB −1 parameterization for simulating λE in the water-unlimited extremes. This is further reflected in the weekly ET bias versus soil moisture and E p /P scatters ( Figure 8) where a negative bias was found for low values of soil moisture and E p /P ratio, conditions that presumably exist in the start of the growing season.
Consistent underestimation of λE at the start of the season (when the latent heat fluxes were generally low) is also associated with the structural uncertainties in SEBS in estimating the relative evaporation (Λ r ) which was estimated by scaling the actual sensible heat flux (H) with the sensible heat fluxes for the driest (H dry ) and wettest (H wet ) conditions [Λ r = 1 − (H − H wet )/(H dry − H wet ); and H dry = Rn − G] [19]. In the start of the growing season (spring), when λE is low, any condition that produces H ≈ H dry would tend to simulate substantially low relative evaporation (Λ r ≈ 0), and λE will be consequently underestimated (Figure 9a). As seen in Figure 9b, the residual daily ET error (∆ ET ) in SEBS had a rather linear relationship with daily Λ r (r = 0.35-0.37, p-value < 0.05) and a systematic underestimation in daily ET (up to − mm day −1 ) was revealed for 0 < Λ r < 0.4.  [19]. In the start of the growing season (spring), when λE is low, any condition that produces H ≈ Hdry would tend to simulate substantially low relative evaporation (Λr ≈ 0), and λE will be consequently underestimated (Figure 9a). As seen in Figure 9b, the residual daily ET error (ΔET) in SEBS had a rather linear relationship with daily Λr (r = 0.35-0.37, p-value < 0.05) and a systematic underestimation in daily ET (up to − mm day −1 ) was revealed for 0 < Λr < 0.4. Uncertainties in the estimation of wet limits of H (i.e., Hwet) and associated aerodynamic conductance in the wet limits (gAwet) are also responsible for the systematic underestimation of λE the relationship between residual error in daily ET (∆ ET ) from SEBS with daily Λ r , which further confirms the systematic underestimation in ET was associated with estimation errors in Λ r , (c) Scatterplots showing the relationship between ∆ ET from SEBS with simulated sensible heat flux for the wet limit (H wet ), (d) Scatterplots showing the relationship between ∆ ET from SEBS with simulated aerodynamic conductance for the wet limit (g Awet ). Scatterplot in the inset shows the relationship between H wet and g Awet .
Uncertainties in the estimation of wet limits of H (i.e., H wet ) and associated aerodynamic conductance in the wet limits (g Awet ) are also responsible for the systematic underestimation of λE and ET in SEBS (Figure 9c,d). A consistent negative bias in daily ET was evident (−1 to −2 mm day −1 ) for high magnitude of H wet (Figure 9c) [r = (−0.29)-(−0.31), p-value < 0.05]. This was mainly due to the inverse exponential relationship between H wet and g Awet (inset of Figure 9d) which was consequently propagated in ET as evident from the scatter between ∆ ET versus g Awet [r = (−0.15)-(−0.19), p-value < 0.05] (Figure 9d).
Like SEBS, g A estimation in SCOPE1.7 is also based on the Monin-Obukhov Similarity Theory and similar empirical sub-models for the roughness lengths (z 0M and z 0H ) [58]. Therefore, the errors in λE and ET in SCOPE1.7 due to large overestimation of g A appear to be similar to SEBS. However, an additional source of uncertainty in SCOPE1.7 is due to a consistent overestimation of g S (Figure 3b) which was consequently propagated into overestimation of λE. In SCOPE1.7, g S parameterization is based on the Ball-Woodrow-Berry (BWB) g S -photosynthesis model [85]. BWB model was originally developed at the leaf-scale and no universally agreed scaling method is available to extrapolate this model for the ecosystem scale [86]. Photosynthesis simulation in SCOPE1.7 depends on the calibration of V cmax (i.e., maximum carboxylation capacity) [87], parameter for dark respiration, Extinction coefficient for vertical profile of V cmax . Uncertainty in photosynthesis simulation in SCOPE1.7 would lead to erroneous g S (due to the g S -photosynthesis feedback in the model) and λE. The major problem of using BWB or BWB-Leuning class of models is that they are valid only for saturated soil water conditions and cannot accommodate the soil drying process or when the soil drying is coupled with high atmospheric vapor pressure deficit. However, there are new generation g S models that could be used to predict stomatal conductance during soil dry-down when accounting for the soil-xylem hydraulics and detailed plant physiological attributes [88][89][90][91]. Assessing the photosynthesis related uncertainty in SCOPE1.7 is beyond the scope of this study, but, the consistent overestimation tendency in ET (Figure 8b) indicates the uncertainty in g S and g A simulation ( Figure 3) under high atmospheric aridity (high E P and low P) (Figure 10 below) to be one of the main sources of errors in SCOPE1.7.

Effects of Ecohydrological Conditions on Conductance Estimation and Implication on Model Performances
To probe into the reasons for the weekly bias in ET versus θ and E p /P ratio (as seen in Figure 8), the effects of ecohydrological conditions on the two biophysical conductance retrievals are presented in Figure 10, which showed the scatters between weekly g S /g A ratio from STIC1.2 and SCOPE1.7 with weekly θ and E p /P ratio. For STIC1.2, a weak and nearly invariant relationship (r = 0.02, p-value > 0.05) was found between g S /g A ratio versus θ (Figure 10a). Since the soil moisture was predominantly high (to the level of saturation) in the UBNP wetland, the conductances retrieved through STIC1.2 appeared to be unaffected due the underlying soil moisture conditions. However, g S /g A ratio from STIC1.2 progressively diminished with increasing E p /P (Figure 10b) with a significant correlation (r = 0.64, p-value < 0.05), which means low g S as compared to g A with increasing atmospheric aridity and evaporative demand. This further emphasizes the uncertainties due to surface moisture characterization in STIC1.2 especially the assumptions associated with the slope of temperature-vapor pressure relationship as described in Section 6.1. For SCOPE1.7, although g S /g A ratio varied between 0.8-1.0 for the entire range of soil moisture, but, there was a significantly positive relationship (r = 0.41, p-value < 0.05) between SCOPE1.7-derived g S /g A ratio with E p /P. This signifies relatively high values of g S (as compared to g A ) with increasing evaporative demand (and atmospheric aridity). The tendency of simulating high g S as compared to g A with increasing evaporative demand (i.e., high E P ) was eventually propagated into an overestimation of λE (and ET) in SCOPE1.7. With the progress of the growing season from spring onwards, an increase in the atmospheric vapor pressure deficit, air temperature and radiative load led to an increase in the evaporative demand of the atmosphere where λE overestimation (due to g S ) was predominant. Water 2018, 10, x FOR PEER REVIEW 20 of 27 (a) (b) Figure 10. (a) Scatterplots of weekly gS/gA ratio simulated from STIC1.2 and SCOPE1.7 (kB −1 for SEBS) versus weekly soil moisture (θ), and (b) weekly climatic aridity index (EP/P ratio). Since SEBS does not simulate gA, the behaviour of kB −1 was assessed with respect to the same ecohydrological conditions (θ and EP/P ratio). This shows that gS/gA ratio simulated through STIC1.2 and SCOPE1.7 was not affected by the soil moisture conditions, but it was strongly affected by EP/P ratio. In SEBS, kB −1 was found to be affected by both θ and EP/P ratio.
As discussed in Section 6.1, errors in the aerodynamic conductance in SEBS appears to be caused due to the empirical parameterizations of z0M and kB −1 that is eventually propagated to z0H [77]. To further understand the role of general ecohydrological conditions on gA simulation and its impact on ET estimation in SEBS (SEBS does not simulate gS), the scatterplots of weekly kB −1 versus weekly soil moisture and Ep/P ratio are also showed in Figure 10. This figure highlights a significant impact of these two ecohydrological factors on the estimation of kB −1 (r = −0.26 to −0.29, p-value < 0.05). The prevailing soil moisture and Ep/P ratio in the UBNP wetland during the start of the spring typically varying between 0.7-0.8 m 3 m −3 and 0-2 and the underestimation tendency of SEBS is maximally noted in this range of soil moisture and Ep/P ratio. Majority of the λE (ET) underestimation was evident when kB −1 values exceeded 6 which is a typical range of kB −1 in the humid ecosystems with low Ep/P ratio (0-2), and the underestimation trend was consistent for EP/P ratio of 0-2. Thus, underestimation of z0M (due to low NDVI in the start of the growing season) in conjunction with an overestimation of kB −1 would lead to an underestimation of z0H [z0H = z0M/exp(kB −1 )]. Underestimation of z0H consequently led to an overestimation of gA and H both for the actual and the wet limits, and a consistent underestimation of λE (ET) was discernible. Previous studies ( [77,92]) also revealed substantial uncertainties in z0M and kB −1 that eventually led to uncertain estimation of ET. Our results demonstrated the critical role of surface to aerodynamic conductance parameterizations on the performance of the SVAT and SEB models.
It is also important to mention that significant differences were noted in the aerodynamic conductances simulated by the three models. The differences between gA from the models were mainly attributed to their structural differences and the nature of the parameterization of gA from all Figure 10. (a) Scatterplots of weekly g S /g A ratio simulated from STIC1.2 and SCOPE1.7 (kB −1 for SEBS) versus weekly soil moisture (θ), and (b) weekly climatic aridity index (E P /P ratio). Since SEBS does not simulate g A , the behaviour of kB −1 was assessed with respect to the same ecohydrological conditions (θ and E P /P ratio). This shows that g S /g A ratio simulated through STIC1.2 and SCOPE1.7 was not affected by the soil moisture conditions, but it was strongly affected by E P /P ratio. In SEBS, kB −1 was found to be affected by both θ and E P /P ratio.
As discussed in Section 6.1, errors in the aerodynamic conductance in SEBS appears to be caused due to the empirical parameterizations of z 0M and kB −1 that is eventually propagated to z 0H [77]. To further understand the role of general ecohydrological conditions on g A simulation and its impact on ET estimation in SEBS (SEBS does not simulate g S ), the scatterplots of weekly kB −1 versus weekly soil moisture and E p /P ratio are also showed in Figure 10. This figure highlights a significant impact of these two ecohydrological factors on the estimation of kB −1 (r = −0.26 to −0.29, p-value < 0.05). The prevailing soil moisture and E p /P ratio in the UBNP wetland during the start of the spring typically varying between 0.7-0.8 m 3 m −3 and 0-2 and the underestimation tendency of SEBS is maximally noted in this range of soil moisture and E p /P ratio. Majority of the λE (ET) underestimation was evident when kB −1 values exceeded 6 which is a typical range of kB −1 in the humid ecosystems with low E p /P ratio (0-2), and the underestimation trend was consistent for E P /P ratio of 0-2. Thus, underestimation of z 0M (due to low NDVI in the start of the growing season) in conjunction with an overestimation of kB −1 would lead to an underestimation of z 0H [z 0H = z 0M /exp(kB −1 )]. Underestimation of z 0H consequently led to an overestimation of g A and H both for the actual and the wet limits, and a consistent underestimation of λE (ET) was discernible. Previous studies ( [77,92]) also revealed substantial uncertainties in z 0M and kB −1 that eventually led to uncertain estimation of ET. Our results demonstrated the critical role of surface to aerodynamic conductance parameterizations on the performance of the SVAT and SEB models.
It is also important to mention that significant differences were noted in the aerodynamic conductances simulated by the three models. The differences between g A from the models were mainly attributed to their structural differences and the nature of the parameterization of g A from all the three models. Therefore, the difference in g A estimates markedly contributed to the different statistical metric between simulated versus observed λE (and ET) from the models. In general, the accuracies in commonly used parametric g A estimates based on wind speed and surface roughness parameters several meters distant from canopy foliage are limited due to the uncertainties related to the attenuation of wind speed close to the vegetation surface [93,94]. The wind speed close to the foliage can be substantially lower than that measured at some reference location above the vegetation canopy [95]. Notwithstanding the inequalities of g A estimated with different methods, inferring the accuracy of the different g A estimates is beyond the scope of the manuscript. However, g A is one of the main anchors in the SVAT and SEB models because it provides feedback to g S [96]. Therefore, ET estimation using SVAT and SEB models are very sensitive to g A and a universally agreed g A parameterization is needed to obtain similar λE results from the models [97,98]. Given the lack of consensus in the community on the "true" g A and from the nature of surface flux validation results (Figure 4), it appears that g A from STIC1.2 tends to be the appropriate aerodynamic conductance that tend to produce the lowest λE error.

Conclusions
This paper demonstrates the role of aerodynamic and canopy-surface (i.e., vegetation-substrate complex) conductance (g A and g S ) parameterizations in estimating ET from SVAT and SEB models. By using high temporal frequency eddy covariance measurements and three different SVAT and SEB models (SCOPE1.7, SEBS, and STIC1.2), we showed the critical role of ecohydrological conditions in influencing the conductance and latent heat flux simulation, and consequently daily ET. Independent validation of the models using observed latent heat flux data from an anomalous and a normal precipitation year (2015 and 2016) from one of the most ecologically important wetlands (Upper Biebrza National Park, Poland) (UBNP) led us to the following conclusions.
(a) Notable differences were found out in the g A and g S estimates from the three models. While SCOPE1.7 revealed substantial overestimation of both g A and g S with respect to the EC tower estimates, STIC1.2 derived g A and g S were within the range of EC tower estimates. SEBS revealed a consistent overestimation of g A during the start of the growing season in spring, and g A estimates were is good agreement with the EC tower during the active vegetative phase in summer. (b) All the models explained significant variability in the observed ET with a root mean square error (RMSE) of 0.4-1 mm day −1 and mean absolute percent error (MAPE) of 16-44%. Model intercomparison showed STIC1.2 to produce the least bias and good agreement with the observations, whereas SEBS and SCOPE1.7 revealed consistent underestimation and overestimation, respectively, in both years. (c) Underestimation of λE (and ET) in SEBS was mainly attributed to the underestimation in the roughness lengths of momentum and heat transfers (z 0M and z 0H ). While the underestimation of z 0M is associated with the empirical modeling structure, the underestimation of z 0H was associated with the overestimation of 'kB −1 -term' under high soil moisture and low atmospheric aridity conditions. Underestimation of both z 0M and z 0H led to an overestimation of the aerodynamic conductance (g A ) and sensible heat flux (H), which was consequently reflected in the underestimation of ET. (d) Although both SEBS and SCOPE1.7 had similar empirical parameterization of g A , a consistent overestimation of λE (and ET) in SCOPE1.7 was associated with the overestimation of the canopy-surface conductance (g S ) under high atmospheric aridity and also presumably due to the g S -photosynthesis modeling uncertainty in SCOPE1.7 under high atmospheric vapor pressure deficit.
(e) Despite all the three model captured substantial variability in λE (and ET), the principal difference between the models appear to be associated with the differences in g A and g S . Different magnitude g A and g S from all the models indicate the critical role of ambiguous parameterizations of these two important conductances for a broad spectrum of ecohydrological conditions. While SEBS require improved roughness length representation for enhancing the performance of g A sub-models under low fractional vegetation cover conditions; SCOPE1.7 requires robust parameterizations for both g A and g S , and default calibration parameters prior to large-scale ET monitoring in the wetlands. (f) The models showed promise as a quick and simple monitoring tool for wetland evapotranspiration. The simplified analytical model STIC1.2, requiring only surface-air temperature, humidity, and radiation data, can produce comparable results to more complex methods like SEBS under fully vegetated conditions and relatively better results under low fractional vegetation cover. Furthermore, this study demonstrated the model's potential for large scale ET mapping in the wetlands to capture the spatio-temporal ET dynamics. A dense network of radiation, temperature and humidity monitoring stations would also help create near-real time ET maps for the eco-hydrological studies in the Upper Biebrza National Park region.