Combining a Two Source Energy Balance Model Driven by MODIS and MSG-SEVIRI Products with an Aggregation Approach to Estimate Turbulent Fluxes over Sparse and Heterogeneous Vegetation in Sahel Region (Niger)

Estimates of turbulent fluxes (i.e., sensible and latent heat fluxes H and LE) over heterogeneous surfaces is not an easy task. The heterogeneity caused by the contrast in vegetation, hydric and soil conditions can generate a large spatial variability in terms of surface–atmosphere interactions. This study considered the issue of using a thermal-based two-source energy model (TSEB) driven by MODIS (Moderate resolution Imaging Spectroradiometer) and MSG (Meteosat Second Generation) observations in conjunction with an aggregation scheme to derive area-averaged H and LE over a heterogeneous watershed in Niamey, Niger (Wankama catchment). Data collected in the context of the African Monsoon Multidisciplinary Analysis (AMMA) program, including a scintillometry campaign, were used to test the proposed approach. The model predictions of area-averaged turbulent fluxes were compared to data acquired by a Large Aperture Scintillometer (LAS) set up over a transect about 3.2 km-long and spanning three vegetation types (millet, fallow and degraded shrubs). First, H and LE fluxes were estimated at the MSG-SEVIRI grid scale by neglecting explicitly the subpixel heterogeneity. Moreover, the impact of upscaling the model’s inputs was investigated using in-situ input data and three aggregation schemes of increasing complexity based on MODIS products: a simple averaging of inputs at the MODIS resolution scale, another simple averaging scheme that considers scintillometer footprint extent, and the weighted average of inputs based on the footprint weighting function. The H and LE simulated using the footprint weighted method were more accurate than for the two other aggregation rules despite the heterogeneity of the landscape. The statistical values are: correlation coefficient (R) = 0.71, root mean square error (RMSE) = 63 W/m2 and mean bias error (MBE) = −23 W/m2 for H and an R = 0.82, RMSE = 88 W/m2 and MBE = 45 W/m2 for LE. This study opens perspectives for the monitoring of convective and evaporative fluxes over heterogeneous landscape based on medium resolution satellite products.


Introduction
The Sahelian region is one of the most vulnerable areas to climate variability and change due to the scarcity of water resources, associated with a high evaporative demand and low and irregular rainfall. The understanding of the interaction between regional climate and the hydrological cycle has substantially advanced through several research projects conducted over this region in recent years (HAPEX-Sahel [1] and AMMA [2]). Among the key surface processes impacting thermo-hydric characteristics of the low atmosphere, the regional exchange of water and energy at the soil-vegetation-atmosphere continuum is of prime importance for weather forecasting and climate studies [3][4][5].
Within this context, several land surface models have been developed to estimate the evapotranspiration (ET). These models are classified into three different categories: (i) conceptual models that are diagnostic and do not incorporate temporal information about the surface state such as the FAO-56 approach [6] which provides only a daily estimation of the evapotranspiration by modeling the evaporative function based on an empirical coefficient; (ii) surface energy balance (SEB) models which calculate ET as a residual term of the surface energy balance (TSEB [7], Mapping EvapoTranspiration at high Resolution with Internalized Calibration (METRIC) [8], and SEBS [9]); and (iii) the mechanistic models which resolve both water and energy budget equations in coupled manner and are often called soil-vegetation-atmosphere transfer models (SVAT models) such as Suivi de l'Etat Hydrique des Sols (SEtHyS) [10], Interaction Soil Biosphere Atmosphere (ISBA) [11], a Simple Soil-Plant-Atmosphere Transfer model (SiSPAT) [12] and the Interactive Canopy Radiation Exchange (ICARE) [13]. The latter can simulate directly the temporal evolution of two surface state variables (soil moisture and land surface temperature) and thus ET. However, contrary to the SEB models, SVAT models need more inputs parameters and data which are not always measurable, especially over heterogeneous surfaces. Especially, a good knowledge of water inputs by precipitation and irrigation is not needed for SEB approaches. This represents a significant advantage in semi-arid areas characterized by sparse meteorological networks and high spatio-temporal variability of precipitation.
In addition to meteorological data, the surface main driver of SEB models is the land surface temperature (Ts). For several applications, in particular large-scale hydrological studies, the spatial heterogeneity of the processes, which results in variability of Ts and other input parameters and variables (such as albedo, LAI, the roughness length (z_0), etc.), needs to be carefully characterized to predict correctly regional H and LE. This is not possible on heterogeneous areas with network of point-sampling devices such as the precision infrared temperature sensor (IRTS-Ps and Apogee) for surface temperature because of the limited representivity of these in-situ measurements. In this regard, remotely sensed data can be a valuable tool to provide representative values of these measurements [14,15]. With the existing sensors, the only Ts products able to sample the high intra-and inter-daily dynamics of Ts are provided either by geostationary or by medium-scale sun-synchronous sensors. Indeed, geostationary sensors such as the Meteosat Second Generation-The Spinning Enhanced Visible and InfraRed Imager (MSG SEVIRI) can provide regional scale of Ts with temporal sampling from 15 min to 1 h with a resolution of about 3 km. The use of MSG SEVIRI data to monitor vegetation or land surface vegetation related processes has been limited until now. Sandholt et al. [16] used SEVIRI data as an input to a regional hydrological model; Fensholt et al. [17] investigated diurnal SEVIRI NDVI for the African continent; and Stisen et al. [18] estimated regional ET from SEVIRI data. More recent studies have demonstrated the potential of SEVIRI data for evapotranspiration modeling [19][20][21] and drought and water stress monitoring [22,23]. In contrast with MSG SEVIRI, sun-synchronous satellites such as the MODIS sensor (Moderate resolution Imaging Spectroradiometer) provide data with a spatial resolution of about 1 km and a daily revisit time. Within this context, the question arises if LE over heterogeneous landscape can be monitored directly at large scale using geostationary data or if sub-pixel heterogeneity should be described, in particular by taking advantage of the finer resolution of sun-synchronous products and aggregation approaches.
In this study, a SEB model driven by MSG SEVIRI and MODIS observations in conjunction with an aggregation scheme is used to estimate H and LE over a transect of about 3.2 km located at the Wankama Catchment (Niger) [24][25][26][27][28][29]. A spatial aggregation approach is conceived as a method which seeks to link the model parameters that control surface exchange (Ts, α, d, z_0, LAI, and hc) on a patch scale with the area-average value of equivalent model parameters applicable on a larger scale or grid-scale, assuming that the same equations are used to describe surface fluxes at both scales [27,[30][31][32][33][34]. In general, aggregation schemes can be applied on the input requirements [27] or directly to predicted fluxes [35]. Nevertheless, the efficiency of this approach lies on the evaluation of their outputs/performances against field observations, which require the development of a measurement network at the local scale such as Eddy covariance systems [36][37][38][39]. Such a system is very costly, and manpower demanding with a competent staff for data processing and maintenance, as well as an important energy supply for operation. To overcome these difficulties, Large Aperture Scintillometer (LAS) is one of the alternative techniques employed to validate the satellite remote sensing H and LE estimates due to the comparable spatial resolutions. It has already been used in various contexts over complex topography and heterogeneous vegetation cover [27,40,41]. Consequently, LAS is becoming popular in hydro-meteorological studies because it is relatively cheap, robust, and easy to operate and maintain. Additionally, the LAS can be potentially used to improve the representation of surface heterogeneity in land-surface-atmosphere models operating at large scales.
In this study, the thermal-based two-source energy balance (TSEB) [7] was chosen as it has been evaluated with success in different contexts including sparse canopy in semi-arid areas [42][43][44]. As far as we know, the TSEB model has never been tested over a heterogeneous Sahelian agro-ecosystem. Therefore, a validation at patch scale was firstly investigated using in-situ measurements from three eddy-covariance stations sampling the dominating canopy (millet, fallow and degraded shrubs) of the studied area. TSEB predictions of the convective fluxes were then evaluated at the scintillometer footprint scale by aggregating the station scale predictions fed by the in-situ observations of albedo, Ts and Leaf area index (LAI). In a second step, remote sensing products were used. A first evaluation was directly conducted at the scale of the scintillometer measurements by using 3 km Ts observations from the MSG SEVIRI instrument. Finally, three aggregation schemes of MODIS products were tested to assess if a better representation of heterogeneity based on the 1-km products improves our large-scale estimates.

Data and Methods
This section describes the experimental data and the satellite products, the two-source energy balance model together with the aggregation approaches. All variables used in this study are listed in Table 1. The region of interest is part of the Wankama catchment, which is located 70 km east of the Niamey city, Niger ( Figure 1). It is situated between an upland at an altitude of 255 m and a pond at 200 m [45]. This semi-arid site belongs to the AMMA-CATCH-Niger observatory [46], one of three meso-sites along the West African latitudinal transect [47]. The climate is typically Sahelian with a short rainy season from June to September with an annual mean of 560 mm for the years 1905-2004, and very strong year-to-year and spatial variability. The site is characterized with high temperatures throughout the year with a daily average which ranges from 24 • C to 35 • C. Additionally, the potential evapotranspiration is about 2500 mm/year [48], greatly exceeding the annual rainfall. The soil is sandy (>90% sand) and poor in nutrients [49]. The study area is covered by: (1) Millet fields cover 58% of the catchment's surface area. This rainfed crop is cultivated traditionally with little chemical fertilizer and no pesticides, being by far the main crop grown. (2) Fallow savannah fields represent 23% of the total surface area and they are no more than five years old and typical of the Niamey region.
(3) Degraded shrubs occupy the remaining area [26]. Vegetation usually emerges after the first rainfall events occurring in June or early July, peaks around the end of the rainy season in early September and then dries out during the senescence phase.

Patch Scale
The experiment was conducted between 23 July and 23 October 2006. The site was instrumented by: three eddy covariance (EC) stations, measuring surface fluxes of sensible heat, latent heat (i.e., evapotranspiration), momentum and CO 2 ; CNR1 radiometer sensors to measure the four components of solar radiation, and soil heat plate fluxes; and Time Domain Reflectometry (TDR) probes, for the soil moisture profiles of the three main land cover types (millet crop, fallow savannah and degraded shrubs) ( Figure 2). Over the millet and fallow sites, belonging to the AMMA-CATCH observatory [46], the EC system consisted of a 3D sonic anemometer (CSAT3, Campbell Scientific Ltd., Loughborough, UK) which measured the fluctuations of the wind velocity components and temperature, and an open-path infra-red gas analyzer (LICOR-7500, Campbell Scientific Ltd.) that measured concentration of water vapor and carbon dioxide. Raw data were sampled at a rate of 20 Hz and processed with the EdiRe software (Version 1.4.3.1167, R. Clement, University of Edinburgh), based on CarboEurope recommendations [50], including despiking, double rotation, cross-correlation for derivation of time lag between the sonic anemometer and the gas analyzer, spectral corrections, Webb corrections and atmospheric stability test, with no gap filling [26]. For the degraded shrubs site, EC system consisted of a CSAT3 and a Krypton hygrometer (KH20, Campbell Scientific Ltd.) and the fluxes were calculated using the ECpack software after performing planar fit corrections [51], correcting the sonic temperature for the presence of humidity [52], frequency response corrections for slow apparatus and path length integration [53], the inclusion of the mean vertical velocity according to Webb et al. [54] and oxygen correction for the Krypton hygrometer which is sensitive to O 2 [55]. Ezzahar et al. [27] compared the Ecpack and EdiRe outputs and they concluded that the difference between the two software packages is not significant (RMSE (root mean square error) = 6.7 W m −2 and R 2 = 0.99).
These automatic measurements have been complemented by manual vegetation characteristics estimates about every 15 days including LAI (using hemispherical photographs) and height for millet and fallow sites [24]. Timings of the season was quite similar with a beginning occurring quite late about end of July. Leaf Area Index values were generally low for both sites with peaks of 0.9 m 2 m −2 at the end of August and of 0.24 m 2 m −2 at mid-September for fallow and millet sites, respectively.

Grid Scale
A Large Aperture Scintillometer (LAS) was set up over a 3.2 km transect spanning the three ecosystem types, i.e., millet and fallow fields and the degraded shrubs. The LAS used in this study was developed and built by the Meteorology and Air Quality Group from the Wageningen University (Netherlands). This instrument was constructed according to the basic design described in Ochs and Wilson [56]. It has an aperture size of 0.15 m and the transmitter operates at a wavelength of 0.94 µm. At the receiver, C 2 n is sampled at 1 Hz and averaged over 1 min intervals by a CR510 datalogger [57]. The transmitter and the receiver were installed on 10 m towers with an altitude difference of approximately 46 m. The receiver was installed at the highest part of the basin (upland) whilst the transmitter was installed at the lowest part of the basin ( Figure 2). The direction of the LAS path was 250 • from the north. At the receiver tower, a 2D anemometer (Wp200, R.M. Young Company, Traverse City, MI, USA) to measure wind speed and direction at 10 m was also installed. Additionally, the air temperature and humidity were measured using Vaisala HMP45C probe.

Scintillometry Theory
The transmitter emits electromagnetic radiation, which is scattered by the turbulent atmosphere, and the resulting variations in signal intensity (scintillations) are recorded by a receiver comprising an identical mirror and a photodiode detector. The intensity fluctuations are related to the path averaged structure parameter of the refractive index of air C 2 n . The scintillations are primarily the result of fluctuations in air temperature and humidity. Strictly speaking, the measured C 2 n is related to the structure parameters of temperature C 2 T , humidity C 2 q , and the covariant term C Tq . For electromagnetic waves in the visible and near-infrared region, however, humidity related scintillations are much smaller than temperature related scintillations. Wesely [58], and more recently Moene [59], demonstrated that for a LAS operating at a near-infrared wavelength, the C 2 T is related to C 2 n as follows: where T a is the air temperature, p is atmospheric pressure and β is the Bowen ratio. The factor involving the Bowen ratio is the correction term for the influence of humidity fluctuations. C 2 n and C 2 T are in (m −2/3 ) and (K 2 m −2/3 ) respectively. Using the Monin-Obukhov Similarity Theory (MOST), the sensible heat flux (H LAS ) can be obtained from C 2 T and an additional wind speed data through the following dimensionless relationship: and U * is friction velocity expressed as: where z LAS is the effective height of the LAS above the surface [60], Ψ is the integrated stability function [61], d is the displacement height, z 0 is the roughness length, k is the von Karman constant, g is the gravitational acceleration, ρ is the density of air and c p is the specific heat of air at constant pressure. During the iteration procedure, the Bowen ratio is evaluated using the H LAS , net radiation (R n ), and soil heat flux Several forms of the function f T have been proposed (e.g., [62][63][64][65]). In this study, we adopted [64]. Finally, the latent heat flux (LE LAS ) from the LAS can be derived by imposing the energy-balance closure [27,30,[66][67][68][69][70][71]].

Footprint Model
The footprint model describes the measured flux as the integral of the contributions from the entire source area. This source can be determined using simple analytical models [72] or more complex Lagrangian models [73]. In this study, the Horst and Weil model, which has been adapted to scintillometer measurements, was used to evaluate the source surface. It has been used successfully in previous studies [27,35].
The footprint function f relates the flux measured at z m , F(x, y, z m ) to the spatial distribution of surface fluxes, F(x,y,z 0 ) = F 0 (x,y), that is: where x and y, respectively, are the upwind and crosswind distances (m) from the point where the measurements are taken. The source surface is calculated by integration of the footprint function.
The footprint function f is calculated using the Horst and Weil model [72]: where z is the mean plume height for diffusion from a surface source (m), z m is the measurement height, and u(z) is the mean wind speed profile. Variables A, b, and c are functions of the gamma function of shape parameter r. Horst and Weil [72] found the following relationship for the diffusion in the lateral direction (which is commonly assumed to be Gaussian): Here, σ y is the standard deviation of the lateral spread, which is related to the plume travel time x/U (where U is the plume advection speed) and the standard deviation of the lateral wind fluctuations σ v : The three-dimensional footprint function for scalar fluxes is obtained: To assess the accuracy of the LAS, a footprint analysis was made. The wind direction pattern showed that only the interval ranging from 70 • to 250 • is considered during the study period. This is because the footprint of the LAS covers the monitored sites where the EC systems are installed. The EC measured fluxes are thus more or less comparable with the LAS fluxes. Contrary to the interval 250 • to 70 • spanning more degraded shrubs fields, therefore, the LAS measures differed strongly to EC measures.

Remote Sensing Data
In the objective of developing a method to characterize convective fluxes at the landscape scale from remote sensing observations only, the main inputs of the TSEB model (see below) were derived from MODIS and MSG SEVIRI products.

MODIS Products
The MODIS instrument operates on both Terra and Aqua spacecraft. It has a viewing swath width of 2330 km and views the entire surface of the Earth every 1-2 days. Its detectors measure 36 spectral bands between 0.405 and 14.385 µm, and it acquires data at three spatial resolutions: 250 m, 500 m, and 1000 m. The collection version 5 of the products was used in the present study. The MOD15A2 LAI is a 1 km global data product updated once every 8 days derived from the MODIS sensor on board TERRA. The MODIS LAI algorithms were developed jointly by personnel at Boston University, the University of Montana SCF and NASA GSFC. The algorithm consists of a main procedure based on the inversion of a 3D radiative transfer model thanks to a look-up-table. This algorithm exploits the spectral information content of MODIS surface reflectances at up to seven spectral bands. Should this main algorithm fail, a back-up algorithm is triggered to estimate LAI empirical relationships with vegetation indices. The LAI products in the collection version 5 are available from 2001 to present.
The MCD43B3 α product is upscaled from MCD43A3 and provides both the white-sky-albedo (α) and the black-sky-α at 1-km resolution every 16 days. Both Terra and Aqua data are used in the generation of this product. We were interested in a single integrated value over the entire solar emission spectrum (0.3-5.0 µm) called shortwave broadband α. The black-sky α (directional hemispherical reflectance) is a directional α corresponding to 100% of direct light; it depends on the solar zenithal angle unlike the white sky α (bihemispherical reflectance) which corresponds to a completely diffuse light. The "blue sky" α used in this study is a weighted average between these two extreme values. The weighting depends essentially on the aerosol content of the atmosphere [74]. In our case, we considered that we have 85% direct light and 15% diffuse, and the final α value is calculated as follows [74]: The MOD11A1 and MYD11A1 at 1 km spatial resolution under clear-sky conditions are the daily land surface temperature (Ts) products, derived from Terra and Aqua, respectively. Ts were derived from the brightness temperature using bands 31 and 32 through a generalized split-window algorithm. The daily level 3 MODIS Ts (collection 5) were used in this study. Table 2 ummarizes the data layer characteristics of MODIS products.
All MODIS products were downloaded from the website (https://search.earthdata.nasa.gov/search). The Spinning Enhanced Visible and InfraRed Imager (SEVIRI) aboard Meteosat Second Generation (MSG) is a multi-spectral sensor, imaging across the visible and near-IR, which provides data with a high temporal and spectral resolution every 15 min in 12 wavebands, and image spatial sampling scales of 3 km for nadir view. The land surface temperature is estimated from top of atmosphere (TOA) brightness temperatures of SEVIRI split-window channels, centered on 10.8 and 12.0 µm, using a generalized split-window (GSW) algorithm [75] with the adoption of SEVIRI data [76]. IR radiance is absorbed and scattered even by thin clouds and aerosols. Therefore, the retrieval of Ts works only for completely cloud-free pixels, and the quality of this product is automatically assessed by means of the accuracy of the parameters used in the Ts algorithm. Ts products used in this study are distributed by EUMETSAT through the Satellite Application Facility on Land Surface Analysis (LSA SAF).

TSEB Model Description and Implementation
The TSEB model was presented and described by Norman et al. [7]. It solves two separate energy balances for the soil and vegetation and estimates evaporation and transpiration as residual terms of the energy balance. The TSEB model adopts the Priestley-Taylor (P-T) parameterization and an iteration procedure estimate the energy partitioning.
In TSEB, the directional radiometric temperature (T rad (θ)) is divided into its soil and vegetation fractions as seen by the radiometer, and is expressed as follows: where T veg and T soil are the vegetation and soil components of temperature (K). The fraction of the field of view of the infrared radiometer occupied by canopy, f(θ), can be calculated by combining the view zenith angle θ and the vegetation cover fraction f c considering vegetation with a spherical distribution of leaf angles [77]: f c is simply f(θ = 0), or The LAI estimates were then used to derive the canopy height by the empirical equation: In this study, the soil emissivity and leaf emissivity were taken as constant values from the literature (0.95 and 0.97, respectively).
The double source energy balance corresponds to the simple source balance shared between the soil and vegetation components: R n,veg = H veg + LE veg (15) with R n,soil being the soil net radiation, H soil the soil sensible heat flux, LE soil the soil latent heat flux, R n,veg the vegetation net radiation, H veg the vegetation sensible heat flux and LE veg the vegetation latent heat flux. The surface sensible heat flux (H) is partitioned between the vegetated canopy and soil using the following relations: H soil = ρc p T soil − T a r ah + r s (17) where ρ (kg m −3 ) is the density of air; c p (J kg −1 K −1 ) is the specific heat capacity of air; T a is the air temperature; T veg and T soil are, respectively, the canopy and soil temperatures; r ah is the aerodynamic resistance to heat transfer across the canopy-surface layer interface; and r s is the resistance to heat flux in the boundary layer immediately above the soil surface. The r ah is calculated from the adiabatically corrected log temperature profile equation [78] expressed as: where z u and z T are the height of wind speed and temperature measurement U a and air temperature measurement T a , respectively; d is the displacement of reference plan (d = 2/3 × h c ; h c , is canopy height); z m is the roughness length for momentum (z m = h c /8); V kar is the von Karman constant taken to be 0.4; and Ψ m and Ψ H are the adiabatic correction factors for momentum and heat, respectively [78]. Although r s is rather complex, as it depends on many factors, Sauer et al. [79] provided a reasonable approximation: where a = 0.004 m s −1 , b = 0.012, and U s is the wind speed at a height above the ground where the effect of soil roughness is minimal (typically between 0.05 and 0.2 m). Using the empirical relations from Campbell and Norman [80]: with a sc being the leaf size expressed as: and U h the wind speed at the top of vegetation canopy given by: The surface soil heat flux is estimated as a fraction of R n,soil : where c g ∼ 0.35 [81]. The latent heat flux from the vegetated canopy is derived from the PT formula: where γ is the psychometric constant (≈67 Pa K −1 ), f g is the fraction of leaf area index (LAI) that is green, ∆ is the slope of the saturation vapor pressure versus temperature curve, and α PT ∼ 1.26 for TSEB model [82]. To solve the system that includes more unknowns than equations, an iterative process is used that closes on the Obukhov length (L). The algorithm starts by calculating the wind and turbulence variables (L, U * , U s , and r ah , expressed by the above equations).
The available energy and heat fluxes are then determined by calculating the net radiation (R n ) and its partition between the vegetation (R n,veg ) and the soil (R n,soil ), as well as the conduction flux in the soil (G). The plant transpiration (LE veg ) is estimated by the PT equation, which allows the residual calculation of the vegetation sensitive heat flux (H veg ) as well as the vegetation temperature (T veg ). The soil temperature (T sol ) is estimated by radiometric temperature and canopy fraction, and is then used to estimate the soil heat flux (H sol ). An energy balance on the soil is used to calculate soil evaporation (LE sol ).
The estimated components of latent heat flux (LE veg and LE soil ) are assumed to be positive, meaning that there is no condensation. If LE soil is negative, then its value is set to zero and a new value of LE veg is calculated as a residual term of the vegetation energy balance. If LE veg is negative, then transpiration is reduced iteratively by decreasing the value of the PT coefficient until the LE veg value becomes positive or zero.
The heat fluxes then allow the recalculation of L. The iterative process is repeated until numerical stability of L is reached.

Aggregation Scheme
Spatial aggregation is a method to bridge the model parameters controlling surface exchange at a patch scale with the areal average value of equivalent model parameters applicable at a larger scale to estimate grid scale surface fluxes using the same equations that govern the patch scale behavior [66]. In general, aggregation can be applied either to the input forcing of the evapotranspiration models, or to the fluxes derived from fine resolution input fields [33]. Here, the first assumption is used to link between local and effective values of land surface parameters [27,66,83,84]. The chosen averaging approach is directly related to the considered variable. For albedo (α) and displacement height d, a simple arithmetic average is used: From the Stefan-Boltzmann law, the aggregated value of the surface temperature T s can be obtained from: where f i is the fraction of the surface covered by the patch i (with ∑ i f i = 1). Finally, the aggregation of the roughness length is obtained from a logarithmic average as follows: Because of the absence of a soil occupation map for the study site, the aggregation scheme using in-situ measurements was done by weighted average of the estimated variables from the data of the 3 EC stations. Following Ezzahar et al. [27], the local values obtained by each EC system are assumed to be representative for the sites where those systems are installed along the basin. Sites millet, fallow, and degraded shrubs represented about 54%, 26%, and 20% of the Wankama catchment, respectively.
Three different aggregation methods were implemented to calculate the inputs of TSEB model from the MODIS products. They differ in the way the pixel scale inputs are used to derive the grid scale ones. The first scheme consists in averaging the MODIS products at their nominative resolution scale on geographic windows of 10 × 8 pixels around the scintillometer transect (i.e., independently of the scintillometer footprint; cf. Figure 3). The second scheme averages the MODIS products over the scintillometer transect without considering the representivity of the corresponding pixels with regards to the footprint (i.e., each pixel belonging to the footprint have the same contribution to the aggregated value). Finally, a weighted average aggregation scheme was implemented. It consists in weighting the MODIS pixel according to their proportion occupied in the footprint. The aggregated values of MODIS products are then used as inputs of the TSEB model at Terra-MODIS and Aqua-MODIS overpass time.

Results and Discussion
As a preliminary step, the quality of experimental data including EC and scintillometer surface fluxes and satellite products were assessed. In a second step, the TSEB predictions were evaluated at the station scale. Finally, the TSEB predictions of the convective fluxes were evaluated at the scale of the scintillometer footprint: (1) by aggregating the station scale predictions using the in-situ observations of albedo, Ts and LAI, which constitutes the ideal case; (2) by using 3 km Ts observations from the MSG SEVIRI instrument directly at the scale of the scintillometer measurements; and (3) by testing three aggregation schemes of MODIS products to assess if a better representation of heterogeneity based on the 1-km products improves our large-scale estimates.

In-Situ Surface Fluxes
The quality of EC measurements is usually assessed by assuming that the energy balance is closed. The comparison of the (R n -G) and the sum of the latent and sensible heat fluxes (LE EC + H EC ), measured independently by the EC systems (data not shown here), showed an underestimation of the turbulent fluxes by about 8%, 17%, and 20% for the millet, fallow and degraded shrubs, respectively. Several practical reasons preclude perfect closure [85], e.g., spatial scale discrepancies between measurements of energy balance components, and spatial heterogeneities; storage between measurement levels (canopy and superficial soil); so-called eddy flux losses (not measured here); and difficulties with ground flux estimation [86,87]. Results reported in other experimental studies have shown that balance shortfalls commonly stand between 10% and 40% of the available energy [85], [88]. As a conclusion, closure is quite satisfactory for our experiment, especially at the millet site which was more homogenous than fallow and degraded shrubs sites. Figure 4 plots the time series of the average between 9:00 a.m. and 5:00 p.m. of H LAS , and LE LAS during the entire study period. Note that the scintillometer only calculates H, and LE is derived as a residual term of the energy balance (LE LAS = R n − G − H LAS ). In this case, the available energy has been aggregated using the values measured at each instrumented field ( R n − G agg = Stages 1, 2 and 3 present, respectively, the beginning of season (1 June-mid-July), the growing stage which took place over two months (mid-July-mid-September) and the senescence stage (from mid-September). Preliminary to the scintillometer data analysis, a comparison between H LAS against weighted EC fluxes was done for the wind direction interval 70 • to 250 • , because the footprint of the LAS covers the monitored sites where the EC systems were installed. The statistical results show a good agreement between the LAS sensible heat fluxes and those derived from the EC systems with a relative error of about 20% [27]. The H LAS and LE LAS dynamics are almost stable throughout the season and LE LAS largely dominates H LAS through the growing stage. Rainfall is rather regular between mid-July and mid-September which ensured enough soil moisture to keep millet growing until mid-September. It can be seen that every rainfall event caused an immediate response with a significant energy shift between H and LE. LE LAS is around 200 W/m 2 and increases after each rain event, and peaks at the end of August with 315 W/m 2 , which corresponds to the LAI peak for Millet crop [24]. In the same period, soil moisture was high and reached its maximum. Some points with very low LE LAS values were recorded during the growing stage which can be explained by the very dry conditions or by a considerable amount of bare soil. A sharp increase of LE LAS was shown after the 40-mm rain event of 20 September. However, LE LAS decreases systematically from the end of September to the end of October, while H LAS reaches its maximum (265 W/m 2 ).

Remote Sensing Products
First, an evaluation of the agreement between MODIS albedo and Ts products and in-situ measurements was required because of their great influence on the predicted fluxes. This comparison was carried out only for clear sky days when MODIS images are available. Figure 5 displays the scatter plots of MODIS derived Ts and albedo against in-situ measurements, for millet and fallow sites, as only one month of measurements were available for the degraded shrub site. The agreement between MODIS and measured albedo is relatively poor (RMSE = 0.05, R = 0.45). This may be attributed to the lack of representivity of the sites covered by relatively dense vegetation with regards to the north from the LAS path mainly composed of bare soil. This lack of representivity is particularly prominent in the Sahel as soils are very clear, especially at the beginning (July) of the season. The dry spell produced a short period of albedo increase. During the rain events, albedo follows a very markedly decreasing general trend, with large fluctuations due to alternating dry and wet spells [26]. This assumption is supported by the observed positive bias of 0.03. The Ts comparison shows a reasonable correlation of R = 0.77 with an RMSE = 3.8 • C, and an underestimation of Ts MODIS by −1.8 • C on average.
The RMSE is on the upper limit of the values reported in the literature that usually range between 2 and 4 K [89]. The observed scattering may also stem from the fact that the in-situ Ts is not representative for the satellite Ts. A finer look at the results highlights also a very contrasted bias between the two sites: over Millet that dominates the Wankama land use with 58% of the area, the bias is logically low (MBE = −0.5 K) while it reaches more than −4 • C over the denser and thus colder vegetation of the fallow site. Finally, Figure 6 shows an intercomparison of MSG SEVIRI and MODIS Ts averaged over a 3 × 3 pixel windows to match the MSG resolution at time of Terra and Aqua MODIS overpasses. The overestimation of MSG SEVIRI Ts with regards to MODIS Ts is prominent along the whole Ts range. This could be linked to the view zenith angle differences between MSG SEVIRI and MODIS satellites. The impact of changes in viewing angle is enhanced with surface heterogeneities. In addition, Ts was retrieved using one MSG SEVIRI window channel (10.8 µm), while the MODIS algorithm is based on two windows channels (11.03 µm and 12.02 µm). These dissimilarities related to the spectral characteristics of the two instruments can induce the observed differences in retrieved Ts. Finally, the positioning error of MSG has been found to reach up to 8 pixels, i.e.,~30 km [90], which could also explain this overestimation, in particular if the pixel is shifted to the north of the area characterized by bare and hot soils.

Multi-Scale Convective and Evaporative Fluxes Predictions
As a preliminary step, the TSEB predictions are evaluated at the station scale both using in-situ and MODIS derived products. Focus is then put on the evaluation of TSEB at the grid scale using in-situ, MODIS and MSG products.

Station Scale
Using In Situ Data The available energy and turbulent fluxes predicted by the TSEB model are compared to those measured over the three eddy covariance stations at half-hourly time step. Table 3 summarizes the statistical results including the number of observations (n), the correlation coefficient (R), the root mean square error (RMSE), and the mean bias error (MBE). The TSEB model estimates correctly the available energy (R n -G) with RMSE of 52, 49, and 57 W/m 2 for millet, fallow and degraded shrubs, respectively. A first remarkable feature is the slight underestimation (−23 W/m 2 ) for the millet site. The second remarkable feature is a stronger dispersion for the degraded shrubs area. This dispersion is mainly related to the estimation of the soil heat flux, with RMSE = 38 W/m 2 and Mean Bias Error MBE = 23 W/m 2 for the millet site and RMSE = 60 W/m 2 and MBE = 54 W/m 2 for degraded shrubs site. Indeed, affected by the complexity of the physical processes that occur in the soil as well as the vegetation cover, soil heat flux is the most difficult parameter to estimate with great precision. Several studies have pointed out this difficulty in estimating G especially on sparse vegetation [27,30,68]. In addition, measuring this parameter at 5 cm and in sparse vegetation represents a major challenge to get representative measurement. The installation of heat flux plates in the soil must consider several criteria: the plates must be totally covered to ensure that it is not directly exposed to the sunlight, especially when the soil is very sandy as in the Wankama catchment where sand represents 80-90% [91]. Additionally, heavy rainfall within a few hours, which is relatively common in the Sahel region, can uncover the plates and thus lead to exposing them to sunlight. Moreover, the use of the Brutsaert formula (developed for clear sky days by Brutsaert [92]) for the incoming longwave radiation (that replaced the measurements when they were not available for this station, to estimate R n ) can create a significant dispersion for low radiation values in cloudy conditions [68].
The agreement between the measured and predicted turbulent fluxes is obviously lower than for the estimated available energy but, despite a significant dispersion, correlation coefficients and RMSEs are encouraging ( Table 3). The observed dispersion can be related to the footprint effect, as the footprint of the EC system is considerably larger than the representivity scale of the measured input variables (mainly α, Ts and LAI). In addition, other parameters such as roughness length (z 0m ) and displacement height (d) were estimated as a fraction of the vegetation height using classical rules of thumb [93,94]. These equations have been established for homogeneous and dense covers while the study sites are heterogeneous and sparse. A further step to minimize dispersions would be to estimate these variables using EC data following Hoedjes et al. [95] and Lagouarde et al. [96]. Finally, several empirical and constant model parameters have been used in the TSEB model, which require calibration according to the specific conditions of the study area. Among them, the Priestley-Taylor parameter (α PT ) directly relates latent heat flux to the available energy at the surface. Most studies driven by TSEB model have used its theoretical value of 1.26 [7,42,97,98]. Other studies have identified that it is variable under different surface and atmospheric conditions (varies in the range 0.5-2). In particular, the α PT was found to be smaller for dry surfaces and higher for humid conditions [99]. In the same context, Ait Hssaine et al. [44] demonstrated, using a new TSEB-based evapotranspiration model (called TSEB-SM), that α PT cannot be considered as a constant because it varies in time according to several factors including LAI, green vegetation cover fraction and soil water availability. In particular, it can take extreme values under dry, water deficit and advective conditions. Interestingly, the partition between latent and sensible heat fluxes is correctly reproduced for the degraded shrub site (mean biases lower than 0.92 W/m 2 and 11 W/m 2 for respectively, H and LE). Indeed, the results of sensible heat flux over the millet crop matches perfectly with the results found by Lhomme et al. [83] (RMSE = 43 W/m 2 ) who used a two layer model to estimate H from Ts. Lhomme et al. [83] used an empirical relationship between Ts and Ta (δT = a(Ts − Ta) m ), m and a being statistically determined by adjusting the modeled to the H observed by the Bowen ratio method. Degraded-shrubs Using MODIS Data An intercomparison between (R n -G) and turbulent fluxes measured by each station and the fluxes estimated by TSEB model using MODIS data was done and the error statistics are reported in Table 3. Overall, the discrepancies between the available energy estimated from TSEB model and measured one, for both sites (millet and fallow) are likely due to greater scatter between MODIS and measured α, and to the difference between MODIS and in-situ Ts. In all cases, LE does largely dominate H through the growing season (July-September), in particular for the fallow field. For the millet site, the late start of the rainy season had significant adverse effects and the crop was not able to take advantage of the abundant rainfall because of insufficient plant development [45]. For both sites, an overestimation of simulated LE compared to measured ones is recorded at the end of the season. Indeed, the saturation of TSEB in the higher range of LE is due to the fixed maximum value for α PT (equal to 1.26) [44]. The structure of the model cannot accommodate large evaporative demand conditions and strong advective conditions [100]. Moreover, the change in wind direction in monsoon periods strengthens advection effects. As a conclusion, TSEB prediction, although slightly biased at the end of the season, reproduces quite well the observations: RMSE/MBE are relatively low for both millet and fallow sites.

Three-Kilometer Grid Scale
Using In-Situ Data Figure 7 displays the comparison between three sites-averaged turbulent fluxes values derived from TSEB model using weighted input data (Ts, α, d, z 0 , LAI, and h c ) compared against the turbulent fluxes measured by the scintillometer. The corresponding statistics are shown in Table 3. TSEB provides satisfying results for the sensible heat flux with a relative error of about 24% (R = 0.87, RMSE = 37 W/m 2 , MBE = −21 W/m 2 ). Indeed, the scintillometer path is very heterogeneous. This can generate a difference in flux even between two adjacent plots. Therefore, the dispersion observed can be related to the differences in the footprint of LAS, as well as to the uncertainties of the similarity stability functions. The surface temperature is a key parameter for turbulent fluxes and for monitoring the energy balance. Weighted Ts over heterogeneous surface using locally measured values can lead to errors in the partition of temperature between soil and vegetation, and consequently to estimate sensible heat flux. A visual assessment of scatter plots in Figure 7b and the statistics ( Table 3) clearly indicates that TSEB overestimates LE fluxes with MBE = 39 W/m 2 , especially at the end of the season. LE measured by the LAS never exceeds 300 W/m 2 while simulations reach 400 W/m 2 , which could be related to the fixed value of α PT to 1.26. Note also that any error in the available energy induces the LE LAS errors.

Using 3-km MSG SEVIRI Ts
Here, MSG SEVIRI Ts is used to simulate H and LE by TSEB directly at the scale of the scintillometer path (3 km) spanning the three vegetation types (only the interval 70 • to 250 • is considered). Those data were used in conjunction with an average of a 3 × 3 pixel MODIS area for α and LAI. TSEB estimates were compared to H LAS and LE LAS . According to precipitated amounts during the growing season, the LE is noticeably higher at the fallow site compared to millet site (patch scale). Note that the MSG SEVIRI pixel covers diverse land cover (millet, fallow, degraded shrubs and bare soil), which explains the discrepancies obtained between MSG SEVIRI-based and measured LE (Figure 8b). By contrast, the overestimation of H at the end of the season is mainly linked to the higher surface temperatures measured by MSG SEVIRI. It is reminded that the scintillometer fetch does not exceed 100 m, while the MSG SEVIRI resolution at the study site is about 3 km.

Using MODIS Products
In this section, the effect of land cover heterogeneity on the estimated turbulent fluxes has been studied by combining TSEB model driven by Terra and Aqua remotely sensed data and an aggregation scheme, at Terra-MODIS and Aqua-MODIS overpass time. Figure 9 displays results for three different aggregation schemes (the simple averaging at the MODIS resolution scale (Scheme 1 (Figure 9a,b)), the simple averaging over the footprint extent (Scheme 2 (Figure 9c,d)), and a weighted average of the footprint extent (Scheme 3 (Figure 9e,f), see Section 2.5). Statistical metrics are reported in Table 3 The simple averaging method using MODIS resolution data in Figure 9a underestimates H with a relative error of about 21% (RMSE = 73 W/m 2 , R = 0.71, MBE = −48 W/m 2 ). The correspondence between simulated and measured fluxes for Scheme 2 is good compared to Scheme 1 (Figure 9c). The relative error was about 13%, the RMSE value was 65 W/m 2 , and the linear regression forced through the origin yielded an R of 0.70 and the mean bias error was reduced to −30 W/m 2 instead of −48 W/m 2 for the first scheme. H simulated by using Scheme 3 in Figure 9e is closer to the 1:1 line, providing a quite significant improvement with regards to the second scheme. The RMSE and the mean bias error between the simulated and measured values were 63 and −23 W/m 2 , respectively. The relative error of 10% indicates that Scheme 3 is very reliable in providing area-averaged sensible heat flux over heterogeneous surfaces.
In a subsequent step, simulated LE was compared to the observed one using the three schemes described before. As shown in Figure 9b, the Scheme 1 overestimates LE with a relative error of about 58% (RMSE = 102 W/m 2 , R = 0.85, MBE = 73 W/m 2 ). That can be explained by the fact that LE is a residual term affected by errors in both available energy and sensible heat flux. Figure 9d,f plots simulated versus observed LE for Schemes 2 and 3. TSEB provides satisfying results with a RMSE of 65 W/m 2 and 63 W/m 2 , and relative error of 39% and 35% for the second and third schemes, respectively. Based on the above results, Scheme 1 could not achieve satisfactorily the upscaling of turbulent fluxes over heterogeneous land surfaces. In fact, the method considers that all pixels belonging to the footprint of LAS contributes 100% of the pixel value, even for pixels of low contribution. This will obviously lead to large differences in Ts, which is a key feature for determining the turbulent fluxes.
To go further into the comparison, the Figure 10 shows the time series of the convective fluxes observed by the LAS (H LAS and LE LAS ) at the time of the satellite overpass and simulated by TSEB based on MODIS products for the footprint weighted method (Hsim and LEsim). According to the availability of MODIS images, and considering the wind direction interval 70 • to 250 • , only LAS data from mid-August to the end of the experiment on 22 October were available. Despite a significant scattering, both observed and simulated fluxes depict a similar behavior with a gradual switching of available energy from latent to sensible heat following the last rainfall around mid-September. An interesting feature is the underestimation of H by TSEB during the senescence while LE is overestimated for the same period. This could be attributed to the constant value of α PT equal to 1.26 used in this study. Indeed, this value is usually accepted for semi-arid to sub-humid agricultural areas under irrigation regime [7,43,101] while it has been shown that lower α PT are obtained for dry areas such as in our region of study after the end of the rainy season. Finally, the obtained significant scattering, in particular after the last rainfall could also be related to the radiative properties of vegetation that are not well represented for dry vegetation.

Summary and Concluding Remarks
The main objective of this work was to estimate the turbulent fluxes over a path covering the three dominant crops (millet, fallow and degraded shrubs) in the Wankama basin by combining TSEB model, an aggregation method and satellite data (MODIS and MSG SEVIRI). The data used to validate this approach were collected within the framework of the AMMA project. Each site was equipped with meteorological stations and Eddy Covariance systems for monitoring the energy and water balance as well as the soil and vegetation characteristics. In addition, a large aperture scintillometer was installed over a transect of about 3.2 km spanning the three dominant vegetation types (with contrasted energy balances and different fluxes behaviors) to derive the integrated H and LE values over this heterogeneous landscape. Our approach was organized into three steps. The experimental data including in-situ surface fluxes and remote sensing products were evaluated. The TSEB model was then evaluated at the station scale for the three study sites (millet, fallow and degraded shrubs) instrumented with EC systems. Finally, the TSEB prediction of surface fluxes were evaluated at the grid-scale based on in-situ data (at the 30 min time), MODIS (at Terra and Aqua overpass time), and MSG products (at MDG-SEVIRI overpass time) and aggregation schemes.
The results obtained at the station scale for H and LE are relevant, especially when using in-situ data given the complexity of the study sites and in the range of what has been reported in the literature in terms of RMSE and correlation coefficient. Some scatter was noted, which is mainly related to the difference in footprints of the measuring instruments (about 100 m for EC systems compared to less than 10 m for inputs variables such as LAI, albedo and Ts), as well as to the fixed maximum value for α PT (equal to 1.26) during the entire season. Despite this dispersion and as already underlined by several authors, TSEB provides with good performances for H and LE prediction at the station scale by using in-situ measured inputs or, to a lesser extent, MODIS derived inputs.
TSEB was then evaluated at the (1 km and 3 km) grid-scale along the scintillometer transect spanning a heterogeneous landscape. To this end, in-situ data were first aggregated to estimate the mean turbulent fluxes. In this ideal case, where the major inputs are well known thanks to in-situ measurements on the dominant vegetation type, the predicted fluxes obviously showed a good agreement with the scintillometer observations (RMSE = 37 and 75 W/m 2 for H and LE, respectively) although some dispersions were observed, which were explained by the difference in the footprints of the measuring instruments as well as because LE LAS was calculated as a residual term of the energy balance using the average available energy values derived from measurements at each site. In a second step, TSEB was validated at the grid-scale but with remotely-sensed inputs. To this objective, MSG SEVIRI and MODIS sensors only offer a revisit time able to sample the high inter-and intra-daily variability of the SEB in semi-arid areas. The first assessment was performed using the MSG SEVIRI Ts at 3 km without explicit representation of the sub-pixel heterogeneity. The agreement between H and LE simulated by TSEB model and scintillometer observations is very poor, in particular in terms of temporal dynamic (R = 0.39 and 0.2 for H and LE, respectively). This dispersion could be attributed to the large bias on MSG SEVIRI Ts linked to the heterogeneity of land cover. In addition, this highlights the need to represent the sub-pixel heterogeneity. This was tested in a final step using 1-km MODIS products as input of the TSEB model and three aggregation schemes. While the results with the simplest aggregation scheme consisting in averaging MODIS inputs independently of the footprint provides poor results (RMSE of 73 W/m 2 and 102W/m 2 for H and LE, respectively), the statistical metrics based on the two other aggregation schemes are encouraging. In particular, Scheme 2, considering the MODIS pixel included in the scintillometer footprint only, and without weighting the values based on the footprint contribution as in Scheme 3, represents a good trade-off between accuracy and ease of implementation to map evapotranspiration over complex landscape.
This work offers perspectives for the spatial mapping of evapotranspiration over poorly gauged areas such as semi-arid areas of the Sahel.