An evapotranspiration model self-calibrated from remotely sensed surface soil moisture, land surface temperature and vegetation cover fraction: application to disaggregated SMOS and MODIS data

Thermal-based two-source energy balance modeling is essential to estimate the land evapotranspiration (ET) in a wide range of spatial and temporal scales. However, the use of thermal-derived land surface temperature (LST) is not sufficient to simultaneously constrain both soil and vegetation flux components. Therefore, assumptions (about either soil or vegetation fluxes) are commonly required. To avoid such assumptions, an energy balance model, TSEB-SM, was recently developed by Ait Hssaine et al. (2018b) in order to consider the microwave-derived near-surface soil moisture (SM), in addition to the thermal-derived LST and vegetation cover fraction (fc) normally used. While TSEB-SM has been successfully tested using in situ measurements, this paper represents its first evaluation in real life using 1 km resolution satellite data, comprised of MODIS (MODerate resolution Imaging Spectroradiometer) for LST and fc data and 1 km resolution SM data disaggregated from SMOS (Soil Moisture and Ocean Salinity) observations. The approach is applied during a 4-year period (2014–2018) over a rainfed wheat field in the Tensift basin, central Morocco. The field used was seeded for the 2014–2015 (S1), 2016–2017 (S2) and 2017–2018 (S3) agricultural seasons, while it remained unploughed (as bare soil) during the 2015–2016 (B1) agricultural season. The classical TSEB model, which is driven only by LST and fc data, significantly overestimates latent heat fluxes (LE) and underestimates sensible heat fluxes (H ) for the four seasons. The overall mean bias values are 119, 94, 128 and 181 W m−2 for LE and −104, −71, −128 and −181 W m−2 for H , for S1, S2, S3 and B1, respectively. Meanwhile, when using TSEB-SM (SM and LST combined data), these errors are significantly reduced, resulting in mean bias values estimated as 39, 4, 7 and 62 W m−2 for LE and −10, 24, 7, and −59 W m−2 for H , for S1, S2, S3 and B1, respectively. Consequently, this finding confirms again the robustness of the TSEB-SM in estimating latent/sensible heat fluxes at a large scale by using readily available satellite data. In addition, the TSEB-SM approach has the original feature to allow for calibration of its main parameters (soil resistance and Priestley–Taylor coefficient) from satellite data uniquely, without relying either on in situ measurements or on a priori parameter values.

Regarding the data availability over extended areas, remote sensing is the only viable technique that can provide representative and multi-resolution measurements of ET. As a consequence, the spatial modeling has become a dominant means of estimating ET fluxes over regional and continental areas (Anderson et al., 2007;Fisher et al., 2017). In this context, numerous models based on land surface temperature (LST) data have been developed, such as (i) residual balance methods that consider ET to be the residual term of the energy balance, like TSEB (two-source energy balance,  and SEBS (surface energy balance system, Su, 2002), (ii) contextual methods that estimate ET as the potential ET times the evaporative efficiency (Moran et al., 1994) or as the available energy times the evaporative fraction (Merlin et al., 2013;Roerink et al., 2000) and (iii) other categories of models that integrate LST into a water balance model  or into the Penman-Monteith energy balance (PMEB) equation to directly estimate ET Mallick et al., 2015Mallick et al., , 2018.
Among well-known temperature-driven energy flux models, the TSEB model proposed by  has been shown to be robust for a wide range of landscapes (Colaizzi et al., 2012;Ait Hssaine et al., 2018a). TSEB has two key input variables, which can be derived from remote sensing data. The first one is the LST and the second is the vegetation cover fraction (f c ). The TSEB model adopts an iterative procedure, in which an initial estimate of the plant transpiration is given by the Priestly-Taylor (PT) formulation (Priestley and Taylor, 1972). This assumption requires few input data and allows a precise estimate of potential ET (Fisher et al., 2008). Nevertheless, several studies (Ait Hssaine et al., 2018b;Fisher et al., 2008;Jin et al., 2011;Yang et al., 2015) have stressed that the PT coefficient cannot be considered a constant value, as it is influenced by several parameters. Other authors (Gonzalez-dugo et al., 2009;Long and Singh, 2012;Morillas et al., 2014) reported that the PT approach may overestimate the canopy ET, especially for low soil wetness, and/or sparse vegetation cover, because it does not include a reasonable reduction of the initial canopy ET under stress conditions. Recently, Boulet et al. (2015) developed the Soil-Plant-Atmosphere and Remote Sensing Evapotranspiration (SPARSE) model, which is similar to the TSEB model in its basic assumption but with additional constraints to improve the ET model performance in heterogeneous vegetation. The former first generates an equilibrium LST from the evaporation efficiency and the transpiration efficiency estimates by assuming that their values are equal to 1. Then, LST is implemented in the SPARSE retrieval mode to circumscribe the output fluxes by both limiting cases (namely the fully stressed and potential conditions). In spite of the good retrieval performances of ET by this model, significant uncertainties are observed during the quasi-senescent vegetation period .
Alternatively to the use of LST as a proxy for ET, numerous studies have stressed that the soil moisture plays a critical role in the partitioning of available energy into latent and sensible heat fluxes and is the prominent controlling factor of actual ET Gokmen et al., 2012;Kustas et al., 1998Li et al., 2006). Several authors have revised the well-known LST-based TSEB model and replaced the LST with microwave-derived surface soil moisture (SM) to estimate daily ET (Bindlish et al., 2001;Kustas et al., 1998Li et al., 2006). Bindlish et al. (2001) found that the impact of SM on surface fluxes is strongly related to the vegetation cover. The impact is high for low fraction cover and relatively weak for the high cover fraction. Moreover, the soil evaporation is constrained by the SM through its soil-texture-dependent coefficients (a rss and b rss ) reported in Sellers et al. (1992). In the same way, Li et al. (2006) indicated that the model performance is sensitive to these two coefficients, and thus they proposed averaging the output of LST-based TSEB and SM-based TSEB models in order to provide more consistent results over a wide range of conditions.
Previous studies, either LST-or SM-based, agree with the view that combining both LST and SM information at a time would enhance the robustness and accuracy of ET estimates in various biomes and climates. Nevertheless, few studies have simultaneously combined both observations in a unique energy balance model. One difficulty lies in developing a consistent representation of the soil evaporation (as constrained by SM, Chanzy and Bruckler, 1993), the total ET (as constrained by LST, Norman et al., 1995) and the plant transpiration (as indirectly constrained by both LST and SM, Ait Hssaine et al., 2018b). Gokmen et al. (2012) explicitly integrated the SM derived from AMSR-E (Advanced Microwave Scanning Radiometer for EOS) data (Owe et al., 2008) into the LST-derived SEBS model via the kB −1 parameter, which plays an important role in the aerodynamic resistance. The updated SEBS model (SEBS-SM) provided a large improvement of sensible heat flux and thus ET estimates under water-limited conditions. The point is that the soil evaporation reduction parameters were calibrated using in situ measurements, which limits the validity of the approach over large areas. In the same vein, Gan and Gao (2015) incorporated a SM-based soil resistance term in the TSEB formalism and calibrated several parameters (including the PT coefficient) using LST data. The obtained results showed that the model calibrated by LST data performed better than the non-calibrated one. Note that the parameters of the soil resistance were set to constant values as in Sellers et al. (1992) and Li et al. (2006). As a further step towards the combination of LST and SM data, Ait Hssaine et al. (2018b) modified the TSEB formalism (Kustas and Norman, 1999; (named TSEB-SM) and proposed a new calibration strategy of the main PT-based TSEB-SM parameters. The TSEB-SM model was tested using in situ measurements and provided an important improvement in terms of latent heat flux/sensible heat flux estimates compared to the classic TSEB all along the agricultural season, especially during the crop emergence and the senescence periods. Such improvements are attributed to stronger constraints exerted on the representation of soil evaporation (via SM data and the calibrated soil parameters) and plant transpiration (via the calibrated daily PT coefficient). It should be noted that only LST and SM data are used for the calibration of yearly a rss and b rss as well as daily α PT , while the flux measurements are needed only for the validation of the TSEB-SM-simulated sensible and latent heat fluxes.
One crucial point is that all the above studies based on remotely sensed SM and LST have neglected the mismatch in the spatial resolutions of readily available SM products. Especially the global-scale SM data sets have a typical resolution of 40-50 km (Entekhabi et al., 2010;Kerr et al., 2010;Njoku et al., 2003). Such spatial resolution is generally unsuitable or even incompatible with many hydrological and agricultural applications. To fill the gap, disaggregation approaches of AMSR-E, SMOS (Soil Moisture and Ocean Salinity) and SMAP (Soil Moisture Active and Passive) like SM data have been developed (Peng et al., 2017) but, to date, there has been no application of SM-based ET models to disaggregate SM data sets. In addition, the use of remote sensing data would be necessary in order to avoid the timeconsuming process of calibrating the TSEB model over each field.
Although TSEB-SM has the capability to calibrate its main parameters from remotely sensed data, the real-life application needs extensive evaluation and testing. The objective of this paper is thus to demonstrate for the first time this capacity using disaggregated SMOS and MODIS (MODerate resolution Imaging Spectroradiometer) data. For this purpose, TSEB-SM is applied to 1 km resolution using MODIS LST/f c data and to SMOS SM data. To make the SMOS data spatially consistent with MODIS data, the SMOS SM is disaggregated at 1 km resolution using the DisPATCh (DISaggregation based on Physical And Theoretical scale Change) algorithm Merlin et al., 2013;Molero et al., 2016). The proposed methodology is evaluated over a rainfed wheat field in the Tensift basin, central Morocco during four agricultural seasons (2014-2018).
2 Data description and methods

Site and in situ data description
The study site is situated in the east (Sidi Rahal) of the Tensift basin in central Morocco (see Fig. 1). The region is characterized by a semi-arid Mediterranean climate, with an average yearly precipitation of about 250 mm and an atmospheric  (Allen et al., 1998;Jarlan et al., 2015). Soil is characterized by a fine texture with 47 % of clay, 33 % of loam and 18.5 % of sand (Er-Raki et al., 2007). The experiment has been set up in a rainfed wheat ("Bour") field since 2013 (Ali Eweys et al., 2017;Amazirh et al., 2018;Merlin et al., 2018). Located within a larger area occupied by rainfed wheat "Bour", this field was chosen to be representative at a scale of 1 km, thus enabling the comparison between 1 km resolution satellite-derived and localized in situ measurements. The field was seeded in September 2014, Septem-  Table 1). The field was instrumented by an eddy covariance (EC) system at a 2 m height. EC tower includes a CSAT3 3-D sonic anemometer that measures the wind and temperature fluctuations and a krypton hygrometer, KH20, that measures the concentration of water vapor. The EC tower is also equipped with a CNR1 radiometer (Kipp and Zonen) to measure the four components of the net radiation (R n ) with several heat flux plates (HFT3-L, Campbell Scientific Ltd) to measure the soil heat flux (G). Energy balance closure analysis indicated that the available energy (R n − G) was generally higher than the EC measurements. The relative closure was about 68 %, 76 %, 79 % and 79 % for S1, S2, S3 and B1, respectively. The sensible and latent heat fluxes (H and LE) were finally corrected to force the closure of the energy balance by the Bowen ratio method (Twine et al., 2000). LST is measured at the EC station by using two Apogee IRTS-P infrared radiometers, oriented downward and measuring the surfaceleaving radiance between 8 and 14 µm, set up at a 2 m height above the ground. An estimate of LST is obtained by averaging both measurements. The soil water content is measured at various depths (5, 10, 20, 30, 50, 70 cm) using Time Domain Reflectometry probes (model CS616) installed in a soil pit at the bottom of the EC tower. A weather station was set up nearby the studied field to measure air temperature, solar radiation, relative humidity, wind speed and rainfall at a 30 min time step.

MODIS
Three products from the MODIS sensor onboard the Terra and Aqua satellites are used in this study: the (1) MODIS Terra/Aqua Land Surface Temperature product (MOD11A1/MYD11A1), the (2) MODIS Terra Vegetation Indices product (MOD13A2) and the (3) MODIS Albedo (combined Terra and Aqua) product (MCD43A3). All products are gridded in the sinusoidal projection.
The MOD11A1 and MYD11A1 provide LST at 1 km spatial resolution under clear-sky conditions, derived from Terra and Aqua, respectively. Brightness temperatures from bands 31 and 32 are used to derive LST through a generalized split-window algorithm. The MOD13A2 provides several vegetation indices at 1 km resolution. One particular vegetation index of interest in this study is NDVI, available at 16day temporal intervals. This product is derived from bands 1 and 2 of the MODIS Terra satellite.
The obtained NDVI is used to derive the leaf area index (LAI) via the following formulation (Wang et al., 2013): The vegetation cover fraction is expressed as (Kustas and Norman, 1997) Finally, the MCD43A3 product provides the surface albedo (α) at 500 m resolution every 16 d. The latter is generated from both the Terra and Aqua products. In this work, the shortwave broadband α is used by integrating its value over the entire solar emission spectrum (0.3-5.0 µm). This value is obtained as a weighted average of the directional hemispherical reflectance (black-sky-α) and the bi-hemispherical reflectance (white sky α) using their two extreme values. The current α (called "blue sky") is a weighted average between these two extreme cases (Lewis and Barnsley, 1994). Herein, the percentages of 85 % and 15 % are used for the direct and diffuse lights, respectively.

SMOS
The SMOS mission measures the natural (passive) microwave radiation around the frequency of 1.4 GHz (L-band). It aims to monitor SM at a depth of about 3-5 cm with a spatial resolution of about 40 km and an accuracy better than 0.04 m 3 m −3 . The revisiting time at the Equator is 3 d for both ascending and descending passes, which are Sun synchronous at 06:00 and 18:00 LT, respectively. The SMOS level-3 1 d global SM product (MIR CLF31A/D) posted on the ∼ 25 km Equal Area Scalable Earth (EASE) version 1.0 grid is used as input to the DisPATCh algorithm.

DisPATCh
The DisPATCh remote sensing algorithm combines the coarse-scale microwave-retrieved SM with high-resolution optical/thermal data within a downscaling relationship to produce SM at a higher spatial resolution. Soil (T s,min , T s,max ) and vegetation (T v,min , T v,max ) temperature endmembers are estimated from the polygon obtained by plotting MODIS LST against MODIS NDVI, where the LST is partitioned into its soil and vegetation components according to the trapezoid method of Moran (Moran et al., 1994). Details on the DisPATCh algorithm and the methodology to determine dry and wet edges can be found in Merlin et al. (2012). The retrieved soil temperature is then used to estimate the soil evaporative efficiency (SEE), which is defined as the ratio of actual to potential soil evaporation. Finally, DisPATCh converts the high-resolution optically derived SEE fields into high-resolution SM fields given a semiempirical SEE model and a first-order Taylor series expansion around the SMOS observation. In our application, we applied DisPATCh to 40 km resolution SMOS level-3 SM and 1 km resolution MODIS optical/thermal data to produce SM at a 1 km resolution (Molero et al., 2016). The input data sets are composed of MODIS LST, MODIS NDVI and the GTOPO digital elevation model (DEM) used to correct LST for topographic effects Merlin et al., 2013).

TSEB-SM
The recently developed TSEB-SM is fully described in Ait Hssaine et al. (2018b). The equations and sub-equations used in TSEB-SM are provided in Table 2; the main equations are given below. The originality of TSEB-SM is to integrate SM observations in addition to LST and vegetation cover fraction data in order to calibrate both the soil resistance to evaporation (constant parameters) and the PT coefficient on a daily basis. The model is based on the original TSEB formalism, meaning that the energy balance for vegetation is the same as in TSEB using the PT formula, although the soil evaporation is estimated as a function of SM using a soil resistance developed by Sellers et al. (1992). The use of the soil resistance formulation is justified by the fact that its main parameters (a rss , b rss ) can be adjusted based on soil texture characteristics  or by combining SM and LST data under bare  or partially covered (Ait Hssaine et al., 2018b) soil conditions.
The surface soil heat flux is estimated as a fraction of R n,soil : where c g ∼ 0.35 (Choudhury et al., 1987).
The vegetation latent heat flux LE veg is estimated via the PT formulation: where α PT is the PT coefficient, f g the fraction of green vegetation, γ the psychometric constant (≈ 67 Pa K −1 ), the slope of the relationship between saturation vapor pressure and air temperature, and R n,veg the vegetation net radiation. Note that f g is set to 1 as the α PT coefficient is varied (through the calibration procedure) to take into account the fraction of transpiring vegetation. The soil latent heat flux is estimated using the resistance formulation: where e s is the saturated vapor pressure at the soil surface, e a the actual air vapor pressure, r ah the aerodynamic resistance calculated from the adiabatically corrected logarithmic temperature profile equation (Brutsaert, 1982) and r s the surface-soil resistance to transport of heat between the soil surface and a height representing the canopy estimated using (Sauer et al., 1995). Both resistances are simulated every 30 min (between 11:00 and 14:00 LT) and at Terra and Aqua overpass times for in situ and satellite data, respectively. r ss is computed as a function of SM and is expressed as (Chirouze et al., 2014;Li et al., 2006;Sellers et al., 1992) with SM being the 0-5 cm SM, a rss and b rss are two empirical parameters (to be calibrated) and SM sat is the SM at saturation expressed as (Cosby et al., 1984) with f sand being the sand percentage of soil. In Ait Hssaine et al. (2018b), an innovative calibration approach of α PT , a rss and b rss is developed from in situ SM and LST data (Ait Hssaine et al., 2018b). The calibration methodology is briefly explained below.

Retrieval and calibration of r ss , a rss and b rss
The r ss is first adjusted by minimizing a cost function defined by with T surf,sim and T surf,mes being the simulated and measured LST, respectively. The T surf,sim was simulated as follows: Cost function for minimizing r ss where T veg and T soil are the vegetation and soil components of temperature (K). The LST is simulated every 30 min (between 11:00 and 14:00 LT) and at Terra and Aqua overpass times for in situ and satellite data, respectively. The LST at the first calibration step is simulated with a constant value of α PT (average value of the α PT retrieved for f c > 0.5).
Then, for the second calibration step, it is simulated using the daily retrieved α PT . The inverted r ss is then correlated with the SM (in situ or DisPATCh) to determine the a rss and b rss parameters by considering that, when f c is lower than a given threshold (f c,thres ), the dynamics of total LE is mainly controlled by the temporal variation of soil evaporation, meaning that both soil parameters are estimated when the PT coefficient can be set to a constant value.

Daily α PT retrieval
Once both parameters a rss and b rss have been estimated, the PT coefficient is retrieved on a daily basis when f c is larger than f c,thres , by minimizing a cost function at the Terra and Aqua-MODIS overpass times: In fact, an iterative loop is run on soil (r ss ) and vegetation (α PT ) parameters to reach convergence of all parameters. LST and SM data are thus used for calibration, while the calibrated TSEB-SM is run on a daily basis using SM data as forcing solely (in addition to vegetation cover fraction data). In this paper, an improvement is made on the former version of TSEB-SM to normalize the output fluxes using the LST-derived available energy. Therefore, the new version of TSEB-SM uses both LST and SM data (in addition to vegetation cover fraction data) as forcing on a daily basis. In practice, the latent and sensible heat fluxes derived from the TSEB-SM model are re-computed using the TSEB-SM derived evaporative fraction (EF, defined as the ratio of latent heat to available energy) and the LST-derived available energy. The rationale is that numerous modeling studies have shown the regularity and constancy of EF during daylight hours in cloud-free days (Gentine et al., 2011;Lhomme and Elguero, 1999;Shuttleworth et al., 1989) and the EF has a strong link with SM availability (Bastiaanssen and Ali, 2003), which is an important factor for estimating latent heat flux. For that purpose, the LST data collected at the Terra and Aqua-MODIS overpass times are used to estimate the instantaneous R n and G. A ratio between the daily (obtained as an average value between Aqua and Terra overpass times) latent heat flux LE daily and the daily available energy (R n,daily − G daily ) is used to calculate an average daily EF: The daily EF and the instantaneous available energy (calculated using Terra and Aqua MODIS LST) are finally used to re-calculate the instantaneous TSEB-SM output of LE and H by the following formulas:

Uncertainty in TSEB-SM input data
The LST collected by MODIS at Terra and Aqua overpass times and the SM product derived at 1 km resolution from the DisPATCh algorithm applied to SMOS data are used as input to TSEB and TSEB-SM models. Validation of TSEB and TSEB-SM input data prior to the evaluation of model output is an important issue, because of the scale discrepancy between the spatial resolution (1 km) of MODIS/DisPATCh data and the footprint of the EC flux measurements that does not exceed 100 m (Schmid, 1994).  order to obtain a better estimate of daily SM at 1 km resolution. They found that the assimilation of DisPATCh data improved quasi-systematically the dynamics of SM. Figure 2 shows the scatterplots of MODIS LST (at Terra and Aqua overpass) versus in situ measurements for the four agricultural seasons separately. The obtained R 2 , root mean square error (RMSE), and mean bias error (MBE) are reported in Table 3. The statistical comparison shows strong linear correlations (0.76 ≤ R 2 ≤ 0.90) for all years. The RMSE is around 4 K for the S2 (2016-2017) and S3 (2017-2018) agricultural seasons, while it reaches 6 K for S1 (2014)(2015) and B1 (2015-2016), respectively. The observed scatter may stem from the fact that the localized (1 or 2 m wide) in situ LST is not fully representative of the 1 km resolution MODIS pixel (Ait Hssaine et al., 2018a;Yu et al., 2017). For all years (S1-3, B1), it can be seen that the MBE is negative. Note that the MBE is the greatest when the temperatures are largest. Such a systematic error is probably due to the nonrepresentativeness of the in situ LST observations when compared to the corresponding scale of MODIS observations. The DisPATCh products have been extensively evaluated, especially over semi-arid areas like the Marrakech region   during four agricultural seasons (S1, B1, S2 and S3), at the site level where the comparison between TSEB and TSEB-SM is undertaken. The statistical results including the coefficient of determination (R 2 ), the RMSE, and the MBE are reported in Table 3. The R 2 ranges from 0.27 to 0.55, the RMSE from 0.04 to 0.09 m 3 m −3 and the MBE from −0.05 to −0.03 m 3 m −3 . These results are encouraging considering the heterogeneous land use composed of rainfed wheat, bare soil, and fallow and farm building (see Fig. 4).
In fact, the localized in situ measurements may not be perfectly representative of the 1 km resolution satellite data. Note that the efficiency of DisPATCh is supposedly higher for low SM values , which is clearly illustrated during the B1 season, while it is lower for high SM values (after rain events). This can be explained by the constraints of atmospheric and vegetation conditions on disaggregation results as well as the saturation of SEE in the higher SM range. Another major issue that can lead to differences between DisPATCh and in situ SM is that the ground SM sensors are buried at a depth of 5 cm, while the penetration of the L-band wave varies between 2 and 5 cm depend-ing on soil conditions (notably SM content, texture). For S2, the SM provided by DisPATCh underestimated field measurements, especially in the higher SM range. This particular behavior could be explained by the particularly low precipitation amount during this year. It is especially possible that the surrounding plots were not sown by neighboring farmers, resulting in a soil that dried quickly compared to our field, which retained the SM for a longer period of time.
Note that despite the relative heterogeneity within the 1 km pixel (characterized by rainfed wheat in addition to bare soil and fallow), the comparison between field measurements and 1 km resolution satellite data reflects acceptable accuracies.

Results and discussion
In this section, the a rss and b rss parameters and the α PT are firstly retrieved by following the two-step calibration based on a threshold of f c (cited in the Methods section). Then, the obtained calibrated values are used to estimate the surface fluxes using TSEB-SM. Finally, TSEB-SM fluxes are evaluated against the eddy covariance measurements, and re- sults are compared with the original TSEB. To facilitate the interpretation of the simulation results using MODIS and SMOS/DisPATCh data as input, the calibration and validation steps are previously tested using in situ (LST and SM) data.

Retrieving a rss and b rss parameters
The soil resistance r ss is inverted for f c ≤ f c,thres between 11:00 and 14:00 LT and at the Terra and Aqua overpass time step for in situ and satellite data, respectively. The result of this inversion is correlated with the actual to saturated soil moisture ratio SM/SM sat to determine a rss and b rss parameters. The calibration process is applied for each season independently. Then a pair (a rss , b rss ) is calculated for the entire study period for in situ and satellite data, respectively. Figure 5a and b plot the ln(r ss ) versus in situ SM/SM sat using in situ and satellite data, respectively. The mean retrieved values (7.62, 2.43) and (7.32, 4.58) for in situ and satellite data, respectively, are close to the values found in Li et al. (2006) (8.2, 4.3) and in Ait Hssaine et al. (2018b) (7.2, 4). However, by comparing both figures ( Fig. 5a and b), one notes that the use of in situ data generates more scatter than with satellite data. The apparent scatter in retrieved r ss could be interpreted by the impact of the daily cycle of meteorological (evaporative demand) conditions or soil property differences (Merlin et al., 2011;Merlin et al., 2016;Merlin et al., 2018). The retrieved soil parameters also vary from year to year: the standard deviation is 0.39 and 1.69 for a rss and b rss , respectively. This can be explained by the compensation effects linking a rss and b rss parameters which prove the empirical nature of the r ss . Another major issue that can lead to these differences is the depth of SM measurements (Merlin et al., 2011). In Sellers et al. (1992), the near-surface soil moisture is defined in the 0-5 cm soil layer, whereas in our field, SM measurements are made at 5 cm depth. Also, the sensing depth of SMOS observations is generally shallower than the in situ surface measurements . Moreover, the variability of a rss and b rss in Fig. 5b using remote sensing data can be linked to the scale difference between DisPATCh SM/MODIS products (1 km) and the field measurements. As shown in Fig. 3, the field is surrounded by trees, buildings and fallows, which causes the spatial heterogeneity within the pixel of 1 km. This heterogeneity can introduce errors into the model inversion. Nevertheless, soil parameters are quite similar for in situ and satellite data sets. Therefore, the heterogeneity issues within the 1 km pixel scale are minor in this study.

Time series of daily retrieved α PT
The second calibration step consists in inverting the daily α PT when vegetation is covering a significant part of soil (f c > f c,thres ), for the three seasons of rainfed wheat (S1-S3), by using in situ data and satellite data, separately. Herein, the calibration of α PT is bounded by minimum (0) and maximum (2) acceptable physical values, in order to avoid unacceptable values of α PT that can be produced because of the uncertainties in daily LST estimates. Such an upper bounding is especially needed when vegetation partially covers the soil.  Figure 6 plots the daily variation of α PT for each season (S1-S3) separately, using in situ data. The mean retrieved values of α PT are 1.26, 1.12 and 1.09 for S1, S2 and S3, respectively. In all cases, the mean α PT is close to the theoretical α PT value (1.26). It is well observed that the retrieved α PT for S1 is slightly larger compared to those obtained for both S2 and S3. This can be explained by the timing and amount of rainfall during each season. Note that unexpected low values of α PT are recorded for S3 during the first few days (25 January-4 March) of the development stage. They may be associated with uncertainties in retrieved α PT , as the impact of soil surface is still significant, as well as with a relatively low evaporative demand especially since this period coincides with cloudy days and abundant precipitations. Indeed, the coupling between transpiration (and hence retrieved α PT ) and LST is expected to be lower under lower atmospheric demand.

Using in situ data
The retrieved α PT is then smoothed as in Ait Hssaine et al. (2018b) to remove outliers and to reduce uncertainties at the daily timescale. The smoothed values of α PT range from 0 to 1.54, 0 to 1.38 and 0.45 to 1.43 for S1, S2 and S3, respectively. The maximum of α PT is close to 1.26 for S2, while it is higher for S1 and S3. This result is in accordance with the total rainfall amounts, which were about 608, 214 and 421 mm for S1, S2 and S3, respectively. Additionally, one can state that the stability of α PT strongly depends on the rainfall distribution along the agricultural season. The daily α PT is more stable for S1 than for S2 and S3. Indeed, the amount of rain during S1 is very important, with two peaks of about 83 mm that occurred at the beginning of the season and during the growing stage. The second one coincides exactly with the maximum value of the retrieved α PT . However, different results are obtained for S2 compared to S1 due to the lowest precipitation amount recorded over that season. As shown in Fig. 6, the amount of rain is concentrated at the beginning of the growing stage (mid December), when the α PT peaks. Afterward, the smoothed α PT tends to decrease because of insufficient soil water reserve in the root zone to enable wheat to continue growing. Rainfall is also significant for S3, and every rainfall event causes an immediate (daily) response of α PT (after 4 March). As mentioned before, the significant error in α PT retrievals for S3 between 25 January and 4 March induces strong uncertainties in the smoothing function estimates.

Using satellite data
Figure 6 also illustrates the daily variation of α PT retrieved from satellite data for each season separately. S1 and S2 have a very similar distribution of the retrieved α PT as compared to the retrieved α PT using in situ data, respectively. For S3, only six retrieved α PT values are available because of the Figure 6. Time series of daily retrieved and smoothed α PT (calibration step 2 -using in situ data and satellite data) collected during S1, S2 and S3. non-availability of MODIS products during cloudy days. For this reason, no information linked to the variability of α PT can be derived during this season. The retrieved values are smoothed and superimposed with the rainfall events. It is clearly shown that the smoothed α PT for S1 and S2 have the same shape with a small variability when compared with the smoothed α PT using in situ data, resulting in an error estimated as the RMSE to the mean α PT ratio, of about 11 % and 19 %, for S1 and S2, respectively. For S1 the maximum of smoothed α PT is reached at the same time as when using the in situ data, with a value of about 1.38, while the maximum for S2 is reached 10 d before the maximum of the α PT derived from in situ data with little response of α PT to rainfall events. These differences may be linked to uncertainties in disaggregated SMOS SM as well as to the weaker availability of satellite data. Because of the small number of data points (retrieved α PT ) during S3, the smoothed α PT remains at a mostly constant value (∼ 0.7) throughout the study period, with a significant relative difference of about 34 % when compared with the α PT retrieved using in situ data.

Interpretation of α PT variabilities
Figure 7 plots variation of calibrated daily α PT , superimposed with NDVI and rainfall events. It is visible that the maximum value of NDVI appears sooner than the maximum value of α PT for both S1 and S3. Such a delay is attributed to the high soil moisture level in the root zone during the maturity stage. Later in the season, α PT decreases as NDVI starts to decline at the onset of senescence. In contrast, the maximum value of NDVI appears later than the maximum value of α PT for S2. This can be explained by the fact that rainfall at the beginning of the development phase satisfies the plant requirements, while the rainfall amount during the development stage is relatively low compared to the crop water needs (Kharrou et al., 2011). Large variations in α PT occur during the agricultural season, as a result of the amount, frequency, and distribution of rainfall along the season. In general, the analysis of the α PT variability using satellite data illustrates the robustness of the proposed approach, which combines microwave and optical/thermal data to retrieve a water stress indicator at the daily timescale.

Surface fluxes
The robustness of TSEB and TSEB-SM for partitioning (R n − G) into H and LE is evaluated using in situ and remotely sensed LST, SM, and NDVI, separately, at Terra and Aqua MODIS overpasses. Figure 7. Time series of calibrated daily α PT (red -using in situ data, green -using satellite data) superimposed with NDVI and the rainfall events during S1, S2 and S3, separately. Figure 8 shows an intercomparison of simulated and observed LE for the four seasons separately. TSEB-SM clearly provides improved results compared to the original TSEB.

Using in situ data
The obtained values of RMSE by TSEB-SM are about 68 and 72 W m −2 for S1 and S2, respectively, which is significantly lower than those revealed by TSEB (109 and 86 W m −2 , respectively) (see Table 4). For B1 (season of bare soil), TSEB largely overestimates LE with a MBE of about 165 W m −2 compared to TSEB-SM, which yields a MBE of 59 W m −2 . This overestimation of TSEB is most probably related to an inadequate value of α PT (= 1.26) for bare soil surfaces. In fact, 1.26 is an optimum value for the potential transpiration rate (Agam et al., 2010;Chirouze et al., 2014). In the case of TSEB-SM, biases are reduced thanks to the calibration of the r ss resistance. Additionally, according to TSEB-SM assumptions, α PT for f c ≤ 0.5 is set to the average value of the α PT retrieved for f c > 0.5. During the B1 season (bare soil conditions), α PT was hence obtained as an average value of the mean α PT retrieved for all seasons S1, S2 and S3 when f c > 0.5 (α PT ∼ 1). However, this value remains relatively high for a bare soil, which yields a slight overestimate of LE measurements (see the B1 case in Fig. 8).
For the S3 season, the error in daily retrieved α PT at the beginning of the development stage has a strong impact on LE predictions and thus yields to greater discrepancies illustrated in Fig. 8. To overcome this error, the threshold on f c to separate calibration steps 1 and 2 was increased to 0.63 (arbitrary value). The TSEB-SM model is then run using the new threshold. The LE simulations are improved, with a RMSE of 73 W m −2 instead of 98 W m −2 and a relative error (estimated as the RMSE divided by the mean observed LE) of about 42 % instead of 58 %. The increase in the threshold is intended to decrease the uncertainties in α PT retrievals when vegetation is not fully covering the soil. It can be concluded that the errors in α PT retrievals have a strong impact on LE estimates.
The ability of TSEB-SM to estimate the sensible heat fluxes is also investigated. Figure 9 displays the comparison between TSEB and TSEB-SM for each season and Table 4 summarizes the different statistical parameters. One can notice that TSEB shows greater discrepancies in H estimation, with a RMSE of about 127, 112 and 103 W m −2 and a MBE of about −41, 1, and −71 W m −2 for S1, S2 and Figure 8. Scatterplot of simulated versus observed LE for (top panels) TSEB and (bottom panels) TSEB-SM models using in situ data collected during S1, B1, S2 and S3, respectively. Table 4. Statistical results (RMSE, R 2 and MBE) between modeled and measured sensible and latent heat fluxes for S1, S2, B1 and S3, and for TSEB and TSEB-SM models, separately (R n and G are forced to their measured value  et al., 2018b). The discrepancies between TSEB-SM and in situ H during S3 are mostly rectified by using the new threshold on f c : the statistical results are improved, the RMSE is about 73 W m −2 and the relative error is 39 % (instead of 52 %). It can be concluded that the uncertainty observed over the α PT during the first few days of the development stage (25 January-4 March) is mainly related to the impact of the soil, which is not negligible during the first weeks of the growing stage. Nevertheless, by considering the overall results obtained for the three seasons, the threshold of f c,thres = 0.5 can be considered an acceptable value to calibrate the soil resistance parameters and the Priestly-Taylor coefficient.
As a further step, the intercomparison between TSEB and TSEB-SM is evaluated by predicting R n and G fluxes instead of forcing them to their measured values. The statistical results of the comparison between simulated and observed R n , G, H and LE are listed in Table 5. The scattering obtained when comparing turbulent flux estimations to measurements is mainly related to the uncertainty in available energy estimates, mainly related to the uncertainty in soil heat flux estimates. Indeed, as reported in Table 5, R n is very well simulated for both TSEB and TSEB-SM. The R 2 between simulated and observed R n is about 0.99 during all seasons. Meanwhile, G shows a poor correlation, with an R 2 varying from 0.05 to 0.45. This is mainly linked to the approach used to estimate G, which requires local calibration. Kustas et al. (1998) hence indicated that the ratio G/R n,soil cannot be considered a constant, because it is affected by different factors such as time of day, moisture conditions and soil texture and structure.

Using satellite data
In order to gain greater insight into how TSEB and TSEB-SM models respond to different surface conditions across a landscape, an analysis of the spatial distributions and the magnitude of the turbulent fluxes using remotely sensing data produced from the two models is conducted. The comparisons between TSEB/TSEB-SM versus observed LE over the four seasons are illustrated in Fig. 10. Figure 10 indicates that TSEB overestimates latent heat flux. The overall MBE values are about 119, 181, 94 and 128 W m −2 for S1, B1, S2 and S3, respectively. The overestimation of LE fluxes can be explained by the fact that α PT is set to 1.26 during the entire agricultural season including stress conditions. This probably causes larger errors in the LE estimation, especially during the growing stage. Indeed, the saturation of TSEB during the senescence period is precisely caused by the PT coefficient fixed to 1.26. The errors are reduced when using TSEB-SM. In fact, the constraint on plant transpira-   tion, while retrieving daily α PT values improves ET estimates, especially for the growing stage. Moreover, during the senescence stage the large positive bias of LE is considerably reduced. In fact, the decrease in calibrated daily α PT is associated with the drop in NDVI during senescence (Ait Hssaine et al., 2018b). Additionally, the constraint on the soil evaporation via the DisPATCh SM clearly reduces the MBE values during the emergence period (f c ≤ f c,thres ). Finally, the constraint applied on TSEB-SM output fluxes using LSTderived available energy and the TSEB-SM-derived evaporative fraction (Eq. 8) improves the LE estimates for the whole study period. The MBE values are about 39, 4, 7 and 62 W m −2 for S1, S2, S3 and B1, respectively.
The intercomparison between TSEB and TSEB-SM is made by forcing the available energy to its measured value. The statistics listed in Table 5 indicate that there are similar differences between modeled versus measured R n using either TSEB or TSEB-SM. Overall, the discrepancies between estimated and measured R n are likely due to a greater scatter between MODIS and in situ measured LST. Note that RMSE values up to 6 K have been noted when comparing LST MODIS with ground-based measurements. These uncertainties are likely to be explained by the huge scale mismatch between the 1 km resolution of MODIS LST and the footprint size (approximately 1 m) of ground-based radiometers. The uncertainties in key input data generate large differences in simulated R n compared to the tower measurements. The greater scatter between modeled and measured G from the two models reflects the fact that there is a major mismatch in scale between the area sampled by the soil heat flux sensors and the 1 km resolution of model inputs. It appears that the LE estimates from TSEB-SM are generally in closer agreement with the measurements than the TSEB model outputs. The RMSE is significantly improved from 103 to 52 W m −2 , from 151 to 30 W m 2 , from 101 to 35 W m −2 and from 83 to 24 W m −2 , during S1, B1, S2 and S3, respectively. For the sensible heat flux H , the difference between TSEB estimates and EC measurements listed in Table 5 indicates a fairly large underestimation, the MBE values varying between −56 and −240 W m −2 . However, the TSEB-SM output provides a quite significant improvement, with an absolute MBE lower than −61 W m −2 during all seasons.

Evaluation of H and LE estimates
In this section, the residual error of the H and LE estimated with the TSEB-SM-retrieved soil/vegetation parameters is analyzed. Figure 12 plots retrieved α PT vs. residual H and LE error. The retrieved α PT is poorly correlated with residual H (R = −0.27) and ET (R = 0.27) errors, especially for seasons S1 and S2. For season S3, few retrieved α PT values were available because of the non-availability of MODIS products during cloudy days. It is shown that the trend between α PT and residual H error is slightly negative for S1, while it is slightly positive for S2. According to these results, no information linked to the variability of α PT versus residual ET and H errors can be derived. Figure 13 plots retrieved r ss vs. residual H and LE error and LST for the four study periods. The retrieved r ss is negatively correlated (R = −0.33) with residual H error (predicted-observed) for the four seasons, while it is positively correlated (R = 0.33) with residual LE error. The residual error covers a wide range (between −150 and 150 W m −2 ) for the lower r ss values, while it is biased for the higher r ss values. Such a result indicates that the r ss formulation as a function of nearsurface SM needs further improvements Merlin et al., 2018) in order to reduce systematic uncertainties in evaporation estimates, especially in dry (moisturelimited) conditions. Consistent with the general decrease in LST with SM, LST is positively correlated with retrieved r ss (R = 0.45). This is very coherent since r ss decreases with the  increase in SM. Regarding the sensitivity analysis of residual H and LE errors to observed SM, Fig. 14 shows that SM is positively correlated with residual H error, while it is negatively correlated with residual LE error for the entire study period. The correlation coefficient is about 0.3 when using DisPATCh SM, while it is about 0.4 when using in situ SM. This difference can be explained by the uncertainty (including spatial representativeness issues at the localized scale) in DisPATCh SM. The positive correlation coefficient between residual LE error and observed SM is likely to be due to the systematic errors in r ss estimates for dry conditions as mentioned previously. For S2, B1 and S3, the residual H error ranges between −150 and 50 W m −2 for SM between 0 and 0.10 m 3 m −3 , while it is slightly overestimated for the higher range of SM. The residual LE error is also found to be influenced by SM, but in the opposite sense.

Conclusions
The microwave-derived near-surface soil moisture (SM) from SMOS and the thermal-derived land surface temperature (LST) from MODIS are integrated simultaneously in the TSEB formalism within a calibration procedure to invert both the soil resistance to evaporation (constant parameters) and the PT coefficient based on a threshold on f c . The TSEB-SM model is applied during a 4-year period (2014-2018) over a rainfed wheat field in the Tensift basin, central Morocco. Significant conclusions are given below.
The constraint applied on the soil evaporation represented explicitly as a function of SM via a soil resistance term reduces the errors when using TSEB-SM instead of TSEB.
-The first step of the TSEB-SM approach is to calibrate r ss for (f c ≤ f c,thres ) at Terra and Aqua overpass times. Despite the scale difference between the MODIS/DisPATCh resolution data and the footprint size of in situ measurements, the parameters (a rss , b rss ) calculated for the entire study period using satellite data are relatively close to those derived from in situ measurements.
-The second step of the TSEB-SM approach is to invert the α PT on a daily basis for f c > f c,thres by using LST and SM data. The maximum values of daily cali-brated α PT are 1.38, 1.25 and 0.87, when using satellite data, for S1, S2 and S3, respectively. Those values are in accordance with the total rainfall amounts, which were about 608, 214 and 421 mm per wheat season for S1, S2 and S3, respectively. S1 and S2 have the same distribution of daily calibrated α PT when compared with the α PT retrieved using in situ data, while the retrieved α PT remains at a mostly constant value (∼ 0.7) throughout the S3 study period because of the non-availability of MODIS products during cloudy days.
-Finally, an analysis of the spatial distributions and the magnitude of the turbulent fluxes using remotely sensing data produced from the two models were conducted. TSEB exhibits larger errors in H and LE estimates. These uncertainties can be linked to the theoretical value of α PT , which is fixed to 1.26 for the whole study period as well as to the scale mismatch between the 1 km resolution of MODIS LST and the footprint size (approximately 1 m) of the ground-based radiometer.
As a short-term prospect, the use of high-resolution products from active sensors (Sentinel-1) would allow application of the TSEB-SM approach at the field scale over heterogeneous (e.g., irrigated) landscapes. Also, the robustness of TSEB-SM in terms of evaporation/transpiration partitioning will be tested by using independent flux measurements derived from lysimeters and sap flow sensors (Rafi et al., 2019). In addition, the evaluation of ET at large scale is missing. Spatialized measurements that could be collected by scintillometers installed at various points in the region would be a solution for that purpose.
Data availability. The data that support the findings of this study are available from the authors, but restrictions apply to the availability of these data, which were carried out within the framework of several projects and so are not publicly available. Data are however available from the authors upon reasonable request and with permission of the co-directors of the LMI TREMA laboratory.
Author contributions. BAH and OM developed the idea of the study, designed the methodology of the work, verified/discussed the results and wrote the original draft of the manuscript. BAH analyzed the data and validated the results. BAH, OM, and JE investigated the experimental protocol and contributed to materials/analysis tools. BAH and NO contributed to remote sensing data processing. OM and SK supervised the work. BAH, OM, JE, SK, and SE contributed to the final manuscript.
Competing interests. The authors declare that they have no conflict of interest.