Adapting the thermal-based two-source energy balance model to estimate energy fluxes in a complex tree-grass ecosystem

The thermal-based Two-Source Energy Balance (TSEB) model has successfully simulated energy fluxes in a wide range of landscapes. However, tree-grass ecosystems (TGE) have notably complex heterogenous vegetation mixtures and 20 dynamic phenological characteristics presenting clear challenges to earth observation and modeling methods. Therefore, the TSEB model was adapted here to consider these significant seasonal changes. To ensure this and understand model dynamics, sensitivity analyses (SA) were conducted on both inputs (local SA) and parameters (global SA). Furthermore, a more physically based wind attenuation sub-model was applied and compared against the classical exponential wind attenuation law. The model was subsequently modified (TSEB-2S) and evaluated against eddy covariance (EC) flux 25 measurements and lysimeters over a TGE experimental site in central Spain. TSEB-2S vastly improved modeled fluxes decreasing the mean bias and RMSD of LE from 34 and 77 W m -2 to 4 and 56 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. The results presented here demonstrate the important role that vegetation, through its structure and phenology, has in controlling ecosystem level energy fluxes, which become important considerations for the modeling procedure. Additionally, 30 TSEB showed to be more sensitive to correctly partitioning incoming radiation, such as characterizing vegetation clumping, compared to accurately modeling the wind profile through the canopy or the aerodynamic roughness. 35 https://doi.org/10.5194/hess-2019-354 Preprint. Discussion started: 29 August 2019 c © Author(s) 2019. CC BY 4.0 License.

instrumental to understand model dynamics by providing insight in how parameters affect model outputs (Song et al., 2015). 70 Model uncertainty stems from three main sources: 1) errors associated with input data, 2) imperfection of model structure, and 3) uncertainty in model´s parameters (Jin et al., 2010;Migliavacca et al., 2012). This is particularly important in overparameterizated models that are unnecessarily complex with large uncertainties in parameters, which is a well-described issue (Beven, 1989;van Griensven et al., 2006). Different SA methods exist which are often distinguished between local and global techniques. Local methods compute the main response (1 st order) of the model with respect to changes in single 75 parameter values while keeping other parameters constant (i.e. one-at-a-time). Global methods evaluate the whole parameter space simultaneously and, thus, compute both the main effect (1 st order) and the interactions between parameters (2 nd order) to obtain the total parameter contribution (total order) to variability in model output. It is now widely argued that local SA techniques are unsuitable for complex non-linear models since there are often strong and significant parameter interactions (Pianosi et al., 2017;Rosolem et al., 2012). While other studies have investigated the sensitivity of specific parameters or 80 inputs within TSEB (e.g. Alfieri et al., 2019;Andreu et al., 2018;Gan and Gao, 2015;Li et al., 2005) or performed a SA to optimize TSEB (e.g. Diarra et al., 2017), a comprehensive SA for TSEB has not been discussed in the literature, especially for complex ecosystems, such as TGEs, where surface heterogeneity may potentially lead to increases in parameterization and complexities.

85
The main objective of this paper is to test whether a simple adaptation to TSEB, namely considering two distinct phenological and modeling periods throughout the year, can produce reliable estimations of energy fluxes for a spatially and temporally complex TGE. Prior to this, a SA is performed on parameters and inputs to quantify and pinpoint the different sources of uncertainties within the modeling procedure. A Sobol' global SA (Saltelli et al., 2010;Sobol′, 2001) is performed on the main parameters within TSEB combined with a local SA of the two main inputs: land surface temperature (LST) and 90 leaf area index (LAI). Two modeling structures, with a varying degree of complexity, are tested: one embedded with the classical Goudriaan (1977) wind attenuation formulation and, another using Massman et al. (2017)'s more physically-based wind attenuation scheme (section 2.2.2). The modified model results are evaluated against three independent eddy covariance (EC) systems, including partitioned LE , and lysimeter measurements located within the Majadas experimental site (Perez-Priego et al., 2017). The primary objective is to evaluate the simulated bulk fluxes (i.e. 95 LE) for the whole vegetation-soil system with a secondary goal to investigate the suitability of the LE partitioning, between LE s and LE c , estimated by TSEB.

Study Site
The TSEB model was applied to estimate energy fluxes in a TGE located in Majadas de Tiétar (39°56′24.68″N, 100 5°46′28.70″W) in central Spain (Casals et al., 2009;El-Madany et al., 2018). These types of ecosystems are prevalent, covering nearly 15% of the total Earth surface (Friedl et al., 2010), and are notably valuable in both an economic (i.e. livestock grazing) and ecological (i.e. biodiversity and carbon sequestration) sense. Majadas de Tiétar is a well-established experimental site composed of scattered oak trees, mostly Holm Oak (Quercus ilex. L.), above an herbaceous vegetation understory or grass layer. Holm Oak trees cover roughly 20% of the total land surface and stand at a mean height of 8m (El-105 Madany et al., 2018). The site is a managed semi-natural agroforested area (Spanish 'dehesa') with low-intensity grazing from livestock (< 0.3 cows ha -1 ). It lies within a continental Mediterranean climate region with mean annual temperature of 16.7 ºC and annual precipitation of about 650mm (with significant inter-annual variability) . The area is characterized by very hot and dry summer periods (June to September), with grass rapidly drying and senescing during these periods. The average grass LAI ranges roughly between 0.3 -3.0 m 2 m -2 throughout the year and can present high variability 110 in spring period (between 0.5 -2.5 m 2 m -2 ) due to its spatial heterogeneity Migliavacca et al., 2017).
Trees have developed extensive root systems enabling them to survive during long drought periods and, thus, have less temporal variability (mean tree LAI ranging between 1.39 -1.75 m 2 m -2 ). As such, the distinct survival strategies between grass and tree species allow for coexistence. Three EC towers are present within this experimental site. They are located relatively close to each other (< 650m, Fig. 1) with similar properties within their footprint, but belong to a large scale 115 manipulation experiment, where nitrogen was added to the northern tower (NT), nitrogen and phosphorus were added to the southern tower (NPT) and the central tower kept as a control (CT) Luo et al., 2018).  Kljun et al. (2015)) in early spring (lower left) and summer 120 (lower right

TSEB Model Overview
The TSEB model was first proposed in Norman et al. (1995), with important adjustments described in Kustas and Norman 125 (1999). Its main inputs are LST, derived from thermal infrared (TIR) radiation, vegetation structural properties (e.g LAI, canopy height) and meteorological forcing. 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 (eq. 1). (1)

130
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 (Norman 135 et al., 1995), and their relationship is not well established (i.e. Colaizzi et al., 2004). TSEB, thus, tackles this 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 eq. 2 taken from Norman et al. (1995): (2) 140 where is the fraction of vegetation observed by the TIR sensor at an angle and is mainly a function of LAI; is the vegetation canopy temperature (K); and is the soil surface temperature (K). Using this scheme, TSEB avoids the use of an empirical method to link radiometric and aerodynamic surface temperature, such as the use of excess resistance in SEBS (Su, 2002) or the use of hot and cold end member pixels in METRIC (Allen et al., 2007) and SEBAL (Bastiaanssen et al., 1998).
Using this two-layer approach, the energy balance is formulated in TSEB for each of the layers separately as follows: 145 (3) where is the net radiation (W m -2 ); is latent heat flux (W m -2 ); is the soil heat flux (W m -2 ); and subscript and refer to soil and vegetation canopy layers, respectively. Note that storage terms as well as the energy used for photosynthesis are neglected here. Radiative transfer and absorption through the canopy ( and ) is simulated through an exponential 150 radiation extinction function as described in chapter 15 of Campbell and Norman (1998), considering spectral differences in https://doi.org/10.5194/hess-2019-354 Preprint. Discussion started: 29 August 2019 c Author(s) 2019. CC BY 4.0 License.
shortwave and longwave radiation and direct and diffuse radiation as incorporated in Kustas and Norman (1999). H is derived using an in 'series' resistance network   (Fig. 2), and applying eq. 4 to 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 0 ); 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 165 m -1 ); is the aerodynamic resistance to heat transfer based on the Monin-Obukhov similarity theory. Refer to appendix A of Norman et al. (1995) for details on the series resistance scheme.
Since eq. 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 ( ) using, in this case, the Priestly-Taylor formulation:  ; is the fraction of vegetation that is green and hence actively transpiring (-); is the slope of the saturation vapor pressure curve at air temperature (kPa K -1 ); and is the psychrometric constant (kPa K -1 ). This initial assumption, where the canopy is transpiring without water stress, permits to solve all the systems of equations presented above (eq. 2 to 7). However, if the vegetation is stressed, the 175 Priestley-Taylor formulation will overestimate the transpiration of the canopy, which, in order to conserve the total surface temperature and energy balance in eq. 2 to 3, will result in unrealistic soil condensation (i.e. negative fluxes). As it is assumed that condensation does not occur during daytime convective conditions, an iteration procedure is applied that reduces until realistic soil fluxes area achieved (i.e. ). A more complete discussion on conditions that reduces is given in Anderson et al. (2005) and Li et al. (2005). For more implementation details, the reader is referred to the 180 source code (https://github.com/hectornieto/pyTSEB) and Norman et al. (1995).
TSEB was originally developed for homogeneous cover types, however, adaptations to the model framework have been implemented to better depict partial canopy cover. These include the addition of vegetation clumping index for radiation interception and transmission (Kustas and Norman, 1999) (section 2.2.1) and incorporating a more physically based withincanopy wind attenuation scheme, incorporating the effect of canopy structure, as proposed by Massman et al. (2017) (section 185 2.2.2). However, as Andreu et al. (2018) demonstrated, further adaptations to model structure and parameters may be necessary to simulate fluxes over complex land surfaces such as TGEs.

Radiation transmission in sparse vegetation
The structure and distribution of foliage in the vegetative layer has a significant impact on the dynamics of radiation interception and transmission through the canopy (Anderson et al., 2005;García et al., 2015;Kustas and Norman, 1999). 190 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 (Anderson et al., 2005). The original TSEB radiative transfer equations assume a randomly distributed homogenous (non-clumped) canopy . However, sparse vegetation is generally clumped and tends to intercept less radiation for the same LAI compared to vegetation randomly dispersed over the surface (Campbell and Norman, 1998;Kustas and Norman, 1999). As such, the clumping index (Ω) 195 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 (1999), TSEB estimates Ω as a function of the difference between the canopy gap fraction compared to a homogenous case using eq. 8 (section 2 in Kustas https://doi.org/10.5194/hess-2019-354 Preprint. Discussion started: 29 August 2019 c Author(s) 2019. CC BY 4.0 License. and Norman, 1999) and estimating the beam extinction coefficient assuming an ellipsoidal leaf angle distribution with eq. 9 (section 15.2 in Campbell and Norman, 1998). 200 Where is the is the clumping index when the vegetation canopy is viewed at nadir (-); is the beam extinction coefficient (-) based on a ellipsoidal leaf angle distribution (LAD) (Campbell and Norman, 1998); is the vegetation 205 fractional cover (-); F is the local LAI (LAI/ ); and is the leaf inclination distribution function chi parameter (-) (Campbell and Norman, 1998). Ω is also dependent on the solar zenith angle ( ) and is estimated using eq. 10 as described in section 15.13 of Campbell and Norman (1998).
where ); and is the vegetation´s width to height ratio (-). 210

Aerodynamic resistances and within-canopy wind speed profile
The semi-empirical derivations of and , as introduced in Norman et al. (1995) and Kustas and Norman (1999), are dependent on wind speeds below and within the vegetative canopy, respectively (eq. 11-12).

215
where is the wind speed just above the surface where the impact of soil roughness is minimal (i.e. ) (ms -1 ); is the wind speed within the canopy-air interspace at the height of momentum source/sink (ms -1 ); 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 (1999) and Norman et al.  (1995), based on the works of , Kondo and Ishida (1997) and McNaughton and Van Den Hurk (1995) (default values are 0.0025, 0.012 and 90 respectively). TSEB determines the two key roughness parameters, the zero-220 plane displacement height (d 0 ) and the aerodynamic roughness length for momentum transfer ( ), 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 (Campbell and Norman, 1998). In the case of the tree canopy layer, a different approach is used, which additionally considers the canopy shape and density´s effect on the roughness parameters. TSEB follows the procedure described in Schaudt and Dickinson (2000), which stems from the work of Raupach (1994) and Lindroth (1993), and is 225 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., 2012). However, Alfieri et al (2019) found that the TSEB modeled fluxes were largely insensitive to differences in estimated d 0 and obtained from various methods. and in eq. 11 and 12 are both derived from an observed wind speed above the canopy that is extrapolated through and below the canopy using a wind profile modeling scheme. The classical modeling approach within TSEB is based on the 230 exponential wind attenuation law (Goudriaan, 1977). However, this approach was largely designed for homogeneous and dense vegetation layers and may not be suitable for sparsely, clumped or non-uniformly structured vegetation that will likely impact the wind attenuation flowing through its porous medium . Recently, Massman et al. (2017) proposed a more physically-based wind attenuation scheme to take into account the effect of canopy shape and structure on wind speed, allowing for a smoother transition between partial to full canopy cover. Their main additional input is the 235 canopy foliage vertical distribution, based on the Plant (leaf + stem) Area Index and the relative canopy shape. A logarithmic profile is dominant near the ground surface and lower part of canopy, while a hyperbolic cosine profile dominates higher parts of the canopy, where its foliage would most affect wind speed . A drag coefficient (c d ) of foliage elements also characterizes the turbid medium canopy, which is usually considered to be 0.2 (Goudriaan, 1977;Massman et al., 2017). Both the Massman et al (2017) and Goudriaan (1977) wind profile schemes are incorporated within TSEB in this 240 study to evaluate on their respective influence on simulated fluxes.

Eddy covariance and biometerological 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) are larger than H due to issues related to the instrumentation (Foken et al., 2011), as done in, among other studies, Guzinski 250 et al. (2014) and Kustas et al. (2012). Measurements from 2015 of the CT tower were selected as the 'main' simulation period. This period was used to perform the SA (section 2.4.2) and to adapt the model for seasonality (section 2.4.4). Other towers (NT and NPT) and years (2016 and 2017) provided independent evaluations of model performance (section 2.4.5).
While TSEB is embedded with different methods to estimate ground heat flux (G) Santanello and Friedl, 2003), this study directly forced G measurements as a model input to limit uncertainty and noise in turbulent flux 255 estimations associated with errors in G. The weighted average of 8 soil heat flux plates represented G at the ecosystem level.
They were located both in open grass and below tree canopy and weighted to consider shadow effects throughout the day.
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 longwave radiation measurements from the 4-component radiometer CNR4 (Kipp & Zonen, Delft, Netherlands) as follows: 260 (13) where and are upwelling and downwelling longwave radiance; is the Stephan-Boltzman constant; is the surface emissivity; and is the vegetation fractional cover. The values of 0.99 and 0.94 in eq. 13 correspond to the emissivity for vegetation and bare soil respectively (Sobrino et al., 2005).

Vegetation biophysical measurements
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 (Fig. 1). The 500m by 500m pixel size was deemed adequate to characterize the tower footprint . The normalized difference vegetation index (NDVI) was derived using band 1 (red: 620nm-670nm) and band 2 (NIR: 841nm-876nm) as follows: 270 (14) An empirical relationship between NDVI and in-situ destructive grass LAI measurements was developed specifically for the Majadas experimental site (Appendix A, Fig.A1) using in-situ grass LAI measurements. These In situ grass LAI measurements were acquired at eleven plots (25mx25m) randomly placed around the CT (Mendiguren et al., 2015).
Destructive samples were collected in two 25cmx25cm quadrants within each plot and LAI was derived using gravitational 275 methods, where green and non-green elements were separated to compute both total and green LAI (as well as ). In order to obtain the ecosystem LAI adjusted for the tree canopy effect within the tower footprint, ground tree LAI measurements were incorporated to achieve a weighted average of total LAI (20% tree and 80% grass). 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 year. Tree LAI ranges between 1.39 and 1.75 and has 280 low inter-annual variability . Similarly, time series were obtained through in-situ destructive measurements from 18 field campaigns during the simulation periods (2015-2017) (appendix A, Fig.A2). A simple linear interpolation between sampling dates was performed to obtain a continuous daily time series.

Model Simulations and Evaluation
Prior to conducting the simulations with TSEB, both a global and local SA was performed on the main parameters and 285 inputs, respectively. Both SAs were conducted using the default TSEB model (section 2.4.1). 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 (2.4.3). The model is then adapted considering two distinct seasonal periods, each having a dominating vegetation layer (section 2.4.4). Parameters related to vegetation structure are adapted for each period in order to characterize the assumed prominent vegetation (i.e. either as tree 290 or grass), as well as taking into account the results from the SA. All simulations are evaluated for daytime fluxes (i.e. when S dn > 25 Wm -2 ) since TSEB is designed to model fluxes for daytime conditions with remote sensing data.

Default TSEB model configuration
In this configuration, the single vegetated layer was parameterized as a mix between tree and grass, where the vegetation cover was dominated by grass (f c = 0.8) but the turbulent resistances were assumed to be largely affected by the sparse tree 295 layer . As such, parameters related to vegetation resistance and roughness were configured using tree characteristics (e.g. h c = 8m, see section 2.2.2). Table 3 summarizes the parameter values for this default model. Two model structures with different wind attenuation schemes were tested: one embedded with Goudriaan (1977)

Sobol´ Global Parameter Sensitivity Analysis
Variance-based SA methods are now more prominently used (Rosolem et al., 2012). They decompose the total variance between simulated and observed data into various parts determining the contribution of the different parameters and their 305 combined interactions on total output variance. Multiple model simulations are required with different parameter sets to quantify the model output variance in regard to the variation in the parameter space.
Among them, the Sobol' SA method is widely preferred as it is able to compute the 1 st , 2 nd and total order sensitivity indices.
The main disadvantage of the Sobol' method is the high computational cost, where many simulations are required to obtain robust results within a sufficient confidence level .To apply this methodology, a selection of model parameters and their 310 respective bounds are configured. Based on this, parameters sets for n simulations are generated using the Sobol sequence (Saltelli et al., 2010;Sobol′, 2001), a quasi-random method with quasi-Monte-Carlo integration, which typically samples bounded space more uniformly than completely random sequences .
where is the total output variance; is the portion of variance contributed by parameter (also known as the first order variance or main effect); and is the portion of variance attributed through the interaction between and and so on. As such, the first-order sensitivity (S i ), second-order sensitivity (S ij ) and total order sensitivity (S Ti ) indices are calculated as: where is the average variance depicted when all parameters except for vary (i.e. is kept fixed); represents the contribution of both the direct (1 st order) and indirect (sum of 2 nd order) effects of on the total variance, .

Input Local Sensitivity Analysis
Errors associated to input data are another important source to the total output uncertainty. Since input data is fluctuating at 345 every time step, a different method is needed compared to the method used for the less time-varying parameters described above Random white noise, mimicking potential errors, over a uniform distribution of +-3K and +-0.35 m 2 m -2 were added to 'observed' LST and LAI, respectively. Sensitivity indices (SI) are then computed based on the partial derivatives of the output results with respect to change in input caused by the implemented random noise, using eq. 21 adapted from van Griensven et al. (2006). where M is model output; refers to the different model inputs; 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 tested with both canopy wind attenuation sub-models (i.e. TSEB G -DF and TSEB M -DF) 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 355 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 hour in which the deviation occurs (i.e. a perturbation in LST during night time will likely impact H differently compared to a daytime 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 are derived for each 360 daily time step using the absolute change in inputs with the resulting absolute change in output.

End Member Simulations
In addition to the default simulations, two 'end member' scenarios were conducted as limiting case studies. 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, clumped evergreen broadleaf forest (hereafter as 365 TSEB tree ). The estimation of the roughness characteristics differed for each scenario as described in section 2.2.2, parameters related to the resistances in TSEB tree were changed to mimic sparsely vegetated and rough soil conditions , where specific parameter values are shown in table 3. These scenarios were performed on the CT tower during 2015.

Two-Season modeling approach 370
We propose here 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 considerable changes observed within the different seasonal periods. The differences in vegetation cover alter the turbulent conditions (i.e. roughness), radiation transmission (i.e. ) and, hence, the surface energy balance. The seasonal transition dates were estimated using an asymmetric gaussian 375 filter over the 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 13 th and ends on October 24 th (refer to appendix A, Fig.A3 for all transitions dates used) The dry period is parameterized as a tree-soil system, assuming the grass has senesced. The is set to 0.2 (Oak trees 380 represent roughly 20% of surface cover) and the roughness and resistances are calculated using tree characteristics (see section 2.2.2). By contrast, the growing period is parameterized ignoring the tree canopy and assuming a grass-soil system since the vegetation is dominated by grass (roughly 80%) during this period (see table 3). This assumption is well supported by the evidences that the understory layer dominates LE in this site . A variable , based on field measurements (appendix A, Fig.A2) is incorporated in the grass-soil configuration during the non-summer/growing 385 periods. During the dry period, is constant and set to 0.9 since the trees have relatively constant foliage throughout the year, and less than 1 to account for branches and other non-green canopy material.

Model evaluation
To validate and benchmark model performance, the seasonally adapted TSEB-2S was independently tested on two nearby EC towers within the study site (NT and NPT). El-Madany et al. (2018) demonstrated that, in the pre-fertilization phase, 390 differences in biophysical surface properties caused differences in sub-hourly fluxes between the three nearby EC towers, identifying some level of spatial variability between the three sites. In the fertilization phase (2015-2017), differences in energy and water partitioning were observed (El-Madany et al., in preparation), and therefore it is interesting to test TSEB-2S against the sites considering that they are structurally similar but show different response in terms of ecosystem functioning. The model was also benchmarked for different years (2016 and 2017) with different intra-annual dynamics. 395 Therefore, an attempt was made to validate the TSEB-2S with independent model simulations for different spatial (three tower sites) and temporal (three different years) characteristics.
Additionally, an evaluation of model LE partitioning, between LE c and LE s , was performed. LE s was evaluated against independent lysimeter LE measurements (LE lys ) at CT during dry summer periods, when 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. (2017). In addition, the TSEB-2S simulated LE c is benchmarked against that derived from a physiologically-based water flux partitioning approach, which allows for the separation of LE measured by EC into transpiration and evaporation components. The method was developed and validated in the same study site .

Sobol´ parameter SA
The parameter global SA showed that has the largest impact on model results, both as the main effect (S i ) and through its interactions with other parameters (S T ), regardless of the wind attenuation scheme applied (Fig. 3). Parameters related to

Input local SA 425
Uncertainties related to LST and LAI showed significant impacts on model results (Fig. 4). Simulated H is particularly sensitive to deviations in LST, where differences of a unit change of LST (∆ 1K) is associated to a median 15.1% and 15.4% change in modelled H for TSEB G -DF and TSEB M -DF, respectively. In terms of relative change, a 1% variation in LST (i.e ∆ 3K for a 300K surface) is associated with a median change of 45.7% and 46.4% in modelled H for TSEB G -DF and TSEB M -respectively. In terms of percent change, a 1% change in LAI is associated with a median change of 0.31% and 0.42% in modeled H for TSEB G -DF and TSEB M -DF, respectively.

TSEB-DF, TSEB grass , TSEB tree
Model results of default TSEB using Goudriaan´s wind attenuation scheme (TSEB G -DF) vastly underestimate H (bias: -39 W m -2 ) and, consequently, overestimate LE (bias: 34 W m -2 ) as compared with observed EC data. As shown in Fig. 5, high 440 errors are observed throughout (RMSD of H: 76 W m -2 ) but errors stem particularly from the very significant underestimations of H during the hot and dry summer period (Fig. 5). The end member simulations demonstrate the 'boundary conditions' when, with the same inputs, the simulations completely ignore one of the vegetation layers throughout the year. Predictably, TSEB grass underestimates H even more significantly than TSEB G -DF, while, by contrast, H is significantly overestimated in TSEB tree . Modeled LE from TSEB tree resulted in many insignificant fluxes (i.e. 0 W m -2 ) even when the tower observed significant LE (Fig. 6). This indicates that, in many cases, an overestimated modeled H left no available energy flux for LE from the energy balance residual. 450

TSEB-2S
To improve the seasonal depiction and related changes to the ecosystem, the TSEB model was adapted using a two-season modeling approach as proposed in section 2.4.4. The global SA results showed that parameters related to the vegetation cover (i.e. and ) and structure (i.e ) had the most influence on model performance. Since these are physical and measureable parameters, these were simply set in each seasonal period according to the assumed dominant vegetation. The 460 resistance coefficients b and c, used in the formulation of R s , were increased for the dry period, as was successfully tested and implemented in Kustas et al. (2016) (Table 4). In Kustas et al. (2016), b and c within TSEB were given a value of 0.0605 and 0.0038, respectively based on literature (Kondo and Ishida, 1997; to consider rough soil and partially vegetated surfaces. Other parameters, such as , which is presumably different for grasses and trees, were not modified since they showed little sensitivity (Fig. 3 Simulations of LE and H for 2015 vastly improve using TSEB-2S with the Goudriaan (1977) wind attenuation scheme compared to TSEB G -DF (Fig. 7). It yields a decrease in RMSD from 77 W m -2 and 76 W m -2 to 56 W m -2 and 52 W m -2 for 470 LE and H, respectively. The seasonal average daily H is well depicted from model results, particularly during the summer/dry period. The winter and spring period (roughly January to May) show a slight systematic H underestimation.
The two different wind attenuation models were applied within TSEB-2S but no significant changes occurred on overall model results when using the different wind profiles (Fig. 7). The magnitude of errors maintains relatively equal, only a slight decrease in RMSD (56 Wm -2 to 54 Wm -2 ) and slightly larger underestimation of H (bias of -3 to -6 W m -2 ) was 475 observed using the Massman et al. (2017) wind profile. Similar outcomes occur when comparing the different wind attenuation schemes across different periods, years and sites (Table 5). The largest difference in model performance occurs during the summer, particularly for NPT in 2015, when the canopy shape of the 'tree-soil' model scheme has a greater impact on the wind attenuation, and thus, H. However, for clarity purposes and since results barely change when using the   Table 5. Errors (RMSD and bias) are expressed in W m -2 for modeled H using different wind attenuation schemes in TSEB-2S at CT, NT and NPT;and years 2015, 2016; divided between growing and dry/summer periods. The 'Change' column refers to the absolute 495 difference in metrics between using TSEB-2S with Goudriaan (1977) compared to Massman et al. (2017) (i.e. a change of -5 to -6 Wm -2 bias is a change of +1)

Year
Season Tower

Spatial and temporal evaluation 500
TSEB-2S was additionally tested for the NT and NPT EC towers in 2016 and 2017. The results from the other sites and years are very satisfactory (Fig. 8), with errors are comparable to the 2015 simulation at the CT. RMSD for H ranges between 50 and 60 W m -2 between all validating simulations. The model demonstrates robustness, able to simulate adequately for different temporal periods and conditions (i.e. nutrient fertilization effects on the NT and NPT sites).

LE partitioning 510
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 . It was additionally compared against an EC LE partitioning method . Results indicate a systematic LE s underestimation during the dry period for the three years analyzed (Fig. 9). The LE s overestimation is particularly significant for 2015 (bias: -28 W m -2 ). In 2017, modeled LE s is most aligned 515 against the fluxes measured by the lysimeter with less, while still significant, errors (RMSD: 25 W m -2 ; bias: -13 W m -2 ).
TSEB-2S, while having a systematic offset, was able to capture the temporal dynamics very well, as illustrated by effectively capturing the evaporation peaks caused by the rare summer precipitation events in 2016 and 2017.  (Fig. 10). Throughout most of the year, LE c,TSEB-2S has a positive bias compared to LE c,PP . However, patterns significantly and drastically change during the dry period, where 525 LE c,TSEB-2S is more in line with LE c,PP (less significantly in 2017). Since LE c,PP contains both tree and grass contribution, the LE c,PP peak is more sustained than LE c,TSEB-2S . This is explained by the successive influence of the grass (DOY~90) followed by the trees in late-May/early-June(DOY~150) and then have a gradual decline during the dry-down as modulated by water stress. By contrast, LE c,TSEB-2S peaks earlier, around late-April, and proceeds to decline very rapidly when summer begins, coinciding with the change in model configuration between the grass-soil system to the tree-soil system. LE c,TSEB-2S also 530 shows a steady increase beginning towards the end of summer (roughly September, DOY~250) until the autumn re-greening (roughly late-October, DOY~300), which is less significantly observed in LE c,PP .

Discussion
TSEB-2S significantly improved model performance in simulating LE and H compared to the default TSEB (Fig. 7). The simple assumption of two separate phenological periods, one dominated by a grass-soil system and the other dominated by a tree-soil system, permits the successful application of a two-layer model in an essentially three (soil-grass-tree) layer 540 ecosystem. The results also demonstrate that significant changes to the vegetated surface have a very important influence on the surface energy balance. This confirms the importance of vegetation characteristics in controlling ecosystem level energy fluxes, where seasonal dynamics and phenology of vegetation are key considerations for land-atmospheric modeling. The TSEB-2S model demonstrated robustness, being able to successfully simulate different intra-annual dynamics for various years and towers with distinct surface conditions from the nutrient fertilization experiment (Fig. 8). Results and the 545 associated magnitudes of errors for all model runs are similar to the error bounds found in other energy balance model studies (e.g. Andreu et al., 2018;Boulet et al., 2015;Gan and Gao, 2015;Gonzalez-Dugo et al., 2009;Guzinski et al., 2014;Kustas et al., 2016;Kustas and Norman, 1999;Timmermans et al., 2007).
The global SA demonstrated that TSEB is most sensitive to . This is largely due to being a key input to estimate the Ω (Kustas and Norman, 1999), which in turn affects the radiation interception and partitioning between the mixed vegetated 550 and soil surfaces. The estimation of Ω is mostly based on the transmission through the vegetated layer where is used to obtain a local LAI (LAI/ ) 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. (2005), who found that the Ω uncertainties had a significant impact on flux outputs from TSEB. They tested incremental values, between 0.1 and 1, that resulted in sharp H changes, particularly between higher values, an indication of a high sensitivity for this parameter. The f g is more sensitive compared to 555 even though both parameters are part of the initial canopy transpiration estimate from the Priestly-Taylor formulation (eq. 7). This is the case since the latter is merely an initial value and is iteratively changed through the modeling scheme, if the resulting soil evaporation flux is unrealistic. The b coefficient was the only more empirically derived parameter that demonstrated relatively large influence. This was more significantly observed within TSEB G -DF since the Massman et al.
(2017) wind profile decreases asymptotically to 0 near the soil therefore the term has less influence than 560 for the estimation (eq. 11) compared to the Goudriaan (1977) wind profile where tends to 0.
The input SA showed that uncertainties in LST and LAI both translate relatively significantly into uncertainties in the modeled H. Gan and Gao (2015) also found TSEB to be quite sensitive to biases in LST, where a 1K change is associated with median daily ~12 -25 W m -2 bias. The sensitivity of modeled H to LAI uncertainties were slightly more pronounced with TSEB M -DF compared to TSEB G -DF (Fig. 4). This is the case since LAI is additionally used in the Massman et al. 565 al. (2005) added +-20% deviation to LAI and investigated the associated relative H change in TSEB. Their associated ~3-8% bias in modeled H is largely in line with results presented here, with a 20% change in LAI associated with a median H bias of 6.1% and 8.4% for TSEB G -DF and TSEB M -DF, respectively.
Overall model results using the different wind attenuation schemes were very similar, without significant change in modeled 570 fluxes. For some years, the use of Massman et al. (2017) marginally decreased the resulting errors, while for other years it marginally increased errors (Table 5). This is in line with Nieto et al. (2019), who investigated the effect of the same two wind attenuation schemes on TSEB over a grapevine in a semi-arid region in California. Very little difference in model performance was found between the annual error statistics. However, Nieto et al. (2019) found a relatively significant reduction of errors in modeled H with Massman et al. (2017) during the spring seasonal period, when the grapevine canopy 575 was most vertically clumped. Similarly, Table 5 shows that Massman et al. (2017) tends to have a more significant effect during the summer period. This is due to the model configuration as the growing period is parameterized as a homogenous grass layer and, thus, the canopy shape should not have a large impact on the wind attenuation. However, errors tend to increase in the summer for Massman et al. (2017), notably in 2015 for NPT. The overall benefit of using Massman et al.
(2017) within TSEB is not conclusive through these results. It does not significantly impact overall model results, unless a 580 very large influence of canopy vertical density is expected. These results, along with the SA, indicate that correctly partitioning incoming radiation, such as characterizing vegetation clumping, has a much stronger influence on model performance compared to accurately modeling the wind profile through the canopy.
As shown, bulk (soil + vegetation) fluxes were well modeled in TSEB-2S for different years and towers. However, the partitioning of LE, between LE s and LE c , has greater uncertainty, with biases observed compared to the lysimeter 585 measurements (Fig. 9) and compared to a physiologically-based EC partitioning method (Fig. 10). Since total LE is well modeled, this indicates errors associated with the partitioning itself. The partition of LE is largely controlled through LAI within TSEB , as evidenced in how the LE c time series (Fig. 9) roughly follows the trend of the LAI input time series shown in appendix A (Fig. A3). The rapid decrease in LAI as the summer begins results in a rapid change in LE c, TSEB-2S , in contrast to the gradual decrease of LE c,PP . This can also be attributed to TSEB-2S´s sudden change in model 590 configuration (i.e. from grass to tree dominated).The transition periods between model configurations had the largest biases (roughly around DOY 120 and DOY 300) while the dry period, when the model considers only the tree-soil system, produced the best simulated LE partitions. TSEB's relatively poor performance in partitioning LE were also observed over a vineyard in Kustas et al. (2019). However, Xu et al. (2016) reported a relatively good performance of TSEB's midday LE partitioning over an irrigated cropland against observations using the isotope approach. Colaizzi et al. (2014) demonstrated 595 that using the Penman-Monteith equation to derive initial canopy transpiration in TSEB resulted in a better agreement with LE s and LE c compared to the Priestley-Taylor approach. Colaizzi et al. (2014) also reported the strong influence of vapour pressure deficit over LE c in an irrigated cotton field within an advective, semi-arid climate, which was captured most https://doi.org/10.5194/hess-2019-354 Preprint. Discussion started: 29 August 2019 c Author(s) 2019. CC BY 4.0 License. successfully using the Penman-Monteith approach. As stated in Kustas et al. (2019), more studies need to evaluate whether the poor partitioning is linked to uncertainties in input values (i.e LAI) or biases caused by the modeling structure itself (i.e. 600 initial potential canopy transpiration, radiation transfer).
As demonstrated, a simple adaptation to the modeling scheme within TSEB, depicting the seasonal changes in the ecosystem, 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 2015 (Fig. 7). This is likely due to the modeling scheme in this period ignoring the effect of the tree canopy layer on the turbulent 605 transport and hence on the calculation of the aerodynamic resistances. Compared to grass layers, 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 . As discussed in El-Madany et al. (2018), 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 aerodynamic resistance 610 characteristics of the ecosystem resulting in H underestimations, in some cases. Further adaptations to TSEB model may be needed for the more operational and larger scale use of this model, notably for different complex systems such as TGEs. The differences between soil, grass and tree layers should inherently be integrated within the modeling structure to robustly consider their different geometric, physical and phenological properties and their resulting effect on energy fluxes. This was also hinted in Andreu et al. (2018), who incorporated different layers in a modified Goudriaan (1977) wind speed profile 615 scheme to consider differences between tree and grass canopy layers for TGEs.

Conclusions
When accounting for different phenological periods, the TSEB model can provide robust LE and H estimations for a threelayered heterogeneous and semi-arid TGE. This 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 620 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, LAI and LST, showed an important influence on output uncertainties.
The more physically-based wind attenuation scheme (Massman et al., 2017) had little impact on the overall modelled fluxes, especially compared to factors affecting the radiation budget between the two sources (i.e vegetation clumping). In addition, 625 the LE partitioning, between canopy and soil, showed larger bias compared to the bulk fluxes, particularly during the transition towards the grass senescence period and the prominent contribution of the tree layer. 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.