Assessment of actual evapotranspiration over a semiarid heterogeneous land surface by means of coupled low-resolution remote sensing data with an energy balance model: comparison to extra-large aperture scintillometer measurements

In semiarid areas, agricultural production is restricted by water availability; hence, efficient agricultural water management is a major issue. The design of tools providing regional estimates of evapotranspiration (ET), one of the most relevant water balance fluxes, may help the sustainable management of water resources.

Remote sensing provides periodic data about actual vegetation temporal dynamics (through the normalized difference vegetation index, NDVI) and water availability under water stress (through the surface temperature T_surf), which are crucial factors controlling ET.

In this study, spatially distributed estimates of ET (or its energy equivalent, the latent heat flux LE) in the Kairouan plain (central Tunisia) were computed by applying the Soil Plant Atmosphere and Remote Sensing Evapotranspiration (SPARSE) model fed by low-resolution remote sensing data (Terra and Aqua MODIS). The work's goal was to assess the operational use of the SPARSE model and the accuracy of the modeled (i) sensible heat flux (H) and (ii) daily ET over a heterogeneous semiarid landscape with complex land cover (i.e., trees, winter cereals, summer vegetables).

SPARSE was run to compute instantaneous estimates of H and LE fluxes at the satellite overpass times. The good correspondence (R² = 0.60 and 0.63 and RMSE = 57.89 and 53.85 W m⁻² for Terra and Aqua, respectively) between instantaneous H estimates and large aperture scintillometer (XLAS) H measurements along a path length of 4 km over the study area showed that the SPARSE model presents satisfactory accuracy. Results showed that, despite the fairly large scatter, the instantaneous LE can be suitably estimated at large scales (RMSE = 47.20 and 43.20 W m⁻² for Terra and Aqua, respectively, and R² = 0.55 for both satellites). Additionally, water stress was investigated by comparing modeled (SPARSE) and observed (XLAS) water stress values; we found that most points were located within a 0.2 confidence interval, thus the general tendencies are well reproduced. Even though extrapolation of instantaneous latent heat flux values to daily totals was less obvious, daily ET estimates are deemed acceptable.

1 Introduction

In water-scarce regions, especially arid and semiarid areas, the sustainable use of water by resource conservation as well as the use of appropriate technologies to do so is a priority for agriculture (Amri et al., 2014; Pereira et al., 2002).

Water use rationalization is especially needed for countries actually suffering from water scarcity, or for countries that probably would suffer from water restrictions according to climate change scenarios. Indeed, the Mediterranean region is one of the most prominent “hot spots” in future climate change projections (Giorgi and Lionello, 2008) due to an expected larger warming than the global average and to a pronounced increase in precipitation interannual variability. The major part of the southern Mediterranean countries, among others Tunisia, already suffer from water scarcity and show a growing water deficit, due to the combined effect of the growth in water needs (soaring demography and irrigated areas extension) and the reduction of resources (temporary drought and/or climate change). This implies that closely monitoring the water budget components is a major issue (Oki and Kanae, 2006).

The estimation of evapotranspiration (ET) is of paramount importance since it represents the preponderant component of the terrestrial water balance; it is the second largest component after precipitation (Glenn et al., 2007); hence ET quantification is a key factor for scarce water resources management. Direct measurement of ET is only possible at local scale (single field) using the eddy covariance method, for example, whereas it is much more difficult at larger scales (irrigated perimeter or watershed) due to the complexity not only of the hydrological processes (Minacapilli et al., 2007) but also of the hydrometeorological processes. Indeed, at landscape scale, surface heterogeneity influences regional and local climate, inducing for example cloudiness, precipitation and temperature pattern differences between areas of higher elevation (hills and mountains surrounding the Kairouan plain) and the plain downstream. Moreover, at these scales, land cover is usually heterogeneous and this affects the land–atmosphere exchanges of heat, water and other constituents (Giorgi and Avissar, 1997). ET estimates for various temporal and spatial scales, from hourly to monthly to seasonal time steps, and from field to global scales, are required for hydrologic applications in water resource management (Anderson et al., 2011). Techniques using remote sensing (RS) information are therefore essential when dealing with processes that cannot be represented by point measurements only (Su, 2002).

In fact, the contribution of RS in vegetation's physical characteristics monitoring in large areas was identified years ago (Tucker, 1978); RS provides periodic data about some major ET drivers, amongst others surface temperature and vegetation properties (e.g., normalized difference vegetation index NDVI and leaf area index LAI) from field to regional scales (Li et al., 2009; Mauser and Schädlich, 1998). Many methods using remotely sensed data to estimate ET are reviewed in Courault et al. (2005). ICARE (Gentine et al., 2007) and SiSPAT (Braud et al., 1995) are examples of complex physically based land surface models (LSMs) using RS data. They include a detailed description of the vegetation water uptake in the root zone, the interactions among groundwater, root zone and surface water. However, the lateral surface and subsurface flows are neglected. This can lead to inaccurate results when applied in areas where such interactions are important (Overgaard et al., 2006).

Moreover, RS can provide estimates of large area fluxes in remote locations, although those estimates are based on the spatial and temporal scales of the measuring systems and thus vary one from another. Hence, one solution is to upscale local micrometeorological measurements to larger spatial scales in order to acquire an optimum representation of land–atmosphere interactions (Samain et al., 2012). However, such an upscaling process is not always possible and results might not be reliable in comparison to the RS distributed products.

Water and energy exchange in the soil–plant–atmosphere continuum have been simulated through several land surface models (Bastiaanssen et al., 2007; Feddes et al., 1978). Among them, two different approaches use remote sensing data to estimate spatially distributed ET (Minacapilli et al., 2009): one that is based on the soil water balance (SWB) and one that solves the surface energy budget (SEB). The SWB approach exploits only visible near-infrared (VIS-NIR) observations to perceive the spatial variability of crop parameters. The SEB modeling approach uses visible (VIS), near-infrared (NIR) and thermal infrared (TIR) data to solve the SEB equation by forcing remotely sensed estimates of the SEB components (mainly the surface temperature T_surf). In fact, there is a strong link between water availability in the soil and surface temperature under water stress; hence, in order to estimate soil moisture status as well as actual ET at relevant space and timescales, information in the TIR domain (8–14 µm) is frequently used (Boulet et al., 2007). The SWB approach has the advantage of high resolution and high-frequency VIS-NIR remote sensing data availability against limited availability of high-resolution thermal imagery for the SEB approach. Indeed, satellite data such as Landsat or Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) data provide field-scale (30–100 m) estimates of ET (Allen et al., 2011), but they have a low temporal resolution (16 days to monthly) (Anderson et al., 2011).

The RS-based SWB models provide estimates of ET, soil water content and irrigation requirements in a continuous way. For instance, at field scale, estimates of seasonal ET and irrigation can be obtained by SWB modeling using high-resolution remote sensing forcing as done in the study with the SAtellite Monitoring of IRrigation (SAMIR) model by Saadi et al. (2015) over the Kairouan plain. However, for an appropriate estimation of ET, the SWB model requires knowledge of the water inputs (precipitation and irrigation) and an assessment of the extractable water from the soil (mostly derived from the soil moisture characteristics: actual available water content in the root zone, wilting point and field capacity), whereas significant biases are found mainly when dealing with large areas and long periods, due to the spatial variability of the water input uncertainties as well as the inaccuracy in estimating other flux components such as the deep drainage (Calera et al., 2017). Hence, the major limitation of the SWB method is the high number of needed inputs whose estimation is highly uncertain, especially over a heterogeneous land surface due to hydrologic processes complexity. Moreover, spatially distributed SWB models, typically those using the Food and Agriculture Organization-FAO guidelines (Allen et al., 1998) for crop ET estimation, generally parameterize the vegetation characteristics on the basis of land use maps (Bounoua et al., 2015; Xie et al., 2008), and different parameters are used for different land use classes. Nevertheless, SWB modelers generally do not have the possibility to carry out RS-based land use change mapping due to time, budget or capacity constraints and often use very generic classes, potentially leading to modeling errors (Hunink et al., 2017). In addition, the lack of data about the soil properties (controlling field capacity, wilting point and the water retention) as well as the actual root depths lead to limited practical use of the SWB models (Calera et al., 2017). The same applies to the soil evaporation whose estimation generally relies on the FAO guidelines approach (Allen et al., 1998). However, it was shown that under high evaporation conditions, the FAO-56 (Allen et al., 1998) daily evaporation computed on the basis of the readily evaporable water (REW) is overestimated at the beginning of the dry down phase (i.e., the period after rain or irrigation where the soil moisture is decreasing due to evapotranspiration and drainage, Mutziger et al., 2005; Torres and Calera, 2010). Hence, to improve its estimation a reduction factor proposed by Torres and Calera (2010) was applied to deal with this problem in several studies (e.g., Odi-Lara et al., 2016; Saadi et al., 2015). Furthermore, SWB models such as SWAP (Kroes et al., 2017, Cropsyst (Stöckle et al., 2003), AquaCrop (Steduto et al., 2009) and SAMIR (Simonneaux et al., 2009) are able to take irrigation into account, either as an estimated amount provided by the farmer (as an input if available) or a predicted amount through a module triggering irrigation according to, say, critical soil moisture levels (as an output). However, the limited knowledge of the actual irrigation scheduling is a critical limitation for the validation protocol of irrigation requirement estimates by SWB modeling. Therefore, SWB modelers must deal with the lack of information about real irrigation, which induces unreliable estimations.

Consequently, ET estimation at regional scale is often achieved using SEB approaches, by combining surface temperature from medium- to low-resolution (kilometer scale) remote sensing data with vegetation parameters and meteorological variables (Liou and Kar, 2014). Recently, many efforts have been made to feed remotely sensed surface temperature into ET modeling platforms in combination with other critical variables, e.g., NDVI and albedo (Kalma et al., 2008; Kustas and Anderson, 2009). A wide range of satellite-based ET models were developed, and these methods are reviewed in Liou and Kar (2014). The majority of SEB-based models are single-source models; their algorithms compute a total latent heat flux as the sum of the evaporation and the transpiration components using a remotely sensed surface temperature. However, separate estimates of evaporation and transpiration make the dual-source models more useful for agrohydrological applications (water stress detection, irrigation monitoring, etc.) (Boulet et al., 2015).

Figure 1 The study area: the downstream Merguellil subbasin is the so-called Kairouan plain; the MODIS grid is the extracted 10 km × 8 km MODIS subimage and in red the scintillometer XLAS transect.

Contrarily to SWB models, most SEB models are run in their most standardized version, using observed RS-based parameters such as albedo in conjunction with a set of input parameters taken from literature or in situ data. However, the SEB model validation with enough data in space and time is difficult to achieve, due to the limited availability of high-resolution thermal images (Chirouze et al., 2014). Therefore, it is usually possible to evaluate SEB models results only at similar scale (km) to medium- or low-resolution images. Indeed, the pixel size of thermal remote sensing images, except for the scarce Landsat7 images (60 m), covers a range of 1000 m (Moderate Sensors Resolution Imaging Spectroradiometer MODIS), to 4000 m (Geostationary Operational Environmental Satellite GEOS). However, direct methods measuring sensible heat fluxes (eddy covariance for example) only provide point measurements with a footprint considerably smaller than a satellite pixel. Therefore, scintillometry techniques have emerged as one of the best tools aiming to quantify averaged fluxes over heterogeneous land surfaces (Brunsell et al., 2011). They provide area-averaged sensible heat flux over areas comparable to those observed by satellites (Hemakumara et al., 2003; Lagouarde et al., 2002). Scintillometry can provide sensible heat using different wavelengths (optical wavelength ranges), aperture sizes (15–30 cm) and configurations (long-path and short-path scintillometry) (Meijninger et al., 2002). The upwind area contributing to the flux (i.e., the flux footprint) varies as wind direction and atmospheric stability and must be estimated for the surface measurements in order to compare them to SEB estimates of the flux which are representative of the pixel (Brunsell et al., 2011). Assessing the upwind area contributing to the flux can be done using several footprint models (Schmid, 2002). Although footprint analysis ensures ad hoc spatial intersecting area between ground measurements and satellite-based surface fluxes, the spatial heterogeneity at subpixel scale should be further considered in validating low-resolution satellite data (Bai et al., 2015). The LAS technique has been validated over heterogeneous landscapes against eddy covariance measurements (Bai et al., 2009; Chehbouni et al., 2000; Hoedjes et al., 2009) and also against modeled fluxes (Marx et al., 2008; Samain et al., 2012; Watts et al., 2000). Few studies dealt with extra-large aperture scintillometer (XLAS) data (Kohsiek et al., 2002, 2006; Moene et al., 2006). Historical survey, theoretical background and recent works in applied research concerning scintillometry are reviewed in De Bruin and Wang (2017). Since the scintillometer provides large-scale area-average sensible heat flux (H_XLAS), the corresponding latent heat flux (LE_XLAS) can then be computed as the energy balance residual term (LE_XLAS = R_n − G − H_XLAS); hence, the estimation of a representative value for the available energy (AE = R_n − G) is always crucial for the accuracy of the retrieved values of LE_XLAS. This assumption is valid only under the similarity hypothesis of Monin–Obukhov (MO) (Monin and Obukhov, 1954), i.e., surface homogeneity and stationary flows. These hypotheses are verified in our study area where the topography is flat, and the landscape is heterogeneous only from an agronomic point of view since we find different land uses (cereals, market gardening and fruit trees, mainly olive trees with considerable spacing of bare soil); however, this heterogeneity in landscape features at field scale is randomly distributed and there is no drastic change in height and density of the vegetation at the scale of the XLAS transect (i.e., little heterogeneity at the kilometer scale, most MODIS pixels have similar NDVI values for instance).

In this study, spatially distributed estimates of surface energy fluxes (sensible heat H and latent heat fluxes LE) over an irrigated area located in the Kairouan plain (central Tunisia) were obtained by the SEB method, using the Soil Plant Atmosphere and Remote Sensing Evapotranspiration (SPARSE) model (Boulet et al., 2015) fed by 1 km thermal data and 1 km NDVI data from MODIS sensors on the Terra and Aqua satellites. The main objective of this paper is to compare the modeled H and LE simulated by the SPARSE model with, respectively, the H measured by the XLAS and the LE reconstructed from the XLAS measurements acquired in the course of 2 years over a large, heterogeneous area. We explore the consistency between the instantaneous H and LE estimates at the satellite overpass times, the water stress estimates and also ET derived at daily time steps from both approaches.

2 Experimental site and datasets

2.1 Study area

The study site is a semiarid region located in central Tunisia, the Kairouan plain (35^∘1^′–35^∘55^′ N, 9^∘23^′–10^∘17^′ E; Fig. 1). The landscape is mainly flat, and the vegetation is dominated by agricultural production (cereals, olive groves, fruit trees, market gardening; Zribi et al., 2011). Water management in the study area is typical of semiarid regions with an upstream subcatchment that transfers surface and subsurface flows collected by a dam (the El Haouareb dam) and a downstream plain (Kairouan plain) supporting irrigated agriculture (Fig. 1). Agriculture consumes more than 80% of the total amount of water extracted each year from the Kairouan aquifer (Poussin et al., 2008). Most farmers in the plain use their own wells to extract water for irrigation (Pradeleix et al., 2015), while a few depends on public irrigation schemes based on collective networks of water distribution pipelines all linked to a main borehole. The crop intensification in recent decades, associated with increasing irrigation, has led to growing water demand, and an over-exploitation of the groundwater (Leduc et al., 2004).

Figure 2 XLAS setup: XLAS transect (white), for which the emitter and the receiver are located at the extremity of each white arrow; half-hourly XLAS footprint for selected typical wind conditions (green); MODIS grid (black); orchards (blue); and the location of the Ben Salem meteorological and flux stations. Background is a three-color (red, green, blue) composite of SPOT5 bands 3 (NIR), 2 (VIS-red) and 1 (VIS-green) acquired on 9 April 2013 and showing in red the cereal plots.

2.2 Experimental setup and remote sensing data

An optical Kipp & Zonen extra-large aperture scintillometer (XLAS) was operated continuously for more than 2 years (1 March 2013 to 3 June 2015) over a relatively flat terrain (maximum difference in elevation of about 18 m). The scintillometer consists of a transmitter and a receiver, both with an aperture diameter of 0.3 m, which allows longer path length. The wavelength of the light beam emitted by the transmitter is 940 nm. The transmitter was located on an eastern water tower (coordinates: 35^∘34^′0.7^′′ N, 9^∘53^′25.19^′′ E; 127 m a.s.l. – above sea level) and the receiver on a western water tower (coordinates: 35^∘34^′17.22^′′ N, 9^∘56^′7.30^′′ E; 145 m a.s.l.) separated by a path length of 4 km (Fig. 2).

The scintillometer transect was above a mixed-vegetation canopy: trees (mainly olive orchards) with some annual crops (cereals and market gardening) and the mean vegetation height is estimated to be about 1.17 m along the transect. Both instruments were installed at 20 m height as recommended in the Kipp & Zonen instruction manual for LAS and XLAS (Kipp & Zonen, 2017). At this height and for a 4 km path length, the devices are high enough to minimize measurement saturation and assumed to be above or close to the blending height where MO applied.

Furthermore, two automatic Campbell Scientific (Logan, USA) eddy covariance (EC) flux stations were also positioned at the same level on the two water-tower-top platforms. Half-hourly turbulent fluxes in the western and the eastern EC stations were measured using a sonic anemometer CSAT3 (Campbell Scientific, USA) at a rate of 20 Hz and a sonic anemometer RM 81000 (Young, USA) at a rate of 10 Hz, respectively. The western station data were more reliable with less measurement errors and gaps; hence, the western EC setup was used to initialize friction velocity u^* values and the Obukhov length L_o in the scintillometer flux computation (Sect. 3.1).

Half-hourly standard meteorological measurements including incoming long wave radiation, i.e., global incoming radiation ($R_{{\mathrm{g}}_{\mathrm{30}}}$), the incoming longwave radiation (i.e., atmospheric radiation, R_atm-30), wind speed (u₃₀), wind direction ($u_{{\mathrm{d}}_{\mathrm{30}}})$, air temperature (T₃₀), relative humidity (RH₃₀) and barometric pressure (P₃₀) were recorded using an automated weather station installed in the study area (Fig. 2), referred to as the Ben Salem meteorological station (35^∘33^′1.44^′′ N, 9^∘55^′18.11^′′ E). Meteorological data were used either to force the SPARSE model or as input data in XLAS-derived sensible and latent heat flux. The global incoming radiation was also used in the extrapolation method to scale instantaneous observed (Sect. 3.3.2) and modeled (Sect. 4.2) available energy as well as modeled sensible heat flux (Sect. 4.2) to daily values.

In addition, an EC flux station, referred to as the Ben Salem flux station (a few tens of meters away from the meteorological station) was installed from November 2012 to June 2013 in an irrigated wheat field (Fig. 2) to measure half-hourly convective fluxes exchanged between the surface and the atmosphere (H_BS-30 and LE_BS-30) combined with measurements of the net radiation $R_{{\mathrm{n}}_{\mathrm{BS}\text{-}\mathrm{30}}}$ and the soil heat flux G_BS-30. Net radiation and soil heat flux measurements were transferred to the meteorological station from June 2013 till June 2015. Since there are no R_n and G measurements in the two water towers EC stations, $R_{{\mathrm{n}}_{\mathrm{BS}}}$ and G_BS measurements were among the input data to derive sensible and latent heat fluxes from the XLAS measurements. In addition, measured available energy (AE_BS = $R_{{\mathrm{n}}_{\mathrm{BS}}}$ − G_BS) and H_BS were used to calibrate the extrapolation relationship of the available energy and the sensible heat flux, respectively (Sects. 3.3.2 and 4.2).

Remotely sensed data were acquired for the study period (1 September 2012 to 30 June 2015) at the resolution of the MODIS sensor at 1 km, embarked on board of the satellites Terra (overpass time around 10:30 local solar time) and Aqua (overpass time around 13:30 local solar time). Downloaded MODIS products were (i) MOD11A1 and MYD11A1 for Terra and Aqua, respectively (surface temperature T_surf, surface emissivity ε_surf and viewing angle ϕ); (ii) MOD13A2 and MYD13A2 for Terra and Aqua, respectively (NDVI); and (iii) MCD43B1, MCD43B2 and MCD43B3 (albedo α). These MODIS data provided in sinusoidal projection were reprojected in the Universal Transverse Mercator coordinate system using the MODIS Reprojection Tool. Then, subimages of 10 km × 8 km centered on the XLAS transect (Fig. 1) were extracted. The daily MODIS T_surf and viewing angle, 8-day MODIS albedo, and 16-day MODIS NDVI contain some missing or unreliable data; hence, days with missing data (35 % of all dates) in MODIS pixels regarding the scintillometer footprint (see later footprint computation in Sect. 3.2) were excluded. Albedo products (MCD43) are available every 8 days; the day of interest is the central date. Both Terra and Aqua data are used in the generation of this product, providing the highest probability for quality input data and designating it as a combined product. Moreover, the 1 km and 16-day NDVI products (MOD13A2/MYD13A2) are available every 16 days and separately for Terra and Aqua. Algorithms generating this product operate on a per-pixel basis and require multiple daily observations to generate a composite NDVI value that will represent the full period (16 days). For both products, data are linearly interpolated over the available dates in order to get daily estimates. For each pixel, the quality index supplied with each product is used to select the best data.

3 Extra-large aperture scintillometer (XLAS): data processing

3.1 Scintillometer-derived fluxes

Scintillometer measurements are based on the scintillation theory; fluxes of sensible heat and momentum cause atmospheric turbulence close to the ground, and create, with surface evaporation, refractive index fluctuations due mainly to air temperature and humidity fluctuations (Hill et al., 1980). The fluctuation intensity of the refractive index is directly linked to sensible and latent heat fluxes. The light beam emitted by the XLAS transmitter towards the receiver is dispersed by the atmospheric turbulence. The scintillations representing the intensity fluctuations are analyzed at the XLAS receiver and are expressed as the structure parameter of the refractive index of air integrated along the optical path $C_{\mathrm{n}}^{\mathrm{2}}$ (m${}^{-\mathrm{2}/\mathrm{3}}$) (Tatarskii, 1961). The sensitivity of the scintillometer to $C_{\mathrm{n}}^{\mathrm{2}}$ along the beam is not uniform and follows a bell-shaped curve due to the symmetry of the devices. This means that the measured flux is more sensitive to sources located towards the transect center and is less affected by those close to the transect extremities.

In order to compute the XLAS sensible heat flux, $C_{\mathrm{n}}^{\mathrm{2}}$ was converted to the structure parameter of temperature turbulence $C_{\mathrm{T}}^{\mathrm{2}}$ (K² m${}^{-\mathrm{2}/\mathrm{3}}$) by introducing the Bowen ratio (ratio of sensible to latent heat fluxes), hereafter referred to as β, which is a temperature–humidity correlation factor. Moreover, the height of the scintillometer beam above the surface varies along the path. In our study site, the terrain is very flat leading to little beam height variation across the landscape, except for what is induced by the different roughness values of the individual fields. Since the interspaces between trees are large, the effective roughness of the orchards is not significantly different from that of annual crop fields. Consequently, $C_{\mathrm{n}}^{\mathrm{2}}$ and therefore $C_{\mathrm{T}}^{\mathrm{2}}$ are not only averaged horizontally but vertically as well.

At visible wavelengths, the refractive index is sensitive to temperature fluctuations. Then, we can relate the $C_{\mathrm{n}}^{\mathrm{2}}$ to $C_{\mathrm{T}}^{\mathrm{2}}$ as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{C}_{\mathrm{n}}\mathrm{2}={\left({\displaystyle \frac{-\mathrm{0.78}\times {\mathrm{10}}^{-\mathrm{6}}\times P}{T^{\mathrm{2}}}}\right)}^{\mathrm{2}}{C}_{\mathrm{T}}^{\mathrm{2}}{\left(\mathrm{1}+{\displaystyle \frac{\mathrm{0.03}}{\mathit{\beta }}}\right)}^{\mathrm{2}},\end{array}$$

with T the air temperature (^∘K) and P the atmospheric pressure (Pa).

Green and Hayashi (1998) proposed another method to compute XLAS sensible heat flux (H_XLAS) assuming full energy budget closure and using an iterative process without the need for β as an input parameter. This method is called the “β-closure method” (BCM, Twine et al., 2000). In the calculation algorithm, β is estimated iteratively with the BCM method, as described in Solignac et al. (2009) with an initial guess using $R_{{\mathrm{n}}_{\mathrm{BS}}}$ and G_BS from the Ben Salem flux station and initial u_* coming from the western water tower EC station.

Then, the similarity relationship proposed by Andreas (1988) is used to relate the $C_{\mathrm{T}}^{\mathrm{2}}$ to the temperature scale T_* in unstable atmospheric conditions as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle \frac{C_{\mathrm{T}}^{\mathrm{2}}{\left(z_{\mathrm{LAS}}-d\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}}{T_{*}^{\mathrm{2}}}}=\mathrm{4.9}\left(\mathrm{1}-\mathrm{6.1}{\left({\displaystyle \frac{z_{\mathrm{LAS}}-d}{L_{\mathrm{O}}}}\right)}^{-\frac{\mathrm{2}}{\mathrm{3}}}\right),\end{array}$$

and for stable atmospheric conditions,

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle \frac{C_{\mathrm{T}}^{\mathrm{2}}{\left(z_{\mathrm{LAS}}-d\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}}{T_{*}^{\mathrm{2}}}}=\mathrm{4.9}\left(\mathrm{1}+\mathrm{2.2}{\left({\displaystyle \frac{z_{\mathrm{LAS}}-d}{L_{\mathrm{O}}}}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}\right),\end{array}$$

where L_O (m) is the Monin–Obukhov length, Z_LAS (m) is the scintillometer height, and d (m) is the displacement height, which corresponds to two-thirds of the averaged vegetation height z_v.

The variable z_v counts for the various heights within the selected footprint using angular zones originating from the center of the transect and supported by high-resolution remote sensing data (see Sect. 4.1).

Furthermore, considering the size of the surface changes in roughness (mean or effective vegetation height ∼ 1.5 m), the XLAS measurement height is assumed to be either close to the blending height or higher. Thus, the fluxes measured by XLAS are area-averaged and the MO similarity hypothesis can be applied in the flux algorithm computation.

Moreover, the fluxes measured with the XLAS were calculated for a stability index (Z_LAS∕L_O) between −2 and 0 (∼ 73 % of the cases). Then, according to Gruber and Fochesatto (2013), the sensitivity of the turbulent flux measurements to uncertainties in source measurements is rather similar between iterative algorithms and analytical solution, in this range of atmospheric stability.

From T_* and the friction velocity u_* (computed based on an iteration approach in the BCM method), the sensible heat flux can be derived as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }H=-\mathit{\rho }{c}_{\mathrm{p}}{T}_{*}{u}_{*},\end{array}$$

where ρ (kg m⁻³) the density of air and c_p (J kg⁻¹ K⁻¹) the specific heat of air at constant pressure.

H_XLAS was computed at a half-hourly time step. Before flux computation, a strict filtering was applied to the XLAS data to remove outliers depending on weak demodulation signal. Negative nighttime data were set to zero and daytime flux missing data (one to three 30 min data) were gap filled using simple interpolation. Furthermore, half-hourly H_XLAS aberrant values due to measurement errors and values higher than 400 W m⁻², arising from measurement saturation, were ruled out (3 % of the total measurement throughout the experiment duration). Finally, daily H_XLAS was computed as the average of the half-hourly H_XLAS.

3.2 XLAS footprint computation

The footprint of a flux measurement defines the spatial context of the measurement and the source area that influences the sensors. In case of inhomogeneous surfaces like patches of various land covers and moisture variability due to irrigation, the measured signal is dependent on the fraction of the surface having the strongest influence on the sensor and thus on the footprint size and location. Footprint models (Horst and Weil, 1992; Leclerc and Thurtell, 1990) have been developed to determine what area is contributing to the heat fluxes as well as the relative weight of each particular cell inside the footprint limits. Contributions of upwind locations to the measured flux depend on the height of the vegetation, height of the instrumentation, wind speed, wind direction and atmospheric stability conditions (Chávez et al., 2005).

According to the model of Horst and Weil (1992), for a one-point measurement system, the footprint function f relates the spatial distribution of surface fluxes, F₀(x, y), to the measured flux at height, z_m, F(x, y, z_m), as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }F\left(x,y,z_{\mathrm{m}}\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{x}{\int }}{F}_{\mathrm{0}}(x^{\prime },y^{\prime })f\left(x-x^{\prime },y-y^{\prime },z_{\mathrm{m}}\right)\mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime }.\end{array}$$

The footprint function f is computed as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }{\stackrel{\mathrm{‾}}{f}}^y\left(x,z_{\mathrm{m}}\right)={\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{z}}{\mathrm{d}x}}{\displaystyle \frac{z_{\mathrm{m}}}{{\stackrel{\mathrm{‾}}{z}}^{\mathrm{2}}}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{m}}\right)}{\stackrel{\mathrm{‾}}{u}\left(c\stackrel{\mathrm{‾}}{z}\right)}}A{e}^{-{\left(z_{\mathrm{m}}/b\stackrel{\mathrm{‾}}{z}\right)}^r},\end{array}$$

where $\stackrel{\mathrm{‾}}{u}\left(z\right)$ is the mean wind speed profile and $\stackrel{\mathrm{‾}}{z}$ is the mean plume height for diffusion from a surface source. The variables A, b and c are scale factors and r is a scale factor of the Gamma function. In the case of a scintillometer measurement, the footprint function has to be combined with the spatial weighting function W(x) of the scintillometer to account for the sensor integration along its path. Thus, the sensible heat flux footprint mainly depends on the scintillometer effective height Z_LAS (Hartogensis et al., 2003), which includes the topography below the path and the transmitter and receiver heights, the wind direction, and the Obukhov length L_O, which characterizes the atmospheric stability (Solignac et al., 2009). In a subsequent step, daily footprints were computed as a weighted sum of the half-hourly footprints by the XLAS sensible heat flux.

In fact, there is an issue with the MODIS pixel heterogeneity and notably the distribution of the land use classes at the intersection between the square pixel and the XLAS footprint (Bai et al., 2015). Hence, in order to provide a first guess on these relative heterogeneities, land use classes within each MODIS pixel of the 10 km × 8 km subimage were studied based on the land use map of the 2013–2014 season (Chahbi et al., 2016). The average footprint of all half-hourly footprints for the whole study period was computed and overlaid on the MODIS grid in order to identify the MODIS pixels partially or totally covered by footprint (Fig. 3).

Figure 3 MODIS pixels partially or totally covered by XLAS source area.

Figure 4 Land use classes' percentage of the MODIS pixels within or outside the footprint.

The percentage of land use classes was computed for (i) the part of each pixel that lies within the footprint and (ii) the complementary part of the pixel located outside of the footprint (Fig. 4). Results show that differences in percentages of each land use classes for the pixel fractions located within or outside the footprint are low, with 1.8, 1.7, 1.0 and 3.5 % for cereals, market gardening, trees and bare soil, respectively. Moreover, the major part of the area above the transect is covered by fallow and orchards. The land use classes' partition inside the 13 MODIS pixels totally covered by the average footprint is comparable.

3.3 XLAS-derived latent heat flux

Instantaneous (LE_residual_XLAS_t-FP) and daily (LE_residual_XLAS_day-FP) XLAS-derived latent heat fluxes (i.e., residual latent heat flux) of the XLAS upwind area were computed using the energy budget closure of the XLAS measured sensible heat flux (H_XLAS) with additional estimations of remotely sensed net surface radiation R_n and soil heat flux G, as available energy (AE = R_n − G), as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }\mathrm{LE}\mathrm{\_}\mathrm{residual}\mathrm{\_}{\mathrm{XLAS}}_{t\text{-}\mathrm{FP}}={\mathrm{AE}}_{t\text{-}\mathrm{FP}}-\mathrm{H}\mathrm{\_}{\mathrm{XLAS}}_t,\text{(8)}& {\displaystyle }& {\displaystyle }\mathrm{LE}\mathrm{\_}\mathrm{residual}\mathrm{\_}{\mathrm{XLAS}}_{\mathrm{day}\text{-}\mathrm{FP}}={\mathrm{AE}}_{\mathrm{day}\text{-}\mathrm{FP}}-\mathrm{H}\mathrm{\_}{\mathrm{XLAS}}_{\mathrm{day}}.\end{array}$$

H_XLAS_t and H_XLAS_day are respectively the instantaneous and daily measured H at the times of the satellite overpass interpolated from the half-hourly fluxes measurements. Daily available energy within the footprint (AE_day-FP) was computed from instantaneous available energy (AE_t-FP) as detailed in Sect. 3.3.1 and Sect. 3.3.2. The subscripts “30”, “day” and “t” refer to half-hourly, daily and instantaneous (at the times of Terra and Aqua overpasses) variables, respectively, while the subscript “FP” means that the footprint is taken into account, i.e., instantaneous or the daily (depending on timescale) footprint was multiplied by the variable.

3.3.1 Instantaneous available energy

Net surface radiation is the balance of energy between incoming and outgoing shortwave and longwave radiation fluxes at the land–atmosphere interface. Remotely sensed surface radiative budget components provide unparalleled spatial and temporal information, and thus several studies have attempted to estimate net radiation by combining remote sensing observations with surface and atmospheric data. Net radiation equation can be written as follows:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{R}_{\mathrm{n}}=(\mathrm{1}-\mathit{\alpha })R_{\mathrm{g}}+{\mathit{\epsilon }}_{\mathrm{surf}}{R}_{\mathrm{atm}}-{\mathit{\epsilon }}_{\mathrm{surf}}\mathit{\sigma }{T}_{\mathrm{surf}}^{\mathrm{4}},\end{array}$$

where R_g is the incoming shortwave radiation (W m⁻²), R_atm is the incoming longwave radiation (W m⁻²), α is the albedo, T_surf is the surface emissivity, ε_surf is the surface temperature (^∘K) and σ is the Stefan–Boltzmann constant (W⁻² K⁴). The soil heat flux G depends on the soil type and water content as well as the vegetation type (Allen et al., 2005). The direct estimation of G by remote sensing data is not possible (Allen et al., 2011); however, empirical relations can estimate the fraction ξ = G∕R_n as a function of soil and vegetation characteristics using satellite image data, such as the LAI, NDVI, α and T_surf. Generally, G represents 5–20 % of R_n during daylight hours (Kalma et al., 2008). In order to estimate the G∕R_n ratio, several methods have been tested for various types of surfaces at different locations. The most common methods parameterize ξ as a constant for the entire day or at satellite overpass time (Ventura et al., 1999), according to NDVI (Jackson et al., 1987; Kustas and Daughtry, 1990), LAI (Choudhury et al., 1987; Kustas et al., 1993; Tasumi et al., 2005), vegetation fraction (f_c) (Su, 2002), T_surf and α (Bastiaanssen, 1995), or only T_surf (Santanello Jr. and Friedl, 2003). These empirical methods are suitable for specific conditions; therefore, estimating G, especially in this type of environment where NDVI values are low and thus G∕R_n values are large, is a critical issue. The approach adopted here was drawn on Danelichen et al. (2014), who evaluated the parameterization of these different models in three sites in Mato Grosso state in Brazil and found that the model proposed by Bastiaanssen (1995) showed the best performance for all sites, followed by the model from Choudhury et al. (1987) and Jackson et al. (1987):

Bastiaanssen (1995):

$$\begin{array}{ll}{\displaystyle }G=& {\displaystyle }{R}_{\mathrm{n}}\left(T_{\mathrm{surf}}-\mathrm{273.16}\right)(\mathrm{0.0038}+\mathrm{0.0074}\mathit{\alpha })\\ \text{(10)}& {\displaystyle }& {\displaystyle }\left(\mathrm{1}-\mathrm{0.98}{\mathrm{NDVI}}^{\mathrm{4}}\right)\end{array}$$

Choudhury et al. (1987):

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }G=\mathrm{0.4}R\mathrm{n}\left(\exp \right(-\mathrm{0.5}\mathrm{LAI}\left)\right)\end{array}$$

Jackson et al. (1987)

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }G=\mathrm{0.583}R\mathrm{n}\left(\exp \right(-\mathrm{2.13}\mathrm{NDVI}\left)\right).\end{array}$$

Hence, these three methods were tested for the Ben Salem flux station measurements, by comparing the measured G_BS-t and the computed G using measured $R_{{\mathrm{n}}_{\mathrm{BS}\text{-}t}}$, T_surf-BS-t, α_BS, NDVI_BS and LAI_BS at Terra and Aqua overpass times (results not shown). The best results are issued from the Bastiaanssen (1995) method with a root mean square error (RMSE) of 0.09 (average value of the two satellite overpass times) followed by Jackson et al. (1987) and Choudhury et al. (1987) with RMSE values of 0.15 and 0.2, respectively. Moreover, daily measured G_BS-day was computed and a G accumulation is generally found as has been already mentioned by Clothier et al. (1986), who showed that G is neither constant nor negligible on diurnal timescales and can constitute as much as 50 % of R_n over sparsely vegetated area. Since G estimation was the most uncertain variable, the three above methods were tested to compute the distributed remotely sensed AE. The Ben Salem meteorological station was used to provide $R_{{\mathrm{g}}_t}$ and R_atm-t. Remote sensing variables α, T_surf, ε_surf and NDVI came from MODIS products. Remotely sensed LAI was computed from the MODIS NDVI using a single equation (Clevers, 1989) for all crops in the study area:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }\mathrm{LAI}=-{\displaystyle \frac{\mathrm{1}}{k}}\ln \left({\displaystyle \frac{{\mathrm{NDVI}}_{\mathrm{\infty }}-\mathrm{NDVI}}{{\mathrm{NDVI}}_{\mathrm{\infty }}-{\mathrm{NDVI}}_{\mathrm{soil}}}}\right).\end{array}$$

The calibration of this relationship was done over the Yaqui irrigated perimeter (Mexico) during the 2007–2008 growing season using hemispherical LAI measured in all the studied fields (Chirouze et al., 2014). Calibration results gave the asymptotical values of NDVI, NDVI_∞ = 0.97 and NDVI_soil = 0.05, as well as the extinction factor k = 1.13. As this relationship was calibrated over a heterogeneous land surface but on herbaceous vegetation only, its relevance for trees was checked. For that purpose, clump-LAI measurements on an olive tree, as well as allometric measurements, i.e., mean distance between trees and mean crown size done using Pleiades satellite data (Mougenot et al., 2014; Touhami, 2013) were obtained. Clump LAI is the value of the LAI of an isolated element of vegetation (tree, shrub, etc.); if this element occupies a fraction cover f and is surrounded by bare soil, then the clump LAI value is equal to the area-average LAI divided by f. Hence, we checked that the pixels with tree-dominant cover show LAI values close to what was expected (of the order of 0.3 to 0.4 given the interrow distance of 12 m on average).

Remote sensed available energy was computed for the 10 km × 8 km MODIS subimages at Terra MODIS and Aqua MODIS overpass times, using the three methods estimating G. Since the measured heat fluxes H_XLAS_t represent only the weighted contribution of the fluxes from the upwind area to the tower (footprint), then instantaneous footprints at the times of Terra and Aqua overpass were selected among the two half-hours preceding and following the satellite's time of overpass (lowest time interval) and then was multiplied by the instantaneous remotely sensed available energy AE_t to get the available energy of the upwind area AE_t-FP.

3.3.2 Daily available energy

Most methods using TIR domain data rely on once-a-day, late morning (such as Terra-MODIS overpass time) or early afternoon (such as Aqua-MODIS overpass time) acquisitions. Thus, they provide a single instantaneous estimate of energy budget components. In order to obtain daily AE from these instantaneous measurements and to reconstruct hourly variations of AE, we considered that its evolution was proportional to another variable whose diurnal evolution can be easily known.

The extrapolation from an instantaneous flux estimate to a daytime flux assumes that the surface energy budget is “self-preserving”, i.e., the relative partitioning among components of the budget remains constant throughout the day. However, many studies (Brutsaert and Sugita, 1992; Gurney and Hsu, 1990; Sugita and Brutsaert, 1990) showed that the self-preservation method gives daytime latent heat estimates that are smaller than observed values by 5–10 %. Moreover, Anderson et al. (1997) found that the evaporative fraction computed from instantaneous measured fluxes tends to underestimate the daytime average by about 10 %; hence, a corrected parameterization was used and a coefficient = 1.1 was applied. Similarly, Delogu et al. (2012) found an overestimation of about 10 % between estimated and measured daily components of the available energy, and thus a coefficient = 0.9 was applied. The corrected parameterization proposed by Delogu et al. (2012) was tested, but this coefficient did not give consistent results; therefore, the extrapolation relationship was calibrated in order to get accurate daily results of AE.

Thereby, the applied extrapolation method was tested using in situ Ben Salem flux station measurements. The incoming short-wavelength radiation was used to scale available energy from instantaneous to daily values, but only for clear-sky days for which MODIS images can be acquired and remote sensing data used to compute AE are available. Clear-sky days were selected based on the ratio of daily measured incoming short-wavelength radiation $R_{{\mathrm{g}}_{\mathrm{day}}}$ to the theoretical clear-sky radiation Rso as proposed by the FAO-56 method (Allen et al., 1998). A day was defined as clear if the measured $R_{{\mathrm{g}}_{\mathrm{day}}}$ is higher than 85 % of the theoretical clear-sky radiation at the satellite overpass times (Delogu et al., 2012).

Daily measured available energy AE_BS-day, computed as the average of half-hourly measured AE_BS-30, was compared to daily available energy (AE_BS-day-Terra and AE_BS-day-Aqua) computed using the extrapolation method from instantaneous measured AE_BS-t-Terra and AE_BS-t-Aqua at Terra and Aqua overpass times, respectively (Eq. 14).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{AE}}_{\mathrm{BS}\text{-}\mathrm{day}\text{-}\mathrm{Terra}}=a_{\mathrm{Terra}}{R}_{{\mathrm{g}}_{\mathrm{day}}}{\displaystyle \frac{{\mathrm{AE}}_{\mathrm{BS}\text{-}t\text{-}\mathrm{Terra}}}{R_{{\mathrm{g}}_{t\text{-}\mathrm{Terra}}}}}+b_{\mathrm{Terra}},\\ \text{(14)}& {\displaystyle }& {\displaystyle }{\mathrm{AE}}_{\mathrm{BS}\text{-}\mathrm{day}\text{-}\mathrm{Aqua}}=a_{\mathrm{Aqua}}{R}_{{\mathrm{g}}_{\mathrm{day}}}{\displaystyle \frac{{\mathrm{AE}}_{\mathrm{BS}\text{-}t\text{-}Aqua}}{R_{{\mathrm{g}}_{t\text{-}\mathrm{Aqua}}}}}+b_{\mathrm{Aqua}},\end{array}$$

where $R_{{\mathrm{g}}_{\mathrm{day}}}$ is the daily measured incoming short-wavelength radiation in the Ben Salem meteorological station; $R_{{\mathrm{g}}_{t\text{-}\mathrm{Terra}}}$ and $R_{{\mathrm{g}}_{t\text{-}\mathrm{Aqua}}}$ are the instantaneous incoming short-wavelength radiation measured at Terra and Aqua overpass time, respectively, and AE_BS-t-Terra and AE_BS-t-Aqua are the instantaneous measured available energy in the Ben Salem flux station, at Terra and Aqua overpass times.

Results gave an overestimation of about 15 %. The corrected parameterizations of AE (Table 1), needed to remove the bias between measured (AE_BS-day) and computed AE (AE_BS-day-Terra and AE_BS-day-Aqua), were applied to compute daily remotely sensed AE (AE_day) from instantaneous AE (AE_t) following the extrapolation method shown in Eq. (14).

Then AE_day was multiplied by the corresponding daily footprint to get the daily available energy of the upwind area AE_day-FP. Finally, estimates of Terra and Aqua observed daily LE (LE_residual_XLAS_day-FP) were obtained based on the three methods used to compute G.

4 SPARSE model

4.1 Energy fluxes derived from SPARSE model

The SPARSE dual-source model solves the energy budgets of the soil and the vegetation. Here we use the “layer approach”, for which the resistance network relating the soil and vegetation heat sources to a main reference level through a common aerodynamic level use a series electrical branching. The main unknowns are the component temperatures, i.e., soil (T_s) and vegetation (T_v) temperatures. Totals at the reference height (the measurement height of the meteorological forcing), as well as the longwave radiation budget, are also solved so that altogether a system of five equations can be built:

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }\left\{\begin{array}{l}H=H_{\mathrm{s}}+H_{\mathrm{v}}\\ \mathrm{LE}={\mathrm{LE}}_{\mathrm{s}}+{\mathrm{LE}}_{\mathrm{v}}\\ R_{\mathrm{ns}}=G+H_{\mathrm{s}}+{\mathrm{LE}}_{\mathrm{s}}\\ R_{\mathrm{nv}}=H_{\mathrm{v}}+{\mathrm{LE}}_{\mathrm{v}}\\ {\mathit{\epsilon }}_{\mathrm{surf}}\mathit{\sigma }{T}_{\mathrm{surf}}^{\mathrm{4}}={\mathit{\epsilon }}_{\mathrm{surf}}{R}_{\mathrm{atm}}-R_{\mathrm{an}},\end{array}\right.\end{array}$$

where R_atm is the atmospheric radiation (W m⁻²), R_an is net longwave radiation which depends on T_s and T_v (W m⁻²), and T_surf and ε_surf are respectively the surface temperature (^∘K) and emissivity as observed by the satellite; indexes “s” and “v” designate the soil and the vegetation, respectively.

The first two (Eq. 15) express the continuity of the latent and sensible heat fluxes from the sources to the aerodynamic level through to the reference level, the third and the fourth (Eq. 15) are the soil and vegetation energy budgets, and the fifth (Eq. 15) relates the surface temperature T_surf derived from observed MODIS surface temperature to T_s and T_v.

Table 1 Corrected parameterizations of available energy for the diurnal reconstitution.

Figure 5 Flowchart of the SPARSE algorithm; T_s, H_ss, H_vs, LE_ss and LE_vs are surface temperature, soil sensible heat flux, vegetation sensible heat flux, soil latent heat flux and vegetation latent heat flux in fully stressed conditions, respectively; T_p, H_sp, H_vp, LE_sp and LE_vp are surface temperature, soil sensible heat flux, vegetation sensible heat flux, soil latent heat flux and vegetation latent heat flux in potential conditions, respectively.

The SPARSE model system of equations is fully described in Boulet et al. (2015). SPARSE is similar to the TSEB model (Kustas and Norman, 1999) but includes the expressions of the aerodynamic resistances of Choudhury and Monteith (1988) and Shuttleworth and Gurney (1990). This system can be solved in a forward mode for which the surface temperature is an output (prescribed conditions), and an inverse mode when the surface temperature is an input derived from satellite observations or in situ measurements in the thermal infrared domain (retrieval conditions). Figure 5 illustrates a diagram showing the flowchart of the model algorithm. System (15) is solved step by step by following similar guidelines as in the TSEB model: the first step assumes that the vegetation transpiration (LE_v) is at its maximum, and evaporation (LE_s) is computed. If this soil latent heat flux (LE_s) is below a minimum positive threshold for vegetation stress detection of 30 W m⁻², the hypothesis that the vegetation is unstressed is no longer valid. In that case, the vegetation is assumed to suffer from water stress and the soil surface is assumed to be already long dry. Then, LE_s is set to 30 W m⁻². This value accounts for the small but nonnegligible vapor flow reaching the surface (Boulet et al., 1997). The system is then solved for vegetation latent heat flux (LE_v). If LE_v is also negative, both LE_s and LE_v values are set to zero, whatever the value of T_surf. The system of equation can also be solved for T_s and T_v only if the efficiencies representing stress levels (dependent on surface soil moisture for the evaporation, and root zone soil moisture for the transpiration) are known. In that case the sole first four equations are solved. This prescribed mode allows all the fluxes in known limiting soil moisture levels to be computed (very dry, e.g., fully stressed, and wet enough, e.g., potential). It limits unrealistically high values of component fluxes, latent heat flux values above the potential rates or sensible heat flux values above that of a nonevaporating surface. The potential evaporation and transpiration rates used later on are computed using this prescribed mode with minimum surface resistance to evaporation and transpiration, respectively.

Some of the model parameters were remotely sensed data while others were taken from the bibliography or measured in situ. Remotely sensed data fed into SPARSE are T_surf, ε_surf, ϕ, NDVI, LAI and α. A grid of the vegetation height (z_v) was also necessary as input in the SPARSE model; for herbaceous crops, vegetation height was interpolated with the help of NDVI time series between fixed minimum (0.05 m) and maximum (0.8 m) values, while for trees, the roughness length (z_om) was linked to the allometric measurements (mentioned before) and computed as a function of canopy area index, drag coefficient and canopy height using the drag partition approach proposed by Raupach (1994) for tall sparse vegetative environments. Then, since SPARSE deals with vegetation height and not roughness length, the same simple rule of the thumb as the one used in SPARSE was used to reconstruct z_v for the tree cover types (z_v = z_om∕0.13). In a final step, to get spatial vegetation height, z_v was averaged over the MODIS pixels. In situ parameters used in SPARSE were mainly meteorological data: R_g, R_atm, T_a, H_a and u. No calibration was performed on the model parameters shown in Table 2.

Table 2 SPARSE parameters.

The retrieval and prescribed modes of the SPARSE model were run for the 10 km × 8 km subimages at the times of Terra and Aqua overpasses, to get instantaneous modeled fluxes H_SPARSE_t, LE_SPARSE_t and AE_SPARSE_t as well as sensible heat flux (H_s-t = H_ss-t + H_vs-t) in fully stressed conditions and latent heat (LE_p-t = LE_sp-t + LE_vp-t) and sensible heat (H_p-t = H_sp-t + H_vp-t) fluxes in potential conditions. Modeled values were then multiplied by the nearest half-hourly footprint to the satellite overpass time, in order to get fluxes corresponding to the upwind area: H_SPARSE_t-FP, LE_SPARSE_t-FP, AE_SPARSE_t-FP, H_s-t-FP, H_p-t-FP and LE_p-t-FP.

In a subsequent step, the retrieval and prescribed modes of SPARSE model was run at a half-hourly time step using the half-hourly meteorological measurements to get half-hourly latent heat flux at potential conditions LE_p-30 and half-hourly modeled available energy AE_SPARSE₃₀. The potential LE weighted by the corresponding half-hourly footprint (LE_p-30-FP) is used later when computing daily LE based on the stress factor method, while the half-hourly AE weighted by the corresponding half-hourly footprint (AE_SPARSE_30-FP) is used to compute daily LE based on the evaporative fraction method (Sect. 4.2).

4.2 Reconstruction of daily modeled ET from instantaneous latent heat flux

Daily ET is usually required for applications in hydrology or agronomy for instance, whereas most SEB methods provide a single instantaneous latent heat flux because the energy budget is only computed at the satellite overpass times (Delogu et al., 2012). In order to scale daily ET from one instantaneous estimate, there are various methods relying on the preservation, during the day, of the ratio of the latent heat flux to a scale factor with known diurnal evolution.

4.2.1 Stress factor (SF) method

The stress factor SF (Eq. 16) is assumed to be invariant during the same day, and the diurnal modeled fluxes are accounted for by recovering the diurnal course of potential ET.

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }\mathrm{SF}=\mathrm{1}-{\displaystyle \frac{\mathrm{LE}\mathrm{\_}{\mathrm{SPARSE}}_{t\text{-}\mathrm{FP}}}{{\mathrm{LE}}_{\mathrm{p}\text{-}t\text{-}\mathrm{FP}}}}\end{array}$$

The daily modeled ET (LE_SPARSE_day-FP) can be expressed as the product of the instantaneous estimate of SF at the satellite overpass times and the daily potential evapotranspiration:

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }\mathrm{LE}\mathrm{\_}{\mathrm{SPARSE}}_{\mathrm{day}\text{-}\mathrm{FP}}=(\mathrm{1}-\mathrm{SF}){\mathrm{LE}}_{\mathrm{p}\text{-}\mathrm{day}\text{-}\mathrm{FP}}.\end{array}$$

LE_p-day-FP was calculated as the sum of the half-hourly modeled latent heat fluxes at potential conditions LE_p-30-FP.

4.2.2 Evaporative fraction method

The evaporative fraction (EF) self-preservation is a valid assumption under dry conditions but no longer under wet conditions (Hoedjes et al., 2008). For these conditions, assuming a constant EF underestimates actual EF and therefore ET (Lhomme and Elguero, 1999). Indeed, according to Gentine et al. (2007), the diurnal shape of EF depends on both atmospheric forcing and surface conditions. Therefore EF was computed every 30 min using the following empirical parameterization (Delogu et al., 2012):

$$\begin{array}{}\text{(18)}& {\displaystyle }{\displaystyle }{\mathrm{EF}}_{\mathrm{30}}=\left[\mathrm{1.2}-\left(\mathrm{0.4}{\displaystyle \frac{R_{{\mathrm{g}}_{\mathrm{30}}}}{\mathrm{1000}}}+\mathrm{0.5}{\displaystyle \frac{{\mathrm{RH}}_{\mathrm{30}}}{\mathrm{100}}}\right)\right]\left({\displaystyle \frac{{\mathrm{EF}}_{\mathrm{SPARSE}\text{-}t}}{{\mathrm{EF}}_{\mathrm{met}\text{-}t}}}\right),\end{array}$$

where $R_{{\mathrm{g}}_{\mathrm{30}}}$ and RH₃₀ are respectively the half-hourly incoming short-wavelength radiation and relative humidity EF_SPARSE-t and EF_met-t are respectively SPARSE EF (Eq. 19) and computed EF using the evaporative fraction method (Eq. 20) at the satellite overpass times.

$$\begin{array}{}\text{(19)}& {\displaystyle }& {\displaystyle }{\mathrm{EF}}_{\mathrm{SPARSE}\text{-}t}={\displaystyle \frac{\mathrm{LE}\mathrm{\_}{\mathrm{SPARSE}}_{t\text{-}\mathrm{FP}}}{\mathrm{AE}\mathrm{\_}{\mathrm{SPARSE}}_{t\text{-}\mathrm{FP}}}},\text{(20)}& {\displaystyle }& {\displaystyle }{\mathrm{EF}}_{\mathrm{met}\text{-}t}=\left[\mathrm{1.2}-\left(\mathrm{0.4}{\displaystyle \frac{R_{{\mathrm{g}}_t}}{\mathrm{1000}}}+\mathrm{0.5}{\displaystyle \frac{{\mathrm{RH}}_t}{\mathrm{100}}}\right)\right],\end{array}$$

where LE_SPARSE_t-FP and AE_SPARSE_t-FP are respectively the latent heat flux and the available energy modeled by SPARSE at the satellite overpass times; $R_{{\mathrm{g}}_t}$ and RH_t are respectively the incoming short-wavelength radiation and relative humidity measured at the times of the satellite overpasses.

The half-hourly modeled ET (LE_SPARSE_30-FP) was computed as the product of the half-hourly EF estimate and the half-hourly modeled available energy AE_SPARSE_30-FP (Eq. 21). AE_SPARSE₃₀ was computed from instantaneous modeled available energy (AE_SPARSE_t) using the same approach detailed in Sect. 3.3.2 and applying Eq. (14) for a half-hourly time step (instead of a daily time step). AE_SPARSE₃₀ was weighted by the corresponding half-hourly footprint to get the modeled AE of the upwind area AE_SPARSE_30-FP. The daily modeled ET (LE_SPARSE_day-FP) was computed as the sum of the half-hourly LE_SPARSE_30-FP (Eq. 22).

$$\begin{array}{}\text{(21)}& {\displaystyle }& {\displaystyle }\mathrm{LE}\mathrm{\_}{\mathrm{SPARSE}}_{\mathrm{30}\text{-}\mathrm{FP}}={\mathrm{EF}}_{\mathrm{30}}\times \mathrm{AE}\mathrm{\_}{\mathrm{SPARSE}}_{\mathrm{30}\text{-}\mathrm{FP}}\text{(22)}& {\displaystyle }& {\displaystyle }\mathrm{LE}\mathrm{\_}{\mathrm{SPARSE}}_{\mathrm{day}\text{-}\mathrm{FP}}=\sum \mathrm{LE}\mathrm{\_}{\mathrm{SPARSE}}_{\mathrm{30}\text{-}\mathrm{FP}}\end{array}$$

4.2.3 Residual method

Besides, daily modeled ET (LE_SPARSE_day-FP) was also estimated as a residual term of the surface energy budget using daily modeled sensible heat flux (H_SPARSE_day-FP) and available energy (AE_SPARSE_day-FP) as follows:

$$\begin{array}{ll}{\displaystyle }\mathrm{LE}\mathrm{\_}{\mathrm{SPARSE}}_{\mathrm{day}\text{-}\mathrm{FP}}& {\displaystyle }=\mathrm{AE}\mathrm{\_}{\mathrm{SPARSE}}_{\mathrm{day}\text{-}\mathrm{FP}}\\ \text{(23)}& {\displaystyle }& {\displaystyle }-\mathrm{H}\mathrm{\_}{\mathrm{SPARSE}}_{\mathrm{day}\text{-}\mathrm{FP}}.\end{array}$$

H_SPARSE_day was computed from modeled sensible heat flux (H_SPARSE_t) following the same extrapolation method used for the available energy (see Sect. 3.3.2). The corrected parameterizations of H were obtained from the comparison of daily measured sensible heat flux H_BS-day computed as the average of half-hourly measured H_BS-30 and daily sensible heat flux (H_BS-day-Terra and H_BS-day-Aqua) computed using the extrapolation method from instantaneous measured H_BS-t-Terra and H_BS-t-Aqua at Terra and Aqua overpass times, respectively (Eq. 24).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{H}_{\mathrm{BS}\text{-}\mathrm{day}\text{-}\mathrm{Terra}}=a_{\mathrm{Terra}}^{\prime }{R}_{{\mathrm{g}}_{\mathrm{day}}}{\displaystyle \frac{H_{\mathrm{BS}\text{-}t\text{-}\mathrm{Terra}}}{R_{{\mathrm{g}}_{t\text{-}\mathrm{Terra}}}}}+b_{\mathrm{Terra}}^{\prime },\\ \text{(24)}& {\displaystyle }& {\displaystyle }{H}_{\mathrm{BS}\text{-}\mathrm{day}\text{-}\mathrm{Aqua}}=a_{\mathrm{Aqua}}^{\prime }{R}_{{\mathrm{g}}_{\mathrm{day}}}{\displaystyle \frac{H_{\mathrm{BS}\text{-}t\text{-}\mathrm{Aqua}}}{R_{{\mathrm{g}}_{t\text{-}\mathrm{Aqua}}}}}+b_{\mathrm{Aqua}}^{\prime },\end{array}$$

where H_BS-t-Terra and H_BS-t-Aqua are the instantaneous measured sensible heat flux in the Ben Salem flux station.

Therefore, the corrected parameterizations of H (Table 3), needed to remove the bias between measured (H_BS-day) and computed H (H_BS-day-Terra and AE_BS-day-Aqua), were applied to compute daily modeled H (H_SPARSE_day) from instantaneous modeled H (H_SPARSE_t) following the extrapolation method shown in Eq. (21). Finally, H_SPARSE_day was weighted by the corresponding daily footprint to get the daily modeled H of the upwind area H_SPARSE_day-FP.

Table 3 Corrected parameterizations of sensible heat flux for the diurnal reconstitution.

5 Water stress estimates

Water stress estimation is crucial to deduce the root zone soil moisture level using remote sensing data (Hain et al., 2009). Water stress results in a drop of actual evapotranspiration below the potential rate. Its intensity is usually represented by a stress factor as defined in Sect. 4.2, ranging between 0 (unstressed surface) and 1 (fully stressed surface).

Modeled values of SF at the times of Terra and Aqua overpasses (SF_mod) have been computed from modeled potential LE (LE_p-t-FP) as follows:

$$\begin{array}{}\text{(25)}& {\displaystyle }{\displaystyle }{\mathrm{SF}}_{\mathrm{mod}}=\mathrm{1}-{\displaystyle \frac{\mathrm{LE}\mathrm{\_}{\mathrm{SPARSE}}_{t\text{-}\mathrm{FP}}}{{\mathrm{LE}}_{\mathrm{p}\text{-}t\text{-}\mathrm{FP}}}}\end{array}$$

where LE_SPARSE_t-FP and LE_p-t-FP are the modeled latent heat fluxes in actual and potential conditions, respectively.

Furthermore, surface water stress factor derived from XLAS measurement, named SF_obs, at the times of Terra and Aqua overpasses was computed as follows (Su, 2002):

$$\begin{array}{}\text{(26)}& {\displaystyle }{\displaystyle }{\mathrm{SF}}_{\mathrm{obs}}={\displaystyle \frac{\mathrm{H}\mathrm{\_}{\mathrm{XLAS}}_t-H_{\mathrm{p}\text{-}t\text{-}\mathrm{FP}}}{H_{\mathrm{s}\text{-}t\text{-}\mathrm{FP}}-H_{\mathrm{p}\text{-}t\text{-}\mathrm{FP}}}},\end{array}$$

where H_s-t-FP and H_p-t-FP are the modeled sensible heat flux in actual and potential conditions, respectively, and H_XLAS_t is the XLAS sensible heat flux at the satellite overpass times.

6 Results and discussion

6.1 XLAS- and model-derived instantaneous sensible heat fluxes

Our primary focus is the comparison between scintillometer measurements and the modeled sensible heat fluxes computed using the Terra and Aqua remotely sensed data. The scintillometer H at the times of the two satellite overpasses (H_XLAS_t) are interpolated from the half-hourly H measurements. Heat flux determination was possible for typically about 87 % of the daytime measurements during the summer, availability of XLAS heat flux values was lower during the cold season due to poor visibility and/or stable stratification.

H_SPARSE was weighted by the XLAS footprint in order to be able to compare the modeled values (H_SPARSE_t-FP) with the XLAS measurements (H_XLAS_t). Therefore, due to XLAS and remote sensing data availability, we got 175 and 118 values for Terra and Aqua respectively. In order to highlight H interseasonality between the drier 2012–2013 and the wetter 2013–2014 seasons, we present an example of 2 days, each in one season: DOY 2013-082 shows a H value ranging between 163 and 342 W m⁻² while DOY 2014-208 shows a H value ranging between 97 and 311 W m⁻² (Fig. 6). The colored area shows the modeled flux and the contours shows the surface source area contributing to the scintillometer measurements. Day 2013-82 (23 March 2013) is chosen in the cold season while day 208-2014 (27 July 2014) is in the warm season to focus on land cover impact on T_surf and thus on modeled H (trees and cereals in winter vs. only irrigated trees and market gardening in summer). Moreover, the first day experiences a strong southern wind while there is a light northern wind during the second day. Generally, a small number of MODIS pixels brings a high contribution to the signal; among these, two are hot pixels (pixel with high T_surf and low NDVI) in which the land use is mainly arboriculture.

Prediction performance is assessed using RMSE and the coefficient of determination (R²). Results for the sensible heat flux are illustrated in Fig. 7 and show good agreement between modeled and measured H at the times of satellite overpasses. This is illustrated by linear regressions of H_SPARSE_t-FP = 1.065 H_XLAS_t − 14.788 (R² = 0.6; RMSE = 57.89 W m⁻²) and H_SPARSE_t-FP = 1.12 H_XLAS_t − 10.57 (R² = 0.63; RMSE = 53.85 W m⁻²) for Terra and Aqua, respectively. This result is of great interest considering that the SPARSE model was run with no prior calibration. However, we noted that bias is a function of the flux level and most outliers are recorded for H greater than 200 W m⁻². This can be explained by (i) the XLAS measurement saturation (according to the “Kipp & Zonen LAS and XLAS instruction manual” (Kipp & Zonen, 2017), for a path length of 4 km and a scintillometer height of 20 m, the saturation measurement problems start from H values higher than 300 W m⁻²), (ii) uncertainties on the correction of stability using the universal stability function, and (iii) potential inconsistencies between the area average MODIS surface temperature and the air temperature measured locally at the meteorological station.

Whereas there are several studies dealing with large aperture scintillometer data whose measurements are compared to modeled fluxes, in the few studies dealing with extra large aperture scintillometer data, the comparison is generally done with eddy covariance station measurements (Kohsiek et al., 2002; Moene et al., 2006). Indeed, our results are in agreement with those found by Marx et al. (2008), who compared LAS-derived and satellite-derived H (SEBAL was applied with NOAA-AVHRR images providing maps of surface energy fluxes at a 1 km × 1 km spatial resolution) and found that modeled H is underestimated with a RMSE of 39 W m⁻² for the site Tamale and 104 W m⁻² for the site Ejura. Moreover, Watts et al. (2000) compared the satellite (AVHRR radiometer) estimates of H to those from LAS over semiarid grassland in north-western Mexico during the summer of 1997. They found RMSE values of 31 and 43 W m⁻² for LAS path lengths of 300 m and 600 m respectively and showed that LAS measurements are less good than those derived from a 3-D sonic anemometer. They also suggested longer LAS path length (greater than 1.1 km) since the LAS is rather insensitive to the surface near the receiver and the emitter.

Figure 6 Model-derived sensible heat fluxes and footprints for (a) DOY 2013-082 at Aqua overpass time and (b) DOY 2014-208 at Terra overpass time. The colored area shows the modeled flux and the contours shows the surface source area contributing to the scintillometer measurements.

Figure 7 Modeled vs. observed sensible heat fluxes at Terra and Aqua overpass times.

Figure 8 Modeled vs. observed latent heat fluxes at Terra and Aqua overpass times.

6.2 XLAS- and model-derived instantaneous latent heat fluxes

In a subsequent step, SPARSE-derived LE (LE_SPARSE_t-FP) was compared to observed LE (LE_residual_XLAS_t-FP). Results are illustrated in Fig. 8, showing a good agreement between modeled and observed LE. However, these results are less good than for the H results, as shown by the following linear regressions: LE_SPARSE_t-FP = 0.94 LE_residual_XLAS_t-FP + 12.47 (RMSE = 47.20 W m⁻²) and LE_SPARSE_t-FP = 0.85 LE_residual_XLAS_t-FP + 11.51 (RMSE = 43.20 W m⁻²) for Terra and Aqua respectively, with an overall R² of 0.55 for both satellites. We note a greater scatter for latent heat flux than for the sensible heat flux (Fig. 7), which can be explained by the fact that LE is here a residual term affected by estimation errors in both AE and H. Despite this moderate discrepancy, the good agreement between both approaches indicates that the methodology adopted in SPARSE for estimating H and AE using MODIS imagery is appropriate for modeling latent heat fluxes.

Figure 9 Modeled vs. XLAS-derived stress index SF at Terra and Aqua overpass times.

Figure 10 Modeled and observed stress index evolution at Terra and Aqua overpass times compared to daily rainfall.

6.3 Water stress

The scattered values of the stress factor as shown in Fig. 9 are consistent with previous studies such as Boulet et al. (2015). SEB retrieval of stress is limited by the scale mismatch between the instantaneous estimate of the surface temperature during the satellite overpass (which can be influenced by high-frequency turbulence) and the aggregated values of other forcing data which are derived from half-hourly averages (Lagouarde et al., 2013, 2015). However, general tendencies are well reproduced, with most points located within a 0.2 confidence interval (illustrated by dotted lines along the 1 : 1 line) as found by Boulet et al. (2015) at field scale, which is encouraging from the perspective of assimilating ET or SF in a water balance model for example. Moreover, it is noted that results include small LE and LE_p values with the same order of magnitude as the measurement uncertainty itself. Most outliers with greater water stress (∼ 1) correspond to high evaporation from bare soil since the dominant land use in the study area is arboriculture, but also, this could be due to saturation of scintillation which led to an underestimation of H XLAS measurements, as pointed out by Frehlich and Ochs (1990) and Kohsiek et al. (2002).

Figure 11 Modeled vs. observed daily latent heat fluxes. Dark grey color shows minimum and maximum daily observed LE. Light grey vertical bars show gaps in XLAS data. Error bars for the modeled ET show the minimum and the maximum daily ET resulting from the three methods used to compute daily ET from instantaneous modeled ET.

Modeled and observed stress indexes at Terra and Aqua overpass times show a consistent evolution with daily rainfall (Fig. 10), although the modeled stress shows a greater dispersion than the observed one. During a rainy episode (or an eventual irrigation period), the surface temperature decreases towards the unstressed surface temperature, thus marking an unstressed state, and SF tends to 0. Conversely, after a long dry down, the water stress appears and the surface temperature increases towards the equilibrium surface temperature computed by SPARSE under stressed conditions, and SF tends towards 1. Besides, it is noted that modeled stress indexes computed on the basis of Aqua MODIS's T_surf are often greater than those computed using Terra MODIS's T_surf due to higher T_surf (higher global solar radiation) at the time of Terra overpass (around midday).

6.4 XLAS- and model-derived daily latent heat fluxes

Daily observed ET, i.e., LE_residual_XLAS_day-FP, was computed using the residual method; hence, six estimates of the daily observed ET were obtained by combining the two satellite datasets and three methods to compute G and thus AE (see Sect. 3.3). Only the residual method was used to estimate daily observed ET for two reasons; on the one hand, to reduce the computations approach since, already, three methods to compute AE have been tested and, on the other hand, because the application of the EF method was not possible since we do not have a measured spatially distributed potential evapotranspiration (only point potential evapotranspiration data at the Ben Salem meteorological station are available). From daily observed ET estimates, minimum and maximum ET were selected for each day and minimum and maximum daily ET time series were interpolated between successive days based on the self preservation of the ratio of AE to R_g as scale factor (Fig. 11).

In addition, three methods were used to compute SPARSE daily ET for the Terra and Aqua overpasses (see Sect. 4.2), providing six estimates of the daily modeled ET. For each day average ET was plotted (260 days) with error bars showing minimum and maximum values, along with precipitation to understand the rainfall impact on the ET evolution (Fig. 11).

Despite the uncertainty in reconstructing the daily ET from instantaneous ET, overall results show a good agreement between XLAS-derived and SPARSE-derived ET values with similar seasonal dynamics. Daily observed and modeled ET over the whole study period were both in the range of 0–4 mm day⁻¹ with an RMSE of 0.7 mm day⁻¹, which is consistent with the land use present in the XLAS path: mainly trees spaced by a considerable fraction of bare soil and fewer herbaceous soil-covering crops (see Sect. 3.2). As expected, ET rates decrease significantly during dry periods (summers) since arid conditions limit the latent heat flux in favor of sensible heat flux and increase immediately after rainfall events due to the high amount of water evaporated from soil. The rainfall peaks that occurred on 3 September 2013 (about 10 mm), 6 October 2013 (about 20 mm), 15 March 2014 (about 100 mm) and 22 April 2014 (about 25 mm) are followed by well-reproduced dry downs.

At seasonal scale, we note a good agreement between modeled and observed daily ET for the 2013–2014 and 2014–2015 seasons, especially when vegetation cover was more developed: from March to July 2014 and from March to May 2015; these periods correspond to cereal vegetation peaks in some plots (March–April) and to market gardening crops (e.g., tomato, water melon, pepper) cultivated generally from spring to the beginning of autumn in the interrow area of trees plots, which is a common farming practice in the Kairouan plain. However, the 2012–2013 season was dry compared with the two other ones, and less accurate results were obtained. Some points with little to null ET were recorded from May to July 2013 which can be explained by the very dry conditions and scattered vegetation cover with a considerable amount of bare soil. This behavior was not observed in the same period of 2014, because 2014 was a rainy year in comparison to 2013; therefore, even supposing that the farmers have the same attitude and cultivate the same crop types between the two years (which is not true in the context of our study area and farmers always change crop types), precipitation favors the growth of spontaneous vegetation over fallows which contribute to ET rise. However, since this year experiences more rain, farmers cultivate a larger part of the land and diversify the crop types; the vegetation cover is denser and contributes to an overall increase in ET. Overall, lower ET values are recorded in autumn (October and November), which correspond to evapotranspiration from trees only, since the latest summer crops (market gardening crops) have been already harvested and the winter crops (mainly cereals) are not yet sown.

Moreover, it can be seen that occasionally SPARSE overestimated ET. For example, three dates can be selected in August 2013 (15, 25 and 29 August 2013) for which modeled ET were 3.30, 3.80 and 2.80 mm while maximum observed ET were 2.0, 2.40 and 1.20 mm, respectively; broader amplitude between modeled (4.00 mm) and observed ET (1.40 mm) was also recorded on 18 May 2013. SPARSE also overestimates ET throughout 10 days in August 2014 with an average difference of 1.1 mm and a maximum difference of 1.60 mm recorded in 23 August 2014. These discrepancies are always recorded under wet conditions (minimum stress factor) which show the difficulty in accurately representing the conditions close to the potential ET. This might be related to the theoretical limit of the model for low vegetation stress, especially when coupled with low evaporation efficiencies (i.e., dry soil surface), as already reported by Boulet et al. (2015) for senescent vegetation. The average difference between SPARSE and XLAS-derived LE estimates when both are available indicate that SPARSE can predict evapotranspiration with accuracies approaching 5 % of that of the XLAS.

7 Conclusions

This study evaluated the performances of the SPARSE model forced by MODIS remote sensing products in an operational context (no model calibration) to estimate instantaneous and daily evapotranspiration. The validation protocol was based on an unprecedented dataset with an extra large aperture scintillometer. Indeed, to our knowledge, this is the first work based on XLAS measurements acquired during the course of more than 2 years, as compared to 3 months in previous works (Kohsiek et al., 2002; Moene et al., 2006). The estimates of the sensible heat flux derived from the SPARSE model are in close agreement with those obtained from the XLAS. These results indicate that the XLAS can be fruitfully used to validate large-scale sensible heat flux derived from remote sensing data (and residual latent heat flux), in particular for the results obtained at the satellite overpass times, providing a feasible alternative to local micrometeorological techniques for measuring the sensible heat flux and validating satellite-derived estimates (i.e., eddy correlation). Furthermore, the extrapolation from instantaneous to daily evapotranspiration is less obvious and three methods were tested based on the stress index, the evaporative fraction and the residual approach. The daily latent heat fluxes derived from the XLAS agreed rather well with those modeled using the SPARSE model, which shows the potential of the SPARSE model in water consumption monitoring over heterogeneous landscapes in semiarid conditions, and especially to locate areas most affected by water stress. However, the precision in ET prediction with the SPARSE model is restricted by several assumptions and uncertainties: for instance, the instantaneous remote sensing data, mainly T_surf which is paramount in stress coefficient computation, are assumed to be reliable. Moreover, there is an issue with the MODIS pixel heterogeneity and notably the distribution of components at the intersection between the square pixel and the XLAS footprint. Uncertainties are also due to half-hourly forcing (meteorological and flux data) and XLAS data as well as to the extrapolation method from instantaneous to daily results. Furthermore, the empirical estimation methods of soil heat flux G (three methods were tested) as well as the possible daily heat accumulation lead to possible errors in available energy estimation and in turn in residual LE estimation.

Even if overall results are encouraging, further work is needed to improve results by (i) being most efficient in the SPARSE model application using calibrated input data specific to our study area, especially input parameters to which the model is particularly sensitive such as the mean leaf width and the minimum stomatal resistance; (ii) taking into account the heterogeneity of the 1 km MODIS pixel by applying MODIS footprint, which is determined by the sensor's observation geometry; and (iii) using a land surface model applied at the field scale (Etchanchu et al., 2017) to analyze the scaling properties from the field to the footprint of the XLAS and the MODIS pixels similarly.

Finally, in a future work, we plan to take advantage of the complementarities between the soil water balance and surface energy balance approaches (i.e., continuous but uncertain estimates using SWB due to poor soil water content control on the one hand and sensitivity of SEB to the actual water stress on the other hand) to implement an assimilation scheme of the remotely sensed surface temperature into land surface models. In fact, in order to provide further information about distributed soil water status over the studied areas, the TIR-derived evapotranspiration products could be assimilated directly either in land surface or hydrological models.