remote Analysis of Multispectral Drought Indices in Central Tunisia

: 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. interpretation N.F., Visualisation, M.L.P.; review G.B., M.L.P. J.C.; ideas 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,


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 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][34][35][36][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 applicablity of SI with routinely available data at a global scale.

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. 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 2 ) 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.

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.
where ρ N IR and ρ 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.

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.
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].

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].
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): where k is an extinction factor, about 1.13, NDV I max = 0.97 and NDV I 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].

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): where 1 {·} is the indicator function, x 1 , . . . , x n are rainfall amounts aggregated at different time scales (years, months or decades) andF n (x) is the UPI value withF n (x) = P(X ≤ x) ∈ [0, 1], with X the random variable representing the rainfall amounts at a given time scale. In contrast, SPI is computed as Φ −1 (F n (x)) where Φ −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.

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 Tsur f by assuming at first that the vegetation is not stressed, that allows a simplification to solve the underdetermination problem. Tsur f 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 pot and LEv pot respectively). The water stress indice 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 (halfhourly).

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 SWG .
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.

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 , . . . , x n ) be the values taken by the index X for a given data set. Then the ECDF is given by : where 1 {·} 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 ∈ [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).

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). 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).

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 objectoriented 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.

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 SWG and SI ERA5 , 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 SWG and SI ERA5 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. 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 ERA5 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 SWG 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 ERA5 , with a Kendall's τ of about 0.64. However, SI ERA5 presents low correlation coefficients with NDVI and SWI.

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 SWG , SI ERA5 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 present 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 SWG , SI ERA5 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 SWG , SI ERA5 and UPI are below the stress threshold. However, each lagged UPI value increase leads to an elevation in TIR indices.  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.  Besides, we observe that SI SWG and SI ERA5 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 SWG . 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.

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 XLAS , 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 SWG , SI ERA5 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 SWG and SI ERA5 succeed to better reproduce SI XLAS fluctuations in Figure 11a. In addition, we observe that these indices, especially SI SWG , follow the SI XLAS 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 XLAS (under the threshold line), observed in October, December and February, are mostly isolated with SI SWG . Moreover, SI SWG succeed in identifying stress conditions defined by SI XLAS , 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 XLAS 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 XLAS . 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 XLAS 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 XLAS in this growth season, particularly using SI SWG . Indeed, SI SWG succeed to identify a water stress period caught by SI XLAS 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.

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 SWG and SI ERA5 . 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).  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 SWG and SI ERA5 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.

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 SWG and SI ERA5 , respectively, are very well correlated with UPI. However, SI ERA5 has a higher Kendall's τ with UPI. Indeed, SI ERA5 is simulated using regional climatic information (ERA5 reanalyses) and local remote sensing information. The use of the large scale meteorological data to simulate SI ERA5 could explain its high correlation with the UPI also provided by gridded data-sets. Moreover, we observe a very strong seasonal correlation between SI SWG and SI ERA5 . 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 SWG 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 ERA5 , 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 SWG , which incorporates information from the local meteorological station close to the XLAS transect, shows better performances than SI ERA5 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 ERA5 . 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 ERA5 . 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.

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 requirement 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 vegetation 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 SWG 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). 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. 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:

Institutional
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)). where, 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): 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).
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. 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.