Sensitivity of Surface Fluxes in the ECMWF Land Surface Model to the Remotely Sensed Leaf Area Index and Root Distribution: Evaluation with Tower Flux Data

: The surface-atmosphere turbulent exchanges couple the water, energy and carbon budgets in the Earth system. The biosphere plays an important role in the evaporation process, and vegetation related parameters such as the leaf area index (LAI), vertical root distribution and stomatal resistance are poorly constrained due to sparse observations at the spatio-temporal scales at which land surface models (LSMs) operate. In this study, we use the Carbon Hydrology Tiled European Center for Medium-Range Weather Forecasts (ECMWF) Scheme for Surface Exchanges over Land (CHTESSEL) model and investigate the sensitivity of the simulated turbulent ﬂuxes to these vegetation related parameters. Observed data from 17 FLUXNET towers were used to force and evaluate model simulations with different vegetation parameter conﬁgurations. The replacement of the current LAI climatology used by CHTESSEL, by a new high-resolution climatology, representative of the station’s location, has a small impact on the simulated ﬂuxes. Instead, a revision of the root proﬁle considering a uniform root distribution reduces the underestimation of evaporation during water stress conditions. Despite the limitations of using only one model and a limited number of stations, our results highlight the relevance of root distribution in controlling soil moisture stress, which is likely to be applicable to other LSMs.


Introduction
The water, energy, and carbon exchanges between the land surface and the atmosphere are key components of the Earth system. These exchanges are crucial for the understanding of the Earth system' response to the increase in greenhouse gas concentration and climate change. The variety and complexity of the underlying processes driving these exchanges (from boundary layer turbulence to plant physiology) are challenging both to observe and to model. From a modelling perspective, these fluxes provide boundary conditions to numerical weather prediction models (NWP) and earth system models (ESM). The longer time scales of some of the land surface processes [1,2], when compared with the atmosphere, are crucial for sub-seasonal to seasonal predictability [3][4][5] and to represent climate feedbacks [6,7].
The turbulent exchanges of latent heat flux (Qle), or evaporation, couples the surface water, energy and carbon budgets. Evaporation requires available water and energy [8]. Additionally, a significant part of evaporation is due to plant transpiration [9][10][11], which is coupled with vegetation photosynthesis and carbon exchanges, with a key role for soil moisture [12]. From an atmospheric perspective, the partition of the net surface energy into Qle and sensible heat flux (Qh) modulates the atmospheric boundary layer evolution with implications ranging from convection and cloud evolution [13] to the amplification of extremes [14,15]. Vegetation transpiration and carbon uptake for photosynthesis are mediated by stomatal opening. Therefore, stomatal conductance is an essential component in the representation of evaporation in land surface models (LSMs). Stomatal conductance is modulated by environmental factors (e.g., water stress) [16], and the leaf-to-canopy upscaling, dependent on the vegetation characteristics [17].
In this study, we focus on the impact of the representation of canopy resistance on the Qle and Qh fluxes in the European Center for Medium-Range Weather Forecasts (ECMWF) land surface model (LSM) Carbon Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (CHTESSEL) [18,19]. This study is motivated by the results of the Land Surface Model Benchmarking Evaluation Project (PLUMBER [20]), in which CHTESSEL participated, and by Ukkola et al. [21], which focused on seasonal-scale evaporative droughts. These studies compared the performance of several LSMs using atmospheric forcing from several FLUXNET towers. Best et al. [20] benchmarked the LSMs against simple physically based models and empirical relationships. The striking results showed that all participating LSMs were outperformed by an out-of-sample linear regression against downward shortwave radiation for the simulation of sensible heat flux and by a three-variable nonlinear regression that uses instantaneous atmospheric humidity and temperature in addition to downward shortwave radiation for the simulation of latent heat flux. Another key finding in PLUMBER was that LSMs had difficulties in simulating latent heat flux at sites with low annual precipitation, suggesting that this could be associated with the representation of water stress conditions in the models. The follow-up study of Ukkola et al. [21] found systematic biases across the LSMs, including CHTESSEL, in the simulation of energy and water fluxes under water-stressed conditions. Despite these limitations, Best et al. [20] also found that LSMs performed well, compared with the empirical benchmarks, when considering metrics for the extremes of the distributions. Therefore, the physical constraints in the LSMs and their continued development are paramount to simulate conditions outside of training conditions, as is the case for climate change scenarios [22].
Canopy resistance in CHTESSEL follows an empirical multiplicative formulation describing the effects of different environmental variables [23], which is a common approach in several LSMs. Another approach is to couple net photosynthesis with canopy resistance [24]. CHTESSEL also has a module to compute natural land carbon exchanges, but the carbon exchanges are decoupled from evaporation due to conflicting results in NWP [25], justifying the focus on Qle and Qh in this study. In the resistance approach used by CHTESSEL, vegetation characteristics are considered via the use of the leaf area index (LAI). The environmental factors account for water stress, which is dependent on soil moisture content varying linearly between unstressed vegetation and the wilting point. In this study, we address these two components of the canopy resistance formulation in CHTESSEL: (i) vegetation characteristics and (ii) water stress.
Vegetation characteristics are investigated by testing the use of a high-resolution LAI climatology based on the Moderate Resolution Imaging Spectroradiometer (MODIS) product. The impact of a high-resolution albedo climatology is also investigated. The use of this high-resolution data is motivated by the increased quality and amount of available Earth observations (EO) [26]. Here, we aim at evaluating the impact of using location specific LAI data when compared with the common LAI data that have been spatially aggregated to a coarse grid for NWP applications. The water stress representation is investigated by testing a uniform root distribution with an associated maximum rooting depth. This is compared with the current formulation in CHTESSEL that assumes an exponential root distribution [27].
The root distribution defines the maximum amount of water that is reachable by the plants, modulating the response of transpiration to moisture stress [12,28].
As in Best et al. [20], observed data from 17 stations of the FLUXNET global network of micrometeorological flux measurements are used to carry out offline simulations with CHTESSEL and evaluate sensible and latent heat fluxes, with particular emphasis on the representation of evaporation and its relation with water stressed conditions. The main objective of the study is to evaluate the representation of canopy resistance by testing (i) the added value of the high-resolution LAI and albedo to represent local station conditions and (ii) the impact of root distribution on soil moisture stress. The data and methods are described in the next section, followed by the presentation of the results. The results are discussed in Section 4 followed by the main conclusion of this study in Section 5.

Tower Flux Data
Observed data from 20 FLUXNET towers were initially considered in this study as in Best et al. [20]. These stations are widely used since they present a variety of vegetation types and climates [20,29,30]. The Spanish stations of ElSaler and ElSaler2 were removed due to the presence of large water bodies nearby and irrigated agricultural fields. The Canadian station of Merbleue was also problematic and thus removed since this tower is in the middle of a wetland, which is difficult to correctly model. This corroborates Haverd et al. [31], who also removed these three stations from their study, resulting in the 17 stations described in Table 1. The data include gap-filled half-hourly driving data for the LSM (air temperature and humidity, surface pressure, wind speed, precipitation and downward solar and thermal radiation) and surface sensible and latent heat flux used for model evaluation. In this study, we use the combined Moderate Resolution Imaging Spectroradiometer (MODIS) LAI and fraction of photosynthetically active radiation (FPAR) product (MCD15A3H.006). This is a 4 day composite dataset with a spatial resolution of 500 m available since July 2002 to the present. This MODIS product has been validated in different studies [32,33] reporting uncertainties in the order of 1 m 2 /m 2 . In addition to the MODIS LAI product, the Copernicus Global Land Service (CGLS) LAI was also used in this study. The product was obtained from SPOT-VGT and PROBA-V satellite observations with 1 km resolution using a neural network algorithm [34] and has been available since 1999. A comparison between MODIS and CGLS (a new version with a 300 m resolution) LAI products with ground based observations over North America reported root mean squared differences of 0.57 for the CGLS product and 0.81 to 0.89 for MODIS [35].
In addition to the LAI product, we also use the Combined MODIS bidirectional reflectance distribution function (BRDF) and Albedo products (MCD43A4.006) [36]. This product has been available since February 2000 with a 500 m resolution. Wang et al. [37] reported root mean squared errors ranging between 0.02 for forest regions and 0.03 for agricultural/grassland regions.
In this study, we apply the most restrictive quality control (QC) criteria in order to select only the best quality data. In the case of the albedo product, the model requires snow-free albedo. Therefore, the MODIS snow products (MOD10A1.006 and MYD10A1.006) were used to mask the albedo data when snow was present. A less restrictive quality filtering was applied for LAI [32], and good quality data were also selected (QC < 64-highest and good quality). A similar quality control data filtering was also applied to the CGLS data (only QC = 0 (land pixel clear observations) or QC = 512 (high latitude (lat) > 55 • and solar zenith angle > 70)).
MODIS LAI and albedo, as well as the CGLS LAI were extracted for the 17 FLUXNET stations (see Table 1) by averaging the data of the central pixel with the four other pixels in the direct vicinity. Monthly means were first computed and the monthly climatology computed excluding months with more than 50% of missing data, for the full available period until 2018 (since 1999 for the CGLS LAI, since 2003 for the MODIS LAI and since 2000 for the MODIS albedo).
The quality based data screening was too restrictive for the northern latitude stations, especially during the winter when satellite products are subjected to cloud contamination. In Hyytiala, Loobos and Palang, the MODIS albedo product has a large amount of cloud contamination, and not enough data are available to calculate the monthly climatology. The quality control data filtering was therefore relaxed and changed not only to the best quality data, but also to good quality data (QC <= 1). After this processing, Hyytiala presented high albedo values in January, associated with snow contamination. January's albedo was replaced by a linear interpolation between December and February. This shows the difficulty in generating snow-free albedo data [38], particularly in regions and during seasons that are predominantly snow (and cloud) covered. However, this is not expected to influence the model simulations as during January, the surface albedo has a small impact due to the reduced amount of available solar radiation in Hyytiala.

CHTESSEL Model
CHTESSEL is the land surface scheme of the ECMWF model [18,25]. In CHTESSEL, each grid-box is divided into up to 6 land fractions representing vegetation, soil, snow and interception. Surface fluxes are calculated separately for the different subgrid surface fraction (or "tile"), leading to a separate solution of the surface energy balance and skin temperature for each of these tiles, which are then aggregated in each grid-box. In each grid-box, two vegetation types can coexist (high and low vegetation) controlled by four parameters: dominant high and low vegetation type and the area fraction for the high and low vegetation. Each vegetation type is characterized by a series of fixed parameters: (i) minimum canopy resistance r smin , (ii) vegetation coverage density C veg , (iii) coefficient GD, for the dependence of the canopy resistance on water vapour pressure deficit, and (iv) the root distribution over the soil layers, specified by an exponential profile modulated by the coefficients a r and b r . The numerical values for the parameters (see Table S1 in the Supplementary Material) are based on the literature [27,[39][40][41][42][43][44] and on expert decisions in the context of the performance of the model in NWP. For example, the introduction of a seasonally variable LAI was accompanied with a revision of the minimum canopy resistance [25]. A detailed description of CHTESSEL can be found in the Integrated Forecasting System (IFS) documentation [45] and here, the details of the evaporation calculation are presented due to its relevance in the results' interpretation and discussion.
For high and low vegetation, the turbulent flux of water is given by: where ρ a is the air density, U L and q L are the wind speed and humidity at the lowest atmospheric model level, q sat (T sk ) the saturation humidity at skin temperature (T sk ) and C H the turbulent exchange coefficient (depending on the atmospheric stability). In addition to the aerodynamic resistance (r a ), the canopy resistance (r c ) [23] is calculated as: where r smin is the minimum stomatal resistance (see Table S1), f 1 is a function of downward short-wave radiation (R s ), f 3 a function of the atmospheric water vapour deficit (D a ) and f 2 the soil moisture resistance given by: where θ pwp and θ cap are the soil moisture at the permanent wilting point and at field capacity, respectively, and θ is a weighted average of the unfrozen soil water computed using the fraction of roots in each layer (R k ) [27] using the a r and b r coefficients (in Table S1): where z k+1/2 is the depth of the bottom layer k in m and z 1/2 is = 0. In CHTESSEL, the soil is discretized into four layers with thicknesses of 0.07, 0.21, 0.72 and 1.89 m, with lower bounds at 0.07, 0.28, 1 and 2.89 m.
In CHTESSEL, the state of vegetation is given by the LAI, entering the canopy resistance calculation normalizing r smin (Equation (2)). In the current operational NWP configuration, a satellite observation based climatology is considered for the representation of LAI. It is based on Collection 5 of MODIS (product MOD15A2). The climatology was derived from 9 years of data (2000 to 2008) and rescaled to a previous LAI static field used at ECMWF before 2012 [19].

Simulation Setup
Offline point simulations were performed using CHTESSEL driven by the meteorological data observed at the towers. The model is initialized with soil moisture at field capacity and runs once for the full length of the available forcing period for each station. The state at the end of the simulation is then used to provide initial conditions to start the main simulation. For most sites, this procedure guarantees enough time for spin-up, in particular for the deeper soil moisture. The exceptions are Palang and Kruger with only two years. In both cases, the top meter soil moisture in the two years of simulation does not show a significant drift. In Palang, there is an increase of the deep soil moisture during the simulation (see Figure S1 in the Supplementary Material), but still within the inter-annual variability. At each site, only one dominant high or low vegetation type is allowed, defined according to the plant functional type (PFT). Howard, Kruger and Mopane stations report PFTs that are not available in CHTESSEL (woody savannah and savannah). These stations were set as tall grass vegetation type (see Table 2), as this is the closest vegetation type in the vicinity of the stations in the global land cover dataset used by the model. Several model configurations were tested and are resumed in Table 3. These configurations aim at investigating the role of (i) high-resolution remote sensing LAI and albedo and (ii) model formulation and parameters. The first configuration, labelled as "CTR" for control, used CHTESSEL original input data with the ECMWF Integrated Forecasting System (IFS). The monthly climatologies of snow-free albedo and LAI were extracted from the nearest grid-point of operational NWP fields with a 9 km resolution. This model configuration serves as a reference against which other model configurations will be compared and was the same as used in [20]. The second experiment, labelled as "MALB ", uses the MODIS extracted high-resolution albedo monthly climatology over the stations. The third configuration, labelled as "MLAI" uses the MODIS extracted high-resolution LAI monthly climatology over the station. In both cases, the albedo and LAI data represent the nearest five MODIS pixels from the station (as described in Section 2.1.2). Table 3. Model simulations acronyms and detailed configuration. CHTESSEL.

CTR
Control simulation with default CHTESSEL parameters and input LAI and Albedo MALB Same as CTR, but replacing the input albedo climatology with the new high-resolution MODIS climatology MLAI Same as CTR, but replacing the input LAI climatology with the new high-resolution MODIS climatology MLAI_NOSMS Same as MLAI, but removing the soil moisture stress function from the canopy resistance (setting f 2 = 1) when the soil moisture is above the wilting point in Equation (3)  MLAI_RSMIN Same as MLAI, but selecting the optimal r smin for each station from a set of simulations with varying r smin between 25 and 500. MLAI_ROOT Same as MLAI, but using a uniform root distribution (Equation (5)) and selecting the optimal R DMAX for each station from a set of simulations with R DMAX of 0.5, 1, 2, and 3 m.
Two additional model configurations were tested to investigate the role of r smin and root distribution, both using the MODIS LAI high-resolution data. For r smin , a set of simulations varying r smin between 25 and 500 with a step of 25 was performed for each site. For the root distribution, a uniform root distribution was adopted up to a maximum rooting depth (R DMAX m) by changing the computation of the root fraction in each layer in Equation (4) to: In this configuration, four maximum rooting depths were tested: 0.5, 1, 2 and 3 m. The CTR root distribution for each vegetation type and for the tested uniform rooting depths are given in Table S1 (in Supplementary Material). The root density distribution and root depth are poorly constrained at large spatial scales as they vary between individual species and climate factors [46][47][48]. Therefore, in this study, a uniform distribution is a simple and straightforward first-order assumption.
An optimal model simulation for each site was then selected and labelled MLAI_RSMINfor the optimal simulations with varying r smin and MLAI_ROOT for the optimal simulations with varying the maximum rooting depth R DMAX . The selection of the optimal configuration is described in the following section. Finally, an idealized experiment removing the soil moisture stress function from the canopy resistance (MLAI_NOSMS, setting f 2 = 1 in Equation (2) when soil moisture is above the wilting point) was performed. This idealized experiment provides an estimate of evaporation in a situation of soil moisture at field capacity.

Evaluation
The simulated Qle and Qh are compared with the tower measurements using 4 metrics (i) the mean bias error (MBE), standard deviation difference (SD), correlation coefficient r , and normalized mean error (NME). The SD is computed as the absolute difference between 1.0 and the ratio between the simulated and observed standard deviation. The NME is the mean absolute error normalized by the mean absolute observed deviations from the mean. These metrics follow Best et al. [20] and are detailed in Appendix A. The calculations were performed on daily simulated and measured fluxes.
For the selection of the optimal simulations with varying r smin or maximum rooting depth, a ranking approach considering the four metrics for both Qh and Qle was considered. For each metric and flux, each simulation was ranked in ascending order from 1 to 20 in the case of r smin and 1 to 4 in the case of R DMAX . This resulted in rankings that were then added, and the simulation with the lowest rank was selected as the optimal simulation. This process was performed independently for each site and assumed that the metrics for Qle and Qh were of equal weight, and it did not distinguish between larger or smaller differences between metrics in the ranking. The use of the rankings of each metric without considering the actual metric value or the differences between simulations is a limitation as simulations ranked differently can be actually very close. This approach can be interpreted as a simple optimization strategy to select the optimal r smin or R DMAX for each site. However, it is not intended to be a proper model calibration as it is site specific and considers the entire set of observations. The search for the optimal simulations at each site by changing r smin or R DMAX was designed to provide an estimate of the upper bound on the best achievable performance of the model. It is a benchmark to determine if the model could be improved just by changing the particular parameter. The reduced number of stations (only 17) and sampling of vegetation types imposed limitations in performing a proper out-of-sample calibration to provide estimates of r smin and R DMAX for each of the vegetation types.

Comparison of LAI and Albedo
In this section, we compare the LAI and high-resolution albedo (station) with the original CHTESSEL input data. The high-resolution MODIS data for each tower are compared with the original CHTESSEL input data in terms of the mean and root mean squared differences (RMSD), which are computed for the mean annual cycle. The MODIS albedo is very close to CTR in all stations (see Table 4) with most of the stations with RMSD around 0.03. There are two exceptions, Fort Peck and Loobos with an RMSD of 0.07 mostly due to lower values in ALB when compared with CTR during winter, which are associated with snow contamination. Despite the importance of albedo in the surface energy balance, these results suggest that for these 17 stations, we should expect a small impact on the CHTESSEL simulation due to the update of the CTR albedo by the high-resolution MODIS albedo. The comparison of CTR and MLAI in Table 5 depicts comparatively large differences in several stations with the RMSD above one in Bugac, Hesse, Howland and Palang. The time series of the climatological LAI of CTR, MODIS and CGLS are shown for all stations in the Supplementary Material ( Figures S1-S17). The large differences in the climatology are mostly in the seasonal cycle amplitude, with MLAI presenting a more pronounced seasonal phenology. These results show some uncertainty in mapping LAI between the different earth observation products. Despite the differences between MODIS and CGLS, our results show a larger discrepancy between the products and the CTR climatology used in CHTESSEL with a median RMSD of 0.86 between CTR and MLAI when compared with a median RMDS of 0.43 between MODIS and CGLS LAI. While the CTR LAI is also based on MODIS data, it was derived from a shorter time period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) and was rescaled to match a previous LAI field used at ECMWF [19]. These LAI differences, both in the mean and in the annual cycle, are expected to impact the model simulations, which are investigated in the following section.

Fluxes Evaluation
The use of the high-resolution MODIS albedo had a negligible impact on model performance (not shown). This was expected due to the similarity between the albedo climatologies. Therefore, the MALB simulations are not shown nor further explored in this study. Figure 1 shows the distribution of the different metrics for Qh and Qle comparing the CTR and MLAI simulations (among other model configurations to be discussed in the next section). In the MLAI simulation, the Qle MBE is reduced from a median value of -3.02 W/m 2 in CTR to -2.19. Similarly, the Qh MBE is reduced from 12.08 W/m 2 in CTR to 8.48 W/m 2 in MLAI. The standard deviation ratio (SD) is also improved with a reduction from 0.26 and 0.17 for Qle and Qh, respectively, in CTR to 0.19 and 0.14 in MLAI. This improved variability is also seen in the temporal correlation with an increase from 0.79 to 0.83 in the Qle from CTR to MLAI and from 0.77 to 0.80 in Qh. Finally, the normalized mean error (NME) has a slight reduction in Qh from 0.78 in CTR to 0.75 in MLAI and slightly increases for Qle (0.62 in CTR and 0.64 in MLAI). We note that the statistical significance of the changes in the model performance were not assessed because of the small sample size. Despite the small sample, only 17 stations, these results suggest a modest, but positive impact of replacing CTR LAI by the high-resolution MODIS data for the station. It was not possible to identify a systematic improvement in model performance associated with the use of the high-resolution LAI dataset. Each station reacts differently, depending on the time period considered. Some stations in some years benefit from the MODIS product, while other stations or time periods show a deterioration of model performance. Moreover, the impacts of the high-resolution LAI dataset on model performance could not be generalized or classified by plant functional type (PFT) as stations within the same PFT do not react similarly to LAI changes. Stations with a good performance in CTR (Fort Peck, Harvard, Howard, Howlandm, Hyytiala, Loobos, University of Michigan) are not impacted by the LAI changes. Additionally, the large differences between the CTR and LAI climatologies are not reflected in the changes in model performance in some stations like Hesse or Howlandm. However, some stations suffering from poorer performances in CTR benefit from MLAI as for example Bugac, Espirra, Palang, Sylvania and Tumba. Considering the four metrics used to characterise model performance, the temporal correlation benefits the most from the use of the high-resolution LAI since this product has a larger inter-seasonal variability than CTR.
LAI enters the canopy resistance formulation normalizing the minimum stomatal resistance (r smin in Equation (2)). It is therefore expected that changes in the LAI climatology would benefit from an adjustment of the r smin parameter. Such an adjustment could be achieved by constraining r smin so that the ratio r smin /LAI would stay unchanged for each PFT. However, due to the small sample of stations for each PFT, such an approach is not feasible. A set of simulations with varying r smin was performed to select an optimal r smin with the best model performance for each station, as described in Section 2.4. The optimal r smin for each site is presented in Table 2. The simulations with the optimal r smin (MLAI_RSMIN) have a negligible impact on model performance (see Figure 1, comparing MLAI and MLAI_RSMIN). The only noteworthy impact is a clear reduction of the SD for Qle (see Figures 1 and S2), with the remaining performance metrics unchanged when compared with MLAI. These results show that r smin acts primarily on the seasonal amplitude of Qle, explaining the large impact on SD, but neutral in the remaining performance metrics for Qle and Qh. In particular, the neutral impact on Qh SD highlights that the improved seasonal variability in Qle does not necessarily lead to an improved Qh seasonality.  There are three sites with a large increase of r smin from CTR to MLAI_RSMIN: Kruger (100 to 300), Mopane (175 to 300) and Sylvania (175 to 400). In all three cases, the mean LAI decreases from CTR to MLAI (see Table 2). This suggests that an adjustment that would keep the ratio r smin /LAI unchanged with changes in LAI would not be optimal or an indication that the default r smin parameters are not optimal. At these three stations, the increased r smin results in a reduction of evaporation, as expected, which is clearly seen in Sylvania (see Figure S3). However, the Qle reduction in Sylvania is associated with an unusually large r smin of 400 sm −1 , suggesting that this might be compensating missing processes in the model and/or errors in the driving/observations data (Ukkola et al. [21] reported that Sylvania was excluded from their study due to precipitation problems). There are three sites with a large decrease of r smin : Fort Peck (100 to 25), Blodgett (250 to 150) and Espirra (240 to 75). At these stations, there is also a reduction of the LAI from CTR to MLAI. From these three stations, Espirra stands out with the largest reduction of r smin , resulting in an increase of Qle and a decrease of Qh, particularly during spring (see Figure S4). We note that Qle in Espirra was underestimated in CTR and was further reduced in MLAI due to the reduction of the LAI (mean LAI changed from 2.38 in CTR to 1.42 in MLAI). The optimal r smin change from 240 to 75 acts to reduce the r smin /LAI ratio from 100 sm −1 to 52 sm −1 . We also note that the other two evergreen broadleaf tree stations, Palang and Tumba, have a mean LAI of about 4.4 m 2 /m 2 with an r smin /LAI ratio of 54 sm −1 . These results suggest that the high-resolution MODIS LAI extracted for Espirra station underestimates the station vegetation LAI [49].
The metric distributions in Figure 1 show some dispersion, associated with the occurrence of a few outliers. The detailed results for each station and score are presented in the Supplementary Material, in Figures S18-S25. It is possible to identify five stations with recurrent problematic results for Qh/Qle in different metrics: Blodgett, Espirra, Mopane, Palang and Kruger. In Kruger, CHTESSEL overestimates substantially the Qle after rainfall events, followed by an underestimation of the Qle and an overestimation of the Qh (see Figure S5). Kruger is a very dry station with soil moisture below the wilting point, and the peak evaporation following rainfall events is mainly driven by bare soil evaporation. The remaining four stations, Blodget, Espirra, Mopane and Palang, share a similar problem of an early evaporative reduction associated with dry soil moisture conditions resulting in the overestimation of Qh. Although not seen clearly in the evaluation metrics, Amplero and Tumba also show similar Qle temporal error patterns. This early evaporation reduction, or evaporative droughts, as proposed by Ukkola et al. [21], is further investigated in the following section.

Soil Moisture Stress
Ukkola et al. [21] identified a systematic overestimation of evaporative droughts in several LSMs, including CHTESSEL, which is consistent with the results in this study mainly in Amplero, Blodget, Espirra, Howard, Mopane, Palang and Tumba. In this study, evaporative drought refers to the regular dry season evaporation reduction, as well as anomalous dry periods as defined by Ukkola et al. [21]. Figure 2 shows the time series of Qle and soil moisture for Amplero, Blodgett and Espirra (time series for all stations are provided in the Supplementary Material). In these three examples, it is possible to identify the early reduction of evaporation associated with the drop of soil moisture (Figure 2: compare Qle reduction in blue with top soil moisture reduction in red). At these stations, the evaporation starts to diverge from the observations when soil moisture in the top meter (where most of the roots are present in CTR; Table S1) falls below field capacity and approaches the wilting point. During this period, the deep layer soil moisture (dashed red line in Figure 2) is above or close to field capacity. This is further illustrated by comparing the MLAI Qle simulations with LAI_NOSMS in Figure 2 (blue vs. dashed-blue line). These idealized simulations, without soil moisture stress, differ from MLAI exactly during the dry-down period, being closer to the observed Qle. The LAI_NOSMS simulations present some unrealistic Qle variations (see for example 2003 in Amplero in Figure 2) associated with the drop of soil moisture below the wilting point that stops evaporation to avoid water conservation problems that could result in numerical instabilities. Despite this limitation, the LAI_NOSMS simulations suggest that the evaporative drought in CHTESSEL is tightly associated with the soil moisture stress formulation. The discrepancy between soil moisture conditions in the top meter versus the bottom layer during the dry-down and associated evaporative drought ( Figure 2) motivated the exploratory revision of the root distribution.
The root distribution with the uniform formulation (Equation (5)) for the different maximum rooting depths (0.5, 1, 2 and 3 m) is presented in Table S1, allowing the comparison with the CTR root fraction for each layer and vegetation type. The main changes of the uniform rooting depth are a decrease of the roots fraction in the top layer(s) and an increase of root fraction in the deeper layers. This will give more weight to the deeper layers when computing the root zone soil moisture. The optimal maximum rooting depth for each station (see Section 2.4) in MLAI_ROOT is shown in Table 2.
Of the three short grass stations, Bugac and Fort Peck show an optimal R DMAX of 0.5 m with negligible differences for MLAI, while Amplero shows a R DMAX of 2 m with an increase of the Qle in late spring and summer. The increased Qle in Amplero in August/September 2005 (Figure 3) partially reduces the evaporative drought in MLAI, but results is an overestimation of the Qle in June/July. A closer investigation of the full time series of simulated fluxes and the high-resolution MODIS LAI ( Figure S6) identifies an LAI drop during the 2005 summer, which cannot be represented when using a climatology. A similar behaviour is also seen in the 2004 summer. Despite the reduction of the evaporative drought in late summer in Amplero with the uniform root distribution, other factors influence the inter-annual and inter-seasonal variability of the Qle and Qh, which are not captured by CHTESSEL, as well as observational errors and scale mismatches. Inter-annual variability is also seen in Howard with a negligible impact of the uniform root distribution in 2005, but with a positive impact in 2002, 2003 and 2004 (see Figure 3 for 2004 and Figure S7 for the full time series). Mopane also displays evaporative drought events, but with very dry soil moisture conditions (see Figure S8), the impact of the uniform root distribution is negligible. In the forest sites, all R DMAX are between 2 and 3 m, with the exception of Harvard with 1 m and Howland and Hyytiala with 0.5 m. These three forest sites, plus Loobos and University of Michigan, are among the best performing in CTR and MLAI. The impact of the uniform root distribution is small since all five sites have abundant precipitation and soil moisture conditions that are close to field capacity. Hesse station is the only case of a consistent deterioration of Qle with the uniform rooting depth and R DMAX of 2 m, mainly due to a considerable overestimation, contrasting with a reduction of the Qh bias. This mixed effect in the Qh and Qle justifies the use of several metrics and joint Qh and Qle in the model evaluation. Other processes and/or driving/observation limitations are likely responsible for the errors. The remaining four forest sites Blodget, Espirra, Palang and Tumba (see Figure 3) all display excessive evaporative drought associated with soil moisture dry-down. At these sites, the experimental formulation of uniform rooting depth with the optimal R DMAX partially reduces the excessive evaporative drought. The largest impact is seen in Palang and Blodgett, which are among the sites with the worst Qle and Qh correlations. From these sites shown in Figure 3, the uniform root distribution has the largest impact on Qle, while r smin and MLAI changes only show minor differences for CTR in Espirra.

Discussion
The systematic underestimation of the Qle and overestimation of the Qh for CTR (see Figure 1) is opposite the findings of Martens et al. (2020) [50], who compared ERA5 data against observations from 143 FLUXNET sites. Similar biases in CHTESSEL considering a subset of 51 FLUXNET stations were also found [51]. The differences of these studies compared to our results are likely associated with the different sampling of stations. It is also worth noting that the energy budget closure of eddy covariance sites can induce systematic errors [52,53]. Therefore, the interpretation of the systematic biases in the simulated Qle and Qh must be cautious. In this study, the use of several evaluation metrics that account for both mean biases and variability, with equal weights, partially mitigates the problems that could arise from systematic observation errors.
Despite the availability of long-term EO records of LAI, in this study, we considered only the high-resolution LAI climatology. This was driven by two reasons: (i) several stations encompass periods before the EO record were available (before 1999 in the case of CGLS or 2002 in the case of MODIS), and (ii) we considered that the first step was to evaluate the impact of the high-resolution station location conditions, neglecting inter-annual variability. However, our results suggest that the lack of inter-annual variability in the vegetation state might be responsible for some of the model problems. Amplero station is an example where we observed some relation between the time-varying high-resolution LAI (see Figure S6, bottom panel, dotted blue line) and the Qle inter-annual variability. There is robust evidence of the added value of assimilating LAI in land-surface models, enhancing their capabilities to monitor vegetation phenology and fluxes [54][55][56][57]. This is particularly relevant in models that have a prognostic evolution of vegetation phenology.
The soil moisture resistance, or water stress factor, in CHTESSEL has two components: (i) the derivation of the root zone soil moisture content (using an exponential root distribution) and (ii) a linear transformation between field capacity (unstressed vegetation) and the wilting point. In this study, we focused only on the derivation of the root zone soil moisture content, testing a new uniform root distribution and the associated maximum rooting depth. There are different approaches to relate root zone soil moisture and the water stress factor, such as curvilinear relationships or the use of soil matric potential [58,59]. Such approaches, combined with the rooting depth distribution, are likely to further enhance the model capability to represent root water uptake. Moreover, Liu et al. [12] reported that a dynamical rooting depth evolution improved the performance of the Noah-MP-Crop model under drought-like situations. Additionally, vapour pressure deficit and stem xylem conductance are also known to influence stomatal regulation [60,61], but were not considered in this study. The presence of a shallow water table could also influence latent heat flux during drought periods [62]. The water table is not represented in CHTESSEL, and we are not aware of such effects at the 17 sites considered.
Finally, the metrics distribution for MLAI_ROOT (in Figure 1) showed more potential to reduce the Qle and Qh bias and NME and increase the correlation than the calibration of r smin alone. This is achieved by partially addressing the excessive evaporative drought at Blodget, Espirra, Palang and Tumba stations. A joint optimization of both r smin and maximum rooting depth is likely to further improve the model simulations. However, it was not the aim of this study to optimize/calibrate CHTESSEL for these stations, but to investigate in the current model formulation which aspects require further attention. The use of the subset of FLUXNET stations to evaluate land surface models' performance is a common practice [20,21,31,63,64].

Conclusions
This study focused on the impact of the representation of canopy resistance on the simulations of latent and sensible heat fluxes by the CHTESSEL model. Observed data from 17 FLUXNET towers were used to carry out offline simulations with a particular emphasis on the representation of evaporation and its relationship with water stress conditions. Three constraints of canopy resistance were evaluated: (i) the role of vegetation characteristics via the use of high spatial resolution LAI representing local station conditions; (ii) the impact of the minimum canopy resistance; and (iii) the impact of a uniform root distribution on soil moisture stress. The replacement of the current LAI climatology used by CHTESSEL, based on an older MODIS climatology with a coarse resolution, by a new high-resolution climatology representative of the stations locations did not significantly affect the simulated surface fluxes. The close relationship between LAI and the minimum canopy resistance was investigated showing some potential to improve the latent heat flux variability via an adjustment of this parameter. However, these changes did not impact the excessive evaporative drought at several stations, as reported by Ukkola et al. [21]. In CHTESSEL, we found that this limitation was tightly associated with the depletion of the top meter soil moisture, while water was still available at deeper layers. However, the current model formulation for root distribution was not capable of reaching this deeper reservoir, resulting in an early reduction of evaporation, when compared with the observations. The replacement of the current exponential roots' profile by a uniform root distribution and associated maximum rooting depth reduced the underestimation of evaporation during water stress conditions. Despite the limitations associated with a reduced number of stations, our results highlight the importance of root distribution in controlling soil moisture resistance in water stress conditions. However, root distributions are weakly constrained observationally on the global scale. Therefore, research is necessary to understand the implications of these changes on the global water and energy budgets, as well as in the coupling with the atmosphere. The proposed uniform root distribution with a single associated parameter, the maximum rooting depth, is also appealing for a parameter optimization. These could be further addressed along with the revision of land cover and vegetation recently proposed for CHTESSEL [65], which identified the necessity to calibrate vegetation related parameters [51]. Considering the coupled nature of the surface water and energy cycles and the relevance of land-atmosphere coupling, a calibration methodology considering multiple observational datasets, covering different water/energy components and temporal time-scales, should be favoured [66], also to possibly include some indirect information of poorly constrained parameters through variables coupled with them.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4433/11/12/1362/s1, Table S1: CHTESSEL vegetation types and associated parameters. R 1 to R 4 denote the root fraction in Layers 1 to 4. The final 4 lines indicate the root fraction when considering a uniform rooting depth with a maximum rooting depth R DMAX of 0.5, 1, 2 and 2.89 m. Figure S1: Time series of Sensible heat flux, latent heat flux, soil moisture and LAI in Palang. The turbulent fluxes' time series compare the observations (grey) with the simulations: CTR (black), LAI (blue), LAI_RSMIN (red), LAI_ROOT (green) and LAI_NOSMS (dashed blue). The soil moisture time series shows the CTR the top 3 layers' meter soil moisture (top meter, solid back) and the bottom layer soil moisture (dashed black), as well as the soil moisture at the wilting point (blue) and field capacity (red). The LAI time series compares CTR (black) with the high-resolution MODIS LAI time series (dotted blue), the high-resolution MODIS climatology (blue), the climatology of MODIS considering the 0.25 • bounding box (dashed blue) and the CGLS LAI climatology (grey). Figure S2: Latent heat flux SD at the 17 stations for each simulation: CTR, MLAI, MLAI_RSMIN and MLAI_ROOT. Figure S3: Same as Figure S1, but for Sylvania station. Figure S4: Same as Figure S1, but for Espirra station. Figure S5: Same as Figure S1, but for Kruger station. Figure S6: Same as Figure S1, but for Amplero station. Figure S7: Same as Figure S1, but for Howard station. Figure S8: Same as Figure S1, but for Mopane station. Figure S9: Same as Figure S1, but for Blodgett station. Figure S10: Same as Figure S1, but for Bugac station. Figure s11: Same as Figure S1, but for Fort Peck station. Figure S12: Same as Figure S1, but for Harvard station. Figure S13: Same as Figure S1, but for Hesse station. Figure S14: Same as Figure S1, but for Howlandm station. Figure S15: Same as Figure S1, but for Hyytiala station. Figure S16: Same as Figure S1, but for Loobos station. Figure S17: Same as Figure S1, but for Tumbarumba station. Figure S18: Same as Figure S1, but for University of Michigan station. Figure