Estimating the Water Deficit Index of a Mediterranean Holm Oak Forest From Landsat Optical/Thermal Data: A Phenomenological Correction for Trees Casting Shadow Effects

Land surface temperature (Ts) integrated in contextual evapotranspiration models can be used to derive the water stress level of forests. One difficulty of such models implemented over forests is related to the impact of trees casting shadows on the remotely sensed Ts, which potentially hide the water stress signature. Until now, there has been no method to correct for this effect at the spatial resolution of current (Landsat) thermal sensors. This study investigates the impact of the solar zenith angle as a proxy of trees casting shadows on the water deficit index (WDI), using Landsat-7 and Landsat-8 data over a 21-km2 area partially covered by a holm oak forest in South-eastern France, for 7 successive summers (2015 to 2021). Results are compared to in situ evaporative fraction (EF) measurements. First, a calibration of the correction of WDI for shadow effects correlates the WDI error (WDI minus 1-EF) with the solar zenith angle. The correction improves the correlation (R) and slope between WDI and 1-EF, from R=0.32 to R=0.62 and slope=0.21 to slope=0.63, for non-corrected and corrected WDI respectively. Second, a calibration of the correction method that does not rely on in situ measurements, evaluates the linear relationship between remotely sensed WDI and the solar zenith angle in dry conditions. The correction improves the correlation and slope between WDI and 1-EF, from R=0.32 to R=0.57 and slope=0.21 to slope=0.50, for non-corrected and corrected WDI respectively.


I. INTRODUCTION
M EDITERRANEAN forests have dealt with high levels of hydric stress for long periods [1], [2]. However, global warming is a threat to most of forests all around the world [3], [4], [5] and particularly in the Mediterranean area, which is bound to be affected by increasingly frequent and severe droughts [6], [7]. Therefore, monitoring the hydric stress of Mediterranean forests is important in the context of climate change. It provides useful information for decision makers and landscape managers The authors are with the Centre Etude Spatial de la BIOsphére, 31000 Toulouse, France (e-mail: victor.penot@iut-tlse3.fr; olivier.merlin1@univ-tlse3.fr).
Digital Object Identifier 10.1109/JSTARS. 2023.3288360 to develop long-term forest/agricultural policies like the choice of species adapted to warmer and/or dryer conditions and mitigation management strategies [8], and also for firefighters to assess locally the short-term water dynamics of fuels [9], [10]. Remote sensing (RS) technology offers a good cost-benefit tradeoff to derive the evapotranspiration (ET) and water stress of ecosystems over wide areas on a daily/weekly basis. Different ET RS methods have been developed for this purpose and are generally based on optical/thermal data. There exists two broad categories: 1) the energy balance methods and 2) the contextual methods. The former are based on solving the physically based energy balance equations [11], [12]. They are generally more complicated to put into operational practice, due to their fine parametrization and their possibly large number of input parameters [13], [14]. The latter are based on semiempirical interpretations of the observed spatial correlations between land surface temperature (Ts) and fractional green vegetation cover (f vg ) [15], [16], [17] and/or albedo [18], [19] data. They may perform similarly to energy balance methods given a good calibration, whereas they are more parsimonious and hence simpler to apply over large areas [20], [21].
Under the assumptions of large gradients of surface hydric and vegetation cover states at the pixel scale and uniform meteorological conditions over the study area, the scatter plot of Ts versus f vg , often called Ts-f vg space, has the shape of a triangle or a trapezoid [15], [16]. The triangle method has, for instance, been used to assess the ability of Ts and f vg to discriminate the vegetation moisture content of a Mediterranean forest in Spain [22]. The water deficit index (WDI), a proxy of surface water stress, is notably determined from the relative position of a point in the Ts-f vg space [15], from a dry edge (fully dry pixels from bare soil to full canopy coverage) and a wet edge (well-watered pixels from bare soil to full canopy coverage). WDI has been widely used to assess the water stress of crops, and has shown good agreement with field measurements [15], [23], [24], [25]. However, to date, thermal RS methods have seldom been tested over Mediterranean forests. In an early study, Vidal and Devaux-Ros [26] used WDI as an index to predict the fire area/occurrence during the summer of 1990 in Southern France. However, the lack of site flux data in that study did not allow a quantitative assessment of the accuracy of the retrieved forest WDI. Whatever thermal-based method is used to estimate the water stress, cast shadows within the canopy and on the ground have a strong influence on the measured Ts [27]. However, it is difficult to model and evaluate accurately the cast shadows impact in complex ecosystems [27]. The use of radiative transfer and/or geometric projection models supports these results, if a mock-up can be designed. The authors in [28], [29], [30], and [31] used these tools to study directional effects in oak woodlands or crops. They especially showed the strong impact of surface component fractions within the satellite field of view on Ts retrievals, and particularly through the contribution of shaded and sunlit soil.
To derive ET or water stress proxies in nonhomogeneous canopies, the authors in [32] and [33] suggested the use of data collected at very high resolution to explicitly take shadow effects into account. Three different methods have been developed from ground-based [33], [34], [35], [36] or aircraft/UAV-based [32], [37], [38], [39] thermal camera. One solution consists removing in the Ts image, the pixels corresponding to shaded soil/leaves to isolate sunlit leaves [32], [37], [38], [39]. It relies on the assumption that sunlit leaves are more likely to be stressed than shaded ones [40] and that they may provide meaningful information about the vegetation water status/flux. Separating shaded from sunlit soil/leaves requires a very high-resolution (centimetric) sensor. Such a classification can be done manually [35], [37], [39] or by automatic detection using supervised classification [39], [41], density-based methods [34], [36], [42], clustering (K-means) and postclassification [38], high resolution elevation thresholds [37], or vegetation index thresholds [32]. However, due to the geometric misregistration of visible/near-infrared and thermal images at such a fine resolution, soil and shaded leaves background may contaminate the thermal data of an area identified as sunny and vegetated [38], [42]. In addition, these methods rely mainly on a priori sunlit leaf properties, which are either physically determined or manually selected from images to build a training database. Therefore, the presence of mixed pixels or sunlit reflectances outside the training database may decrease the accuracy of the segmentation [41].
Another method to take shadow effects into account, consists in simply setting a ground-based camera in a position that avoids as much as possible shaded leaves/soil in the field of view of the thermal radiometer. The camera is hence generally set in a nadir/solar noon position so that only the sunlit part of the tree is observed without soil background [33], [34]. Note that this method strongly depends on the solar zenith angle and on the remaining ground and shaded leaves contamination [33], [34].
Rather than extracting sunlit leaves temperature from the thermal image or avoiding shadow effect in the sensor field of view, Luquet et al. [31] suggested a third option, that is to correct the WDI method to take into account directional effects. This correction involved redefining the wet and dry edges of the Ts-f vg space by using a 3-D model that was run in each solar and observation angle. However this method uses realistic and efficient 3-D representation of every plant present within the study area. It requires both extensive and exhaustive ground measurements and a relatively small study area due to the high computation cost.
The above literature review indicates that until now no operational solution has been proposed to correct the remotely sensed water stress proxies for the cast shadow/directional effects at the scale of current (Landsat and MODIS/Sentinel-3) spaceborne thermal sensors. This is also true for near-future thermal missions which will provide Ts data at 57-and 50-m spatial resolution from 2025 and 2028 for TRISHNA and LSTM, respectively. In this context, the objective of this study is to propose an operational method for correcting the satellite-derived WDI for trees casting shadow effects, by relying on the solar zenith angle (θ S ). Specifically, this article aims at 1) quantitatively assessing the impact of θ S on WDI and 2) developing a self-calibration correction of cast shadow effects.
To reach these goals, the study focuses on a common holm oak forest at Puechabon site, in Southern France, where in situ flux measurements are available between 2015 and 2021. Section II describes the study area, the in situ and Landsat data used and the new method to correct the Landsat-derived WDI for cast shadow effects. Section III presents the results obtained by calibrating the correction method using in situ data (site calibration) or satellite data solely (self-calibration) and discusses them in the prospect of future regional applications. Finally, Section IV concludes this article.  [43]. This site is representative of Mediterranean evergreen broadleaf forests with dense coppice. The holm oak (Quercus Ilex L.) is the dominant tree species, with a mean tree height of 5.5 m. The understory is composed of sparse shrubs of 2-m height [2]. The site was built in 1984 to study Mediterranean ecosystems and specially the holm oaks response to climate change and to severe droughts [1], [43].

A. Puechabon Site and in situ Data
The mean annual temperature over the 2008-2021 period is 14.5 • C. The mean annual precipitation (2008-2021) is 877 mm with strong intra annual variations. Summer is dry and hot, while heavy precipitations occur primarily in fall but also in spring. For illustration, Fig. 1 plots the monthly cumulative rain and the monthly mean air temperature between 2015 and 2021 at the Puechabon site. This region is subject to a cool and dry north-west wind, known as "Mistral." Table I summarizes intra and interannual differences of cumulative rain and air temperature during the 2015-2021 study period. Three particular groups of years can be highlighted. 2015 was a particularly wet year during summer: it observed 267-mm cumulative rain between June and September, and 166-mm cumulative rain between July and August. 2017 was a particularly dry and hot year with only 74-mm cumulative rain between June and September, and 13-mm cumulative rain between July and August with five days over 35 • C between June and September. 2016 also experienced a very dry summer (17 mm-cumulative rain between July and August), while 2018, 2019 and 2020 were hot years as regards the number of days above 30 • C between June and September (38, 51, 41, respectively) and/or the number  of days above 35 • C between June and September (5, 3, 1, respectively). It is important to note that in 2016, 2018, and 2019, the month of October was the rainiest of the year, suggesting that cloud cover during that month is likely to weaken the use of optical RS.
The soil at Puechabon site is silty claim loam with a rock and stone volumetric content of 75% in the first 50 cm and 90% in the whole profile [44]. It is poorly developed and does not store water well [2]. However, it was shown that Quercus Illex was able to extract water deeper than 4 m during drought [44], [45].
2) Flux and Meteorological Data: The station of Puechabon is equipped with an eddy covariance flux measurement system that measures meteorological data and energy fluxes. The flux tower is 10-m height and equipped with a Sonic anemometer (Solent R3A Gill) and an infra red gas analyzer (IRGA, LI-6262,Li-COR) for latent and sensible heat flux measurements. The sampling is done at a rate of 21 Hz and recorded every 30 min. Allard et al. [45] studied the measurement footprint and its influence on the flux data in all weather conditions. It was shown that the measured fluxes come from a dense Quercus Illex coppice, surrounded by a less dense Quercus Illex dominated vegetation and that there is no seasonal pattern of footprint location. A temperature and humidity transmitter (MP100, Rotronic) is used for air temperature and rainfall measurement [45], [46]. Records are done every 30 min. Flux and meteorological data from 2015 to 2021 were downloaded from the European Fluxes Database Cluster.
Regarding the monitoring and modeling of ET/hydric stress from thermal data, few studies have been undertaken at the Puechabon site. The authors in [47], [48], and [49] evaluated the 1-km resolution eight-day ET MODIS product. More recently, Ollivier et al. [2] developed an ET model driven by the MODIS EVI (enhanced vegetation index) for assessing the groundwater resource in the surrounding Karst watershed.
3) Evaporative Fraction: The evaporative fraction (EF) is defined as the ratio of latent heat to the available energy [16] with LE being the latent heat (W/m 2 ), Rn the net radiation (W/m 2 ) and G the ground conduction (W/m 2 ). G is not available at Puechabon site. So, under the assumption of energy balance closure, EF is expressed as with H being the sensible heat (W/m 2 ). Allard et al. [45] showed that the energy closure (77%) at this site is reasonable and consistent with other similar sites.
The property of EF to remain constant during daytime hours [17], [50], [51] allows to check the quality of instantaneous EF estimates by means of the Bowen ratio β. The Bowen ratio is the ratio of the sum of H from 8 A.M. to 3 P.M. by the sum of LE from 8 A.M. to 3 P.M. [17]. Hence the daytime EF (EFd) can be estimated by means of the Bowen ratio A statistical analysis during the study period shows a good agreement between instantaneous EF at 10:30 A.M. and the daily EFd (R 2 = 0.89, slope = 0.97). Hence the instantaneous EF is kept as a reference to evaluate the satellite-derived WDI. In practice, due to the inverse change direction of EF and WDI, the reference in situ hydric stress is estimated as 1-EF.

1) Remote Sensing Data:
Landsat Enhanced Thematic Mapper Landsat-7 (ETM+) and Landsat-8 Operational Land Imager (OLI) and Thermal Infrared (TIR) both provide reflectance and Ts data with a 16-day revisit period, with an 8-day shift between Landsat-7 and Landsat-8. Landsat data for Path/Row 197/030 were downloaded from USGS Earth Explorer platform/Landsat Collection 2 Level-2 Science Products from 2015 to 2021. Shortwave optical data were corrected to be bottom of atmosphere reflectances and provided at their native resolution of 30 m. Landsat-7 data were processed with the landsat ecosystem disturbance adaptive processing system (LEDAPS) algorithm (Version 3.4.0) [52] and Landsat-8 data with the land surface reflectance code (LaSRC) algorithm (Version 1.5.0) [53].
After 2017, Landsat-7's orbit drifted from its original schedule to an earlier overpass time. Changes in reflectance and thermal data are limited until 2020 [54]. To ensure consistency between Landsat-7 and Landsat-8 data, we chose to keep only Landsat-7 acquisitions after 10:00 A.M. in the dataset.
f vg is estimated using the expression of [55] f vg = NDVI − NDVI min NDVI max − NDVI min 2 (4) with NDVI being the normalized difference vegetation index defined as the difference of near-infrared and red reflectances divided by their sum. NDVI min = 0.19 (bare soil) corresponds to the quantile 0.01 of NDVI over time (2015-2021) and space in the study area, and NDVI max = 0.81 (fully green vegetation cover) corresponds to the quantile 0.97 of NDVI over time (2015-2021) and space in the study area.
Landsat-7 and Landsat-8 thermal data have a native resolution of 60 and 100 m, respectively. Both datasets were downsampled at 30-m resolution and provided over the same 30-m resolution grid of Landsat reflectances. They were processed with the Landsat surface temperature algorithm (Version 1.3.0) to produce Ts maps from raw thermal radiances [56], [57].
2) Study Area Extent and Description: Contextual ET methods based on the Ts-f vg space require the largest gradients of vegetation coverage and hydric state [58], [59]. Meteorological conditions must also be relatively homogeneous within the study area [14]. Thus the study area extent is defined to include the Puechabon site and to observe the largest gradient of fraction cover, from bare soil (quarry and vineyard) to fully covering forests, and also the largest gradient of hydric state (from the Herault river banks to bare dry soils), while minimizing the gradient of elevation. The 21-km 2 selected area is presented in Fig. 2. On its northern part, it is delimited by the narrow Herault river and a plateau mainly covered by closed deciduous forests, were the flux site is set. Note that according to the data provider definition (Institut Géographique National, IGN), a closed canopy forest is covered by more than 40% of trees that reach a height of 5 m and more. Vineyards and some orchards (olive trees) occupy most of the Southern area. No town or major built area is included. The northern plateau and southern agricultural area are delimited by a topographic break [see Fig. 2

C. WDI Method and Image Selection
The methodology to set the dry and wet edges of the Tsf vg space (Section II-C1 and II-C2) and to estimate the WDI (Section II-C3) is first presented. Then, a quality control of the input data is undertaken to automatically filter out the dates when meteorological conditions weaken the WDI application (Section II-C4).
1) Dry Edge Estimation: Dry conditions are generally encountered in the study area, due to the climatic context during late spring, summer, and early autumn. Therefore, the dry edge is defined as the linear regression of the 99% quantiles of Ts evaluated in equal-length f vg bins against f vg , as in [25] (see Fig. 3).
2) Wet Edge Estimation: Wet conditions are seldom met in the study area due to its climatic context. Hence, a physically based assumption must be set to solve this situation. The wet edge corresponds to pixels where LE is maximal and H minimal whatever f vg . Thus when LE is maximal, H is negligible and can be set equal to 0 [17]. According to [60], Ts is approximatively equal to aerodynamic temperature, so H is as a first guess proportional to the difference between Ts and air temperature (Tair) [15]. Following this reasoning, the wet edge in this study is set constant according to the equation Ts = Tair(10:30) at the satellite overpass hour, i.e., 10:30 A.M. (see Fig. 3). This assumption has been successfully assessed in other previous studies undertaken in semiarid environments [19], [61], [62].
3) WDI Calculation: Fig. 3 illustrates how WDI is computed in the Ts-f vg space from dry edge and wet edge distances [15] WDI = BC AB (5) with C being the point of interest, A the point in the dry edge line of same abscissa, and B the point of the wet edge line of same abscissa. By definition, WDI is the complement in 1 of EF, i.e., WDI = 1 − EF. 4) Date Selection: First, the use of thermal RS to evaluate hydric stress requires water-limited conditions meaning that Ts is driven mainly by water availability rather than by the energy available at the surface [14], [59]. Hence, the study period is  Second, the presence of clouds, cloud shadows, or recent rainfalls is likely to decrease Ts under the wet edge, defined as Tair(10:30), as previously mentioned. Moreover, unstable/changing meteorological conditions may result in heterogeneous solar radiation, wind speed, and Tair(10:30) within the area, which do not allow the use of contextual ET methods [14]. Therefore, an approach is proposed to automatically remove such dates. Only the pixels qualified as clear in the quality masks provided with Landsat-7/8 reflectance products are used during this filtering step. Zhang et al. [61] suggested that pixels with a Ts lower than Tair(10:30) are contaminated by clouds or located over a water body. Accordingly, a date is removed from the dataset if min(Ts) ≥ Tair(10:30) − 1, considering that the image is contaminated by clouds. In addition, a minimum number of pixels within the image must be valid (qualified as clear in Landsat-7/8 masks) to get a representative view of all the required conditions in the study area. A threshold of 85% of clear pixels within the study area is chosen to filter out cloudy conditions. This limit quantifies the clear or low cloudy conditions that are commonly used in previous studies [15], [19], [63], [64].

D. Correction Method of WDI for Solar Zenith Angle
1) General Approach: The correction method proposed in this article relies on the assumption that cast shadows effects are strongly linked to solar zenith angle (θ S ) through cast shadow geometries, although this link is not explicit in complex geometries, such as tree canopies. For clarity, the framework of the correction strategies is summarized in the Fig. 4.
Once the WDI is classically estimated (as described in Section II-C3, see Fig. 4 upper part), a relation is built between WDI and θ S with ΔWDI being the error on WDI attributed to cast shadows and a and b two empirical parameters. Two approaches are then proposed to calibrate a and b 1) if flux data are available, (a site , b site ) is determined from the linear regression of ΔWDI = WDI − EF + 1 with θ S in particularly dry conditions. This correction strategy is called site-calibrated correction (see Fig. 4 left part); 2) if flux data are unavailable, which is the case in most operational situations, a self can be estimated by interpreting the regression of WDI with θ S in particularly dry conditions, and b self is set to the minimum of θ S in the May-September period, namely, θ S,Min . This correction strategy is called self-calibrated correction (see Fig. 4 right part). For both correction strategies, ΔWDI is subtracted from WDI. Finally, the corrected WDI (WDI c ) is compared to in situ 1-EF (see Fig. 4 lower part) to evaluate the correction performance. Each component of the correction method is described as follows.
2) Identification of Dry and Very Dry Conditions: Hydric stress and shadow effects are generally mixed, and the main challenge herein is to disentangle both effects from WDI and θ S data solely. Especially, the shadow effects on WDI may overwhelm the signature of hydric stress. During droughts, or during long periods without significant rains and high Tair, hydric stress grows to its maximum value, as long as trees and shrubs sustain the lack of water. In such extreme situations, the variations over time of WDI are assigned to variations of cast shadows geometry, which are mainly linked to θ S . Hence, to disentangle cast shadow effects from hydric stress, only the driest dates of a time series should be used. We define the Very Dry dataset, as the dates of our dataset whose 15-day cumulative rain (P 15 d ) is under the first quartile of the 15-day cumulative rain of all days during the 2015-2021 study period (q 0.25 (P 15 d ) = 7 mm). This dataset will be used for both site and self calibration strategies (see Section II-D3 and II-D4). The 15-day time depth is chosen on the basis of the regression quality of site calibration. We also define the Dry dataset (Wet dataset), as the dates of our dataset whose 15-day cumulative rain (P 15 d ) is under (over) the median of the 15-day cumulative rain of all days during the 2015-2021 study period (q 0.5 (P 15 d ) = 23mm). These datasets will be used for evaluation purposes (see Section II-D5). The whole dataset is called All dataset.
3) Site-Calibrated Correction: When and where flux data are available, the difference between WDI and 1-EF regressed over θ S can be used to correct the WDI [see (6)]. Only the Very Dry dataset is used for this regression. On such dates, the evolution of ΔWDI may be attributed to cast shadow effects rather than to hydric stress. This regression is useful to quantitatively assess the impact of θ S on WDI at Puechabon site and to verify the underlying assumption of the correction methodology. The sitecalibrated correction equation is

4) Self-Calibrated Correction:
In practice, flux data are rarely available and in general, they cannot be used for calibration purposes of (6) due to the strong heterogeneity of land surfaces. Therefore, a self-calibration method must be imagined to operationally apply our proposed correction method, i.e., first by calibrating the slope a of (6) and then by setting the value of b.
A two-step procedure is proposed to isolate the dates when the WDI evolution is strongly influenced by cast shadow effects and to calibrate the slope a = a self of (6). First, only the Very Dry dataset is used, i.e., rainfall conditions when WDI is likely to be close to 1. Note that selecting the driest dates using the cumulative rain (instead of using large WDI values directly) allows for filtering out the dates when WDI is much affected by cast shadow effects. Then, WDI is plotted against θ S and the upper edge of this space is linearly interpolated. The points located on the upper edge correspond to the pixels with an actual hydric stress index close to 1, so their WDI evolution with θ S is assumed to be the result of cast shadow effects only. Hence the slope of the linear regression, called a self quantitatively expresses the evolution of WDI with θ S during droughts or very dry periods. While θ S increases, cast shadow effect increases, so we expect that a self < 0.
The intercept b self = b is set to the minimum of θ S (θ S,Min ) reached in late June at the site/pixel location. For a given canopy geometry, the impact of shaded areas at this time should be minimal when the Landsat viewing angle is near nadir [28], [29], which is the case in the Puechabon study area (< 1.5 • ). Finally, we assume that the effect of θ S on WDI is negligible when θ S = θ S,Min and that it can only decrease WDI when θ S increases. Therefore, the self-calibrated correction equation is

5) Evaluation of Calibration and Correction Strategies:
To assess the quality and the improvements brought by both (siteand self-calibrated) corrections, the corrected WDI is evaluated against the site measurements of 1-EF for All dataset, Wet dataset, or Dry dataset separately. In each case, the root mean squared error (RMSE), mean absolute error (MAE), coefficient of correlation (R), slope and bias of the linear relationship between satellite WDI, and in situ 1-EF are calculated as metrics of quality and accuracy.

III. RESULTS AND DISCUSSION
To assess the impact of θ S on WDI, Section III-A focuses on the results of the site-calibrated correction. In Section III-B, the results the self-calibrated correction are presented and discussed to assess the correction method in an operational context. In Section III-C, a spatial analysis of the self-calibrated correction method is proposed.

A. Site-Calibrated Correction of WDI for Solar Zenith Angle
1) Site Calibration: As mentioned previously, (1) is calibrated with the Very Dry dataset (including the 25% driest dates) with the objective of identifying the dates when θ S significantly drives the error ΔWDI = (WDI − EF + 1). To assess the usefulness of this assumption, Fig. 5(a) and (b) plot the error ΔWDI versus θ S for All dataset and Very Dry dataset, respectively.   when All dataset and Very Dry dataset is used for calibration separately. Fig. 6(a) and (b) plot the noncorrected and site-calibrated corrected WDI (using Very Dry dataset for calibration) versus 1-EF, respectively. Both site calibration strategies improve the quality of the corrected WDI compared to the noncorrected WDI [see Table II and Fig. 6(a) and (b)]. The slope and correlation coefficient increases from 0.21 and 0.32 (no correction) to 0.50 and 0.57 for the All dataset and 0.63 and 0.62 for the Very Dry dataset used for calibration, respectively. The bias and RMSE also significantly decrease from bias = −0.10 and RMSE = 0.17 (no correction) to bias = 0.00 and RMSE = 0.12, and bias=0.02 and RMSE=0.12 for the All dataset and Very Dry dataset-based calibration corrections, respectively. Note that the RMSE remains in the standards of other studies in similar context [64], [65], [66].
Even if the calibration regression using All dataset is significant (determination coefficient R 2 = 0.25 and pvalue< 5%), its quality is clearly improved when only the Very Dry dataset is used for calibration [R 2 = 0.59 and pvalue< 5%, see Fig. 5(a)  and (b)]. These results suggest that θ S significantly drives the evolution of ΔWDI when only very dry conditions are encountered, i.e., when cast shadow effect is the main driver of WDI. In fact, when the actual hydric stress index is close to 1 and when the contrast between sunny and shaded areas is large, the impact of hydric stress on Ts is expected to be small compared to cast shadow effects. Based on the above analysis, the Very Dry dataset is henceforth used for calibrating the correction method.
The slope a site in Fig. 5 is negative, which is consistent with the fact that the cast shadow effect on WDI increases with θ S (initial assumption). The site calibration provides a correction slope a site = −0.020 [see Fig. 5   negligible when θ S is minimum and that ΔWDI becomes more and more negative when θ S increases. This is consistent with a minimum cast shadow effect on Ts and hence on WDI in late June, while the Landsat viewing angle is near nadir (< 1.5 • ) at the Puechabon site [28], [29]. In brief, the site calibration provides an equation that efficiently corrects for the cast shadow effects over the full θ S range (from May to September). The value of the correction is null when θ S is minimum in late June, it linearly increases with θ S and reaches its maximum value in late September when θ S is maximum.
2) Evaluation of the Site-Calibrated Correction: Table III summarizes the statistics obtained for different calibration strategies and evaluation datasets. The site-calibrated corrected WDI outperforms the noncorrected WDI when evaluated on the Dry dataset (37 dates). When evaluated on the Dry dataset, the correction improves the slope of the linear regression (and correlation coefficient) from 0.07 to 0.84 (from 0.08 to 0.59) for the noncorrected and site-calibrated corrected WDI, respectively. In parallel the bias decreases from −0.15 to 0.00 for noncorrected and site-calibrated corrected WDI, respectively. The good performance of the correction method for the Dry dataset is explained by strong cast shadow effects when the contrast between sunny areas and shaded areas is maximum.
However, the site-calibrated correction is not so successful when it is evaluated on the Wet dataset (20 dates). Even if R and slope still increase from R=0.21/slope=0.14 (no correction) to R = 0.36/slope = 0.37 (site calibrated correction), the bias significantly increases from bias = −0.02 (no correction) to bias = 0.08 (site-calibrated correction) (see Table III). These results support the fact that wet dates/areas are less likely to be affected by shadow effects. Consequently, a correction equation that is calibrated over Very Dry dataset tends to overcorrect WDI for the Wet dataset.
Note that the land-atmosphere feedback on air temperature [67] was neglected in this study. Especially, due to the general absence of wet conditions in the study area, the wet edge was set to a constant value based on the equation Ts=Tair (measured at Puechabon site at 10:30 A.M.), regardless of the stress level of the underlying surface. This strategy is likely to overestimate the temperature of the wet edge and hence to underestimate the satellite-derived stress index in dry conditions. Disentangling both (feedback on air temperature and trees casting shadow) effects on the satellite-derived stress index should be investigated in the future.
3) Yearly Analysis of Site-Calibrated Correction Results: Fig. 7(a) and (b) plots noncorrected and site-calibrated corrected WDI versus 1-EF per year. The correction not only improves the 1-EF estimation over the whole study period but also for almost every year separately (2016, 2017, 2018, 2019, 2020, 2021) but 2015.
2015 is a particularly wet year that illustrates the limitation of the correction methods: the proposed correction for trees casting shadow effects does not bring any improvement to the classical WDI without correction [see Fig. 7(a) and (b)]. The slope decreases from 0.30 to 0.23 and the correlation coefficient from 0.53 to 0.28 for noncorrected WDI and site-calibrated corrected WDI, respectively. Particularly, the correction for the Wet dataset tends to overcorrect WDI. This year is the rainiest of the seven studied years during summer, with 267-mm cumulative rain between June and September, after the driest beginning of year (177 mm of cumulative rain between February and May). Locally, Limousin et al. [1] actually measured in 2015 the maximum of minimum predawn summer water potentials-a These results suggest that when water is available, cast shadows is not the main driver of WDI variations. However, the relatively small number (6-12) of observations available for each year separately, does not allow a deeper evaluation of the correction method on an interannual basis.

B. Self-Calibrated Correction of WDI for Solar Zenith Angle 1) Self-Calibration Equation:
For operational needs, a selfcalibrated correction method is developed by plotting WDI versus θ S , without relying on flux data. The calibration of the slope a self of (6) is undertaken using both the All dataset and the Very Dry dataset separately to assess the performance of each calibration strategy. As for the site calibration case, selecting the Very Dry dataset allows to better identify the impact of cast shadow effects. Therefore, the calibration of a self using the Very Dry dataset provides a more efficient correction than when it is calibrated with the All dataset. In Table II, we see that with the All dataset, the slope (and R) only increases from slope = 0.21 (0.32) to 0.31 (0.44) for the noncorrected and self-calibrated corrected WDI, respectively. In contrast, using the Very Dry dataset the slope (and R) significantly increases from 0.21 (0.32) to 0.50 (0.57) for the noncorrected and self-calibrated corrected WDI, respectively. In parallel the bias decreases (in absolute value) from -0.10 (no correction) to −0.07 (All dataset-based calibration correction) and to −0.02 (Very Dry dataset-based calibration correction). These results suggest that isolating particularly dry conditions during the calibration step allows to get a more efficient correction on WDI [see Fig. 5(c) and (d)], and are consistent with the same analysis undertaken for the site calibration. We now consider the WDI correction results obtained from the calibration of a self using the Very Dry dataset only.
The negative value of the slope a self = −0.013 is consistent with the increase of cast shadow effects with θ S , as already mentioned for the site-calibrated correction case [see Fig. 5(b) and (d)]. Even if a self is lower (in absolute value) than a site , the previous discussion in Section III-A1 remains valid.
Only seven years of data have been used for the selfcalibration, meaning that the calibration strategy is likely to be sensitive to outliers. Even if this dataset covers a wide range of meteorological conditions, including particularly wet (2015) and dry (2017) years, a longer calibration period is likely to provide more statistically significant results. In fact, the upper edge of the WDI-θ S space is strongly influenced by extreme points that may be encountered during dry and hot years. In other words, drier conditions might change the calibrated value of the self-calibration slope a self .
To define the final calibration (8), the intercept is set to the minimum θ S between May and September-b self = θ S,Min = 25.7 • -which is close to the site-calibrated b site = 25.6 • . It is suggested that cast shadow effects are negligible in mid-June, when the difference between the viewing zenith angle of Landsat (near nadir) and θ S is minimum. Despite the good agreement between b site and θ S,Min at the Puechabon site, further researches should be conducted to assess the variability of b parameter in time and space.
2) Evaluation of the Self-Calibrated Correction: Results presented in Fig. 6 and Table III indicate that the self-calibrated corrected WDI outperforms the noncorrected WDI when evaluated on the Dry dataset. In particular, the self-calibrated correction improves the slope and correlation coefficient (see Table III), which increase from 0.07 to 0.59 and from 0.08 to 0.48, respectively. In parallel the bias decreases in absolute value from 0.15 to 0.05 for noncorrected and self-calibrated corrected WDI, respectively. However, the self-calibrated corrected WDI is not so successful when evaluated on the Wet dataset. Even though both R and slope are slightly improved from 0.21 to 0.35 and from 0.14 to 0.29, respectively, the absolute bias increases from 0.02 (no correction) to 0.04 (site-calibrated correction) (see Table III).
These results are fully consistent with those obtained for the site-calibrated correction case discussed previously. The behavior of the self-calibrated corrected WDI is very similar to that of the site-calibrated corrected WDI, as illustrated in the per year analysis of Fig. 7(a) and (c). Because a site < a self , the self-calibrated correction is slightly less effective to correct WDI for Dry dataset, and consequently limits the overcorrection of WDI for Wet dataset. The self-and site-calibrated corrections are generally of equivalent quality. This result confirms that the slope of the linear upper edge of the WDI-θ S space in particularly dry conditions (Very Dry dataset) is a good approximation of the a site parameter in (8).

C. Spatial Analysis of the Self-Calibrated Correction
The self-calibrated correction of WDI for solar zenith angle has been assessed in terms of WDI accuracy at the Puechabon site. It is also interesting to investigate the spatial behavior of the corrected WDI and the correction parameters (a self ) at the pixel scale within the study area. The study area [see Fig. 2(a)] is characterized by a large variability of land use and land covers from bare soil, organized vineyards and orchards in the south, to wild Mediterranean evergreen forests in the north. To assess spatially the self-calibrated correction in such media, the method is applied at the 30-m resolution over the entire study area and during the entire study period (2015-2021). The satellite-derived WDI is smoothed with a 5 × 5 pixels split window to take into account the actual spatial resolution of thermal Landsat data (60 and 100 m for Landsat-7 and Landsat-8, respectively). θ S is supposed to be uniform withinover the whole study area on each Landsat overpass. Note that for few pixels, the retrieved correction slope a self is positive; in this case it is set to 0. The minimum θ S is set to 25.7 • at each pixel of the study area. The self-calibrated correction equation can finally be applied to correct WDI for cast shadow effects over the entire study area. Fig. 8 shows maps and density histograms of WDI and WDI c,self , respectively, for data on the 22th, July 2016. The summer 2016 is particularly dry and hot at the Puechabon Site. Cumulative rainfall during the 15 previous days is only 0.2 mm. The measured 1-EF at the flux site is 0.88, site-calibrated corrected WDI c,self is 0.81 while noncorrected WDI is 0.67. It corresponds to an error reduction of 77% by correcting the WDI for trees casting shadow effects.
In Fig. 8(a), the noncorrected WDI ranges from 0.2 to 1. In fact, the difference between the minimum Ts (299.7 K) and Tair at 10:30 A.M. (297.7 K, wet edge value) is 2 K, which implies strictly positive WDI values. The spatial distribution of WDI generally reflects the variability in land cover/use. On the one hand, the southern part of the study area is covered with vines and orchards, and shows a mean WDI of 0.84. This is consistent with both large bare soil areas of intercrop (the intercrop herbaceous species are senescent at this time) and the relatively large size (compared to the interrow spacing of vineyard for instance) of a Landsat-7/8 thermal pixels. On the other hand, the northern part of the study area is mainly covered by holm oak dominated forest with a mean WDI of 0.66. As holm oak is able to efficiently deal with low levels of soil water availability, a lower value of WDI lower than in vineyard/orchards would be expected [1]. Further north, the slope along the Herault river witnesses particularly low values of WDI. Fig. 8(b) presents the self-calibrated corrected WDI. Vineyards and orchards are not particularly affected by the correction-the (mean WDI c,self is 0.87-due to the near-zero slope correction a self (see next paragraph). In contrast, WDI c,self is particularly contrasted and enhanced in the forest area. The mean WDI c,self now reaches 0.74 (increase of 12%), and very  few values are below 0.4. This is particularly true in the denser part of the forested area and along the Herault river.
The correction slope a self is mapped in Fig. 9(a) and density histograms of a self values per dominant land use/cover class are provided in Fig. 9(b). Values of a self are strongly related to the land cover map [see Fig. 2(a)] and to the elevation map [ Fig. 2(b)]. In the southern part of the study area, vineyards and orchards (olive trees) show close to zero correction, with a mean slope correction of -0.003 [see Fig. 9(a) and (b)]. These land covers are characterized by large intercrops including bare soil areas or fast drying herbaceous species. Note that on the southern part, along the Herault river, quarries provide hight value of a self . On the northern plateau where the flux tower is located, the Quercus Illex dominated forest is less open and witnesses more complex landscapes than orchards/vineyards, which implies a significant decrease in the mean a self equal to -0.008 [see Fig. 9(a) and (b)]. The northern valley along the Herault river is characterized by very low a self values, consistent with topographic effects. Note that Landsat-7 artefacts (strips) are clearly visible. This is due to the impact of a key date (in particularly dry conditions) of the Landsat-7 time series used in the calibration of a self .

IV. CONCLUSION
Until now, no method has allowed to correct satellite-derived WDI for trees casting shadow effects over large areas. To fill the gap, this study is based on solar zenith angle (θ S ) to predict in time the effect of cast shadows on WDI.
The objectives are: 1) to quantitatively assess the impact of θ S on WDI during dry periods (high values of hydric stress), 2) to build a self-calibration method that corrects WDI on a pixel basis from θ S and satellite data solely for operational use, and 3) to evaluate calibration strategies against in situ 1-EF used as ground truth. The classical WDI method (without correction from cast shadow effects) is used as benchmark to evaluate the improvement of the corrections. The approach is tested using in situ data collected between 2015 and 2021 at the Puechabon site located in a holm-oak forest in Southern-east of France.
As a first assessment of the effect of θ S on WDI, the error in WDI at the Puechabon site [WDI -(1-EF)] is correlated to θ S (R 2 = 0.59). By applying the site calibrated correction of WDI for cast shadow, the slope of the linear regression between remotely sensed WDI and in situ 1-EF is much closer to 1: it increases from 0.21 to 0.63. In addition, the RMSE decreases from 0.17 to 0.12 and the correlation coefficient (R) increases from 0.32 to 0.62 for noncorrected and corrected WDI, respectively.
As a further step, a self-calibrated correction method is developed at the pixel scale, without relying on in situ data, by correlating high WDI values to θ S on particularly dry dates. This correction method still succeeds in improving WDI estimates: the R, the slope of the linear regression, and the RMSE are 0.57 (0.32), 0.50 (0.21), and 0.12 (0.17) for the corrected (noncorrected) WDI, respectively. Therefore, the self-calibrated method provides an operational workflow for practical use.
When comparing the correction results in dry and wet conditions, relatively low values of correlation (0.36 and 0.35) and slope of the linear regression (0.37 and 0.29) between satellite WDI and in situ 1-EF are obtained for the Wet dataset for self-and site-calibrated corrected WDI, respectively. In fact, the calibration strategy that relies on dry dates is not suitable for wet conditions, as the Ts contrast between sunny and shaded areas decreases in cool/wet conditions. It is also assumed that other physical factors than θ S can affect the thermal behavior of forests, especially when the hydric stress is relatively low (e.g., thermal inertia [68], [69]).
Although quite promising results were obtained at the Puechabon site using our simple θ S -based correction method, the approach should be tested to other sites including a variety of tree species and landscapes. In particular, we used the cumulative rain with a 15-day time depth to extract the Very Dry dataset for calibration purposes. The genericity of such an assumption should be assessed at other sites as this value might be linked to both soil water availability and tree physiological functioning. In addition, we set the intercept b self to the minimum of θ S , while Landsat viewing zenith is near nadir in the Puechabon site area. Note that it is not general as the Landsat-8 viewing angle can reach 8 • , which suggests nonnegligible directional effects at the edges of the swath.
Several of the limitations mentioned above could be removed by implementing a physically based 3-D radiative transfer model like DART [70]. However, this would require the combined used of energy balance models to precisely solve the temperature of soil and vegetation components [71], [72]. Last but not least, this would also require a precise and exhaustive 3-D digital mock-up of the study area, which remains challenging due to both the lack of detailed knowledge of complex land surfaces and the constraints associated with calculation times. The TRISHNA mission, to be launched in 2025, will provide thermal data at high spatio-temporal resolution. The high temporal frequency of TRISHNA data is related to the relatively large swath width of the instrument, which will imply significant viewing angle directional effects [73], [74]. Thus correcting for directional effects is especially relevant for the future operational use of these valuable data.