<figure><img src="image_1.png" /></figure>

HAL
open science

# Analysis of Multispectral Drought Indices in Central Tunisia

Nesrine Farhani, Julie Carreau, Zeineb Kassouk, Michel Le Page, Zohra Lili

Chabaane, Gilles Boulet

► To cite this version:

Nesrine Farhani, Julie Carreau, Zeineb Kassouk, Michel Le Page, Zohra Lili Chabaane, et al.. Analysis of Multispectral Drought Indices in Central Tunisia. Remote Sensing, 2022, 14 (8), pp.1813. 10.3390/rs14081813. [hal-03821759](https://hal.science/hal-03821759)

HAL Id: hal-03821759

https://hal.science/hal-03821759

Submitted on 20 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Article

# Analysis of Multispectral Drought Indices in Central Tunisia

Nesrine Farhani<sup>1,2,*</sup> and Julie Carreau<sup>3</sup>, Zeineb Kassouk<sup>1</sup>, Michel Le Page<sup>2</sup>, Zohra Lili Chabaane<sup>1</sup>

<figure><img src="image_4.png" /></figure>

check for
updates

**Citation:** Farhani, N.; Carreau, J.; Kassouk, Z. ; Le Page, M.; Lili Chabaane, Z.; Boulet, G. Analysis of Multispectral Drought Indices in Central Tunisia. *Remote Sens.* **2022**, *14*, 1813. https://doi.org/10.3390/rs14081813

Academic Editors: Alban Kuriqi and
Luis Garrote

Received: 8 March 2022

Accepted: 6 April 2022

Published: 9 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_5.png" /></figure>

**Copyright:** © 2022 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

**Abstract:** Surface water stress remote sensing indices can be very helpful to monitor the impact of drought on agro-ecosystems, and serve as early warning indicators to avoid further damages to the crop productivity. In this study, we compare indices from three different spectral domains: the plant water use derived from evapotranspiration retrieved using data from the thermal infrared domain, the root zone soil moisture at low resolution derived from the microwave domain using the Soil Water Index (SWI), and the active vegetation fraction cover deduced from the Normalized Difference Vegetation Index (NDVI) time series. The thermal stress index is computed from a dual-source model Soil Plant Atmosphere and Remote Evapotranspiration (SPARSE) that relies on meteorological variables and remote sensing data. In order to extend in time the available meteorological series, we compare the use of a statistical downscaling method applied to reanalysis data with the use of the unprocessed reanalysis data. Our study shows that thermal indices show comparable performance overall compared to the SWI at better resolution. However, thermal indices are more sensitive for a drought period and tend to react quickly to water stress.

**Keywords:** indicators; droughts; evapotranspiration; reanalyses

## 1. Introduction

Droughts are a recurring natural climate event that result from a reduction in precipitation amount received over an extended period of time, such as a season or a year [1]. In arid and semi-arid areas characterized by a significant temporal and spatial climatic variability, the vulnerability advanced by recurrent droughts is considerable as it makes serious threats to agroecosystem health and productivity [2]. More specifically, in several Mediterranean countries, we observe during the last decade a significant warming trend, more pronounced in summer, and a decrease in rainfall during the wet season [3]. This could have a strong impact on water resources and constitutes the main driver of agricultural droughts which appear as a result of a long-term period of precipitation deficiency and lead to lower soil moisture. Droughts at a sensitive stage of crop development (emergence, pollination, and grain filling) can lead to significant damages [4,5] specifically on rainfed agriculture. Thus, an important issue for rainfed agriculture is to receive adequate rainfall at the appropriate timing. In contrast, drought responses over irrigated areas are related to surface or groundwater resource availability which can also suffer a severe drop during prolonged periods of drought. Therefore, it is crucial to improve the monitoring of agricultural droughts and the prediction of their occurrence in the future [6]. In this work, we are interested in agronomic drought identification.

Several drought indices have been developed to quantify agronomic drought periods. Initial drought indices rely essentially on meteorological variables. Most meteorological

<footnote n="1">¹ Institut National Agronomique de Tunisie, Université de Carthage, LR17AGR01 (Lr GREEN-TEAM), Tunis 1082, Tunisia; zeineb.kassouk@inat.u-carthage.tn (Z.K.); zohra.lilichabaane@inat.u-carthage.tn (Z.L.C.)</footnote>

<footnote n="2">² Centre d'Études Spatiales de la Biosphère, Université de Toulouse, CNRS, CNES, IRD, UPS, INRAE, 31400 Toulouse, France; michel.le_page@ird.fr (M.L.P.); gilles.boulet@ird.fr (G.B.)</footnote>

<footnote n="3">³ Department of Mathematics and Industrial Engineering, Polytechnique, Montreal, QC 3453, Canada; julie.carreau@polymtl.ca</footnote>

<footnote>* Correspondence: nesrine.farhani@ird.fr; Tel.: +216-20-212-738</footnote>
indices use precipitation either on its own or in combination with other meteorological vari-
ables [6]. The most commonly used meteorological index is the Standardized Precipitation
Index (SPI) [7], which classifies precipitation intensities according to a probability distri-
bution derived from historical records, over different durations (1 month SPI, 3 months
SPI, etc.). SPI is particularly useful for drought identification and frequency analysis [7–9].
However, meteorological indices bring information on meteorological droughts defined
only in terms of rainfall shortage, but not specifically on their consequences on vegetation
for instance, especially for agronomic droughts. On the other hand, remote sensing (RS)
is an effective tool to provide information on whether vegetation functioning is optimal
or sub-optimal. Indeed, RS has shown high potential to retrieve land surface properties
related to plant functionality, such as vegetation cover, from the visible/near infrared
domain [10,11], surface soil moisture from the microwave remote sensing using passive
or active radiometers [12–15], or surface temperature from the thermal infra-red (TIR)
data that exploits surface energy balance components to identify land surface energy
interactions [16].

RS observations of green biomass vegetation indices are relevant indicators of year-to-
year drought conditions and vegetation greenness, especially over the past few decades
thanks to the increasing number of multi and hyperspectral sensors [11,17], such as
Advanced Very High Resolution Radiometer (AVHRR), SPOT-Vegetation and Moderate-
Resolution Imaging Spectroradiometer (MODIS). These sensors provide spatial resolution
ranging from 250 to 1000 m. Historically, the Normalized Difference Vegetation Index
(NDVI) [18], built from the normalized difference between near-infrared (NIR) and RED
reflectances, has been used for decades to provide information on vegetation health due to
its strong correlation with the photosynthetic activity of plant cover. A suite of drought
indices based on NDVI were subsequently developed in order to provide more focused
information, such as the vegetation condition index (VCI) [19]. VCI uses the NDVI infor-
mation provided from each pixel normalized with the maximum and minimum statistical
range over the available time series of images. It has been tested in several regions of
the world, showing a relevant efficiency to characterize periods of droughts [20,21]. The
vegetation anomaly index (VAI) developed by [22], also forms a standardization of the
NDVI time series. Although RS-based vegetation observation indices have proven their
efficiency to assess drought impacts on vegetation growth, they are not able to detect
incipient stress because of the time lag needed for the plant to reflect drought-induced
physiological changes. Therefore, a delayed response to prolonged stress will be provided
using these indices [23]. Response of vegetation to drought is a gradual process altering
soil moisture, that leads to a time lag effect [24].

Surface soil moisture (SSM), which is inferred from microwave RS, is another efficient
key parameter for agricultural planning and water resource management [12,25]. Neverthe-
less, for drought monitoring and agricultural modeling, a representation of root-zone soil
moisture is needed [26]. A number of soil moisture-based indices have been developed and
widely used for drought monitoring, such as the Soil Water Index (SWI) [27]. This index is
retrieved from the Advanced SCATterometer (ASCAT) SSM estimates in C band, which
was found to be more suitable in bare soils or sparse vegetation cover [28], to determine soil
moisture (SM) distributions within the soil profile using an exponential relationship that
links surface and subsurface SM via a characteristic time length [29]. SWI is also used on
analysis of statistical anomalies over several years to derive the Moisture Anomaly Index
(MAI) [15], that has shown a good performance in providing quantitative information on
drought periods over the central region of Tunisia. The Soil Moisture Deficit Index (SMDI),
developed by [30], was helpful for planning to mitigate drought impacts. This index is
calculated from weekly percentage SM deficit or excess using long-term median, maximum
and minimum soil water [30]. The main limitation of drought indices based on passive
microwave data is their low resolution. Moreover, quantitative soil moisture estimation is
still a challenge nowadays, especially under vegetation cover [12,31].
RS maps of land surface temperature (LST) are very informative of water availability and form a good indicator of incipient droughts [32]. Indeed, water stress induces a stomatal closure that generates in turn an elevated canopy temperature [33–37]. Thermal channels of Landsat, AVHRR or MODIS were exploited to retrieve remotely sensed LST estimates. LST can be used as a signature of the land surface energy budget and, in particular, to determine whether the dissipation of available energy is more into sensible (dry conditions) and/or latent heat (wet conditions) [38]. It also enables the computation of evapotranspiration (ET) from latent heat, computed as the residual of the surface energy budget [39,40]. ET is particularly informative about water stress conditions [40,41]. Indeed, the crop water requirements are adjusted to match the water losses from actual evapotranspiration (ETa) [42]. Consequently, in water stress conditions, the plant reduces its transpiration in comparison with unstressed vegetation under the same atmospheric conditions [43]. Thus, the water stress intensity may be inferred from the ratio of actual and potential (unstressed) ET. Several drought indicators that integrate ET have been explored and showed improved drought monitoring. In [44], the authors proposed the so-called Water Deficit Index (WDI) based on the approach Vegetation Index/Temperature Trapezoid (VIT). WDI is estimated for each pixel to retrieve dry/wet conditions based on the fractional vegetation cover and the difference between LST and air temperature (Ta). A good example of TIR indices is the evaporative stress index (ESI) [35] that describes anomalies in the actual and the reference ET ratio. ESI was shown to provide early warning particularly for flash droughts that could produce damaging impacts over short time periods [5]. The stress index (SI) is defined as the ratio between actual and potential evapotranspiration rates.

In summary, several drought indices are commonly used for drought monitoring and impact assessment in agriculture. However, it is hard to assess their respective performance in deciphering stress/non stress and drought/non drought situations. An inter-comparison between the performance of these different drought indicators provided from different wave lengths, in terms of consistency, reliability and ability to detect incipient water plant stress, is thus relevant. Kogan [19] evaluated VCI and the Temperature Condition Index (TCI) in different regions of the world. The analyses show that VCI has an excellent ability to detect each period of stress, the drought onset, intensity, duration, and impact on vegetation. The TCI provides additional information about vegetation stress through LST. In Kogan [19], the authors used VCI and TCI to develop the so-called Vegetation Health Index (VHI), a well known combined stress index that is widely used for drought detection and assessment of drought severity and duration [45]. However, as a simple averaging is performed, there is no real theoretical support for the way to weight the relative impact of stress (through LST) and vegetation development (through NDVI). The relationship between the changes in canopy temperature and the soil water supply in the fields is widely applied as a drought indicator [46]. The Vegetation Supply Water Index (VSWI), which is the ratio between the LST and NDVI, is found reliable to detect vegetation stress, moisture and drought-affected areas [46,47]. However, in this work, we suggest to compare different indices provided from individual biophysical variables (1) the “Normalized difference vegetation index” (NDVI) from the solar (Visible/Near InfraRed spectrum), (2) “Soil Water Index” (SWI) from the microwave domains and, (3) a “Thermal stress Index” (SI) from Thermal InfraRed retrieved from a surface energy budget model. Analyses are carried out in the Kairouan area in central Tunisia which is subject to a semi-arid climate. The thermal stress index used is based on ET simulated from a dual source energy balance model that provides robust estimates of ET when meteorological forcing and vegetation cover are accurately known [48]. A statistical downscaling method [49] combined with reanalysis data is used to generate surrogate series in order to extend the observation period. Then, simulated ET series from the downscaling method or from the unprocessed reanalyses are used to constrain the dual-source model Soil Plant Atmosphere and Remote Evapotranspiration (SPARSE) [50]. In this paper, we address three objectives:

(1) We first assess at a regional scale, the consistency between those indices and a precipitation-based drought index called the “Uniformized Precipitation Index” (UPI)
which is a variant of the SPI without the transformation into Normal quantiles. We assume
that a robustness of the indices at very low resolution (12 km) translates into a good
accuracy at higher resolution (kilometric).

(2) We assess the reliability of very low (UPI, SWI) and low (SI, NDVI) resolution indices at higher resolution (kilometric) at local scale, by focusing on drought monitoring for pixels either with a high proportion of rainfed wheat or that correspond to the extra large aperture scintillometer (XLAS) transect;

(3) We compare SI when the energy balance model is either forced with the downscaled meteorological data or with the unprocessed reanalysis data, in order to test the applicability of SI with routinely available data at a global scale.

## 2. Study Area

The study site, the Merguellil catchment (Figure 1), is located in central Tunisia, which is typical of semi-arid environments. The study area is influenced both by the Mediterranean climate (dry subhumid) and the pre-Saharan climate (arid). It is characterised by the inter-annual irregularity of precipitation, with an average annual rainfall of about 300 mm per year, with a short rainy season mostly between September and May [51] and with a high evaporative demand of about 1600 mm per year. In our study region, the annual mean air temperature did not show considerable variation during the study period. In Figure 2, we observe that this variable is mainly steady during the study period, ranging from about 0 °C and 45 °C with a median of about 20 °C.

<figure><img src="image_6.png" /><figcaption>**Figure 1.** Merguellil catchment in central Tunisia with seasonal (winter–spring) land cover map of 2019–2020.</figcaption></figure>

The upstream and the downstream sub-catchments are separated by El Haouareb dam that protects the village from inundations and provides surface irrigation water for the plain. However, most of the water used for irrigation is extracted from ground water. Merguellil is emblematic of hydrological processes which have been profoundly modified by human activities since the end of the 1960s, by the inclusion of soil and water conservation works, the construction of large and small dams, and intensification of irrigated farming [51,52]. As expected, agriculture is the main water consumer of available water resources in this region (around 80 of the available water%) [53]. Therefore, the aquifer system shows a considerable decrease over time and space due to over-exploitation [52]. This makes
the Merguellil catchment an interesting case study to investigate the ability to monitor agricultural water management in the recent years. In this work, we are interested only by the lower sub-catchment (3000 km²) which forms an alluvial plain mainly flat with altitude between 50 and 200 m. The plain incorporates mainly small cultivated areas [54], with characteristic vegetation of semi-arid regions: rainfed agriculture (olive trees and cereals) and vegetables (melons, chilis and tomatoes), see Figure 1. The gauged network in the Merguellil plain catchment has three stations (Figure 1), installed since 2012 for the earliest ones.

<figure><img src="image_7.png" /><figcaption>Figure 2. Annual air temperature (using ERA5 data) variation over the study period (2000–2019).</figcaption></figure>

## 3. Data and Methodology

### 3.1. Biophysical Indices Derivation

#### 3.1.1. NDVI

We used the MOD13 version 6 product from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor on board the NASA Terra and Aqua Earth Observation System satellites that provides each 16-day composite at 1 km of spatial resolution. Time series are retrieved between 2000 and 2019 over our study area (see Figure 1). NDVI uses the normalized difference between near-infrared (NIR) and (RED) reflectances.

$$NDVI = \frac{\rho_{NIR} - \rho_{RED}}{\rho_{NIR} + \rho_{RED}} \qquad (1)$$

where $\rho_{NIR}$ and $\rho_{RED}$ are the near-infrared (NIR) and red reflectances, respectively. NDVI increases with the amount of healthy green photosynthetically active vegetation. Uniformized NDVI measures reflect the “current” vegetation conditions according to a longer-term historical average NDVI value over the same study area.

#### 3.1.2. SWI

The soil water index (SWI) data are derived from the surface soil moisture (SSM) using an infiltration model to describe the relation between surface soil moisture and profile soil moisture as a function of time.

$$SWI(t) = \frac{\sum SSM(t_i)e^{-(t-t_i)/T}}{\sum e^{-(t-t_i)/T}}, \qquad (2)$$

where SSM is the surface soil moisture estimate from the ASCAT at time $t_i$. The parameter T, called the characteristic time length, represents the time period to integrate SSM data. T is the most important input parameter to derive the profile soil moisture content from remotely sensed surface time series: a high T value describes a deeper soil layer if the soil diffusivity is constant, and the same T value describes different depths for different soils (texture) [27]. In Ceballos et al. [55], the authors studied a semi-arid region characterized with a similar heterogeneous texture distribution in the soil profiles as our study region, having surface soils of sandy texture and a higher clay concentration to a depth of 2 m [56], and found that T = 60 days is the best choice for the 50–100 cm soil profile, which corresponds to the root layer according to the predominant agricultural use of our study region
(olive trees and cereals). Indeed, the upper part of the profiles promotes high infiltration rates and low water retention forces [55].

The SSM is retrieved from scatterometer observations from the ASCAT instruments on board the MetOP satellites, which measures radar backscatter at C-band in VV polarization, with an initial spatial resolution of 25 km, re-sampled at 12.5 km [57]. SWI product is retrieved from the Copernicus Global Land service and provides global daily information since January 2007 to present [58].

### 3.1.3. Thermal Infrared Stress Index

We used a stress index (SI) obtained from TIR. The surface temperature (LST) provides indirect estimates of water stress since it enters in the surface energy balance equations and it is thus related to the evapotranspiration rate flux. Daily ET simulation is detailed in Appendix A section. The stress index is defined as the ratio between actual ($ET_a$) and potential evaporation rates ($ET_p$) [40].

$$SI = \frac{ET_a}{ET_p} \qquad (3)$$

If the actual evaporation value is close to the potential value, the stress index takes values close to one that reflect unstressed conditions. However, if the actual evaporation is low compared to its potential value, the stress index values may reach zero, which represents fully stressed conditions. $ET_a$ and $ET_p$ are simulated from the energy balance model SPARSE [50] using meteorological observations and remote sensing variables.

LST, viewing angle and emissivity data are retrieved from the daily 1-kilometer LST product (MOD11A1) from MODIS VI products on board Terra and Aqua sensors. Quality control (QC) assessment for LST and emissivity product is also provided from the LST product. For stress index retrieval, the Leaf Area Index (LAI) and albedo products are also needed. We used the 8-day of albedo series (MCD43A3) from MODIS at 500 m, resampled to a spatial resolution of 1 km. Finally, the LAI is retrieved from the NDVI time series using an empirical equation [59], defined in the following Equation (4):

$$LAI = -\frac{1}{k} \ln \frac{NDVI_{max} - NDVI_{i}}{NDVI_{max} - NDVI_{min}} \qquad (4)$$

where k is an extinction factor, about 1.13, $NDVI_{max}$ = 0.97 and $NDVI_{min}$ = 0.05 corresponding to the NDVI value of bare soil [48]. The region is characterized by tree-dominant cover with interrow distance of 12 m corresponding to low LAI values (about 0.3 to 0.4) on an image [60].

### 3.1.4. UPI

We used the Climate Hazards group Infrared Precipitation with Stations (CHIRPS) data which are available from 1981 till now. This product incorporates three types of information: global climatologies, satellite estimates and in situ observations [61]. CHIRPS is available daily and at a spatial resolution of about 5 km. The CHIRPS precipitation data are highly reliable to identify wet and dry periods in our study region [62].

The CHIRPS precipitation time series is used to compute the Uniformized Precipitation Index (UPI), which is constructed following the same steps as the SPI except for the last step that transforms into Normal standard quantiles [7]. SPI maps precipitation amounts to the [0,1] interval by using either the empirical or a parametric cumulative function. Then, a final transformation is applied to obtain values that are normally distributed. UPI is computed with a single transformation: the empirical cumulative function is used to map precipitation amounts derived from historical records, to uniform quantiles in the [0,1] interval, as expressed in Equation (5):

$$\hat{F}_n(x) = \frac{1}{n} \sum_{i=1}^{n} 1_{\{x_i \le x\}}, \qquad (5)$$
where $1_{\{\cdot\}}$ is the indicator function, $x_1, \dots, x_n$ are rainfall amounts aggregated at different time scales (years, months or decades) and $\hat{F}_n(x)$ is the UPI value with $\hat{F}_n(x) = P(X \le x) \in [0, 1]$, with X the random variable representing the rainfall amounts at a given time scale. In contrast, SPI is computed as $\Phi^{-1}(\hat{F}_n(x))$ where $\Phi^{-1}$ is the quantile function (the inverse of the cumulative function) of the standard normal distribution.

We believe the interpretation in terms of uniform quantiles, which are simply probabilities of non-exceedance, is more straightforward than normal quantiles. For example, if UPI equals 0.1 on a given year, it means that 10% of the yearly precipitation amounts are inferior or equal to the amount observed that year. In other words, that year is among the 10% driest years. The SPI value corresponding to 0.1 is approximately -1.28.

## 3.2. Thermal Stress Index Derived from Energy Balance Model

### 3.2.1. SPARSE Model

We rely on the dual source energy balance model SPARSE [50] to simulate evapotranspiration. It estimates evaporation (E) and transpiration (T) separately. This is particularly relevant for arid and semi-arid areas which are characterized by sparse crop canopy and by an uneven relative contribution of evaporation and transpiration [63].

The model can be run in two modes: retrieval and prescribed mode. In retrieval mode, the respective stress levels (between non evaporating/transpiring and fully evaporating/-transpiring, i.e., potential rates) are two unknowns which are determined from the single piece of information provided by *Tsurf* by assuming at first that the vegetation is not stressed, that allows a simplification to solve the underdetermination problem. *Tsurf* is used to estimate the latent heat component from the soil (*LEs*), corresponding to the soil evaporation E. If the vegetation is suffering from water stress, the resulting *LEs* will decrease to unrealistic levels (negative values). In that case, we assume that the soil surface is stressed and *LEs* is set to a minimum value close to zero. Then the energy budget equation is solved for the vegetation component of the latent heat flux (*LEv*), which corresponds to the transpiration T. If *LEv* is also negative, fully stressed conditions are imposed for both soil and surface components [60]. The prescribed mode provides an estimate of the potential latent flux for the soil and the vegetation (*LEs*<sub>pot</sub> and *LEv*<sub>pot</sub> respectively). The water stress index SI is then defined from the actual and unstressed evapotranspiration rates (potential) at the time of the satellite overpass.

The SPARSE model is forced by a series of meteorological observations (air temperature, relative humidity, global radiation and wind speed) and remote sensing variables (NDVI, LAI, albedo and LST). Outputs are derived at meteorological time steps (half-hourly).

### 3.2.2. Meteorological Forcing

It is important to consider a long meteorological series in order to perform a robust statistical characterization of drought periods. In addition, owing to the temporal and spatial scales of climate variability, evapotranspiration and water stress index must be monitored at high temporal (subdaily to daily scales) and spatial resolution. However, available gauged stations (see Section 2) are scarce with short observation periods and numerous gaps, insufficient to perform a robust statistical analyses and characterization of drought periods. We rely, therefore, on surrogate meteorological series that extend in time the original series in order to be used to constrain the SPARSE energy balance model. In this work, we used two types of surrogate meteorological data:

1. Unprocessed reanalysis data ERA5 extracted at the grid cell closest to the region of interest (see Figure 1): ERA5 reanalyses are available at a 31 km spatial resolution [64]) from 1950 to present at an hourly temporal scale [64]. Reanalysis series that correspond to the four meteorological variables required for the energy balance model to simulate the corresponding index $SI_{ERA5}$ are: the incoming global solar radiation at the surface (bottom of atmosphere), wind speed at 10 m, air temperature at 2 m and the relative humidity that was derived from 2 m air temperature and 2 m dewpoint temperature
ERA5 products, according to the procedures defined in [65]. The specific aim of using unprocessed reanalysis data for our study, is to assess its performance to constrain an energy balance model for regions with no gauged stations.

2. Simulated series from a Stochastic Weather Generator (SWG) called "MetGen" [49]: Its implementation is publicly and freely available as an R library. MetGen generates scenarios of meteorological variables at sub-daily temporal resolution in order to extend local observations in the past. It relies on low resolution ERA5 reanalysis data and exploits observations provided by three gauged stations located in our study region (see Figure 1) to simulate regional climatic information. The corresponding index simulated using the SWG meteorological to constrain SPARSE model, is denoted SI<sub>SWG</sub>.

Inter-comparison between indices will be computed at a daily time scale for a long period in the past from 2000 to 2019. Characteristics of indices used are as follows, see Table 1.

<table><caption>**Table 1.** Characteristics of indices used.</caption><thead><tr><th>Wave Lengths</th><th>Visible/Near Infrared</th><th>Infrared + Visible + Meteo.</th><th>Microwave</th><th>Meteo. + Satellites Data</th></tr></thead><tbody><tr><td>Indices</td><td>NDVI</td><td>SI<sub>SWG</sub></td><td>SI<sub>ERA5</sub></td><td>SWI</td><td>UPI</td></tr><tr><td>Satellites</td><td>MODIS</td><td>MODIS</td><td>MODIS</td><td>ASCAT</td><td>CHIRPS</td></tr><tr><td>Model used</td><td>×</td><td>SPARSE</td><td>SPARSE</td><td>×</td><td>×</td></tr><tr><td>Spatial resolution</td><td>1 km</td><td>1 km</td><td>1 km</td><td>12.5 km</td><td>5 km</td></tr><tr><td>Temporal resolution</td><td>daily</td><td>daily</td><td>daily</td><td>daily</td><td>daily</td></tr><tr><td>Temporal availability</td><td>since 2000</td><td>since 2000</td><td>since 2000</td><td>since 2007</td><td>since 1981</td></tr></tbody></table>

### 3.3. Indices Standardization

For added interpretability, the different indices are usually re-scaled or standardized in order (1) to be expressed as frequencies that can be easily interpreted [66] and (2) to statistically characterize the deviation (anomaly) of the index values according to a long series recorded over the study period. For this aim, we standardize each index X with its empirical cumulative density function (ECDF). Let $(x_1, \dots, x_n)$ be the values taken by the index X for a given data set. Then the ECDF is given by:

$$ \hat{F}_n(x) = \frac{1}{n} \sum_{i=1}^{n} 1_{\{x_i \le x\}}, \qquad (6) $$

where $1_{\{\cdot\}}$ is the indicator function, $X_i$ is the index value for time step i and $Y_i$ is the standardized value of the index x. The ECDF standardizes by mapping a value x taken by X to its empirical non exceedance frequency, i.e., the percentage of occurrences of values lesser than or equal to x [67]. This is illustrated in Figure 3: the ECDF maps the values $x \in [0.1, 0.6]$ to $[0, 1]$. In particular, the value $x_i = 0.31$ is mapped to $y_i = 0.8$ which means that 80% of the values of the index X are expected to be below 0.31.

The advantages of the proposed standardization approach with the ECDF are that (1) it does not rely on the assumption that the raw indices values follow a normal distribution which is implicit in standardization approaches that rely on centering and scaling, (2) it allows for flexibility in the choice of reference sample (depending on seasons for example), and therefore can be used for short scale monitoring, and (3) it produces results that are easily interpreted in terms of non-exceedance probabilities, as for UPI, especially when simultaneous variables are analyzed [67]. Standardization using the empirical cumulative
distribution function is performed over the whole period (20 years). At an annual scale, standardization is performed over each year according to available values observed during the whole study period. At a decadal scale, standardization is performed according to season (winter, autumn, summer and spring).

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Empirical cumulative density function illustration.</figcaption></figure>

### 3.4. Evaluation of the Different Drought Indices Performance

#### 3.4.1. At Regional Scale

For comparison at a regional scale, all indices are aggregated at a very low spatial resolution (about 12 km). For this comparison, the CHIRPS precipitation time series is used as a reference to identify water stress conditions at different time scales: annual and decadal scales. In order to assess robustness and consistency between indices, we used for instance time series for visual interpretation and Kendall's τ, which is a rank-based correlation coefficient that measures the strength of the relationship between two indices.

To focus on consistency between the precipitation index and other indices, we assess the ability of the different indices to correctly identify the stress periods defined according to UPI. We define the stress periods based on UPI as follows. We assume that a drought event is a stress period with a probability of occurrence of approximately 20% and that it lasts at least two decades. The threshold of 20% in terms of UPI corresponds to a threshold of about -0.84 in terms of SPI. The drought classes defined in [7] begin when SPI values fall below 0 (this corresponds to 50% in terms of UPI). In our study region, the threshold of 20% corresponds to precipitation amounts inferior to 200 mm, which is approximately the precipitation amount average at a growth season. We transform the decadal standardized indices into Boolean values according to the stress threshold: 0 to denote an unstressed situation (>0.2) and 1 for stress conditions (<0.2). The second step is to check for each time step (t) how many stress periods are recorded within the period between (t - 1) and (t + 1), one decade being considered as a buffer (see Figure 4).

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Identification of stress periods using UPI index transformed into Boolean values.</figcaption></figure>

This method allows us to take into account the memory effects for three decades. Moreover, the use of a moving average of UPI allows us to take further consideration of memory effects and delayed responses of the different indices. We consider, therefore, the use of a UPI moving average lagged up to 2, 3 and 4 decades. If the number of stress periods based on UPI or lagged moving averages UPI exceeds two of the three decades
considered, we assign a drought condition for this time step. Then we assess the response stress/unstress of the other drought indices ($SI_{SWG}$, $SI_{ERA5}$, SWI and NDVI).

### 3.4.2. At Local Scale

For comparison at finer spatial resolution, a kilometric scale, we use local precise information based on:

* XLAS in-situ measurements: Sensible heat flux measurements using an extra-large aperture scintillometer (XLAS) are provided as part of the work of [60], for the period ranging between March 2013 and June 2015. The scintillometer (XLAS) was installed close to the Ben Salem village over a 4 km transect above a mixed vegetation canopy: trees (mainly olive orchards) with some annual crops (cereals and market gardening) [60]. For our analyses, pixels enclosed in the mean XLAS are selected in order to compare the different drought indices with the stress index derived from the sensible heat flux measurements, denoted $SI_{XLAS}$.

* Historical rainfed areas selection: Rainfed crops are more sensitive to rainfall depletion and thus to drought. For our analyses, we identify historical rainfed wheat areas relying on a non-irrigated cereal mask (see Figure 5a), computed for the agricultural year 2011–2012, as part of the work computed by [68]. It is computed using an object-oriented classification technique basing on the Spot image of 31 March 2012. We generate the percentage of non-irrigated cereal fields for this year, over each MODIS pixel, (see Figure 5b). Then, we select pixels that contain more than 40% of rainfed cereal cover. Rainfed cereal pixels selected are used as reference to locate non-irrigated cereal fields in precedent years, in order to assess the response of the different indices over a dry and a wet year.

In order to facilitate interpretation, we define four stress classes that characterize the severity of water stress identified using the different drought indices. The four stress classes are listed in Table 2, “High stress”, “Stress”, “Moderate stress” and “No stress” according to the values of a given index.

<table><caption>**Table 2.** Drought indices classification.</caption><thead><tr><td>Index Values</td><td>Drought Class</td></tr></thead><tbody><tr><td>0–0.15</td><td>High stress</td></tr><tr><td>0.15–0.3</td><td>Stress</td></tr><tr><td>0.3–0.7</td><td>Moderate stress</td></tr><tr><td>0.7–1</td><td>No stress</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 5.** *Cont.*</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 5.** Identification of rainfed cereals in MODIS pixels over the 2011–2012 agricultural season. (a) Rainfed cereal mask, produced for the 2011–2012 agricultural season [68]. (b) Percentages of rainfed cereal in each Modis pixel of the 2011–2012 agricultural season.</figcaption></figure>

# 4. Results

## 4.1. Drought Indices Inter-Comparison at Regional Scale

### 4.1.1. Annual Scale

Figure 6 presents the cumulative rainfall during the growth season (September to May) with the corresponding NDVI variation for each year (Figure 6a), in order to visually characterize the inter-seasonal variation over the study period. Rainfall data are provided from CHIRPS precipitation time series. The red line in Figure 6a corresponds to 200 mm/growth season of precipitation amount. We assume that a dry year presents a cumulative rainfall below 200 mm during a growth season. The Figure 6b shows the time series of different standardized drought indices: the TIR SI simulated when the SPARSE model is constrained by the SWG or ERA5 meteorological data *SI*<sub>SWG</sub> and *SI*<sub>ERA5</sub>, the standardized NDVI index, the UPI, and the standardized SWI which is available since 2007. The red line in Figure 6b indicates the drought threshold, fixed in this study at 0.2 which corresponds to a drought event with an approximately 20% non exceedance frequency. The stress threshold is necessary to discuss drought events, as well as their onset, duration and intensity. The farther below the red line a standardized index is, the worse the drought is.

Overall, we see that annual indices show a similar variation. *SI*<sub>SWG</sub> and *SI*<sub>ERA5</sub> show similar inter-annual variability. Moreover, we observe that they reproduce UPI variability and succeed in identifying wet and dry conditions simultaneously. Dry periods, which correspond to abnormally dry situations that fall below the drought threshold, are well depicted also by a cumulative rainfall less than 200 mm during each growing season, and a low NDVI level. However, we observe in some instances that UPI is not very well correlated with other standardized indices. This is the case for 2004–2005 or 2008–2009 where we observe that UPI is either above the unusual unstressed condition or below the drought threshold while other indices present similar stress level. To better explore the allocation of precipitation intensity along these years, we present in Figure 7, the monthly cumulative precipitation and monthly NDVI variation over these agronomic years. In addition, in
order to facilitate interpretation, we present in Figure 7c, the monthly rainfall and NDVI average over the whole study period between 2000 and 2019. There are in Figure 7a, low cumulative rainfalls in winter and spring in comparison with the monthly rainfall amount mostly observed in our study region during this period (see Figure 7c). However, almost all the largest cumulative rainfall amounts are observed from September to December. Rainfall during this period was sufficient, and coincided with the beginning of the growth season. For this reason, we observe a high level of NDVI in March and May that reflect a good development of vegetation, particularly cereal crops. Consequently, unstressed conditions were also shown by the rest of the standardized indices. We also observe an abnormally unstressed situation depicted by UPI during the agronomic year 2008–2009. However, the other standardized indices indicate a medium stress condition. Indeed, we observe very low rainfall amounts during the whole year, especially in November and December that correspond to seeding and early growth. Then, we observe an abnormally high rainfall amount that exceeds 100 mm/month in January. This high cumulative amount increases the annual cumulative rainfall but it is late for vegetation growth. Indeed, an adequate response of the stress is always dependent on receiving sufficient rainfall at the appropriate moments, with sufficient amounts. In contrast, the agronomic year 2017–2018 shows very low cumulative seasonal rainfall. The TIR SI values succeed to detect the stress period, by presenting low values. However, standardized SWI and NDVI show a medium stress condition. SWI and NDVI did not succeed to detect the drought intensity for this year.

<figure><img src="image_12.png" /><figcaption>**Figure 6.** Time-series of seasonal variation of standardized drought indices accumulated in growth season. (a) Seasonal cumulative rainfall (September to May) with the corresponding NDVI variation. The red line indicates the 200 mm precipitation amount. (b) Drought indices.</figcaption></figure>

To assess the consistency between the various indices at an annual scale, we use Kendall's τ. Kendall's τ is suitable for non-Gaussian distributions, as opposed to the Pearson correlation coefficient [69]. Positive values indicate that both indices tend to increase or decrease simultaneously, while negative values indicate that they tend to vary in an opposite manner. A value near zero signals a lack of dependence. In Figure 8, we see
that all indices present positive correlation coefficients. All indices present high Kendall's coefficients with UPI that exceed 0.35, particularly SI<sub>ERA5</sub> which presents a Kendall's τ of about 0.46. However, SWI presents a very low correlation coefficient with UPI. This is could be explained by the different responses of these indices to a rainfall supply. UPI presents an instantaneous response to rainfall variation. In contrast, SWI which presents the root zone soil moisture, presents a delayed response to rainfall supply. On the other hand, we observe that SWI presents the highest τ value with NDVI, which reflects a very strong seasonal correlation between them. Furthermore, SI<sub>SWG</sub> is highly correlated with other indices. It presents a Kendall's τ of about 0.58 with NDVI and with SWI and the highest correlation with SI<sub>ERA5</sub>, with a Kendall's τ of about 0.64. However, SI<sub>ERA5</sub> presents low correlation coefficients with NDVI and SWI.

<figure><img src="image_13.png" /><figcaption>**Figure 7.** Precipitation and NDVI monthly variation over agronomic years (2004–2005) and (2008–2009). (a) 2004–2005. (b) 2008–2009. (c) Monthly rainfall and NDVI over the study period 2000–2019.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 8.** Seasonal Kendall's τ correlation between drought indices.</figcaption></figure>

### 4.1.2. Decadal Scale

We focus on a finer temporal scale, a decadal scale, in order to assess short term variations. Indeed, a short time scale (monthly or decadal) is more appropriate to provide information on the growing season development over a few months and, therefore, on the water supply, unlike the annual time scale [6].

In Figure 9, we propose a comparison analysis carried out at a decadal scale of the cereal growing season (from September to May) selected from some dry and wet years. We used the values of precipitation index UPI lagged up to five decades rather than the instantaneous value, in order to be more correlated with the aggregated response of the different indices. The agronomic year 2012–2013 (see Figure 9a) is considered as a dry year, based on the cumulative rainfall over the corresponding growth season (Figure 6a). We observe that the standardized NDVI index presents moderate drought conditions (50% of probability of occurrence) at the beginning. Then, it decreases in December to 0.4 and then it reaches 0.2 in the end of the growth season. However, the standardized TIR indices show mainly similar variations, more correlated with lagged UPI variations. UPI shows considerable fluctuations along the whole period. We observe very low values during several decades in September and October and in January until March, that lead to important drops of SI<sub>SWG</sub>, SI<sub>ERA5</sub> and extended low values of SWI which fluctuate around 0.2 in the beginning of the growth season.

The period 2013-2014 represents a typical wet growth season. In Figure 9b, we observe that the different drought indices are above the drought threshold over the whole period. The standardized NDVI index presents high values in spring that reflect a good vegetation health and a maximum crop development during this period. High NDVI values are observed after an important water supply provided during the growing season. UPI presents abnormally high values that reach 100% probability of non exceedance for several decades. SI<sub>SWG</sub>, SI<sub>ERA5</sub> succeed to reproduce the main fluctuation in the lagged UPI: relevant increases in UPI values lead to an important elevation of both the TIR indices and a relevant drop in UPI which falls below the stress threshold mainly in January and February and induces an important rapid decline of both TIR SI values. SWI presents, however, a more delayed response to UPI variations. In fact, we observe a depression in SWI in autumn and spring due to successive drop fluctuations of UPI during these periods.

In Figure 9c, we present another dry growth season. The lack of precipitation observed mainly at the beginning of the growth season induces a decrease in the standardized NDVI and SWI indices at the end of the growing season. The rainfall shortage during this season triggers a relevant depletion of the TIR indices, especially in winter where we observe SI<sub>SWG</sub>, SI<sub>ERA5</sub> and UPI are below the stress threshold. However, each lagged UPI value increase leads to an elevation in TIR indices.
<figure><img src="image_15.png" /><figcaption>**Figure 9.** Decadal times series of standardized stress indices derived from different wavelengths over the growing season of two dry years (2012–2013 and 2017–2018), and a wet year (2013–2014). (a) 2012–2013. (b) 2013–2014. (c) 2017–2018.</figcaption></figure>

A further analysis is carried out to assess more closely the consistency between the different drought indices in identifying stress periods at the decadal timescale according to UPI. As explained in Section 3.4.1, we define stress periods according to the different indices. Besides, we identify three cases: (1) the stress indices correctly reproduce the drought condition identified with UPI (stress period or not stress period), flagged as “Congruent with UPI”, (2) the stress indices identify a stress period not confirmed by UPI, called a “Not confirmed by UPI”, and the last case, (3) the stress indices do not succeed to identify a stress period detected by the UPI, called a “missed stress”. The first case estimates the consistency between the different drought indices and UPI to identify the same stress conditions. Analyses are performed according to seasons, in order to assess drought index responses through different seasons.

In Figure 10, we present the probability to identify these three different cases, according to seasons and UPI values in Figure 10a and lagged moving average UPI to 2, 3 and 4 decades, respectively, for Figure 10b–d. Overall, we observe that the different indices show a good performance to identify the same stress condition in autumn, winter and especially in spring. However, in summer, the probability to detect the same stress conditions using the different stress indices is decreased at the expense of missed stress. Indeed, in summer the probability of rainfall events is extremely low in our study region. Therefore, drought periods identified by UPI increase, but they are not necessarily identified by other drought indices due to the supply of water provided from irrigation as an alternative water supply for this period. The NDVI shows high probability to detect stress status congruent with UPI, especially in spring where we have the maximum of vegetation development. The case of “congruent with UPI” could be related to a stressed or an unstressed period. In this
case, the high level of the “congruent with UPI” performed by NDVI is mainly related to a
no stress period over the rainy season and particularly in spring.

<figure><img src="image_16.png" /><figcaption>**Figure 10.** Probability of drought detection based on UPI and lagged UPI during the period between 2007 and 2019. (a) Stress identification basing on UPI. (b) Stress identification based on a UPI 2 decade moving average. (c) Stress identification basing on UPI 3 decade moving average. (d) Stress identification basing on UPI 4 decade moving average.</figcaption></figure>
Besides, we observe that *SI*<sub>SWG</sub> and *SI*<sub>ERA5</sub> show a good performance to correctly detect the stress conditions determined from the UPI and it presents a low probability of identifying stress periods not confirmed by UPI, especially using *SI*<sub>SWG</sub>. This performance is maintained during different seasons and even enhanced with the comparison with the lagged moving average of the UPI. On the other hand, we observe that SWI shows a better performance in winter and spring. Indeed, during this season, vegetation is well developed with deeper root zones, which could explain the high performance of SWI to identify a stress period or unstressed period according to UPI. Furthermore, this performance is improved incrementally when we increase the lagged moving average value of UPI. The application of a moving average of 4 decades seems to be the best delay period for adequate responses for all drought indices. Beyond this value, too many memory effects are lost.

## 4.2. Drought Indices Inter-Comparison at Local Scale

### 4.2.1. Evaluation with XLAS In-Situ Measurements

In Figure 11, we present a decadal time series of the different standardized indices during two growing seasons, where XLAS measurements are available at both local (Figure 11a,c) and regional scale (Figure 11b,d). As mentioned above, we used for this comparison the lagged UPI (up to 5 decades) and *SI*<sub>XLAS</sub>, respectively, as regional and local references. The aim of the comparison of indices at a local and a regional scale is to assess if an index found to be robust at low resolution (12 km) maintains its performance at higher resolution (kilometric). For interpretation, we consider two groups of indices according to their spatial resolution: a group with low resolution (kilometric) (*SI*<sub>SWG</sub>, *SI*<sub>ERA5</sub> and NDVI), and a second group with very low spatial resolution, at dozens of kilometers (SWI and UPI).

In Figure 11a,b, we present a wet year (2013–2014). We see that indices *SI*<sub>SWG</sub> and *SI*<sub>ERA5</sub> succeed to better reproduce *SI*<sub>XLAS</sub> fluctuations in Figure 11a. In addition, we observe that these indices, especially *SI*<sub>SWG</sub>, follow the *SI*<sub>XLAS</sub> variations more closely, particularly during the rainfall season (from October to February) and the last decades of this growth season. Besides, we observe that stress periods defined by the *SI*<sub>XLAS</sub> (under the threshold line), observed in October, December and February, are mostly isolated with *SI*<sub>SWG</sub>. Moreover, *SI*<sub>SWG</sub> succeed in identifying stress conditions defined by *SI*<sub>XLAS</sub>, even using its aggregated values at the regional scale (see Figure 11b). The NDVI, in spite of being defined among higher resolution indices, did not succeed to reproduce *SI*<sub>XLAS</sub> decadal variations. However, it brings, as expected, an overall delayed information about stress conditions. We observe, for example, in Figure 11a, which presents a wet year, an elevation in NDVI values in winter and spring that also correspond to several elevated values observed by *SI*<sub>XLAS</sub>. On the other hand, we observe a different behavior of the standardized SWI according to the different spatial scale. At a local scale (Figure 11a), SWI did not succeed in correctly reproducing *SI*<sub>XLAS</sub> variations. However, at a regional scale (Figure 11b), SWI succeeded in identifying dry or wet conditions, according to the lagged UPI fluctuations. Dry conditions were observed at the beginning of the growth season and in February. Wet conditions were observed in November–December and the end of the growth season.

The 2014–2015 period represents a relatively dry year. We observe that NDVI presents mainly moderate stress values in Figure 11d. In Figure 11c, we observe that TIR stress indices succeed to better reproduce *SI*<sub>XLAS</sub> in this growth season, particularly using *SI*<sub>SWG</sub>. Indeed, *SI*<sub>SWG</sub> succeed to identify a water stress period caught by *SI*<sub>XLAS</sub> which persists for several decades during March, as well as an unstressed period observed later in April and May, visible in the same figure. At a local scale (Figure 11c), SWI shows low values observed below the stress threshold during several decades in Autumn and March, which amply testify an overall dry situation during this growth season. Even at a regional scale (Figure 11d), SWI shows moderate stress values. However, SWI was not able to accurately identify stress conditions in comparison with in situ measurements. This could be explained by its very low spatial resolution, inadequate for local stress detection.
<figure><img src="image_17.png" /><figcaption>**Figure 11.** Decadal times series of standardized stress indices derived from different wavelengths over the growing season of two years (2013–2014) and (2014–2015) at a local scale (a,c), in comparison with *SI*<sub>XLAS</sub> derived from XLAS measurements over pixels containing the scintillometer transect and at a regional scale (b,d). (a) Local scale (2013–2014). (b) Regional scale (2013–2014). (c) Local scale (2014–2015). (d) Regional scale (2014–2015).</figcaption></figure>

4.2.2. Drought Detection in Rainfed Areas

In Figure 12, we present stress class maps of the different indices at the seasonal scale (growth season) of the year 2011–2012: Figure 12a for the standardized thermal index constrained by meteorological SWG data, Figure 12b for the standardized thermal index constrained by meteorological ERA5 data, Figure 12c for standardized NDVI index, Figure 12d for standardized SWI index and Figure 12e for standardized precipitation index. In addition, we present the percentage of stress classes identified by each index using extracted rainfed cereal pixels, as explained in Section 3.4.2 (see Figure 12f). Combined results presented in Figure 12f show that most indices depicted a moderate stress class for this year with the highest percentage. However, we observe that the TIR stress indices and SWI show an elevated percentages of stress class and even a very low percentage of a
high stress class for SI<sub>SWG</sub> and SI<sub>ERA5</sub>. This finding is in agreement with the hydrological year characteristics presented in Figure 6a. Indeed, this year characterized by a cumulative seasonal rainfall about 228 mm did not belong to very wet or very dry year. Stress classes depicted by thermal indices and SWI could be more precise than NDVI characterized by a delayed response and than the UPI which is less informative when aggregated at a low temporal scale (seasonal scale).

<figure><img src="image_18.png" /><figcaption>**Figure 12.** Stress classes for each drought index used over the agronomic year 2011–2012, (a) for the standardized thermal index constrained by meteorological SWG data, (b) for the standardized thermal index constrained by meteorological ERA5 data, (c) for standardized NDVI index (d) for standardized SWI index, (e) for standardized precipitation index and (f) for percentage of the stress class identified by each index.</figcaption></figure>

Then, the same analysis was carried out for a dry year, 2009–2010. In Figure 13, we observe that the rainfed cereal areas selected correspond mainly to high stress or stress classes. More than 80% of rainfed cereal pixels indicate a high stress according to SI<sub>SWG</sub> and SI<sub>ERA5</sub> maps. Using the SWI information, we also identify a high stress of more than 50%. However, using the NDVI, we mainly identify stress class and moderate stress, as well as the response of the UPI. Expected results are obtained mostly by the TIR indices and SWI. In fact, these indices are related directly with vegetation water requirement, especially
for the TIR that also depends on atmospheric demand. Precise information concerning vegetation and surface stress should be obtained from these indices. Due to its low spatial resolution, SWI might be less accurate for local analyses.

<figure><img src="image_19.png" /><figcaption>**Figure 13.** Stress classes for each drought index used over the agronomic year 2009–2010, (a) for the standardized thermal index constrained by meteorological SWG data, (b) for the standardized thermal index constrained by meteorological ERA5 data, (c) for standardized NDVI index (d) for standardized SWI index, (e) for standardized precipitation index and (f) for percentage of the stress class identified by each index.</figcaption></figure>

**5. Discussion**

We focused on growing seasons to perform our analysis (from September to May). Indeed, droughts could be damaging for crop health in these periods of the year. Analyses were carried out at different spatial scales: regional and local scale. At the regional scale, UPI is considered as the reference index to characterize the year-to-year water statement, and to assess the response of the aggregated indices used at a large scale (12 km). Overall, we find that drought indices from the different wave lengths show mostly the same variations at the annual time scale. All the indices succeed to reproduce the UPI variations. We also find that thermal stress indices provided by the SPARSE model, when constrained either from SWG or from ERA5 meteorological data, SI<sub>SWG</sub> and SI<sub>ERA5</sub>, respectively, are very well
correlated with UPI. However, *SI*<sub>ERA5</sub> has a higher Kendall's τ with UPI. Indeed, *SI*<sub>ERA5</sub> is simulated using regional climatic information (ERA5 reanalyses) and local remote sensing information. The use of the large scale meteorological data to simulate *SI*<sub>ERA5</sub> could explain its high correlation with the UPI also provided by gridded data-sets. Moreover, we observe a very strong seasonal correlation between *SI*<sub>SWG</sub> and *SI*<sub>ERA5</sub>. This correlation, measured by Kendall's coefficient, was also conserved at the decadal scale (see Figure 9), due to their high spatial resolution. SWI and NDVI also present high correlation due to their aggregated response to a stress condition. Indeed, the NDVI presents a delayed response in comparison with the other indices, and shows less variability because of the time lag needed for vegetation to respond for a supply or a lack of water. On the other hand, SWI, due to the time needed for the water supply to be infiltrated to the root zone soil, also shows a delayed response to water supply in comparison to the other indices.

We assess the consistency between drought indices and UPI as a reference stress index for regional analyses. We find that there is more consistency between UPI and SI<sub>SWG</sub> to identify the stress condition, especially in autumn, winter and spring. Moreover, this index presents the smallest probability to identify a stress period not confirmed by UPI. SWI shows more consistency with UPI in spring. Indeed, during this season, vegetation presents deeper root zones, which could explain the high performance of this index to identify a stress period or unstressed period according to UPI. Our analyses show also that using a moving average of UPI improves the adequation of drought indices to the hydric status. In fact, all indices need a time lag to respond to a water supply or shortage. This temporal lag depends obviously on the index used.

We can conclude, therefore, that TIR stress indices are more sensitive for a drought period. Indeed, a shortage of water supply will affect the soil (evaporation) and the vegetation (root zone), inducing stomatal closure and a surface temperature elevation detected by thermal sensors. However, SWI, which carries information on the soil moisture in the root zone, did not show immediate response as a result of a lack of precipitation. Indeed, the availability of some water storage in the profile soil after surface soil evaporation could delay the response of SWI to a water stress. For this aim, we need further analyses to identify the optimal parameter for the time infiltration (T) according to the characteristics of our study region. Besides, satellite measurements to estimate soil moisture may not be sufficiently precise because of the low spatial resolution of the SSM product used (see Section 3.1.2). However, *SI*<sub>ERA5</sub>, which is derived from very low resolution unprocessed ERA5 reanalyses, shows a good performance in identifying water stress conditions, since it is simulated at kilometric scale with precise satellite information. Analyses at the local scale strengthen our findings. Indeed, comparison of the different drought indices with the index derived from the XLAS measurements over the scintillometer transect show that the thermal stress indices are promising tools for drought identification at low spatial resolution scale (kilometric). In this case, *SI*<sub>SWG</sub>, which incorporates information from the local meteorological station close to the XLAS transect, shows better performances than *SI*<sub>ERA5</sub> and allows an accurate identification of the stress condition. Moreover, the identification of drought classes in non-irrigated cereal pixels shows, therefore, an ambiguous situation corresponding to *SI*<sub>ERA5</sub>. This could be explained by the large spatial resolution of meteorological data (31 km [64]) used to constrain the SPARSE model in order to simulate *SI*<sub>ERA5</sub>. For this task, although the SPARSE model allows the retrieval of separate estimates of evaporation and transpiration, vegetation analyses seem to be insufficient to identify plant stress periods, owing to the sparsity of the vegetation cover and the lack of consistent information about rainfed areas.

## 6. Conclusions

Semi-arid areas are characterized by their high exposure to extreme climatic variability. The occurrence of hot temperatures along with a deficit of rainfall leads to droughts and impacts agricultural production. These trends in drought occurrence could be quantified by a long time series of historical indicators. We rely on RS drought indicators provided by
different wave lengths: the visible/near infrared, the thermal infrared and the microwave
domains. These indicators provide information about vegetation health, water require-
ment and soil moisture that can help us to identify the plant hydric status. We use the
SPARSE model, a dual source energy balance model [50], in order to retrieve estimates of
evapotranspiration and water stress indices from the thermal infrared domain. SPARSE
relies on satellite information and meteorological observation series to characterize vege-
tation cover and atmospheric demand. As far as we know, not many studies of drought
index comparison have been proposed in the literature. Our study shows that NDVI is
very informative on the year-to-year water stress conditions. However, in spite of its low
resolution (kilometric), its delayed response to water stress forms a major disadvantage.
Thermal indices and SWI show consistent information at an annual scale. Nevertheless,
SI<sub>SWG</sub> seems to be more precise for relevant water stress identification, especially at a local
scale. Indeed, SWI is not sufficiently reliable to accurately identify stress intensity at finer
time scales due to it is large spatial scale. On the other hand, this work highlights the
performance of large scale meteorological variables (reanalyses in our case) to identify
periods of droughts. We proposed an efficient alternative when local meteorological data
are not available, particularly at regional analyses. Perspectives for this work include
the simulation of a drought index simulated from transpiration only in order to focus on
agronomic droughts. Indeed, information about vegetation water requirements rather than
information provided from the "vegetation + soil" composite could be an efficient tool
for drought management in semi-arid regions characterized by crop disparity. Moreover,
adding mesoscale numerical weather prediction could provide efficient drought prediction.
Lastly, the SWG-SPARSE tool will be tested for similar applications but in different climatic
conditions (e.g., coastal areas).

**Author Contributions:** writing original draft preparation, N.F.; data processing, N.F. and M.L.P.; Software, N.F., M.L.P. and J.C.; data analysis, N.F. and G.B.; results interpretation N.F., G.B., Z.L.C. and J.C.; Visualisation, G.B. and M.L.P.; writing review and editing, G.B., M.L.P. and J.C.; ideas and discussions, J.C., M.L.P. and G.B.; data curation, Z.L.C., Z.K. and G.B.; funding acquisition, J.C., Z.L.C. and G.B.; Conceptualisation, N.F. and G.B.; methodology, G.B. and N.F.; supervision, G.B.; validation, G.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work has been supported by the NAILA International Laboratory through a PhD grant, as well as the PHC Utique French-Tunisian bilateral "Amande" project, the H2020 RISE mobility "ACCWA" project, the WaterWorks2017 "FLUXMED", PRIMA 2018 "ALTOS" and ERANETMED "CHAAMS" ERA-NET co-fund projects.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are available in Environment Information System of CESBIO site: https://sie.cesbio.omp.eu (generated on 11 October 2021).

**Acknowledgments:** We extend our warm thanks to all the technical teams of the IRD, GREEN-TEAM, CESBIO and INAT for their strong collaboration and support.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Appendix A. Post-Processing of Instantaneous Evapotranspiration Estimates**

Evapotranspiration is simulated at the satellite overpass time and only during days
with clear sky. Some post-processing must be performed to reconstruct seasonal evapo-
transpiration useful for hydrological studies and particularly in agricultural applications.
The post-processing involves three steps: fusion, extrapolation and interpolation, that are
described below.
## Appendix A.1. Extrapolation

ET derived from TIR RS data rely on once-a-day acquisitions. An extrapolation algorithm is used to reconstruct its sub-daily variations using a method based on the evaporative fraction (EF), well described in [70]. EF at the satellite overpass t is defined as the ratio between the latent heat flux (LE) and the available energy (AE), which is derived from the net radiation (Rn) and soil heat flux (G), as defined in the following equation:

$$EF_{obs} = \frac{LE}{AE_t} = \frac{LE}{Rn - G} \qquad (A1)$$

According to [71], EF is relatively stable during the daylight hours on clear days. This hypothesis allows the reconstruction of the diurnal behavior (hourly or half-hourly time scale), see Equation (A2), based on the empirical equation of EF defined by [72] (denoted $EF_{sim}$), see Equation (A3), and observed EF at satellite overpass (Equation (A1)).

$$EF_{(30 \text{ min})} = EF_{sim}(30 \text{ min}) * \frac{EF_{obs}}{EF_{sim}} \qquad (A2)$$

where,

$$EF_{sim} = [1.2 - (0.4 * \frac{GR_t}{1000} + 0.5 * \frac{RH_t}{100})] \qquad (A3)$$

GR and RH are, respectively, the global radiation and relative humidity for a time scale t (semi-hourly or satellite overpass).

Finally, daily evaporation ($E_d$) can be simply obtained using the Equation (A4):

$$E_d = EF * AE_d, \qquad (A4)$$

where $AE_d$ is computed from an equation suggested by [73]: where AE is assumed to present the same diurnal variation as the global radiation, see (A5).

$$AE_d = R_g d * \frac{AE_t}{R_g t}, \qquad (A5)$$

## Appendix A.2. TERRA ET and AQUA ET Merge

ET simulations are derived from TERRA and AQUA satellites. We choose to combine these two sources of data in order to increase the overall temporal availability. Results presented by [60] show that ET derived from these different sensors is potentially biased because of their different overpass time (10:30 a.m. and 13:30 a.m. for Terra and Aqua, respectively) and that TERRA simulations are more correlated with ET derived from measured sensible heat flux (H) from an extra large aperture scintillometer (XLAS). Therefore, in order to take into consideration systematic differences in distributional properties, we perform a bias correction approach for the AQUA simulations (available since 2002) before merging. We considered CDF-t [74], a univariate bias correction method which allows non-linear corrections to reproduce the statistical distributions of the reference time series. Our aim is to keep ET simulations derived from TERRA RS (available since 2000) and to add bias corrected ET derived from AQUA only in days with no TERRA acquisition.

## Appendix A.3. Interpolation

to derive seasonal evapotranspiration, we have to fill the gaps between satellite acquisitions in order to reconstruct days with missing ET data. Methods that rely on self preservation or known diurnal shape of the ratio of evapotranspiration to a scale factor are usually used for this aim [70]. Indeed, this scale factor will monitor interpolation between two acquisition dates of successive images. In this work, we used the global radiation as an interpolation reference quantity. Delogu et al. [41] have shown that it is the most robust and best performing scale factor for seasonal timescales.
## Appendix B. Daily ET Simulations

Figure A1 shows time series of the ET estimates derived from merging ET Terra and ET Aqua over the agronomic year 2013–2014 (from September to August), for which the residual ET derived from scintillometer (XLAS) measurements are available. We present ET simulations separately when the SPARSE model is constrained by the SWG meteorological data (see Figure A1a) and when the SPARSE model is constrained by the ERA5 reanalyses meteorological data (see Figure A1b). Overall, we observe in both figures that ET constrained from Aqua RS data (ET Aqua) is always biased in comparison with ET constrained from Terra RS data (ET Terra), and that ET Terra present higher values. In Figure A1a, we observe that, in almost all cases, ET derived from XLAS measurements shows similar variation to ET Terra, excepted for summer periods. In addition, ET provided by the merging of Terra and Aqua reproduces more closely ET XLAS and succeeds in removing biases observed originally between ET Terra and Aqua. Using unprocessed ERA5 reanalyses (Figure A1b), ET from merging succeeds as well to reduce biases observed initially between ET Terra and Aqua. However, the three ET series simulated from ERA5 meteorological data are not as close to the ET XLAS variations as the ET series simulated from the SWG meteorological data. Indeed, ERA5 series are always biased in comparison with observations data.

<figure><img src="image_20.png" /></figure>

**Figure A1.** Weekly times series of evapotranspiration simulated initially from TERRA and AQUA satellite data superposed with a time series computed from the merging (Merge) of TERRA and bias corrected AQUA series and ET from XLAS measurements for the agronomic year 2013–2014. (a) [ET constrained by SWG meteorological surrogate series. (b) ET constrained by ERA5 meteorological surrogate series.

## References

1. Wilhite, D.A.; Svoboda, M.D. Drought early warning systems in the context of drought preparedness and mitigation. In *Early Warning Systems for Drought Preparedness and Drought Management*; World Meteorological Organization: Geneva, Switzerland, 2000; pp. 1–21.

2. Schilling, J.; Hertig, E.; Tramblay, Y.; Scheffran, J. Climate change vulnerability, water resources and social implications in North Africa. *Reg. Environ. Chang.* **2020**, *20*, 1–12. [CrossRef]

3. MassonDelmotte, V.; Zhai, P.; Pirani, A.; Connors, S.L.; Péan, C.; Berger, S.; Caud, N.; Chen, Y.; Goldfarb, L.; Gomis, M.I.; et al. (Eds.) Summary for Policymakers. In *Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change*; Cambridge University Press: Cambridge, UK, 2021.

4. Anderson, M.C.; Hain, C.; Wardlow, B.; Pimstein, A.; Mecikalski, J.R.; Kustas, W.P. Evaluation of drought indices based on thermal remote sensing of evapotranspiration over the continental United States. *J. Clim.* **2011**, *24*, 2025–2044. [CrossRef]

5. Otkin, J.A.; Anderson, M.C.; Hain, C.; Mladenova, I.E.; Basara, J.B.; Svoboda, M. Examining rapid onset drought development using the thermal infrared-based evaporative stress index. *J. Hydrometeorol.* **2013**, *14*, 1057–1074. [CrossRef]
6. Mishra, A.K.; Singh, V.P. A review of drought concepts. *J. Hydrol.* **2010**, *391*, 202–216. [CrossRef]

7. McKee, T.B.; Doesken, N.J.; Kleist, J. The relationship of drought frequency and duration to time scales. In Proceedings of the 8th Conference on Applied Climatology, Boston, MA, USA, 17-22 January 1993; Volume 17, pp. 179–183.

8. Mishra, A.; Desai, V. Drought forecasting using stochastic models. *Stoch. Environ. Res. Risk Assess.* **2005**, *19*, 326–339. [CrossRef]

9. Mishra, A.; Singh, V.P. Analysis of drought severity-area-frequency curves using a general circulation model and scenario uncertainty. *J. Geophys. Res. Atmos.* **2009**, *114*, D06120. [CrossRef]

10. Zhang, J.; Xu, Y.; Yao, F.; Wang, P.; Guo, W.; Li, L.; Yang, L. Advances in estimation methods of vegetation water content based on optical remote sensing techniques. *Sci. China Technol. Sci.* **2010**, *53*, 1159–1167. [CrossRef]

11. AghaKouchak, A. Advancements in Satellite Remote Sensing for Drought Monitoring. In *Drought and Water Crises*; CRC Press: Boca Raton, FL, USA, 2017; pp. 225–258.

12. Dorigo, W.; Wagner, W.; Albergel, C.; Albrecht, F.; Balsamo, G.; Brocca, L.; Chung, D.; Ertl, M.; Forkel, M.; Gruber, A.; et al. ESA CCI Soil Moisture for improved Earth system understanding: State-of-the art and future directions. *Remote Sens. Environ.* **2017**, *203*, 185–215. [CrossRef]

13. Wigneron, J.P.; Chanzy, A.; Calvet, J.C.; Bruguier, N. A simple algorithm to retrieve soil moisture and vegetation biomass using passive microwave measurements over crop fields. *Remote Sens. Environ.* **1995**, *51*, 331–341. [CrossRef]

14. Zribi, M.; Nativel, S.; Le Page, M. Analysis of Agronomic Drought in a Highly Anthropogenic Context Based on Satellite Monitoring of Vegetation and Soil Moisture. *Remote Sens.* **2021**, *13*, 2698. [CrossRef]

15. Amri, R.; Zribi, M.; Lili-Chabaane, Z.; Wagner, W.; Hasenauer, S. Analysis of C-band scatterometer moisture estimations derived over a semiarid region. *IEEE Trans. Geosci. Remote Sens.* **2012**, *50*, 2630–2638. [CrossRef]

16. Zhang, D.; Tang, R.; Zhao, W.; Tang, B.; Wu, H.; Shao, K.; Li, Z.L. Surface soil water content estimation from thermal remote sensing based on the temporal variation of land surface temperature. *Remote Sens.* **2014**, *6*, 3170–3187. [CrossRef]

17. Jiao, W.; Wang, L.; McCabe, M.F. Multi-sensor remote sensing for drought characterization: Current status, opportunities and a roadmap for the future. *Remote Sens. Environ.* **2021**, *256*, 112313. [CrossRef]

18. Tucker, C.J. Red and photographic infrared linear combinations for monitoring vegetation. *Remote Sens. Environ.* **1979**, *8*, 127–150. [CrossRef]

19. Kogan, F.N. Application of vegetation index and brightness temperature for drought detection. *Adv. Space Res.* **1995**, *15*, 91–100. [CrossRef]

20. Quiring, S.M.; Ganesh, S. Evaluating the utility of the Vegetation Condition Index (VCI) for monitoring meteorological drought in Texas. *Agric. For. Meteorol.* **2010**, *150*, 330–339. [CrossRef]

21. Kogan, F.N. Operational space technology for global vegetation assessment. *Bull. Am. Meteorol. Soc.* **2001**, *82*, 1949–1964. [CrossRef]

22. Amri, R.; Zribi, M.; Lili-Chabaane, Z.; Duchemin, B.; Gruhier, C.; Chehbouni, A. Analysis of vegetation behavior in a North African semi-arid region, using SPOT-VEGETATION NDVI data. *Remote Sens.* **2011**, *3*, 2568–2590. [CrossRef]

23. Farrar, T.; Nicholson, S.; Lare, A. The influence of soil type on the relationships between NDVI, rainfall, and soil moisture in semiarid Botswana. II. NDVI response to soil moisture. *Remote Sens. Environ.* **1994**, *50*, 121–133. [CrossRef]

24. Ozelkan, E.; Chen, G.; Ustundag, B.B. Multiscale object-based drought monitoring and comparison in rainfed and irrigated agriculture from Landsat 8 OLI imagery. *Int. J. Appl. Earth Obs. Geoinf.* **2016**, *44*, 159–170. [CrossRef]

25. AghaKouchak, A.; Farahmand, A.; Melton, F.; Teixeira, J.; Anderson, M.; Wardlow, B.D.; Hain, C. Remote sensing of drought: Progress, challenges and opportunities. *Rev. Geophys.* **2015**, *53*, 452–480. [CrossRef]

26. Paulik, C.; Dorigo, W.; Wagner, W.; Kidd, R. Validation of the ASCAT Soil Water Index using in situ data from the International Soil Moisture Network. *Int. J. Appl. Earth Obs. Geoinf.* **2014**, *30*, 1–8. [CrossRef]

27. Paulik, C.; Naeimi, V.; Dorigo, W.; Wagner, W.; Kidd, R. A global validation of the ASCAT Soil Water Index (SWI) with in situ data from the International Soil Moisture Network. In Proceedings of the EGU General Assembly Conference Abstracts, Vienna, Austria, 22–27 April 2012; p. 10189.

28. Baghdadadi, N.; Cresson, R.; El Hajj, M.; Ludwig, R.; La Jeunesse, I. Soil parameters estimation over bare agriculture areas from C-band polarimetric SAR data using neural networks. *Hydrol. Earth Syst. Sci. Discuss.* **2012**, *9*, 2897–2933. [CrossRef]

29. Babaeian, E.; Sadeghi, M.; Jones, S.B.; Montzka, C.; Vereecken, H.; Tuller, M. Ground, proximal, and satellite remote sensing of soil moisture. *Rev. Geophys.* **2019**, *57*, 530–616. [CrossRef]

30. Narasimhan, B.; Srinivasan, R. Development and evaluation of Soil Moisture Deficit Index (SMDI) and Evapotranspiration Deficit Index (ETDI) for agricultural drought monitoring. *Agric. For. Meteorol.* **2005**, *133*, 69–88. [CrossRef]

31. Zhang, N.; Hong, Y.; Qin, Q.; Liu, L. VSDI: A visible and shortwave infrared drought index for monitoring soil and vegetation moisture based on optical remote sensing. *Int. J. Remote Sens.* **2013**, *34*, 4585–4609. [CrossRef]

32. Jackson, R.D.; Idso, S.; Reginato, R.; Pinter, P., Jr. Canopy temperature as a crop water stress indicator. *Water Resour. Res.* **1981**, *17*, 1133–1138. [CrossRef]

33. Wang, X.; Zhao, C.; Guo, N.; Li, Y.; Jian, S.; Yu, K. Determining the canopy water stress for spring wheat using canopy hyperspectral reflectance data in loess plateau semiarid regions. *Spectrosc. Lett.* **2015**, *48*, 492–498. [CrossRef]

34. Jones, H.G.; Serraj, R.; Loveys, B.R.; Xiong, L.; Wheaton, A.; Price, A.H. Thermal infrared imaging of crop canopies for the remote diagnosis and quantification of plant responses to water stress in the field. *Funct. Plant Biol.* **2009**, *36*, 978–989. [CrossRef]
60. Saadi, S.; Boulet, G.; Bahir, M.; Brut, A.; Delogu, E.; Fanise, P.; Mougenot, B.; Simonneaux, V.; Lili Chabaane, Z. Assessment of actual evapotranspiration over a semi arid heterogeneous land surface by means of coupled low-resolution remote sensing data with an energy balance model: Comparison to extra-large aperture scintillometer measurements. *Hydrol. Earth Syst. Sci.* **2018**, *22*, 2187–2209. [CrossRef]

61. Funk, C.; Peterson, P.; Landsfeld, M.; Pedreros, D.; Verdin, J.; Shukla, S.; Husak, G.; Rowland, J.; Harrison, L.; Hoell, A.; et al. The climate hazards infrared precipitation with stations—a new environmental record for monitoring extremes. *Sci. Data* **2015**, *2*, 1–21. [CrossRef]

62. Bouaziz, M.; Medhioub, E.; Csaplovisc, E. A machine learning model for drought tracking and forecasting using remote precipitation data and a standardized precipitation index from arid regions. *J. Arid. Environ.* **2021**, *189*, 104478. [CrossRef]

63. Massman, W. A surface energy balance method for partitioning evapotranspiration data into plant and soil components for a surface with partial canopy cover. *Water Resour. Res.* **1992**, *28*, 1723–1732. [CrossRef]

64. Hersbach, H.; Bell, B.; Berrisford, P.; Hirahara, S.; Horányi, A.; Muñoz-Sabater, J.; Nicolas, J.; Peubey, C.; Radu, R.; Schepers, D.; et al. The ERA5 global reanalysis. *Q. J. R. Meteorol. Soc.* **2020**, *146*, 1999–2049. [CrossRef]

65. Allen, R.G.; Pereira, L.S.; Raes, D.; Smith, M. Crop evapotranspiration-Guidelines for computing crop water requirements-FAO Irrigation and drainage paper 56. *Fao Rome* **1998**, *300*, D05109.

66. Mega, N.; Medjerab, A. Statistical comparison between the standardized precipitation index and the standardized precipitation drought index. *Model. Earth Syst. Environ.* **2021**, *7*, 373–388. [CrossRef]

67. Hoffman, R.N.; Boukabara, S.A.; Kumar, V.K.; Garrett, K.; Casey, S.P.; Atlas, R. An empirical cumulative density function approach to defining summary NWP forecast assessment metrics. *Mon. Weather. Rev.* **2017**, *145*, 1427–1435. [CrossRef]

68. Chahbi Bellakanji, A.; Zribi, M.; Lili-Chabaane, Z.; Mougenot, B. Forecasting of cereal yields in a semi-arid area using the simple algorithm for yield estimation (SAFY) agro-meteorological model combined with optical SPOT/HRV images. *Sensors* **2018**, *18*, 2138. [CrossRef]

69. Joe, H. *Multivariate Models and Multivariate Dependence Concepts*; CRC Press: Boca Raton, FL, USA, 1997.

70. Delogu, E.; Boulet, G.; Olioso, A.; Coudert, B.; Chirouze, J.; Ceschia, E.; Le Dantec, V.; Marloie, O.; Chehbouni, G.; Lagouarde, J.P. Reconstruction of temporal variations of evapotranspiration using instantaneous estimates at the time of satellite overpass. *Hydrol. Earth Syst. Sci.* **2012**, *16*, 2995–3010. [CrossRef]

71. Lhomme, J.P.; Elguero, E. Examination of evaporative fraction diurnal behaviour using a soil-vegetation model coupled with a mixed-layer model. *Hydrol. Earth Syst. Sci.* **1999**, *3*, 259–270. [CrossRef]

72. Hoedjes, J.; Chehbouni, A.; Jacob, F.; Ezzahar, J.; Boulet, G. Deriving daily evapotranspiration from remotely sensed instantaneous evaporative fraction over olive orchard in semi-arid Morocco. *J. Hydrol.* **2008**, *354*, 53–64. [CrossRef]

73. Jackson, R.D.; Hatfield, J.L.; Reginato, R.; Idso, S.; Pinter, P., Jr. Estimation of daily evapotranspiration from one time-of-day measurements. *Agric. Water Manag.* **1983**, *7*, 351–362. [CrossRef]

74. Michelangeli, P.A.; Vrac, M.; Loukos, H. Probabilistic downscaling approaches: Application to wind cumulative distribution functions. *Geophys. Res. Lett.* **2009**, *36*, L11708. [CrossRef]