Seasonal Adaptation of the Thermal-Based Two-Source Energy Balance Model for Estimating Evapotranspiration in a Semiarid Tree-Grass Ecosystem

The thermal-based two-source energy balance (TSEB) model has accurately simulated energy fluxes in a wide range of landscapes with both remote and proximal sensing data. However, tree-grass ecosystems (TGE) have notably complex heterogeneous vegetation mixtures and dynamic phenological characteristics presenting clear challenges to earth observation and modeling methods. Particularly, the TSEB modeling structure assumes a single vegetation source, making it difficult to represent the multiple vegetation layers present in TGEs (i.e., trees and grasses) which have different phenological and structural characteristics. This study evaluates the implementation of TSEB in a TGE located in central Spain and proposes a new strategy to consider the spatial and temporal complexities observed. This was based on sensitivity analyses (SA) conducted on both primary remote sensing inputs (local SA) and model parameters (global SA). The model was subsequently modified considering phenological dynamics in semi-arid TGEs and assuming a dominant vegetation structure and cover (i.e., either grassland or broadleaved trees) for different seasons (TSEB-2S). The adaptation was compared against the default model and evaluated against eddy covariance (EC) flux measurements and lysimeters over the experimental site. TSEB-2S vastly improved over the default TSEB performance decreasing the mean bias and root-mean-square-deviation (RMSD) of latent heat (LE) from 40 and 82 W m−2 to −4 and 59 W m−2, respectively during 2015. TSEB-2S was further validated for two other EC towers and for different years (2015, 2016 and 2017) obtaining similar error statistics with RMSD of LE ranging between 57 and 63 W m−2. The results presented here demonstrate a relatively simple strategy to improve water and energy flux monitoring over a complex and vulnerable landscape, which are often poorly represented through remote sensing models.


Introduction
Land surface models, mathematical representations of surface-atmospheric exchanges, are important tools to understand fluxes of energy and mass, which drive climatic and Earth system processes [1]. These models provide vital information to understand the response of ecosystems to climate and environmental changes, and monitor Earth system dynamics [2,3]. Latent heat flux (LE), the aggregated water flux consisting of evaporation from the soil and other wet surfaces (LE s ) and plant transpiration (LE c ), has recently been the subject of extensive research [4] due to its importance in evaluating ecosystem functional properties. It is a key process that interlinks the water and energy budget along with carbon cycling through the processes of transpiration and photosynthesis [5]. However, LE is a complex and difficult to estimate phenomena, as it is affected by numerous variables related to the characteristics of the soil-surface, atmosphere and vegetation [6]. As a result, LE is highly variable and dynamic, making remote sensing techniques notably useful to predict it at different temporal and spatial scales.
Different remote sensing methods have been proposed to estimate LE, which range from empirically-based to more process-based modeling schemes [7][8][9]. Among these, surface energy balance (SEB) models estimate LE as the residual of the energy balance, which exploit radiometric land surface temperature (LST) from thermal remote sensing as a key boundary condition [10][11][12][13]. There are two main types of SEB models: one-source models that do not discriminate between vegetation and soil components and dual-source models, which explicitly separate the temperature and energy exchange between the vegetation and soil sources considering the directional effects in the remotely sensed radiometric temperature. The two-source energy balance (TSEB) model [10,14], has been applied in a variety of landscapes obtaining reliable estimates of LE [15][16][17][18][19][20], including for water stressed conditions [21][22][23]. TSEB was originally developed for homogeneous cover types, however, adaptations to the model framework have been implemented to better depict partial canopy cover [14]. Past studies demonstrated that TSEB better accounted for the effects of partial vegetation cover compared to one-source models [15,23]. However, in heterogeneous and complex surfaces such as tree-grass ecosystems (TGE), Earth observation and modeling methods tend to have greater uncertainties, e.g., [24,25].
TGEs, a prevalent savanna-like landscape covering nearly 15% of the total Earth surface [26], are frequently located in semi-arid regions with limited and highly seasonal water availability [27], where a greater urgency is needed to monitor the scarce water resources. Since TSEB treats the vegetated layer as a unique layer, the parameterization and application of this model poses greater difficulty in a TGE landscape where two (sometimes three) distinct vegetation covers are present (i.e., isolated trees over a grassland and/or shrubland), making it difficult to represent the vegetation according to a unique set of parameters. Additionally, in these semi-arid ecosystems, the degree in which each type of vegetation (i.e., tree and grass) influence land-surface interactions changes throughout the year depending on their differentiated phenological stages [28]. During the growing season, trees and the grass understory, along with underlying soil, all interact to contribute to the radiative transfer and turbulent exchanges [29]. However, during the dry summer periods, the grass layer senesces due to meteorological and soil conditions (i.e., water availability, air temperature, vapor pressure deficit), converting the system into (nearly) bare, rather rough, soil with scattered trees, substantially changing land-atmospheric dynamics [30]. As Andreu et al. [31] discussed, further adaptations to the TSEB model structure and parameters may be necessary to accurately simulate energy fluxes over TGEs due to the multiple vegetation layers present. In this regard, remote sensing data and algorithms can play a pivotal role in improving land-surface modelling by fully exploiting information from different sensors (e.g., synthetic-aperture radar (SAR) and light-detection and ranging (LiDAR)) to better estimate auxiliary variables (e.g., canopy height) that are normally assumed or static. To account for these spatial and temporal complexities, TSEB was adapted here to consider two major phenological periods, based on when the grass layer is active during the growing season and not active (i.e., senesced) during the dry summer period.

TSEB Model Overview
The TSEB model was first proposed in Norman et al. [10], with important adjustments described in Kustas and Norman [14]. Its main inputs are LST, derived from thermal infrared (TIR) radiation, vegetation structural properties (e.g., LAI, canopy height), and meteorological forcing (irradiance, air temperature and wind speed). The principle source of uncertainty within TSEB lies in the estimation of the sensible heat flux (H), which is calculated through the heat transport equation (Equation (1)). (1) where H is sensible heat flux (W m -2 ); is the volumetric heat capacity of air (J m -3 K -1 ); is the aerodynamic temperature of the surface (K); is the air temperature at a reference/measurement height (K); and is the aerodynamic resistance to heat transport (s m -1 ). The heat transport equation is satisfied when using aerodynamic surface temperature (i.e., surface temperature at the canopy source-sink height), however, LST obtained from TIR remote sensing (i.e., skin radiometric surface temperature) can differ up to several degrees compared to the aerodynamic surface temperature [10], and their relationship is not well established (i.e., [47]). TSEB thus tackles this issue by assuming that the total blackbody thermal radiance that is emitted by the bulk surface is weighted by the fraction of vegetation observed by the sensor and the emission of both soil and vegetation surfaces, as expressed in Equation (2) taken from Norman et al. [10]: where is the fraction of vegetation observed by the TIR sensor at zenith angle θ and is mainly a function of LAI; is the vegetation canopy temperature (K); and is the soil surface temperature

TSEB Model Overview
The TSEB model was first proposed in Norman et al. [10], with important adjustments described in Kustas and Norman [14]. Its main inputs are LST, derived from thermal infrared (TIR) radiation, vegetation structural properties (e.g., LAI, canopy height), and meteorological forcing (irradiance, air temperature and wind speed). The principle source of uncertainty within TSEB lies in the estimation of the sensible heat flux (H), which is calculated through the heat transport equation (Equation (1)).
where H is sensible heat flux (W m −2 ); ρC p is the volumetric heat capacity of air (J m −3 K −1 ); T o is the aerodynamic temperature of the surface (K); T A is the air temperature at a reference/measurement height (K); and R H is the aerodynamic resistance to heat transport (s m −1 ). The heat transport equation is satisfied when using aerodynamic surface temperature (i.e., surface temperature at the canopy source-sink height), however, LST obtained from TIR remote sensing (i.e., skin radiometric surface temperature) can differ up to several degrees compared to the aerodynamic surface temperature [10], and their relationship is not well established (i.e., [47]). TSEB thus tackles this issue by assuming that the total blackbody thermal radiance that is emitted by the bulk surface is weighted by the fraction of vegetation observed by the sensor and the emission of both soil and vegetation surfaces, as expressed in Equation (2) taken from Norman et al. [10]: Remote Sens. 2020, 12, 904 5 of 29 where f (θ) is the fraction of vegetation observed by the TIR sensor at zenith angle θ and is mainly a function of LAI; T c is the vegetation canopy temperature (K); and T s is the soil surface temperature (K). Using this two-layer approach, the energy balance is formulated in TSEB for each of the layers separately as follows: where R N is the net radiation (W m −2 ); LE is latent heat flux (W m −2 ); G is the soil heat flux (W m −2 ); and subscripts s and c refer to soil and vegetation canopy sources, respectively. Note that horizontal heat advection, canopy heat storage as well as the energy used for CO 2 fixation, are neglected in Equation (3). Radiative transfer and absorption through the canopy (R N,c and R N,s ) are simulated through a radiative transfer model considering spectral differences in the visible, near-infrared and longwave radiation, multiple scattering between soil and canopy, and direct and diffuse irradiance, as described in chapter 15 of Campbell and Norman [48] and incorporated in TSEB by Kustas and Norman [14]. H is derived using an in 'series' resistance network [10] (Figure 2), and applying Equations (4)- (6), which allows for heat turbulent interchange between the canopy and soil layers: where T ac is the air temperature in the canopy space (K) and is equivalent to the aerodynamic temperature (T o ); R S is the resistance to heat transfer in the boundary layer above soil layer (s m −1 ); R X is the bulk canopy resistance to heat transfer (s m −1 ); R A is the aerodynamic resistance to heat transfer based on the Monin-Obukhov similarity theory. Refer to Appendix A of Norman et al. [10] for details on the series resistance scheme.
Remote Sens. 2020, 12, 904 5 of 31 (K). Using this two-layer approach, the energy balance is formulated in TSEB for each of the layers separately as follows: where is the net radiation (W m -2 ); LE is latent heat flux (W m -2 ); G is the soil heat flux (W m -2 ); and subscripts s and c refer to soil and vegetation canopy sources, respectively. Note that horizontal heat advection, canopy heat storage as well as the energy used for CO2 fixation, are neglected in Equation (3). Radiative transfer and absorption through the canopy ( , and , ) are simulated through a radiative transfer model considering spectral differences in the visible, near-infrared and longwave radiation, multiple scattering between soil and canopy, and direct and diffuse irradiance, as described in chapter 15 of Campbell and Norman [48] and incorporated in TSEB by Kustas and Norman [14]. H is derived using an in 'series' resistance network [10] (Figure 2), and applying Equations (4)- (6), which allows for heat turbulent interchange between the canopy and soil layers: where is the air temperature in the canopy space (K) and is equivalent to the aerodynamic temperature ( ); is the resistance to heat transfer in the boundary layer above soil layer (s m -1 ); is the bulk canopy resistance to heat transfer (s m -1 ); is the aerodynamic resistance to heat transfer based on the Monin-Obukhov similarity theory. Refer to Appendix A of Norman et al. [10] for details on the series resistance scheme. Since Equation (2) has two unknowns ( and ), the canopy layer is assumed, as a first estimate, to be initially transpiring at a potential rate ( ) using the Priestley-Taylor formulation: where is the initial canopy transpiration estimate (W m -2 ); is the Priestley-Taylor coefficient (default is 1.26) (-), defined in this case only for the canopy component [9]; is the fraction of vegetation that is green and hence actively transpiring (-); ∆ is the slope of the saturation vapor Since Equation (2) has two unknowns (T c and T s ), the canopy layer is assumed, as a first estimate, to be initially transpiring at a potential rate (LE ci ) using the Priestley-Taylor formulation: Remote Sens. 2020, 12, 904 where LE ci is the initial canopy transpiration estimate (W m −2 ); α PT is the Priestley-Taylor coefficient (default is 1.26) (-), defined in this case only for the canopy component [9]; f g is the fraction of vegetation that is green and hence actively transpiring (-); ∆ is the slope of the saturation vapor pressure curve at air temperature T a (kPa K −1 ); and γ is the psychrometric constant (kPa K −1 ). This initial assumption, where the canopy is transpiring without water stress (i.e., α PT = 1.26), is used to initialize the model, which allows to solve all the systems of equations presented above (Equations (2) to (7)). However, if the vegetation is stressed, the Priestley-Taylor formulation will overestimate the transpiration of the canopy and hence underestimate canopy temperature, which, in order to conserve the total surface temperature in Equation (2), will result in an overestimation of both soil temperature and soil sensible heat flux (i.e., Equation (4)). At this point, in order to preserve the energy balance in the soil (Equation (3)), LE s would likely yield negative values, meaning soil condensation, which is unrealistic under daytime conditions. As it is assumed that condensation does not occur during daytime convective conditions, an iteration procedure is applied that reduces α PT in small steps until realistic soil fluxes area achieved (i.e., LE s ≥ 0). A more complete discussion on conditions that reduces α PT is given in Anderson et al. [49] and Li et al. [39]. For more implementation details, the reader is referred to the source code (https://github.com/hectornieto/pyTSEB) and Norman et al. [10].

Radiation Transmission in Sparse Vegetation
The structure and distribution of foliage in the vegetative layer has an important impact on the dynamics of radiation interception and transmission through the canopy [14,49,50]. This in turn has a very important implication for radiometric temperature partitioning between the soil and vegetation components and their resulting contribution to heat fluxes [49]. The original TSEB radiative transfer equations assume a randomly distributed homogenous (i.e., non-clumped) canopy [10]. However, sparse vegetation is generally clumped and tends to intercept less radiation for the same foliage density compared to vegetation randomly dispersed over the surface [14,48]. As such, the clumping index (Ω) quantifies the spatial distribution of foliage to account for non-randomness in vegetated structures. LAI is multiplied by the clumping factor to obtain effective LAI (ΩLAI). As incorporated in Kustas and Norman [14]. TSEB estimates Ω as a function of the difference between the canopy gap fraction compared to a homogenous case using Equation (8) (Section 2 in Kustas and Norman [14]) and estimating the beam extinction coefficient assuming an ellipsoidal leaf angle distribution (LAD) with Equation (9) (Section 15.2 in Campbell and Norman [48]).
where Ω(0) is the clumping index when the vegetation canopy is viewed at nadir (-); k be is the beam extinction coefficient (-) based on a ellipsoidal LAD function (Campbell and Norman [48]); f c is the vegetation fractional cover (-); F is the local LAI (F = LAI/ f c ); and χ LAD is the leaf inclination distribution function chi parameter (-) [48]. Ω is also dependent on the solar zenith angle (θ S ) and is estimated using Equation (10), as described in Section 15.13 of Campbell and Norman [48].

Resistances within and below Canopy
The semi-empirical derivations of R S and R X , as introduced in Norman et al. [10] and Kustas and Norman [14], are dependent on wind speeds below and within the vegetative canopy, respectively (Equations (11)-(12)).
where u s is the wind speed just above the surface where the impact of soil roughness is minimal (i.e., z 0,soil (ms −1 ); u d o +z 0m is the wind speed within the canopy-air interspace at the height of momentum source/sink (ms −1 ); l w is the effective leaf width size (m); and c (m s −1 K −1/3 )), b (-) and C' (s 1/2 m −1 ) are coefficients taken from Kustas and Norman [14] and Norman et al. [10] based on the works of Sauer and Norman [51], Kondo and Ishida [52], and McNaughton and Van Den Hurk [53] (default values are 0.0025 m s −1 K −1/3 , 0.012 and 90 s 1/2 m −1 respectively). u s and u d o +z 0m in Equations (11) and (12) are both derived from an observed wind speed above the canopy that is extrapolated through and below the canopy based on the exponential wind attenuation law [54]. The two key roughness parameters, the zero-plane displacement height (d 0 ) and the aerodynamic roughness length for momentum transfer (z 0m ) are estimated from the vegetation structure. When considering the grass canopy layer, the traditional fixed ratio to canopy height method is used to estimate d 0 and z 0m [48].
In the case of the tree canopy layer, a different approach is used (see Appendix B), which additionally considers the canopy shape and density´s effect on the roughness parameters. This study follows the procedure described in Schaudt and Dickinson [55], which stems from the work of Raupach [56] and Lindroth [57], and is generally viewed as more suitable for tall wooded vegetation. These methods do not apply corrections associated to the roughness sub-layer (e.g., Weligepolage et al. [58]). However, Alfieri et al. [37] found that the TSEB modelled fluxes were largely insensitive to differences in estimated d 0 and z 0m obtained from various methods in a vineyard. This is due to the fact that rougher surfaces, such as tall canopies, are better coupled to the atmosphere and, hence, the transport of heat and water vapour is more dominated by the canopy´s physiological mechanisms rather than aerodynamic roughness [59].

Eddy Covariance and Bio-Meteorological Measurements
The three EC towers (CT-FLUXNET identifier ES-LMa, NT-FLUXNET ID ES-LM1, and NPT-FLUXNET ID ES-LM2), simultaneously operating within the Majadas study site, provided all necessary inputs to run TSEB (Table 1), except for continuous vegetation LAI estimates (Section 2.3.2). H and LE estimations from the EC system served to benchmark model performance. Details on the data processing are found in El-Madany et al. [43]. Data were obtained between the period of 2015-01-01 and 2017-12-31. Data which did not satisfy the assumptions of the EC methods [60] or with low turbulent conditions [61] were filtered from the dataset. Since TSEB forces the closure of energy balance (via Equation (3)), energy balance closure of the EC flux measurements was ensured by allocating the residuals to the observed LE, assuming that errors in LE are larger than H due to issues related to the instrumentation (i.e., Foken et al. [62]), as in previous studies (e.g., Guzinski et al. [63]; Kustas et al. [64]; Twine et al. [65]; Niaghi et al. [66]). Measurements from 2015 of the CT tower were selected as the 'main' simulation site/period and used to perform the SA (Section 2.4.2) and the main model runs (Section 2.4). Other towers (NT and NPT) and years (2016 and 2017) provided independent evaluations of model performance (Section 2.4.5).
While TSEB has used different methods to estimate ground heat flux (G) (e.g., [10,67]), this study directly forced in situ G measurements as a model input to limit uncertainty and noise in turbulent Remote Sens. 2020, 12, 904 8 of 29 flux estimations associated with errors in G. However, for the validation model runs, G was modeled with the time difference approach of Santanello and Friedl (Section 2.4.5). At each tower, the weighted average of eight soil heat flux plates, all buried at 5cm depth, represented G at the ecosystem level. Four of these were located in open grass and the other four plates were placed below the tree canopy. The weighted contributions of the plates below the tree and in open grass were estimated based on the solar zenith angle (SZA) throughout the day. At solar noon (i.e., SZA is 0 degrees), it is estimated that the tree canopy cover is 20% and that the plates in the shade below the trees contribute this same percentage (and the open grass plates contribute 80%). At sunset/sunrise (i.e., SZA > 85 degrees), the open grass area (i.e., non-shaded area) does not contribute to the ecosystem G. As such, the SZA is normalized between these two limits and used to scale the relative contribution between the soil heat flux plates buried below the tree canopy and those in the open grass. Note that corrections related to heat storage above the soil heat flux plates were not applied and, as such, G is likely slightly underestimated. The radiometric LST was estimated using proximal sensing from the CNR4 (Kipp & Zonen, Delft, Netherlands) longwave radiometers installed in the EC system (Equation (13)). These radiometers cover a footprint representing roughly 19-25% of tree cover for the three towers (i.e., CT, NT and NPT) used. These are very similar to the average tree contribution to the EC footprint reported in El-Madany et al. [43].
where L out and L in are upwelling and downwelling longwave radiance; σ is the Stephan-Boltzman constant; ε sur f = 0.99 f c + 0.94(1 − f c ) is the surface emissivity where leaf and bare soil broadband emissivity were taken as 0.99 and 0.94, respectively [68,69]. These broadband emissivity values were supported by computing broadband emissivity from typical soils and vegetation from the ECOSTRESS spectral library (https://speclib.jpl.nasa.gov/).

Vegetation Biophysical Measurements
In situ grass LAI measurements were acquired in the context of the FLUXPEC project (http: //www.lineas.cchs.csic.es/fluxpec/) at eleven plots (25 m × 25 m) randomly placed around the CT in nineteen field campaigns between 2013-10-31 and 2016-11-16 [71]. Destructive samples were collected in two 25 cm × 25 cm quadrants within each plot and LAI was derived using gravitational methods, where green and non-green elements were separated to compute both total and green LAI. Time series of reflectance factors from the MODIS/Terra and Aqua Nadir BRDF-adjusted Reflectance Daily L3 500m v006 (MCD43A4) product were acquired for pixels centered in each of the three towers ( Figure 1) [72]. The 500m by 500m pixel size was considered adequate to characterize the tower footprint. The normalized difference vegetation index (NDVI) was derived using band 1 (red: 620 nm-670 nm) and band 2 (NIR: 841 nm-876 nm) as follows: An empirical relationship between NDVI and in-situ destructive grass LAI measurements was developed specifically for the Majadas experimental site (Appendix A, Figure A1). In order to obtain the ecosystem LAI adjusted for the tree canopy effect within the tower footprint ( Figure 3), in-situ tree LAI measurements were incorporated to achieve a weighted average between tree (20%) and grass (80%) ( Figure A2). The average tree LAI acquired using the LAI-2200 plant canopy analyzer (LAI-2200) (LICOR Bioscience USA, 2011) during five field campaigns at different seasonal periods between 2017-2018 was used as a reference. Local tree LAI ranges between 1.39 and 1.75 m 2 m −2 (about 0.3 m 2 m −2 effective/landscape LAI) and has low inter-annual variability [28]. Linear interpolation between sampling dates was performed to obtain a continuous daily time series ( Figure A2).

Model Simulations and Evaluation
Global and local SAs were performed on the main parameters and inputs, respectively using the default TSEB model (sections 2.4.1 and 2.4.2). In addition, two end member simulations, where either a pure grassland or broadleaf forest was assumed for the entire year, were performed to better understand the effect of the vegetation canopy on simulated fluxes for extreme cases (section 2.4.3). The model was then adapted considering two distinct seasonal periods, each having a dominating vegetation layer (section 2.4.4). All simulations were evaluated for daytime fluxes (i.e., when Sdn > 50 W m -2 ) since TSEB is designed to model fluxes for daytime conditions with remote sensing data

Model Simulations and Evaluation
Global and local SAs were performed on the main parameters and inputs, respectively using the default TSEB model (Sections 2.4.1 and 2.4.2). In addition, two end member simulations, where either a pure grassland or broadleaf forest was assumed for the entire year, were performed to better understand the effect of the vegetation canopy on simulated fluxes for extreme cases (Section 2.4.3). The model was then adapted considering two distinct seasonal periods, each having a dominating vegetation layer (Section 2.4.4). All simulations were evaluated for daytime fluxes (i.e., when S dn > 50 W m −2 ) since TSEB is designed to model fluxes for daytime conditions with remote sensing data (Section 2.4.5).

Default TSEB Model Configuration
In this default configuration (hereafter as TSEB-DF), the vegetated layer was parameterized attempting to depict the mix between tree and grass observed over tower footprint, where the vegetation cover was dominated by grass, but the turbulent resistances were assumed to be largely affected by the sparse tree layer [43]. As such, ecosystem LAI (weighted average of grass and tree LAI, Figures 3 and A2) was used as input and vegetation resistances and roughness were derived using the weighted average of grass and tree canopy height (h c = 0.8 × 0.5m + 0.2 × 8m = 2m). The roughness parameters were estimated following the procedure as described in Schaudt and Dickinson [55]. Table 2 summarizes the parameter values used for this default model configuration. The green fraction ( f g ) was set to 0.7 being roughly the average value observed from grass field measurements. The model simulations were performed at the sub-hourly (30-min) time step for 2015 over the CT EC footprint and benchmarked against EC derived LE and H. Table 2. Parameter values of the TSEB-DF model configuration.

Sensitivity Analysis Sobol´Global Parameter Sensitivity Analysis
Variance-based SA methods decompose the total variance between simulated and observed data into various parts, determining the contribution of the different parameters and their combined interactions on total output variance. In this work, the Sobol' SA method was used as it is able to compute the 1st, 2nd and total order sensitivity indices. To apply this methodology, a selection of model parameters and their respective bounds were configured (Table 3). Based on this, parameters sets for n simulations were generated using the Sobol´sequence [40,41], a quasi-random method with quasi-Monte-Carlo integration.
where y is an indicator of the model performance and {θ 1 , θ 2 , θ 3 . . . , θ n } is the model parameter set which controls the model behavior and, thus, its performance, y. The Sobol' approach decomposes the variance observed in y through the variance in all θ factors (Equation (15)) by separating parameters into terms of increasing dimensionality, with each successive dimension representing greater interaction of parameters, represented as: where Var(y) is the total output variance; V i is the portion of variance contributed by parameter θ 1 (also known as the first order variance or main effect); and V ij is the portion of variance attributed through the interaction between θ i and θ j and so on. The Equation (16) holds if all the parameter factors are independent [40], which is the case for this study. Parameter values do not correlate and are independently assigned over a uniform distribution between the specified bounds (Table 3). Aerodynamic resistances (Rx) [53] The first-order sensitivity (S i ), second-order sensitivity (S ij ) and total order sensitivity (S Ti ) indices are calculated using Equations (17) to (19) as: where V ∼i is the average variance depicted when all parameters except for θ 1 vary (i.e., θ 1 is kept fixed); S Ti represents the contribution of both the direct (1st order) and indirect (sum of 2nd order) effects of θ i on the total variance, Var(y).
Model variance was calculated based on H since LE is computed through the residual of the energy balance in TSEB, where errors in H are transposed into LE. Root-Mean-Square Deviation (RMSD) was used as the objective function to quantify model output variance for the entire simulation year, where tower-based EC H provided the observed time series.
where H TSEB 30min is modeled H at a 30 minutes time step; H EC 30min is the observed EC H; and N is the total number of observations used. The SA was undertaken with TSEB-DF, considering the 11 parameters within TSEB (i.e., second group of variables in Table 1), which are used in the sub-modules of radiation transfer between vegetation and soil ( f c , w c and X LAD ) roughness and resistance scheme (h c , z0 soil , l w , b, c, C ), and the initial canopy transpiration estimate from the Priestley-Taylor formulation (α PT and f g ). Their respective bounds were based on literature or the physical limits of the parameter in question (Table 3). We note that the lower bound for α PT = 1.26 is merely used as a first estimate for the model initialization. TSEB has an internal procedure to reduce α PT due to conditions of stress, and, since no a priori information on stress is known, it is not recommended to use a value lower than 1.26 for model initialization [20]. A total of 45,500 simulations were used for this SA. The Python package SALib (Sensitivity Analysis Library in Python) [73] was integrated with the Python implementation of TSEB (pyTSEB, https://github.com/hectornieto/pyTSEB) to perform these analyses.

Input Local Sensitivity Analysis
Errors associated to the remote sensing derived input data (i.e., LST and LAI) are another important source to the total output uncertainty. Since input data are dynamic, a different method was needed compared to the static parameters described above for the global SA. Random white noise, mimicking potential errors of each input variable, over a uniform distribution of +-3K and +-0.4 m 2 m −2 were added to the 'observed' LST and LAI, respectively. The range of errors in LST was chosen to fit within the typical uncertainties associated with remote sensing based LST retrievals over vegetation (e.g., Sobrino et al. [74]). The range in LAI perturbation was based on the errors associated with the LAI empirical model used in this study (Appendix A, Figure A1). The sensitivity indices (SI) were then computed based on the partial derivative of the output result with respect to change in input caused by the implemented random noise, using Equation (21) adapted from van Griensven et al. [34].
where M is model output; ρ refers to the different model inputs; ∆ρ i refers to the pertubation of the input caused by artificial noise. The SI is based on the change observed for modeled H. The input SA was implemented for 365 days (during 2015) with a 30-min time step. Since TSEB incorporates LAI at the daily time scale, the SI is derived based on the daily aggregated absolute change in H in relation to the absolute change in LAI for that day. For the LST analysis, the effect of a deviation in LST on the resulting H will likely be impacted by the time in which the deviation occurs (i.e., a perturbation in LST during morning time will likely impact H differently compared to a midday perturbation). To normalize the time of perturbation, for each simulation day, only the mean net perturbation of LST during midday (between 11:00 and 13:00 UTC) was used to analyze the net effect on H during that same time period. Hence, SIs for both LAI and LST were derived for each daily time step using the absolute change in inputs with the resulting absolute change in output (i.e., H).

End Member Simulations
In addition to the TSEB-DF simulations, two end member scenarios were conducted as limiting case studies. These type of end member configurations are typical for models applied to global remote sensing products, which often assign model parameter values (e.g., h c or l w which are typically more difficult to estimate directly) based on a specific or assumed land/vegetation type (e.g., [75]). However, errors in global land cover maps, which tend to be more common for TGEs (e.g., [24]), will cause for a less accurate parameterization and depiction of the surface conditions. For example, the Climate Change Initiative (CCI) land cover product [76] designates the Majadas study site as either 'mosaic cropland' or 'mosaic herbaceous cover'. Therefore, these end member scenarios represent a similar or expected parameterization to the configurations applied for global scale energy flux products (e.g., [75,77]).
The first case assumes a pure grassland, ignoring the tree components, for the entire year (hereafter as TSEB grass ). The second case ignores the grass layer, simulating the vegetated component as a scattered evergreen broadleaf forest (hereafter as TSEB tree ). The estimation of the roughness characteristics differed for each scenario as described in Section 2.2.2 (i.e., fixed ratio for grass layer and the method of Schaudt and Dickinson [55] when considering tree canopy) and the most influential parameters (as derived from the SA) were modified to characterize the assumed dominant vegetation with standard values (see Section 3.1.3). In addition, instead of ecosystem LAI, each end member scenario uses the LAI corresponding to the assumed vegetation layer (i.e., grass or tree LAI as shown in Figure A2). These end member scenarios were performed on the CT tower during 2015.

Two-Season Modeling Approach
We propose here an adaptation of the TSEB-DF by using a two-season modeling approach (hereafter as TSEB-2S) that divides the annual simulation into two main phenological periods: a grass dominated growing period (i.e., grass-soil system) and a tree dominated summer/dry period (i.e., tree-soil system). This modeling scheme attempts to depict the vast changes observed in TGE surface condition due to phenology. Phenological processes are important in this ecosystem [28] and these changes in vegetation cover alter the turbulent conditions (i.e., roughness), radiation transmission and, hence, the surface energy balance [29]. In semi-arid TGEs, during the summer period, the grass understory rapidly transforms into a dry and rough layer with minimal transpiration from vegetation. Therefore, the TSEB-2S separates the simulation period for these two major conditions observed: 1) the grass is active and transpiring (growing season) and 2) the grass species are senesced and not transpiring (summer), neglecting any transition period between seasons due to the rapid grass senescence period that characterize semi-arid TGEs [28]. The transition dates were estimated using an asymmetric gaussian filter over the MODIS (MCD43A4 product) derived NDVI time series, where the dry period was assumed to begin when vegetation (i.e., NDVI) begins to decay (downward inflection point) and the dry period ends when vegetation (i.e., NDVI) begins to re-green (upward inflection point). For instance, for the simulation year of 2015, the dry period begins on May 13th and ends on October 24th (refer to Figure 3 for all transitions dates used). As such, this method can be a general and automatic approach to derive the seasonal transition dates at the global scale, which may be applied to other TGE or semi-arid sites with marked seasonality.

Model Evaluation
TSEB-2S was evaluated using independent model simulations for different spatial (three tower sites) and temporal (three different years) characteristics. The benchmark simulation run (i.e., TSEB-DF) and the SA were forced with observed G to limit the noise and uncertainties related to the errors in estimated G. This way we can more accurately pinpoint the uncertainties related to the H (and, by consequence, LE) retrievals. However, the model validation runs were performed by modeling G using the approach of Santanello and Friedl instead of using observed G. This was implemented to better evaluate the model in a more operational setting, where remote sensing-based retrievals may not have in-situ G measurements available. Model performance was evaluated against the observed EC tower measurements in question and quantified with RMSD (Equation (20)), mean bias (Equation (22)) and the Pearson´s correlation coefficient (r).
Additionally, an evaluation of the modeled LE partitioning, between LE c and LE s , was performed. Since TSEB-2S assumes a different and unique vegetation structure and cover for different seasonal periods, it was interesting to evaluate if the model can also obtain reliable estimates of the LE partitioning along with the bulk fluxes. Modeled LE s was evaluated against independent lysimeter LE measurements (LE lys ) of the understory during the dry summer periods at CT, when the grass is senesced and, thus, LE lys is assumed to be equivalent to soil evaporation (i.e., LE s ). For a detailed description of the lysimeter set-up in the Majadas experimental site, refer to Perez-Priego et al. [42].

Global Sobol´Parameter SA
The parameter global SA showed that f c had the largest impact on model results, both as the main effect (S i ) and through its interactions with other parameters (S T ) (Figure 4). In general, parameters related to vegetated structure (h c ) and cover ( f c and f g ) demonstrated the largest sensitivity. Notably, many parameters had low first order sensitivities (close to 0) but relatively important total order sensitivities, indicative of large interactions between parameters. The b parameter, used in the computation of R S S , was the only more empirically derived parameter that demonstrated a relatively high influence on modelled output.

Local Input SA
Uncertainties related to LST and LAI showed important impacts on model results ( Figure 5). Simulated H is particularly sensitive to deviations in LST, where differences of a unit change of LST (∆ 1K) is associated to a median 17.3% change in modelled H. The uncertainty in LAI demonstrated, compared to LST, relatively less influence on H. A 0.1 m 2 /m 2 change in LAI is associated with a 2.9% change in simulated H with TSEB-DF.

TSEB grass , TSEB tree and TSEB-2S Model Configuration
The global SA results illustrated in Figure 4 showed that parameters related to the vegetation cover (i.e., f c and f g ) and structure (i.e., h c ) had the most influence on model performance. Since these are physical and measureable parameters, these were set for each end-member scenario (Section 2.4.3) according to the assumed dominant vegetation as described in Table 4. The f g was set to 0.7 for the grass configuration according average field observations, while it was set to 0.9 for the evergreen tree cover to account for the presence of non-photosynthetic material such as branches. The h c for the tree layer was set to 8 m according to the mean height of trees within the study site [43] and its f c was set 0.2, this being about the average tree canopy fraction observed for the tower footprints [43]. Additionally, the l w was increased to 0.05 m for the tree cover to better represent the larger leaves of the oak species (compared to the grass species) within the study site. The resistance b coefficient, used in the formulation of R S , was increased for the scattered tree cover, as was successfully tested and implemented in Kustas et al. [21] to better simulate rough surfaces and partially vegetated surfaces. The b parameter was given a value of 0.034 for the tree-soil summer period, based on the range reported in Sauer et al. [70] for partial vegetation conditions [21,70]. Other parameters, such as χ LAD , which is presumably different between grasses and trees, were not modified since they showed little sensitivity ( Figure 4).  To improve the seasonal depiction and related changes to the ecosystem, the TSEB model was adapted using a two-season modelling approach (TSEB-2S) as described in Section 2.4.4. The simulation was divided into two main phenological periods where a dominant vegetation is assumed in each case and TSEB is parameterized according to the respective end-member parameter values shown in Table 4 (i.e., TSEB grass during the growing season and TSEB tree during the dry summer).

Model Evaluations for Main Simulation Period (2015 CT)
TSEB-DF, the end-member configurations (TSEB grass , TSEB tree ), and TSEB-2S were all first evaluated against EC measurements from CT for the entire year of 2015. Figure 6 shows the TSEB-DF modelled vs. observed half hourly turbulent fluxes as well as the annual trend of daytime average sensible heat flux, both modeled and observed. Model results of TSEB-DF vastly underestimate H (bias: -45 W m −2 ) and, consequently, overestimate LE (bias: 40 W m −2 ) as compared with observed EC data. As shown in Figure 6 (left panel), high errors are observed throughout (RMSD of H: 82 W m −2 ) but errors stem particularly from the very large underestimations of H during the hot and dry summer period ( Figure 6, right panel). TSEB-DF was attempted to be parameterized as a mixed layer between tree and grass vegetation. However, it was not able to accurately reproduce LE (and H) largely because its parameterization was not able to consider the drastic change occurring during the summer, where roughly 80% of surface cover (i.e., grass understory) becomes non-physiologically active and, thus, non-transpiring. The end member simulations demonstrate the 'boundary conditions' when the simulations completely ignored one of the vegetation layers throughout the year (Figure 7). Predictably, TSEBgrass, underestimates H even more than TSEB-DF (bias: -54 W m -2 ), while, by contrast, LE is largely underestimated in TSEBtree (bias: -23 W m -2 ), with large errors (RMSD: 82 W m -2 ). Both TSEBgrass and TSEBtree performed poorly because they do not adequately depict the surface conditions for the different phenological periods with one of the two major vegetation covers being ignored throughout the year. TSEBtree performed slightly better compared to TSEBgrass because it better captured the conditions and fluxes for the summer drought season, which was the period with the largest biases and uncertainties.  The end member simulations demonstrate the 'boundary conditions' when the simulations completely ignored one of the vegetation layers throughout the year (Figure 7). Predictably, TSEB grass , underestimates H even more than TSEB-DF (bias: −54 W m −2 ), while, by contrast, LE is largely underestimated in TSEB tree (bias: −23 W m −2 ), with large errors (RMSD: 82 W m −2 ). Both TSEB grass and TSEB tree performed poorly because they do not adequately depict the surface conditions for the different phenological periods with one of the two major vegetation covers being ignored throughout the year. TSEB tree performed slightly better compared to TSEB grass because it better captured the conditions and fluxes for the summer drought season, which was the period with the largest biases and uncertainties.
The performance of the proposed TSEB-2S parameterization is shown in Figure 8. The simulation of LE and H for 2015 improved considerably when using TSEB-2S compared to TSEB-DF (and end-members). Errors with TSEB-2S decreased substantially (decrease in RMSD from 82 W m −2 to 59 W m −2 and 55 W m −2 for LE and H, respectively) and the seasonal average daily H was much improved, particularly during the summer/dry period (RMSD of daily mean H decreased from 59 W m −2 to 26 W m −2 ). TSEB-2S allowed different parameterizations to be applied according to the assumed 'dominant' vegetation structure and cover in contrasting phenological periods observed in TGEs. All metrics for model performance improved with the use of TSEB-2S. RMSD and bias decreased substantially compared to TSEB-DF and the correlation (i.e., r) between modeled and observed data increased. However, the winter-autumn period (roughly between DOY 25 and 125) still showed a slight systematic H underestimation with TSEB-2S similar to TSEB-DF. and time series of simulated (red) and observed (black) daytime daily mean H (right) for TSEB-DF in 2015 at CT.
The end member simulations demonstrate the 'boundary conditions' when the simulations completely ignored one of the vegetation layers throughout the year (Figure 7). Predictably, TSEBgrass, underestimates H even more than TSEB-DF (bias: -54 W m -2 ), while, by contrast, LE is largely underestimated in TSEBtree (bias: -23 W m -2 ), with large errors (RMSD: 82 W m -2 ). Both TSEBgrass and TSEBtree performed poorly because they do not adequately depict the surface conditions for the different phenological periods with one of the two major vegetation covers being ignored throughout the year. TSEBtree performed slightly better compared to TSEBgrass because it better captured the conditions and fluxes for the summer drought season, which was the period with the largest biases and uncertainties. The performance of the proposed TSEB-2S parameterization is shown in Figure 8. The simulation of LE and H for 2015 improved considerably when using TSEB-2S compared to TSEB-DF (and end-members). Errors with TSEB-2S decreased substantially (decrease in RMSD from 82 W m -2 to 59 W m -2 and 55 W m -2 for LE and H, respectively) and the seasonal average daily H was much improved, particularly during the summer/dry period (RMSD of daily mean H decreased from 59 W m -2 to 26 W m -2 ). TSEB-2S allowed different parameterizations to be applied according to the assumed 'dominant' vegetation structure and cover in contrasting phenological periods observed in TGEs. All metrics for model performance improved with the use of TSEB-2S. RMSD and bias decreased substantially compared to TSEB-DF and the correlation (i.e., r) between modeled and observed data increased. However, the winter-autumn period (roughly between DOY 25 and 125) still showed a slight systematic H underestimation with TSEB-2S similar to TSEB-DF.

TSEB-2S Validation
TSEB-2S substantially improved estimated fluxes compared to the default implementation of TSEB (TSEB-DF). To further validate this modeling strategy, TSEB-2S was tested independently for different years (2016 and 2017) with very different intra-annual dynamics, and against two other EC towers (NT and NPT) within the experimental site (using the same meteorological forcing). In addition to evaluating the bulk fluxes, the LE partitioning between soil evaporation and transpiration was evaluated using lysimeter measurements at the CT.

Independent Evaluations for 2016 and 2017 over CT, NT and NPT Towers
TSEB-2S was tested for all EC towers (CT, NT and NPT) in 2016 and 2017. The results from these other towers and years ( Figure 9) were within similar error ranges to the benchmark TSEB-2S 2015 simulation at the CT. RMSD and bias for LE, in all validating simulations, range between 57 and 63 W m -2 and -2 and -12 W m -2 , respectively. This range in error is very much in line to the benchmark simulation period even though modeled, rather than observed, G was implemented. The available energy (AE), Rn -G, was also generally well modeled for all model runs ( Figure A3).

TSEB-2S Validation
TSEB-2S substantially improved estimated fluxes compared to the default implementation of TSEB (TSEB-DF). To further validate this modeling strategy, TSEB-2S was tested independently for different years (2016 and 2017) with very different intra-annual dynamics, and against two other EC towers (NT and NPT) within the experimental site (using the same meteorological forcing). In addition to evaluating the bulk fluxes, the LE partitioning between soil evaporation and transpiration was evaluated using lysimeter measurements at the CT.

Independent Evaluations for 2016 and 2017 over CT, NT and NPT Towers
TSEB-2S was tested for all EC towers (CT, NT and NPT) in 2016 and 2017. The results from these other towers and years ( Figure 9) were within similar error ranges to the benchmark TSEB-2S 2015 simulation at the CT. RMSD and bias for LE, in all validating simulations, range between 57 and 63 W m −2 and −2 and −12 W m −2 , respectively. This range in error is very much in line to the benchmark simulation period even though modeled, rather than observed, G was implemented. The available energy (AE), Rn -G, was also generally well modeled for all model runs ( Figure A3).

LE Partitioning
The LE partitioning in TSEB-2S was evaluated against independent estimates from the lysimeter at the CT footprint during the summer drought, when there is no grass transpiration and, thus, LE measured by the lysimeter (LElys) corresponds to soil evaporation [30,42]. Results indicate a general, slight LEs underestimation during the dry period for the three years analyzed ( Figure 10). Additionally, important errors stem from the greater short-term variability of the modeled fluxes compared to the relatively stable lysimeter measurements. The LEs underestimation (bias: -18 W m 2 ) and errors (RMSD: 41 W m 2 ) were largest for 2016. The larger errors during this year were partly explained by the greater underestimation from the dry-down period (roughly DOY170-185), which is also apparent, while less substantial, for 2015. As shown in Figure 10, modeled LEs had a much more rapid decrease in magnitude compared to the gradual decrease observed from the lysimeter during this dry-down period. TSEB-2S was able to capture the temporal dynamics and evaporation peaks relatively well but the dry-down after the rare rainfall events was also more gradual in the observed lysimeter compared to the modeled fluxes. This may be explained by the evaporation

LE Partitioning
The LE partitioning in TSEB-2S was evaluated against independent estimates from the lysimeter at the CT footprint during the summer drought, when there is no grass transpiration and, thus, LE measured by the lysimeter (LE lys ) corresponds to soil evaporation [30,42]. Results indicate a general, slight LE s underestimation during the dry period for the three years analyzed ( Figure 10). Additionally, important errors stem from the greater short-term variability of the modeled fluxes compared to the relatively stable lysimeter measurements. The LE s underestimation (bias: −18 W m 2 ) and errors (RMSD: 41 W m 2 ) were largest for 2016. The larger errors during this year were partly explained by the greater underestimation from the dry-down period (roughly DOY170-185), which is also apparent, while less substantial, for 2015. As shown in Figure 10, modeled LE s had a much more rapid decrease in magnitude compared to the gradual decrease observed from the lysimeter during this dry-down period. TSEB-2S was able to capture the temporal dynamics and evaporation peaks relatively well but the dry-down after the rare rainfall events was also more gradual in the observed lysimeter compared to the modeled fluxes. This may be explained by the evaporation being sustained from the deeper soil, which is more difficult to capture from the model since it is forced with the 'skin' radiometric surface temperature (i.e., LST).
Remote Sens. 2020, 12, 904 20 of 31 being sustained from the deeper soil, which is more difficult to capture from the model since it is forced with the 'skin' radiometric surface temperature (i.e., LST).

Discussion
The proposed TSEB-2S vastly improved model performance in simulating LE and H compared to the default TSEB configuration (RMSD and bias of modelled H decreases from 82 and -45 W m -2 to 55 and -5 W m -2 , respectively). The simple assumption of two separate phenological and modelling periods, one dominated by a grass-soil system and the other dominated by a tree-soil system, allowed for a two-layer model to accurately simulate turbulent energy fluxes in an essentially three (tree-grass-soil) layer ecosystem. TSEB assumes only one vegetated layer, being more or less photosynthetically active, over a non-photosynthetically active layer (i.e., bare soil or similar). Therefore, it becomes difficult to properly parameterize the model when numerous vegetation covers are simultaneously present. Additionally, as in the case of many TGEs, the influence of the grass understory changes throughout the year, where it becomes largely non-photosynthetically active (i.e., not transpiring) during the summer. Therefore, these vastly different seasonal conditions and characteristics were better captured using TSEB-2S, by assuming a 'dominant' vegetation cover for different phenological periods, compared to the default, single mixed vegetation, configuration of TSEB (TSEB-DF). The relatively simple and automatic separation of seasons from a MODIS-NDVI time series can easily be extrapolated to other TGEs or similar temporally dynamic sites, as these data are globally available. The results also demonstrate that large changes to the vegetated surface have a considerable influence on the surface energy balance. This confirms the importance of Figure 10. Time series of simulated (blue) and observed (black) daytime daily mean LE s at CT during the peak summer in 2015, 2016 and 2017 using TSEB-2S. LE measured by the lysimeter is assumed to be equivalent to LE s during the summer period due to grass senescence.

Discussion
The proposed TSEB-2S vastly improved model performance in simulating LE and H compared to the default TSEB configuration (RMSD and bias of modelled H decreases from 82 and −45 W m −2 to 55 and −5 W m −2 , respectively). The simple assumption of two separate phenological and modelling periods, one dominated by a grass-soil system and the other dominated by a tree-soil system, allowed for a two-layer model to accurately simulate turbulent energy fluxes in an essentially three (tree-grass-soil) layer ecosystem. TSEB assumes only one vegetated layer, being more or less photosynthetically active, over a non-photosynthetically active layer (i.e., bare soil or similar). Therefore, it becomes difficult to properly parameterize the model when numerous vegetation covers are simultaneously present. Additionally, as in the case of many TGEs, the influence of the grass understory changes throughout the year, where it becomes largely non-photosynthetically active (i.e., not transpiring) during the summer. Therefore, these vastly different seasonal conditions and characteristics were better captured using TSEB-2S, by assuming a 'dominant' vegetation cover for different phenological periods, compared to the default, single mixed vegetation, configuration of TSEB (TSEB-DF). The relatively simple and automatic separation of seasons from a MODIS-NDVI time series can easily be extrapolated to other TGEs or similar temporally dynamic sites, as these data are globally available. The results also demonstrate that large changes to the vegetated surface have a considerable influence on the surface energy balance. This confirms the importance of vegetation characteristics in controlling and mediating ecosystem level energy fluxes, where seasonal dynamics and phenology of vegetation are key considerations for land-atmospheric modelling. The findings presented here describe a relatively simple and general approach to account for changes in land-atmospheric exchanges due to phenology such as those occurring in semi-arid TGEs. In this regard, a similar strategy can be implemented in other global soil-vegetation-atmosphere-transfer (SVAT) or prognostic models to better accommodate seasonality for ecosystems with comparable characteristics.
The TSEB-2S model demonstrated robustness, being able to accurately simulate different intra-annual dynamics for various years (i.e., different phenological timings) and sites with differing surface conditions, in this case, derived from a nutrient fertilization experiment (Figure 9). Results and the associated magnitudes of errors for all TSEB-2S model runs are similar to the error bounds found in other energy balance model studies (e.g., [14,15,21,22,38,63,78]) and close to the typical uncertainty of surface turbulent flux measurement systems (i.e.,~50 W m −2 ; [79]). For instance, RMSD of LE between 62 and 70 W m −2 were achieved by Timmermans et al. [23], who compared the use of TSEB against a one-source energy balance model for a sparsely vegetated grassland and rangeland. In the work of Boulet et al. [78], different dual-source model schemes were tested and obtained an RMSD between 53 and 73 W m −2 for midday instantaneous LE for both irrigated and rainfed wheat fields. Therefore, the results presented here, are in line with past studies related to LE retrievals, with these considering much more homogeneous land cover types, which better fit the assumptions inherent in SEB models compared to the multiple and more structurally complex vegetation cover in TGEs. Andreu et al. [31] applied various modified versions of TSEB, notably with different wind profiles and roughness schemes, in similarly complex TGEs and reported errors between 44 and 60 W m −2 for simulated LE. These are comparable to the error bounds presented here even though quite different approaches were used. In our study, the basic TSEB model structure was not modified, instead opting to alter the model parameterization depending on the phenological period using typical land cover characteristics (i.e., grasslands and evergreen broadleaved trees). As such, the findings here may complement the methods used by Andreu et al. [31], who also reported the largest errors during the dry, summer period.
The global SA demonstrated that TSEB was most sensitive to f c . This is largely due to f c being a key input to estimate the Ω (Kustas and Norman [14]), which in turn affects the radiation interception and partitioning between the mixed vegetated and soil surfaces. The estimation of Ω is mostly based on the transmission through the vegetated layer where f c is used to obtain a local LAI (LAI/ f c ) and as an important weighting factor for the gap fraction estimation (Section 2.2.1). These results are largely in line with results from Li et al. [39], who found that the Ω uncertainties had a large impact on flux outputs from TSEB. They tested incremental f c values, between 0.1 and 1, that resulted in sharp H changes, particularly between higher f c values, an indication of a high sensitivity for this parameter. The f g was found to be more sensitive compared to α PT even though both parameters are part of the initial canopy transpiration estimate from the Priestly-Taylor formulation (Equation (7)). This is the case since α PT is merely used as a priori estimate for initializing the model and is gradually decreased when there are conditions of water stress, until the energy balance is achieved with realistic daytime fluxes using LST as a boundary condition (as previously discussed in [10,78]). Figure A4 shows the daily average trend for the retrieved effective α PT with TSEB-2S at CT in 2015. Note in Figure A4 that α PT (defined for only the canopy source i.e., Kustas and Anderson [9]) maintains closer to the initial value of 1.26 during the dry summer period since TSEB-2S is simulating the canopy transpiration as a scattered broadleaved evergreen tree cover, which has extensive root system able to withdraw water even under drought conditions. The b coefficient, used to estimate the soil resistance to heat transport (i.e., R S ), was the only more empirically derived parameter that demonstrated a relatively large influence on model results. However, no significant relationship ( Figure A5A, r = 0.01) was obtained between R S S and errors (i.e., residuals) in modelled H and a very small relationship between R x (canopy boundary resistance) and errors in H was found ( Figure A5B, r = −0.19). The latter is likely caused by uncertainties in tree LAI, as R x is inversely proportional to LAI since the canopy is considered as a set of single-leaf resistors placed in parallel.
The input SA showed that uncertainties in LST and LAI both translate into uncertainties in modeled H. However, model results were found to be much more sensitive to LST compared to LAI. Gan and Gao [38] also found TSEB to be sensitive to biases in LST, where a 1K change is associated with median daily~12-25 W m −2 bias. A local LAI SA in Li et al. [39] added +-20% deviation to LAI and investigated the associated relative H change in TSEB. Their associated~3-8% bias in modeled H is similar to the results presented here, with a 20% change in LAI being associated with a median H bias of 6.1% with TSEB-DF. As such, the more comprehensive SA analysis presented here, considering parameter interaction in the global SA combined with a local input SA, were largely in line with results presented in different local and parameter specific SAs from the literature. These results will be useful for future studies, as we quantified the relative influence of the different modeling components in TSEB, which may be similar for other thermal-based energy balance models.
This work additionally investigated the partitioning of LE, which is a complex process to separate in ecosystems that have multiple vegetation layers. As shown, bulk (soil + vegetation) fluxes were well modeled in TSEB-2S for different years and towers. However, the partitioning of LE had greater uncertainty, with biases observed comparing modeled LE s to the lysimeter measurements ( Figure 10). Since total LE was well modeled, this indicates errors associated with the partitioning itself. This may be due to TSEB interpreting moderately stressed vegetation with a moist soil as fully transpiring vegetation with a dry soil. However, the partition of LE is largely controlled through LAI within TSEB [17]. The poorer performance in the LE partitioning could be due to the assumed vegetation cover of the two different seasons not properly depicting the complex vegetation characteristics observed, affecting net radiation partitioning, notably during the growing season when the tree layer is neglected. TSEB's relatively poor performance in partitioning LE was also observed over a vineyard in Kustas et al. [17]. As stated in Kustas et al. [17], more studies need to evaluate whether the poor partitioning is linked to uncertainties in input values (i.e., LAI) or biases caused by the modelling structure itself (i.e., initial potential canopy transpiration, radiation transfer).
As demonstrated, a simple adaptation to the modelling scheme in TSEB, depicting the phenological change in the vegetation source, were able to successfully simulate LE and H in a complex ecosystem. However, certain limitations are still present including the slight systematic underestimation of H during the grass growing season, particularly visible in the 2015 daily H time series (Figure 8). This is likely due to the modeling scheme in this period ignoring the effect of the tree canopy layer on the turbulent transport and hence on the calculation of the aerodynamic resistances. Compared to grasslands, tree canopies are more aerodynamically coupled to the atmosphere and hence their aerodynamic resistance is lower, resulting in that trees can dissipate heat more efficiently [43]. As discussed by El-Madany et al. [43], tree canopies in the Majadas experimental site were an additional H source, even though they tend to have a lower LST than the grass layer. As such, ignoring the tree canopy during the growing periods may not adequately represent the turbulent transport characteristics of the ecosystem, resulting in H underestimations, in some cases. On the other hand, neglecting the tree canopy during the growing season was supported by the evidence that the understory layer dominates LE in this site (as shown by Perez-Priego [42]). This may explain why TSEB-2S was largely capable to accurately model LE during this seasonal period for the different years assessed. Further adaptations to the TSEB model scheme may be needed for the more operational and larger scale use of this model, notably for similarly complex ecosystems such as TGEs. The differences between soil, grass, and tree layers should inherently be integrated within the modelling structure to robustly consider their different geometric, aerodynamic, and phenological properties and their resulting effect on energy fluxes.

Conclusions
When accounting for different phenological periods through an automatic vegetation seasonal change detection using a MODIS time series analysis, the TSEB model provided robust LE and H estimations for a three-layered heterogeneous and semi-arid TGE using a combination of satellite (LAI) and proximal (LST) sensors. This, in conjunction with the SA, confirms the important role that vegetation characteristics, notably its structure (i.e., h c ) and cover (i.e., f c and f g ), have on ecosystem level energy fluxes. The f c was the single most influential parameter on model performance, largely due to its role in characterizing vegetation clumping and how this, in turn, interacts significantly with other parameters. In addition, the uncertainties related to the traditionally remotely sensed derived inputs, notably LST, showed an important influence on output uncertainties.
The LE partitioning, between canopy and soil, showed a larger bias compared to the bulk fluxes. Based on this, further research should focus on the understanding of the radiation partitioning between the canopy and soil layers, and particularly, inherently accounting for the important differences between soil, grass, and tree layers of TGEs within the modeling structure. However, the relative simplicity of TSEB-2S has the potential to be applied with both satellite and/or airborne data over similarly complex landscapes in different regions. This may provide a simple solution to improve remote sensing-based flux products over semi-arid TGEs and other savanna-like ecosystems, which often assume a homogeneous vegetation layer.   Figure A2. Grass, tree, and ecosystem LAI time series for the CT tower for 2015. Note that tree LAI shown here refers to the effective (landscape) LAI, meaning that the tree-based measurement was weighted by its fractional cover (i.e., 0.2). Figure A2. Grass, tree, and ecosystem LAI time series for the CT tower for 2015. Note that tree LAI shown here refers to the effective (landscape) LAI, meaning that the tree-based measurement was weighted by its fractional cover (i.e., 0.2).

Appendix D
The roughness parameters (i.e., z0m and d0) for the tree canopy are estimated following the procedure described by Schaudt and Dickinson [55]. The reader is also referred to the source code of pyTSEB (https://github.com/hectornieto/pyTSEB/). The roughness parameters are estimated as follows: where λ frontal area index f /w , is the vegetation fractional cover (-), is the canopy width to height ratio (-).
The Lindroth [57] correction factors, fz and fd for z0m and d0m respectively, are then implemented as: