The SPARSE model for the prediction of water stress and evapotranspiration components from thermal infra-red data and its evaluation over irrigated and rainfed wheat

Evapotranspiration is an important component of the water cycle, especially in semi-arid lands. A way to quantify the spatial distribution of evapotranspiration and water stress from remote-sensing data is to exploit the available surface temperature as a signature of the surface energy balance. Remotely sensed energy balance models enable one to estimate stress levels and, in turn, the water status of continental surfaces. Dual-source models are particularly useful since they allow derivation of a rough estimate of the water stress of the vegetation instead of that of a soil–vegetation composite. They either assume that the soil and the vegetation interact almost independently with the atmosphere (patch approach corresponding to a parallel resistance scheme) or are tightly coupled (layer approach corresponding to a series resistance scheme). The water status of both sources is solved simultaneously from a single surface temperature observation based on a realistic underlying assumption which states that, in most cases, the vegetation is unstressed, and that if the vegetation is stressed, evaporation is negligible. In the latter case, if the vegetation stress is not properly accounted for, the resulting evaporation will decrease to unrealistic levels (negative fluxes) in order to maintain the same total surface temperature. This work assesses the retrieval performances of total and component evapotranspiration as well as surface and plant water stress levels by (1) proposing a new dual-source model named Soil Plant Atmosphere and Remote Sensing Evapotranspiration (SPARSE) in two versions (parallel and series resistance networks) based on the TSEB (Two-Source Energy Balance model, Norman et al., 1995) model rationale as well as state-of-the-art formulations of turbulent and radiative exchange, (2) challenging the limits of the underlying hypothesis for those two versions through a synthetic retrieval test and (3) testing the water stress retrievals (vegetation water stress and moisture-limited soil evaporation) against in situ data over contrasted test sites (irrigated and rainfed wheat). We demonstrated with those two data sets that the SPARSE series model is more robust to component stress retrieval for this cover type, that its performance increases by using bounding relationships based on potential conditions (root mean square error lowered by up to 11 W m from values of the order of 50–80 W m), and that soil evaporation retrieval is generally consistent with an independent estimate from observed soil moisture evolution. Published by Copernicus Publications on behalf of the European Geosciences Union. 4654 G. Boulet et al.: SPARSE model for the prediction of evapotranspiration from TIR data


Abstract.
Evapotranspiration is an important component of the water cycle, especially in semi-arid lands. A way to quantify the spatial distribution of evapotranspiration and water stress from remote-sensing data is to exploit the available surface temperature as a signature of the surface energy balance. Remotely sensed energy balance models enable one to estimate stress levels and, in turn, the water status of continental surfaces. Dual-source models are particularly useful since they allow derivation of a rough estimate of the water stress of the vegetation instead of that of a soil-vegetation composite. They either assume that the soil and the vegetation interact almost independently with the atmosphere (patch approach corresponding to a parallel resistance scheme) or are tightly coupled (layer approach corresponding to a series resistance scheme). The water status of both sources is solved simultaneously from a single surface temperature observation based on a realistic underlying assumption which states that, in most cases, the vegetation is unstressed, and that if the vegetation is stressed, evaporation is negligible.
In the latter case, if the vegetation stress is not properly accounted for, the resulting evaporation will decrease to unrealistic levels (negative fluxes) in order to maintain the same total surface temperature. This work assesses the retrieval performances of total and component evapotranspiration as well as surface and plant water stress levels by (1) proposing a new dual-source model named Soil Plant Atmosphere and Remote Sensing Evapotranspiration (SPARSE) in two versions (parallel and series resistance networks) based on the TSEB (Two-Source Energy Balance model, Norman et al., 1995) model rationale as well as state-of-the-art formulations of turbulent and radiative exchange, (2) challenging the limits of the underlying hypothesis for those two versions through a synthetic retrieval test and (3) testing the water stress retrievals (vegetation water stress and moisture-limited soil evaporation) against in situ data over contrasted test sites (irrigated and rainfed wheat). We demonstrated with those two data sets that the SPARSE series model is more robust to component stress retrieval for this cover type, that its performance increases by using bounding relationships based on potential conditions (root mean square error lowered by up to 11 W m −2 from values of the order of 50-80 W m −2 ), and that soil evaporation retrieval is generally consistent with an independent estimate from observed soil moisture evolution.

Introduction
Evapotranspiration is an important, yet difficult to estimate (Jasechko et al., 2013), component of the water cycle, especially in semi-arid lands. Its quantification is crucial for a sustainable management of scarce water resources. The recent development of remote-sensing products and data assimilation methods has led to a new era in the use of remote-sensing data in the various spectral domains to derive more robust estimates of evapotranspiration at various spatial scales (Crow et al., 2008;Olioso et al., 2005). Amongst those products, surface temperature provides access to a rough estimate of water stress. Indeed, moisture-limited evapotranspiration triggers an increase in surface temperature above a theoretical equilibrium value in unstressed conditions (Amano and Salvucci, 1997;Boulet et al., 2007). Most algorithms based on the use of a remotely sensed surface temperature evaluate a total latent heat flux corresponding to the sum of the evaporation and the transpiration components: they are named "single-source models". Total latent heat flux representing the whole surface is derived as the residual term of the surface energy balance at the time of satellite overpass (Kalma et al., 2008). Single-source models require a method to relate the temperature at the aerodynamic level and the surface temperature obtained by remote sensing (Matsushima, 2005). It is very often based on an additional resistance term or kB −1 (Carlson et al., 1995;Verhoef, 1997) that is heavily parameterized. Even though the use of single-source models is widespread, dual-source models are particularly useful because they allow retrieval of separate estimates of evaporation and transpiration. Those components are particularly needed for ecohydrological or agrohydrological applications (irrigation management, water stress detection, etc.). Moreover, dual-source models provide a more realistic description of the main water and heat fluxes, even if the vegetation is seen as a single "big leaf" and the soil as a single "big pore" (Kustas et al., 1996). This is especially true for sparse vegetation, when commonly used scalar profiles within the canopy no longer apply. It also avoids the use of a parameterized kB −1 (Kustas and Anderson, 2009).
Beyond evapotranspiration, estimating water stress is also important for inferring the surface water status and the root zone soil moisture level (Hain et al., 2009). Water stress can be obtained for the surface as a whole by combining the simulated latent heat flux and the potential latent heat flux, i.e. the theoretical value of the latent heat in current climatic conditions if the surface was still undergoing stage one (unstressed) evapotranspiration (Lhomme, 1997). Dualsource energy balance models allow derivation of a rough estimate of the water stress, but of the vegetation instead of a soil-vegetation composite. They also provide an estimate of the climate-controlled and moisture-limited soil evaporation rates. Such frameworks use as input data either the component surface temperatures (e.g. soil and vegetation components retrieved from directional surface temperature data; Jia et al., 2003, or Colaizzi et al., 2012 or a single soilvegetation composite surface skin temperature. For the former, there is no current operational satellite that offers estimates of temperatures at two contrasted view angles with a very small interval between both acquisitions, even though the soon to be launched Sentinel-3 mission will have such capability (Donlon et al., 2012). For the latter, the TSEB model proposes a realistic underlying assumption to downsize the number of unknowns from two (evaporation E and transpiration T ) to one (E or T ;Norman et al., 1995). The TSEB model assumes that in most eco-or agro-systems vegetation has access to enough water in the root zone to transpire at a potential rate, so that a modelled potential transpiration rate is a valid first guess estimate for T . This assumption implies that, if vegetation stress is not properly taken into account, the resulting evaporation will decrease to unrealistic levels (negative fluxes) in order to maintain the same total surface temperature, so that a retrieved negative evaporation is a good witness of plant water stress. This assumption is sometimes misleading, and we propose studying its limits.
The original version of TSEB (Norman et al., 1995) provides two algorithms to describe the soil-vegetationatmosphere interactions, representing, respectively, the "patch" and "layer" approaches following the terminology proposed by Lhomme et al. (2012). In the "layer" approach, one assumes that the air is well mixed within the canopy space so that air temperature at the aerodynamic level is rather homogeneous. The vegetation layer completely covers the ground and prevents the soil from interacting directly (in terms of radiation and turbulent heat transfer) with the atmospheric reference level: soil and vegetation heat sources are fully coupled through a resistance network organized in series ( Fig. 1). In the "patch" approach, soil and canopy sources are located side by side, and the soil interacts directly with the air above the canopy. There is a possible lateral gradient in air temperature around the aerodynamic level even though heat transfer around the canopy is associated with the same momentum transfer: soil and vegetation heat sources are thermally uncoupled and fluxes are computed with two parallel resistance schemes. In the original TSEB version, total net radiation is split into soil and vegetation components according to a simple Beer-Lambert law. Several improvements have been proposed later on and implemented in various TSEB versions. Amongst them, one can mention the development of a more complex net radiation scheme, with an initialization of soil and vegetation temperatures in separate formulations of the net radiation of the soil and the canopy or the use of an incremental decrease in a transpiration efficiency (Kustas and Norman, 1999; it corresponds roughly to the ratio between the actual and potential transpiration rates and matches the definition of the efficiency used in the present work). The TSEB rationale has been translated into several algorithms, with the possibility of using directional radiative temperatures (Kustas and Norman, 1997), day-night temperature difference (Guzinski et al., 2013;Nor- man et al., 2000), correcting for clumping effects in sparsely vegetated areas (Anderson et al., 2005), and finally by taking into account a Penman-Monteith formulation for potential transpiration (Colaizzi et al., 2012).
Here, we propose revisiting the "layer/series" and "patch/parallel" formulations in order to build a new model based on the same rationale that provides the foundation for all TSEB model versions.
First, we build on the statement by Colaizzi et al. (2012) that, in semi-arid lands, it is more relevant to use a resistance scheme based on a Penman-Monteith expression instead of the Priestley-Taylor equation, so that adiabatic exchanges are explicitly described. The most common value of the Priestley-Taylor coefficient (close to 1.3) has indeed been challenged for natural vegetation and sites with strong vapour pressure deficit values where root zone moisture does not limit transpiration (Agam et al., 2010). According to Colaizzi et al. (2014), potential transpiration using the Penman-Monteith equation showed better performances compared to the Priestley-Taylor equation. In particular, these authors showed a consistent underestimation of T and overestimation of E when using a Priestley-Taylor formulation with the classical 1.3 coefficient, even if total evapotranspiration was similar for both models.
Second, since in the layer approach the vegetation is a semi-infinite cover overlaying the ground, it appears more consistent that this version of the model takes into account not only the soil-vegetation interactions of the turbulent fluxes, but also of the radiative fluxes. Conversely, in the patch approach there is no radiation exchange between the soil and the vegetation patches. This is achieved for the series model through a multiple-reflection description between the soil and the overlaying vegetation cover in order to stick more closely to the patch and layer representations schematized in Fig. 1.
Based on those studies, we propose a generalization of the TSEB model (named SPARSE: Soil Plant Atmosphere and Remote Sensing Evapotranspiration model) as a linearization of the full set of energy budget equations and the Choudhury and Monteith (1988) and Shuttleworth and Gurney (1990) expressions of the aerodynamic resistances. The series model is very close to the soil-plant-atmosphere interface of the SiSPAT model (Braud et al., 1995). The full set of equations can be solved either in prescribed conditions (for example, in fully stressed or potential conditions) to compute transpiration and evaporation rates for given stress levels, or in retrieval mode, identically to TSEB. In that case, stress levels are deduced from a known (observed) surface temperature. We propose a third improvement to the existing TSEB model versions, which is similar to what is done in a postprocessing step in the single-source SEBS model (Su, 2002). It consists in bounding each retrieved individual flux component (T , E) by its corresponding potential level deduced from running the model in prescribed potential conditions. Indeed, transpiration can be above its potential level when there is a strong coupling between the soil and the vegetation through conditions at aerodynamic level (stability correction notably): maximum transpiration for a plant surrounded by very dry bare soil is increased above the potential transpiration rate as computed in a fully wet environment. This coupling might be excessive and a potential transpiration of a wet environment is an interesting baseline to assess excess in this coupling.
The main objective of the paper is twofold.
1. Describe the SPARSE model, evaluate it against in situ data and relate its performance to those of the "patch/parallel" and "layer/series" TSEB model formulations, with a focus on the potential gain in robustness obtained when limiting evaporation and transpiration outputs by their corresponding potential rates derived from SPARSE.
2. Test the retrieval capacities of both "patch/parallel" and "layer/series" versions of the model, not only for total evapotranspiration as well as its components (soil evaporation and transpiration), but also for water stress, first with synthetic data (simulation experiment) and second with in situ data collected over two wheat fields in a semi-arid climate, one irrigated and one rainfed. The purpose of the simulation experiment is specifically to test the limits of the underlying first guess assumptions of SPARSE, which are identical to those used in most TSEB versions.

SPARSE system of equations
The SPARSE model computes the equilibrium surface temperatures of the soil (T s ) and the vegetation (T v ) at the meteorological time step as a signature of the energy budget equations of each source. Five main equations are solved simultaneously. The first two express the continuity (series version) or the summation (parallel version) of the latent and sensible heat fluxes from the soil and the canopy to the aerodynamic level and above, the third and the fourth represent the energy budget of the soil and the vegetation, and the fifth describes the link between the radiative surface temperature T rad and its two component temperature sources (soil T s and vegetation T v ). Two versions are derived, which can be regarded as fully coupled (series) and fully uncoupled (parallel) soilvegetation-air exchanges (Fig. 1). This corresponds to (respectively) the "layer" and "patch" approaches described in Lhomme et al. (2012). However, the interpretation of the situations for which one or the other approach is valid differs between TSEB and Lhomme et al. (2012). In TSEB, both soil and vegetation patches share a common surface boundary layer (and therefore the same aerodynamic resistance from the aerodynamic level to the reference level), but the patch representation allows definition of different aerodynamic temperatures at the aerodynamic level over the soil and the vegetation. As pointed out by Lhomme et al. (2012), the patch representation should in theory only apply to patches large enough to develop different surface boundary layers, e.g. fallow fields amongst wetter and taller vegetated areas rather than bare soil patches even a few metres large. Here, we keep the TSEB assumption for our parallel version and assume that the wind profiles above the aerodynamic level in the canopy and above the soil surface are identical in both versions.
The various aerodynamic resistances are computed according to Choudhury and Monteith (1988), Shuttleworth and Wallace (1985) and Shuttleworth and Gurney (1990), while the stomatal resistance is modelled according to Braud et al. (1995) for all environmental control factors except water stress, which is replaced by a transpiration efficiency β v , and the moisture-limited evaporation which is governed by an evaporation efficiency β s (Mahfouf and Noilhan, 1991). Definitions of β s and β v are given just below.

The series model version
In the series model the latent heat flux components for the soil (LE s ) and the vegetation (LE v ) are representative averages for the surface as a whole: where ρc p is the product of air density and specific heat, γ the psychrometric constant, r as the soil to aerodynamic level resistance and r vv the minimum total resistance for latent heat exchange between the vegetation and the aerodynamic level (see Appendix A); e sat (T x ) is the saturated vapour pressure at temperature T x (x refers to "s" for soil, "v" for vegetation) and e 0 is the partial pressure of vapour at the aerodynamic level; T s and T v are the soil and vegetation temperatures, respectively. This formulation is different from that of the most common TSEB algorithms which use the Priestley-Taylor relationship to derive a first estimate of LE v . Efficiencies β x are functionally equivalent to surface resistances (again, x referring to "s" for soil, "v" for vegetation, and left blank for the total evapotranspiration flux). Their range of validity is [0, 1]: if β v = 1, then the vegetation transpires at the potential rate, and if β s = 1, the soil evaporation rate is that of a saturated surface, while β v = 0 or β s = 0 corresponds to a nontranspiring or non-evaporating surface, respectively. Scaling between those extremes depends on the soil moisture content around the root zone (for β v ) or in the top few centimetres (for β s ). Here, r vv β v represents a total canopy resistance including stomatal processes, while r as β s corresponds to a total soil evaporation resistance, both in actual conditions. There is no minimum resistance to vapour extraction from the soil porous medium; therefore, resistances above the soil are the same for sensible and latent heat transfers.
In order to reduce the computational cost of solving the system for all unknown variables including T s and T v , all non-linear expressions are linearized though Taylor expansion around air temperature so that the model can be solved through a simple matrix inversion. This is a requirement if one wants to run the model for a large number of pixels. Equations (1) and (2) are converted to Eqs. (3) and (4): where is the slope of the saturation vapour curve at air temperature T a .
The only non-linear term that is kept in either version is the dependence of the aerodynamic resistance to the stability correction. The latter depends on the difference between the aerodynamic temperature and the reference air temperature (Richardson number; cf. Appendix A). Aerodynamic temperature is updated iteratively until convergence.
According to the layer representation in Fig. 1, total fluxes (net radiation, sensible heat flux, latent heat flux, soil heat flux) are computed as the sum of the soil and vegetation components. The continuity of the latent heat flux below and above the aerodynamic level implies that where LE s is expressed in Eq.

Continuity of the sensible heat reads as
where T 0 is the aerodynamic temperature and H s = ρc p T s − T 0 r as and (7) H r a and r av are the aerodynamic level to reference level and vegetation to aerodynamic level aerodynamic resistances, respectively; see Appendix A for their complete expression. Net radiation depends on the grey-body emissions of the soil and vegetation surfaces at temperature T s and T v . Taylor expansion for those emission terms in the net radiation estimates leads to where σ is the Stefan-Boltzman constant and r rad represents a "radiative resistance". Net radiation is computed according to the radiative transfer scheme of Merlin and Chehbouni (2004) which takes into account the multiple reflections between the soil and the vegetation layer in the shortwave and longwave domains. Application of Eq. (9) to the various equations of this scheme leads to a forcing term depending on the incoming shortwave and longwave radiations, A x , and a linear expression of the unknown surface temperatures T s and T v divided by the appropriate radiative resistances r radx (for the expression of those terms, see Appendix A2). For the soil, this leads to and for the canopy, The total flux is The soil heat flux G is a fraction ξ of the net radiation available for the whole of the soil surface (G = R ns ). If the model is run at the same time of the day, for instance with surface temperatures acquired with a sun-synchronous satellite, ξ depends mostly on the bare soil fraction cover. For diurnal variations of G, a time-dependent expression (e.g. Santanello and Friedl, 2003) should be preferred.
The resulting energy balance for the soil (R ns − G = H s + LE s ) and the canopy (R nv = H v + LE v ) for the series model can be written as follows: for the soil and for the vegetation. Finally, the link between the radiative surface temperature T rad and the net longwave radiation components is where R atm is the incoming atmospheric radiation and R an is the net longwave radiation of the whole surface, which depends on T s and T v and can be expressed as follows: The forcing term for the net longwave radiation A atm is given in Appendix A1. The equation relating the radiative surface temperature T rad and the surface temperatures T s and T v is thus

The parallel model version
For the parallel model, all fluxes are representative of each patch (Fig. 1). The total resistance is the sum of the aerodynamic resistance r a and the surface resistances r as (for the soil) or r vv (for the canopy). The transpiration rate of the vegetated subpixel (in W m −2 ) is thus while for the separate patch of bare soil the evaporation rate is where D a = e sat (T a ) − e a is the vapour pressure deficit at reference level. For the parallel model, the sensible heat flux rate above each patch is for the soil, and for the vegetation.
The value of the leaf area index used for the parallel model is a "clump LAI" obtained by dividing the total LAI by the fraction cover of the vegetation f c (Lhomme and Chehbouni, 1999). Total fluxes are the sum of the soil and vegetation components also weighted by their relative contribution, f c for the vegetation and 1-f c for the soil: where LE s is expressed according to Eq. (20) and LE v to Eq. (21), and where H s is expressed according to Eq. (22) and H v according to Eq. (23). The stability correction for the aerodynamic resistance r a depends on an average aerodynamic temperature computed from the total sensible heat flux H : For the parallel model, incoming solar and atmospheric radiations are fully available for each source. The net radiation components are solved independently and, like the turbulent fluxes, summed according to their respective cover fraction. The radiative transfer scheme is simpler than for the series model. The Taylor expansion of the net radiation expression for the soil is written as and, for the vegetation, where A s and A v are the radiation forcing terms for the soil and the vegetation, respectively (see Appendix A1 for their numerical expression). The total flux is The soil heat flux G is a fraction ξ of the net radiation available on the bare soil patch (G = (1 − f c ) ξ R ns ). Finally, the respective energy balance equations for the soil and the vegetation patches of the parallel model are and For the parallel version, the net longwave radiation also has a simpler expression than for the series model: The equation relating the radiative surface temperature T rad and the surface temperatures T s and T v is thus

"Prescribed" and "retrieval" modes
The system of five equations to be solved simultaneously consists of Eqs. (5), (6), (13), (14) and (17) for the series model, and Eqs. (24), (25), (30), (31) and (33) for the parallel model. This system can be solved in a forward mode for which the surface temperature is an output, and an inverse mode when the surface temperature is an input. The SPARSE model combines both modes (cf. Fig. 2). If the soil and the vegetation efficiencies are known (for example through an ancillary two-compartment water budget model), then the model is run in a forward mode from prescribed water stress conditions (from fully stressed to potential). In that case the system is solved for the following unknowns: T rad , T s , T v , e 0 and T 0 . T rad in this prescribed mode is then an output of the system computed from Eqs. (17) and (33) after solving for T s , T v , e 0 and T 0 in the other four equations. This mode has two direct applications. It can be used independently of the retrieval mode to generate an equilibrium surface temperature at the time of the satellite overpass in order to assimilate surface temperature measurements Hydrol. Earth Syst. Sci., 19, 4653-4672 from known β s and β v values computed at the daily or subdaily time steps from a hydrological model (e.g. Er-Raki et al., 2008). It is also implemented as a final step in the retrieval mode to provide theoretical limits corresponding to maximum reachable levels of sensible heat (fully stressed conditions) or latent heat (potential conditions) for each component (the soil and the vegetation). Output fluxes from the retrieval run are bounded by those limiting cases. In full potential conditions, β s = β v = 1, while in fully stressed conditions, β s = β v = 0. In retrieval conditions (inverse mode), T rad is known and is derived from satellite observations or in situ measurements in the thermal infra-red domain. In order to compute the various fluxes of the energy balance, the full set of five equations must be solved simultaneously by inverting the same matrix corresponding to Eqs. (5), (6), (13), (14) and (17) for the series model and Eqs. (24), (25), (30), (31) and (33) for the parallel model. In that case, however, contrarily to the prescribed mode, the problem is initially ill-posed since the system contains six unknowns: evaporation LE s and transpiration LE v , surface temperature components T s and T v , and aerodynamic level conditions e 0 and T 0 . LE s and LE v values are directly converted into stress levels β s and β v using Eqs. (3) and (4) (series model) or (20) and (21) (parallel model). In order to downsize the number of unknowns, SPARSE carries out the same rationale as the TSEB model: as a first guess, the vegetation is supposed to transpire at potential rate; therefore, β v is set to 1, and the system is solved for unknown LE s (thus β s ), T s , T v , e 0 and T 0 . If a negative LE s is obtained, then the assumption of an unstressed canopy proves to be inconsistent with the observed surface temperature level. In that case, one assumes that the vegetation is suffering from water stress. This means that root zone soil moisture is depleted under critical levels, and that, most probably, the soil surface is already long dry. Therefore, β s is set to 0 and the system is solved for LE v (thus β v ) instead of LE s . Finally, if LE v is negative, fully stressed conditions are imposed for both the soil and the vegetation independently of T rad . Of course, inconsistent positive values of LE s corresponding to slightly stressed vegetation conditions can occur when one assumes that the vegetation is unstressed, but in that case the model will not be able to detect this inconsistency. The limit of this hypothesis will be assessed in Sect. 3 through a synthetic case study.
Finally, in order to ensure that LE x outputs are within realistic bounds, LE x values obtained by running SPARSE in "retrieval" conditions are limited by the evapotranspiration components in potential conditions LE x (β s =1, β v =1) computed by SPARSE in prescribed potential conditions (Fig. 2). This procedure is the dual-source equivalent of what is done in the SEBS single-source model (Su, 2002). For consistency, if LE x is limited by LE x (β s =1, β v =1), all fluxes of the corresponding component energy balance (Rn x , H x and G) are set to their values obtained by the "prescribed" mode in potential conditions, i.e. Rn x (β s = 1, β v = 1), H x (β s = 1, β v = 1) and G(β s = 1, β v = 1). The impact of limiting LE x outputs on the model performance will be assessed in Sect. 4.
Also, an arbitrary minimum positive value of LE s = 30 W m −2 is used as the threshold for vegetation stress detection instead of 0, in order to take into account the contribution of vapour transfer from within the topsoil porous network (Boulet et al., 1997).
3 Assessing the retrieval properties of SPARSE through a synthetic case study

Principles of the simulation experiment
The strong underlying assumptions behind SPARSE are that (i) in a first guess the vegetation is supposed to be unstressed, and that (ii) water stress of the vegetation is always concomitant to a non-evaporative soil. This simplification of the soilvegetation-atmosphere continuum impacts not only the total evapotranspiration retrieval, but also its resulting partition between transpiration and soil evaporation. It is thus important to assess the limits of both assumptions. To do so, a synthetic simulation experiment is proposed. The rationale of the synthetic test is as follows: for each combination of known water stress levels affecting either the transpiration or the evaporation of the soil, one can simulate through the energy budgets of the soil and the vegetation the resulting component temperatures T s and T v and the surface temperature of the whole surface (synthetic T rad ). If one assumes that the satellite is actually measuring this temperature, it can be used as input data to get back to the soil evaporation and transpiration levels and their corresponding efficiencies through the retrieval mode. If there was a unique bijective relationship between the component temperatures and the temperature of the whole surface, the retrieved stress lev-els would correspond to the exact combination of the stress levels used to generate the synthetic T rad . Of course this is not the case, and many different combinations of soil and vegetation efficiency values will correspond to the same equilibrium surface temperature. However, one expects that the whole surface energy balance will be well constrained by the knowledge of T rad , i.e. that each value of T rad will correspond to only one surface stress level (or total efficiency). In other words, we expect that SPARSE will not always partition accurately total total evapotranspiration ET in E and T , but will retrieve the ET value relatively satisfactorily.
The objective of the synthetic stress is to assess the inconsistencies of the decision tree that distributes acceptable stress values between the soil and the vegetation, as well as its impact on the component and total evapotranspiration retrieval performances.

Set-up of the synthetic test
In this simulation experiment, the SPARSE model is run sequentially in its two operating modes: the "prescribed" or "forward" mode to generate an estimate of the radiative surface temperature from prescribed β s and β v efficiencies, and the "retrieval" or "inverse" mode to retrieve β s and β v efficiencies using as input data the surface temperature obtained previously through the "prescribed" mode ("synthetic test" branch of Fig. 2). The test consists therefore in computing a mixed surface radiative temperature (T rad ), soil evaporation (LE s ), transpiration (LE v ) and evapotranspiration (LE) for each possible combination of soil evaporation (β s ∈ [0, 1]) and transpiration (β v ∈ [0, 1]) efficiencies in 0.1 increments with the SPARSE model in prescribed mode, then in forcing the SPARSE model with T rad to retrieve new LE s , LE v and total evapotranspiration LE values as well as the corresponding efficiencies (β s , β v and β for the total). β is deduced as the ratio between two total evapotranspiration estimates: one with actual β s and β v and one with β s = β v = 1. In order to assess the limits of the model assumptions for each version, the prescribed and retrieval modes are run for the same version (series or parallel): the surface temperature obtained by each combination of β s and β v for the series model (or the parallel model) in prescribed conditions is used as input for the series model in retrieval mode (or the parallel model). The retrieval performance is then assessed by comparing these new retrieved β s , β v and β values and the ones used to generate T rad . If the retrieval is fully consistent, those efficiencies must match. The test is carried out for average dry climate conditions (R g = 800 W m −2 , RH = 50 %, u a = 2 m s −1 , T a = 25 • C) and a leaf area index characteristic of the maximum development stage of a cereal cover in dry climates (LAI = 3).

Results
Results for the total evapotranspiration efficiency retrieval are illustrated in Fig. 3. One expects rather good performances (albeit some bias) close to the first guess assumptions (transpiration close to potential conditions, i.e. β v ∼ = 1, and low soil evaporation, i.e. β s ∼ = 0) with a degradation when soil evaporation is high and transpiration is low. In Fig. 3, retrieved total efficiency is compared to the prescribed total efficiency for various incremental values of β v for two discrete levels of β s (0.6 and 0.2, top plots), and for incremental values of β s for two discrete levels of β v (0.8 and 0.4, bottom plots).
Total evapotranspiration and its corresponding β efficiency value are well retrieved for each [β s , β v ] combination for the series model formulation (blue points all aligned along the [1 : 1] line), while for the parallel model β is reasonably well retrieved for situations close to the model assumptions, i.e. a low β s and a high β v . For extreme stress values when the assumption underlying SPARSE algorithms is challenged (low transpiration and non-negligible soil evaporation), the parallel model tends to overestimate β.
In Fig. 4, the performances of transpiration (top plots) and evaporation (bottom plots) efficiency retrievals are assessed separately. Since the first guess of SPARSE is that the vegetation is unstressed, the model will tend to overestimate β v . This is the case for all transpiration efficiency values, with, as expected, a larger difference close to a fully transpiring canopy when the inconsistency in β s retrieval is not yet detected. Indeed, for β v values close to 1, the initial guess of an unstressed canopy leads one to assign a fixed value of 1 to β v . The vegetation temperature is therefore underestimated, and the soil temperature that matches the total surface radiative temperature is overestimated. In turn, sensible heat over the soil is overestimated, the soil net radiation is underestimated, and the resulting soil evaporation computed as a residual term is underestimated. As long as this underestimation does not lead to a negative value of β s , the model does not detect the discrepancy. Consequently, especially for a wet soil (top plot on the left-hand side, β s = 0.6), β v retrievals match poorly the prescribed values, and β v values cling to the unstressed boundary, except for very high prescribed stress levels (β v below 0.4 for the series model, 0.2 for the parallel one).
Despite this overestimation, β v retrievals are relatively consistent if the soil is very dry (top plot on the right-hand side, β s = 0.2). Once again, β v retrievals by the series model are closer to the prescribed values than those of the parallel model. Conversely, soil evaporation retrievals (bottom plots) show, as expected, a slight underestimation when the vegetation is close to unstressed (left-hand plot, β v = 0.8). Its amplitude is fairly constant and mirrors the overestimation of the transpiration efficiencies when the soil is dry. In that case, blue dots (series) and red squares (parallel) of the retrievals are close to the [1 : 1] line for all β s levels.   For conditions far from the initial assumption, e.g. low transpiration efficiencies, soil evaporation is largely underestimated. One must note that this is the case for both models and all β s values. Again, moderately stressed vegetation and a low-level soil evaporation rate will always be interpreted in terms of composite surface temperature as a dry soil and fully transpiring vegetation. As a consequence, very small rain events on an otherwise dry soil will most probably be interpreted as a dry soil surface with slightly stressed vegetation. Those cases, not very frequent but not rare either, must be treated with care from a data assimilation perspective.
All those biases should be kept in mind when interpreting results from all dual-source models based on the same rationale: the fact that the total flux is well simulated does not always means that the component fluxes are consistent, let alone realistic. This has been shown for this particular synthetic data set. This test has been carried out using SPARSE due to the possibility the model offers of combining both modes in a consistent synthetic experiment. Its outcomes are illustrated for this model and a single set of vegetation and climatic conditions. We do not claim that those differences between series and parallel retrieval capacities also fully apply to TSEB, but since they share the same strong underlying assumptions and differ mostly by their parameterization of the fluxes, we are convinced that similar differences would be found with TSEB if TSEB could be run in a prescribed mode.

Data sets
Two data sets were used to assess the performance of the series and parallel versions of the SPARSE model over a whole growing season. The first experimental data set was collected over a rainfed wheat with green leaf area index values up to 2 and the second over an irrigated wheat with green LAI values up to 4. Both have been grown in a semiarid climate (central Tunisia and Morocco). Surface temperature data were acquired with a nadir-looking Apogee ther-moradiometer, while energy fluxes were measured according to classical FLUXNET recommendations (Baldocchi et al., 2001) with Campbell ™ CSAT sonic anemometers and Krypton fast-response hygrometers. Observed and simulated latent heat flux values (half-hourly averages in W m −2 ) are compared at midday (local standard time) in all sky conditions. For the rainfed wheat site, there was clearly a problem with the fast-response psychrometer with an energy balance closure of 60 %. Thus for that site the closure was forced and the corrected LE was computed as Rn-H-G. For the irrigated site, the half-hourly closure was of the order of 80 %. For this site closure was achieved with the conservation of the Bowen ratio H / LE; thus, the corrected LE was computed as (Rn-G)/(1 +H / LE). Data for the irrigated wheat site were acquired during the 2004 growing season (B124 site, Boulet et al., 2012), while the experiment for the rainfed wheat took place in 2012.
Leaf area index was estimated with hemispherical photography every 2 to 3 weeks depending on the phenological cycle, validated by destructive measurements during key stages (growth and full cover). Vegetation height was measured at the same dates. Temporal interpolation of leaf area index for both sites is shown in Fig. 5.

Evapotranspiration estimates
Two sets of SPARSE simulations are derived for each model version (series or parallel): in the set most faithful to the original TSEB, outputs are not limited by potential heat flux values; in the second set, outputs are, like in SEBS, bounded by the potential and fully stressed flux rates considered as absolute maximum and minimum reachable values for evaporation as well as transpiration, whatever the "oasis" or microadvection heat transfer might be. Again, this is legitimate for the parallel version, but for the series version one must inquire whether local advection effects do not enhance latent heat flux values over the total potential value of a uniformly wet surface. No calibration is performed, the minimum stomatal resistance value is arbitrarily set to a realistic level for herbaceous vegetation (100 s m −1 , Gentine et al., 2007) and the G / Rn s ratio ξ is set to 40 % (a value often encountered  around midday for bare soils in arid climates). This is consistent with the potential use of this model which is designed to estimate ET routinely from remote-sensing data, based on surface properties derived per land use type in a similar way to most soil-vegetation-atmosphere transfer models applied to continental scales. Those values are of course less sensitive than the uncertainty in the input variable T rad (not shown). In order to relate those first guess results to those obtained by the series and parallel Kustas et al. (1999) TSEB versions, TSEB is also applied with a default value for the Priestley-Taylor coefficient (1.26).
Total flux values are shown in Figs. 6 and 7 for the bounded sets and RMSE values for both bounded and unbounded sets are reported in Table 1. In both cases (series and parallel versions) the RMSE values are of similar order of magnitude and consistent with values found in the literature (cf. Li et al., 2005). The bounded series outputs display the best performances, with RMSE values lowered by 4 to more than 10 W m −2 . Without bounding, values of evapora- tion and transpiration above potential levels are obtained for the series version during vegetation growth, and some negative values of transpiration are found during late maturity and beginning of senescence.
RMSE values for the parallel TSEB version of Kustas et al. (1999) are very close to that of the SPARSE parallel version, while RMSE values for the TSEB series model are similar to the RMSE values displayed by both parallel versions.
Retrieval performances of the other energy balance components in the bounded case have also been assessed. Statistics are shown in Table 2. The series model shows slightly better retrieval performances for soil heat flux for both sites, but only for net radiation for the irrigated wheat and for sensible heat for the rainfed wheat site. This is consistent with Li et al. (2005) and Morillas et al. (2013), who showed that the series TSEB version was more robust than the parallel version, even though, their relative performances were close.

Water stress estimates
Low RMSE values for the total latent heat flux do not guarantee that total water stress is correctly simulated. Indeed,  if moisture availability in the root zone is large enough to maintain ET at potential levels, the prescribed model in potential conditions can already explain a very large amount of the information content within the observed time series, and the added value of thermal infra-red (TIR) data might be limited. It is thus important to assess the amount of information introduced by the surface temperature itself, i.e. information on moisture-limited evaporation and transpiration rates (i.e. second-stage evaporation; cf. Boulet et al., 2004). Water stress is usually defined as the complementary part to 1 of the ratio between the actual and potential evapotranspiration rates. It is expected to scale between 0 (unstressed surface) and 1 (fully stressed surface). Retrieved and observed surface water stress values have been estimated from potential evapotranspiration rates generated with the SPARSE model in prescribed conditions (β s = β v = 1). Simulated and observed water stress values are computed as 1−LE / LE p and 1−LE obs / LE p , respectively, where LE obs is the instantaneous observed latent heat flux, while LE and LE p are the simulated latent heat fluxes in actual and potential conditions, respectively. Total stress is thus functionally equivalent to 1-β. Results are shown in Figs. 8 and 9. As expected, surface stress is much higher for the rainfed than for the irrigated wheat field. The scatter is quite large, therefore showing the intrinsic limit of stress retrieval from naturally noisy TIR data, as already pointed out by numerous studies (Gentine et al., 2010;Katul et al., 1998;Lagouarde et al., 2013Lagouarde et al., , 2015. However, broad tendencies are well reproduced, with most points located within a confidence interval of 0.2 indicated by dotted lines along the 1 : 1 line. This is encouraging in a data assimilation perspective. One must also note that it includes small LE and LE p values for which measurement uncertainty can be as large as the flux itself. To scale those stress values back to potential evapotranspiration, the LE p order of magnitude is indicated as a marker size in Figs. 8 and 9. Most outliers have smaller LE p values, while the points with the largest LE p fall within the space delimited by the two dotted lines of the confidence interval. Some points with little to no evaporation attest to the difficulty in representing accurately the conditions close to the potential levels and might be related to the theoretical limit of the model for small vegetation stress values illustrated in Fig. 3, especially at low evaporation efficiencies.

Soil evaporation efficiency
As shown in the previous sections as well as many previous studies on soil-vegetation-atmosphere interactions in the literature Morillas et al., 2013), series and parallel versions have fairly similar performances in total flux retrieval even though the series version shows slightly better values for the selected statistical criteria. However, as illustrated with the synthetic case, it might not be the case for component flux retrieval. In order to check the consistency of component flux retrieval, one needs a measurement of either soil evaporation or transpiration. At neither site have transpiration data been collected: measuring transpiration for a cereal cover is quite challenging. On the other hand, surface soil moisture data (at a depth of around 5 cm) are available at both sites. Of course, soil moisture at 5 cm does not always react to small rainfall events, but it is a good driver of soil evaporation despite its influence by shallow roots.
We therefore decided to compare the retrieved soil evaporation efficiency to a fairly independent evaluation noted β s_e derived from the observed time series of soil moisture in the top 5 cm (θ 0−5 cm ) instead of using TIR data. We used the efficiency model of Merlin et al. (2011) to derive β s_e : where θ sat is the in situ water content at saturation (0.30 for the rainfed site and 0.48 for the irrigated wheat) and p is fixed to 1 for the loamy site (rainfed wheat) and 0.5 for the clay site (irrigated wheat) according to 1-LE / LE p observations at the beginning and the end of the growing season when the soil is almost bare. Since the surface temperature (and thus the partition between LE s and LE v ) reacts immediately to atmospheric turbulence (Lagouarde et al., 2015) or very small rainfall events, β s instantaneous retrievals by SPARSE show larger fluctuations than β s_e . Indeed, the latter reacts mostly to the largest rainfall events (wetting of the entire 5 cm topsoil). Meteoro-logical forcing can vary quickly and impact the potential soil evaporation rate LE sp , but the latter is less sensitive to turbulence than T rad . In order to smooth out the quick fluctuations of β s retrievals by SPARSE, we compare 5-day running averages of β s and β s_e .
The resulting β s and β s_e evaporation efficiencies are shown on Fig. 10 (rainfed wheat) and 11 (irrigated wheat). For both sites, increasing and decreasing trends of β s and β s_e are mostly synchronous, although their amplitude varies throughout the growing season. Due to irrigation, β s values are on average higher for the irrigated than the rainfed wheat site.
For the rainfed site, both models simulate fairly large values of β s compared to β s_e at the beginning of the season. The parallel model agrees well with β s_e towards the end of the growing stage , while the series model matches very closely β s_e at maximum cover and early senescence (reduction of β s from DOY 70 to DOY 100). Both models agree well with β s_e at the end of the season (DOY 120-170) except for the last 10 days. The small rainfall event around DOY 125 is not sufficient to impact β s_e , but affects β s in both model versions, whereas the soil moisture increase around DOY 105 is mostly missed out by either version.
For the irrigated wheat, soil evaporation is mostly in the energy-limited stage for the first half of the observation period, and β s remains close to 1. This is due to the complement irrigation up to the middle of the maturation phase. The magnitude of both drying events around DOY 40 and DOY 100 is very well retrieved by the series model and somewhat less by the parallel model. Again, β s reacts more strongly to the small rainfall event around DOY 90 than what is indicated from soil moisture.
At the very end of the season both model versions differ greatly from the β s_e estimates and remain close to the potential rate for both sites.

Discussion and conclusion
A new model based on the TSEB rationale, SPARSE, has been presented. Innovation lies mostly in the formulation of the energy balance equations and the use of complementary modes (prescribed and retrieval) which allow one to bound the outputs by realistic limiting flux values which ensure increased robustness. We demonstrated with two data sets that using bounding relationships based on potential conditions decreases the root mean square error by up to 11 W m −2 from values of the order of 50-80 W m −2 . Theoretical limitations of the performance of the evapotranspiration component (evaporation and transpiration) retrievals from a single radiative surface temperature have been inferred over rainfed and irrigated wheat fields at seasonal scales, as well as through a theoretical simulation exercise. According to results obtained in Sect. 3, it is almost impossible to retrieve a non-zero soil evaporation at medium to large LAI values for very high vegetation stress levels. Also, and by construction, transpiration tends to be overestimated in most ranges, but specifically when only slightly stressed. Within these limits, the SPARSE model shows good retrieval performances of evapotranspiration compared to the original TSEB. This comparison must be treated with special care since both models are run with no prior calibration of the poorly known parameters such as the minimum stomatal resistance (for SPARSE) or the Priestley-Taylor coefficient (for TSEB). If a value of r stmin = 50 s m −1 is used, a value also reported for wheat crops in more temperate regions, RMSE on latent heat flux increases by 4 W m −2 in bounded conditions for the rainfed wheat site (62 W m −2 ) and 13 W m −2 for the irrigated wheat site (66 W m −2 ) for the series version. For the parallel model it increases by 12 W m −2 (82 W m −2 ) and 8 W m −2 (74 W m −2 ), respectively.
As expected for cereal covers whose homogeneity is usually well represented by a "layer" approach, the series version provides in general better estimates of latent heat flux values in both real and synthetic cases tested. Those cases are representative of cereals typically grown in semi-arid lands in irrigated and non-irrigated areas. Both models should be tested for other conditions of heterogeneity (sparse crops, orchards, row crops) whose geometrical features are closer to the "patch" description.
Estimates of water stress have also been looked at. Water stress is an interesting variable that can be assimilated in all hydrological or SVAT models in order to compute moisturelimited evapotranspiration rates. Even if the points in the simulated vs. observed scatterplots have a significant number of outliers, i.e. points outside the 0.2 range along the 1 : 1 line in Figs. 8 and 9, the results indicate that the information re- trieved from TIR data is useful from a data-assimilation perspective since the broad tendencies are well reproduced.
Estimates of soil evaporation efficiency have been evaluated against a reconstructed time series relying on observed soil moisture at the soil surface and therefore independent of any surface temperature measurement. This reconstruction is of course model-dependent (Merlin et al., 2011, in our case) and must be considered with care, but despite this we found that both efficiency values are consistent, except at the beginning and the end of the season, partly due to very small rainfall events, but also probably to the poor understanding of turbulence processes over low or senescent vegetation. It seems that the transpiration of the quasi-senescent vegetation encountered at this period of the year is not always well simulated by the model even if total and green LAI values seem realistic. This could be related to the change in soilvegetation radiation exchange and drag partition in a drying vegetation with shrinking leaves and standing straw. In order to smooth out the scale differences between the information provided by soil moisture (a time-continuous variable) and that of surface temperature (influenced by highfrequency turbulent fluctuations), we compared 5-day moving averages. This is consistent with the potential data assimilation method of β or LE estimated from TIR data that one could use in a SVAT model for example: a smoother is more likely to outperform a sequential assimilation algorithm for short observation windows since the former will naturally smooth out the high-order fluctuations due to high-order fluctuations of T rad . Simpler models would perhaps provide similar performances of soil evaporation efficiencies, for instance in rainfed agriculture where surface soil moisture is well constrained by rainfall, but in irrigated areas it is interesting to get proper timing of water inputs, and this can be achieved with relatively good confidence with this model provided that TIR information is available frequently enough. Future work will assess the potential use of microwave data (radar) to infer topsoil moisture and constrain the inversion procedure using a first guess efficiency value generated from topsoil moisture estimates. Current work is directed towards assessing the model performance over other crops, including orchards, and other climates. SPARSE needs more input data than TSEB, for instance relative humidity. The impact of uncertainty on available meteorological data (reanalysis or remote-sensing meteorological products vs. local meteorological stations network) on SPARSE model performance will also be assessed in the future. Coefficient of the stability function n sw Coefficient in r av r a Aerodynamic resistance between the aerodynamic level and the reference level (s m −1 ) R an Longwave net radiation (W m −2 ) r as Aerodynamic resistance between the soil and the aerodynamic level (s m −1 ) R atm Incoming atmospheric radiation (W m −2 ) r av Aerodynamic resistance between the vegetation and the aerodynamic level (s m −1 ) R g Incoming solar radiation (W m −2 ) Ri Richardson number R n Total net radiation (W m −2 ) R ns Net radiation over the soil (W m −2 ) R nv Net radiation over the canopy (W m −2 ) r rad Radiative resistance (s m −1 ) r rads Soil radiative resistance for the parallel model (s m −1 ) r radv Canopy radiative resistance for the parallel model (s m −1 ) r radss Soil radiative resistance for the soil net radiation in the series model (s m −1 ) r radsv Canopy radiative resistance for the soil net radiation in the series model (s m −1 ) r radvs Soil radiative resistance for the vegetation net radiation in the series model (s m −1 ) r radvv Canopy radiative resistance for the vegetation net radiation in the series model (s m −1 ) r stmin Minimum stomatal resistance (s m −1 ) r vv Surface resistance between the aerodynamic level and the reference level (s m −1 ) T 0 Aerodynamic temperature (K) T a Air temperature at reference level (K) T rad Radiative surface temperature (K)  Vegetation emissivity Slope of the vapour pressure deficit at T a (Pa K −1 ) γ Psychrometric constant (Pa K −1 ) ρ Air density (kg m −3 ) σ Stefan-Boltzmann constant (W m −2 K −4 ) θ 0−5 cm Integrated volumetric soil moisture in the top 5 cm θ sat Volumetric soil moisture at saturation φ View zenith angle (rad)