Estimating the water decit index of a Mediterranean holm oak forest from Landsat optical/thermal data: a phenomenological correction for trees casting shadow effects

Victor Penot, Olivier Merlin

Abstract-The water stress level of forests, as a measure of how trees react to heat waves and droughts, provides crucial information for forest management and climate change mitigation policies. Thermal remote sensing data can be used to derive the water stress of ecosystems. Especially the contextual evapotranspiration models, which rely on spatial correlations between land surface temperature (Ts) and vegetation index data, have strong potential for operational applications. However, very few studies have tested such remote sensing methods over Mediterranean forests. One difculty 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. To ll the gap, this study investigates the impact of the solar zenith angle as a proxy of trees casting shadows on the satellite-retrieved water decit index (WDI). The WDI method is implemented using Landsat-7 and Landsat-8 data over a 21 km 2 area partially covered by a holm oak forest in South-eastern France. The study period extends from May to September for 7 successive years from 2015 to 2021 and the results are evaluated using the evaporative fraction (EF) measured in situ at the Puechabon site. As a rst step, a correction method of WDI for shadow effects is developed at the Puechabon site by correlating the WDI error (WDI minus in situ 1-EF) with the solar zenith angle. The correction signicantly improves the correlation between satellite WDI and in situ 1-EF from 0.32 to 0.62 and increases the slope of the linear regression from 0.21 to 0.63. As a second step, a calibration approach of the correction method is proposed without relying on in situ measurements, by evaluating the linear relationship (on a pixelby-pixel basis) between remotely sensed WDI and the solar zenith angle in dry conditions. The spatially-corrected WDI still shows a signicant improvement compared to the non-corrected WDI with a correlation with in situ 1-EF increasing from 0.32 to 0.57 and a slope of the linear regression increasing from 0.21 to 0.50. This approach offers a nonparametric while efcient technique to account for trees casting shadow effects on the remotely sensed stress indices of forests.

Index Terms-thermal data, hydric stress, forest, evapotranspiration, Mediterranean ecosystems

I. INTRODUCTION

M EDITERRANEAN forests have dealt with high levels of hydric stress for long periods [START_REF] Limousin | Drought acclimation of quercus ilex leaves improves tolerance to moderate drought but not resistance to severe water stress[END_REF], [START_REF] Ollivier | An evapotranspiration model driven by remote sensing data for assessing groundwater resource in karst watershed[END_REF]. However, global warming is a threat to most of forests all around the world [START_REF] Seidl | Increasing forest disturbances in europe and their impact on carbon storage[END_REF]- [START_REF] Peters | Effects of global warming on forests[END_REF] and particularly in the Mediterranean area, Manuscript created January, 2023; Victor Penot, Olivier Merlin, are members of Centre Etude Spatial de la BIOsphère, 31000 Toulouse France, e-mail: victor.penot@iut-tlse3.fr, olivier.merlin1@univ-tlse3.fr which is bound to be affected by increasingly frequent and severe droughts [START_REF]Climate Change, Climate Change 2013 -The Physical Science Basis: Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF], [START_REF] Tramblay | Challenges for drought assessment in the mediterranean region under future climate scenarios[END_REF]. 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 to develop long-term forest/agricultural policies like the choice of species adapted to warmer and/or dryer conditions and mitigation management strategies [START_REF] González-Dugo | Long-term water stress and drought assessment of mediterranean oak savanna vegetation using thermal remote sensing[END_REF], and also for reghters to assess locally the short-term water dynamics of fuels [START_REF] Nolan | Linking forest ammability and plant vulnerability to drought[END_REF], [START_REF] Littell | A review of the relationships between drought and forest re in the united states[END_REF].

Remote sensing (RS) technology offers a good cost-benet trade-off to derive the evapotranspiration (ET) and water stress of ecosystems over wide areas on a daily/weekly basis. Different ET remote sensing methods have been developed for this purpose and are generally based on optical/thermal data. There exists two broad categories: the energy balance methods and the contextual methods. The former are based on solving the physically-based energy balance equations [START_REF] Boulet | The sparse model for the prediction of water stress and evapotranspiration components from thermal infra-red data and its evaluation over irrigated and rainfed wheat[END_REF], [START_REF] Su | The surface energy balance system (sebs) for estimation of turbulent heat uxes[END_REF]. They are generally more complicated to put into operational practice, due to their ne parametrization and their possibly large number of input parameters [START_REF] Delogu | Evaluation of the sparse dual-source model for predicting water stress and evapotranspiration from thermal infrared data over multiple crops and climates[END_REF], [START_REF] Chirouze | Intercomparison of four remote-sensing-based energy balance methods to retrieve surface evapotranspiration and water stress of irrigated elds in semi-arid climate[END_REF]. The latter are based on semi-empirical interpretations of the observed spatial correlations between land surface temperature (Ts) and fractional green vegetation cover (f vg ) [START_REF] Moran | Estimating crop water decit using the relation between surface-air temperature and spectral vegetation index[END_REF]- [START_REF] Jiang | A satellite-based daily actual evapotranspiration estimation algorithm over south orida[END_REF] and/or albedo [START_REF] Roerink | S-sebi: A simple remote sensing algorithm to estimate the surface energy balance[END_REF], [START_REF] Merlin | An original interpretation of the wet edge of the surface temperature-albedo space to estimate crop evapotranspiration (seb-1s), and its validation over an irrigated area in northwestern mexico[END_REF] 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 [START_REF] Tang | An application of the ts-vi triangle method with enhanced edges determination for evapotranspiration estimation from modis data in arid and semiarid regions: Implementation and validation[END_REF], [START_REF] Wang | Estimation of evaporative fraction from a combination of day and night land surface temperatures and ndvi: A new method to determine the priestley-taylor parameter[END_REF].

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 [START_REF] Moran | Estimating crop water decit using the relation between surface-air temperature and spectral vegetation index[END_REF], [START_REF] Jiang | A methodology for estimation of surface evapotranspiration over large areas using remote sensing observations[END_REF]. 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 [START_REF] Chuvieco | Combining ndvi and surface temperature for the estimation of live fuel moisture content in forest re danger rating[END_REF]. The Water Decit Index (WDI), a proxy of surface water stress, is notably determined from the relative position of a point in the Ts-f vg space [START_REF] Moran | Estimating crop water decit using the relation between surface-air temperature and spectral vegetation index[END_REF], 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 eld measurements [START_REF] Moran | Estimating crop water decit using the relation between surface-air temperature and spectral vegetation index[END_REF], [START_REF] Collins | Temporal and spatial changes in grassland transpiration detected using landsat tm and etm + imagery[END_REF]- [START_REF] Virlet | Stress indicators based on airborne thermal imagery for eld phenotyping a heterogeneous tree population for response to water constraints[END_REF]. However, to date, thermal remote sensing methods have seldom been tested over Mediterranean forests. In an early study, [START_REF] Vidal | Evaluating forest re hazard with a landsat tm derived water stress index[END_REF] used WDI as an index to predict the re area/occurrence during the summer of 1990 in Southern France. However, the lack of site ux 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 inuence on the measured Ts [START_REF] Prata | Thermal remote sensing of land surface temperature from satellites: Current status and future prospects[END_REF]. However, it is difcult to model and evaluate accurately the cast shadows impact in complex ecosystems [START_REF] Prata | Thermal remote sensing of land surface temperature from satellites: Current status and future prospects[END_REF]. The use of radiative transfer and/or geometric projection models supports these results, if a mock-up can be designed. [START_REF] Guillevic | Directional viewing effects on satellite land surface temperature products over sparse vegetation canopies-a multisensor analysis[END_REF]- [START_REF] Luquet | Using directional tir measurements and 3d simulations to assess the limitations and opportunities of water stress indices[END_REF] 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 eld of view on Ts retrievals, and particularly through the contribution of shaded and sunlit soil.

To derive ET or water stress proxies in non-homogeneous canopies, [START_REF] Sepulcre-Cantó | Detection of water stress in an olive orchard with thermal remote sensing imagery[END_REF], [START_REF] Testi | Crop water stress index is a sensitive water stress indicator in pistachio trees[END_REF] 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 [START_REF] Testi | Crop water stress index is a sensitive water stress indicator in pistachio trees[END_REF]- [START_REF] Egea | Assessing a crop water stress index derived from aerial thermal imaging and infrared thermometry in super-high density olive orchards[END_REF] or aircraft/UAV-based [START_REF] Sepulcre-Cantó | Detection of water stress in an olive orchard with thermal remote sensing imagery[END_REF], [START_REF] Matese | Estimation of water stress in grapevines using proximal and remote sensing methods[END_REF]- [START_REF] Lu | Effect of the shadow pixels on evapotranspiration inversion of vineyard: A high-resolution uav-based and ground-based remote sensing measurements[END_REF] thermal camera. One solution consists in removing in the Ts image, the pixels corresponding to shaded soil/leaves to isolate sunlit leaves [START_REF] Sepulcre-Cantó | Detection of water stress in an olive orchard with thermal remote sensing imagery[END_REF], [START_REF] Matese | Estimation of water stress in grapevines using proximal and remote sensing methods[END_REF]- [START_REF] Lu | Effect of the shadow pixels on evapotranspiration inversion of vineyard: A high-resolution uav-based and ground-based remote sensing measurements[END_REF]. It relies on the assumption that sunlit leaves are more likely to be stressed than shaded ones [START_REF] Gardner | Infrared thermometry and the crop water stress index. ii. sampling procedures and interpretation[END_REF] and that they may provide meaningful information about the vegetation water status/ux. Separating shaded from sunlit soil/leaves requires a very high-resolution (centimetric) sensor. Such a classication can be done manually [START_REF] García-Tejero | Thermal imaging to monitor the crop-water status in almonds by using the non-water stress baselines[END_REF], [START_REF] Matese | Estimation of water stress in grapevines using proximal and remote sensing methods[END_REF], [START_REF] Lu | Effect of the shadow pixels on evapotranspiration inversion of vineyard: A high-resolution uav-based and ground-based remote sensing measurements[END_REF] or by automatic detection using supervised classication [START_REF] Lu | Effect of the shadow pixels on evapotranspiration inversion of vineyard: A high-resolution uav-based and ground-based remote sensing measurements[END_REF], [START_REF] Leinonen | Combining thermal and visible imagery for estimating canopy temperature and identifying plant stress[END_REF], density based methods [START_REF] Liu | Thermal remote sensing of plant water stress in natural ecosystems[END_REF], [START_REF] Egea | Assessing a crop water stress index derived from aerial thermal imaging and infrared thermometry in super-high density olive orchards[END_REF], [START_REF] Wang | Automated canopy temperature estimation via infrared thermography: A rst step towards automated plant water stress monitoring[END_REF], clustering (Kmeans) and post-classication [START_REF] Poblete | Automatic coregistration algorithm to remove canopy shaded pixels in uav-borne thermal images to improve the estimation of crop water stress index of a drip-irrigated cabernet sauvignon vineyard[END_REF], high resolution elevation thresholds [START_REF] Matese | Estimation of water stress in grapevines using proximal and remote sensing methods[END_REF] or vegetation index thresholds [START_REF] Sepulcre-Cantó | Detection of water stress in an olive orchard with thermal remote sensing imagery[END_REF]. However, due to the geometric misregistration of visible/near-infrared and thermal images at such a ne resolution, soil and shaded leaves background may contaminate the thermal data of an area identied as sunny and vegetated [START_REF] Poblete | Automatic coregistration algorithm to remove canopy shaded pixels in uav-borne thermal images to improve the estimation of crop water stress index of a drip-irrigated cabernet sauvignon vineyard[END_REF], [START_REF] Wang | Automated canopy temperature estimation via infrared thermography: A rst step towards automated plant water stress monitoring[END_REF]. 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 reectances outside the training database may decrease the accuracy of the segmentation [START_REF] Leinonen | Combining thermal and visible imagery for estimating canopy temperature and identifying plant stress[END_REF].

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 eld 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 [START_REF] Testi | Crop water stress index is a sensitive water stress indicator in pistachio trees[END_REF], [START_REF] Liu | Thermal remote sensing of plant water stress in natural ecosystems[END_REF]. Note that this method strongly depends on the solar zenith angle and on the remaining ground and shaded leaves contamination [START_REF] Testi | Crop water stress index is a sensitive water stress indicator in pistachio trees[END_REF], [START_REF] Liu | Thermal remote sensing of plant water stress in natural ecosystems[END_REF].

Rather than extracting sunlit leaves temperature from the thermal image or avoiding shadow effect in the sensor eld of view, [START_REF] Luquet | Using directional tir measurements and 3d simulations to assess the limitations and opportunities of water stress indices[END_REF] suggested a third option, that is to correct the WDI method to take into account directional effects. This correction involved redening the wet and dry edges of the Ts-f vg space by using a 3D model that was run in each solar and observation angle. However this method uses realistic and efcient 3D 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 m 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 satellitederived WDI for trees casting shadow effects, by relying on the solar zenith angle (θ S ). Specically, this paper aims at (i) quantitatively assessing the impact of θ S on WDI and (ii) 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 ux 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.

II. MATERIAL AND METHODS

A. Puechabon site and in situ data 1) Puechabon site:

The Puechabon site (43 • 44'29"N;3 • 35'46"E, 270m a.s.l.) is located in Southern France, 35 km northwest of Montpellier [START_REF] Limousin | Long-term transpiration change with rainfall decline in a mediterranean quercus ilex forest[END_REF]. This site is representative of Mediterranean evergreen broadleaf forests with dense coppice. The holm oak (Quercus Ilex L.) is the dominant tree specie, with a mean tree height of 5.5 m. The understory is composed of sparse shrubs of 2-m height [START_REF] Ollivier | An evapotranspiration model driven by remote sensing data for assessing groundwater resource in karst watershed[END_REF]. The site was built in 1984 to study Mediterranean ecosystems and specially the holm oaks response to climate change and to severe droughts [START_REF] Limousin | Drought acclimation of quercus ilex leaves improves tolerance to moderate drought but not resistance to severe water stress[END_REF], [START_REF] Limousin | Long-term transpiration change with rainfall decline in a mediterranean quercus ilex forest[END_REF].

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, Figure 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 inter-annual 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 5 days over 35 • c between June and September. 2016 also experienced a very dry summer The soil at Puechabon site is silty claim loam with a rock and stone volumetric content of 75% in the rst 50 cm and 90% in the whole prole [START_REF] Rambal | Drought controls over conductance and assimilation of a mediterranean evergreen ecosystem: scaling from leaf to canopy[END_REF]. It is poorly developed and does not store water well [START_REF] Ollivier | An evapotranspiration model driven by remote sensing data for assessing groundwater resource in karst watershed[END_REF]. However, it was shown that Quercus Illex was able to extract water deeper than 4 m during drought [START_REF] Rambal | Drought controls over conductance and assimilation of a mediterranean evergreen ecosystem: scaling from leaf to canopy[END_REF], [START_REF] Allard | Seasonal and annual variation of carbon exchange in an evergreen mediterranean forest in southern france[END_REF].

2) Flux and meteorological data: The station of Puechabon is equipped with an eddy covariance ux measurement system that measures meteorological data and energy uxes. The ux tower is 10-m height and equipped with a Sonic anemometer (Solent R3A Gill) and an infra red gas analyser (IRGA, LI-6262,Li-COR) for latent and sensible heat ux measurements. The sampling is done at a rate of 21 Hz and recorded every 30 minutes. [START_REF] Allard | Seasonal and annual variation of carbon exchange in an evergreen mediterranean forest in southern france[END_REF] studied the measurement footprint and its inuence on the ux data in all weather conditions. It was shown that the measured uxes 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 [START_REF] Allard | Seasonal and annual variation of carbon exchange in an evergreen mediterranean forest in southern france[END_REF], [START_REF] Maruyama | Reciprocal analysis of sensible and latent heat uxes in a forest region using single height temperature and humidity based on the bowen ratio concept[END_REF]. Records are done every 30 minutes. Flux and meteorological data from 2015 to 2021 were downloaded from the European Fluxes Database Cluster.

Regarding the monitoring and modelling of evapotranspiration/hydric stress from thermal data, few studies have been undertaken at the Puechabon site. [START_REF] Jiang | Multi-scale evaluation of global gross primary productivity and evapotranspiration products derived from breathing earth system simulator (bess)[END_REF]- [START_REF] Ryu | On the temporal upscaling of evapotranspiration from instantaneous remote sensing measurements to 8-day mean daily-sums[END_REF] evaluated the 1km resolution 8-day ET MODIS product. More recently, [START_REF] Ollivier | An evapotranspiration model driven by remote sensing data for assessing groundwater resource in karst watershed[END_REF] 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 dened as the ratio of latent heat to the available energy [START_REF] Jiang | A methodology for estimation of surface evapotranspiration over large areas using remote sensing observations[END_REF]:

EF = LE Rn -G (1)
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:

EF = LE LE + H (2)
with H being the sensible heat (W/m 2 ). [START_REF] Allard | Seasonal and annual variation of carbon exchange in an evergreen mediterranean forest in southern france[END_REF] 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 [START_REF] Jiang | A satellite-based daily actual evapotranspiration estimation algorithm over south orida[END_REF], [START_REF] Cragoa | Daytime evaporation and the selfpreservation of the evaporative fraction and the bowen ratio[END_REF], [START_REF] Jiang | Uncertainties in latent heat ux measurement and estimation: implications for using a simplied approach with remote sensing data[END_REF] 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. [START_REF] Jiang | A satellite-based daily actual evapotranspiration estimation algorithm over south orida[END_REF]. Hence the daytime EF (EFd) can be estimated by means of the Bowen Ratio :

EFd = 1 1 + β (3)
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.

B. Remote sensing data and study area extent 1) Remote sensing data: Landsat Enhanced Thematic Mapper Landsat-7 (ETM+) and Landsat-8 Operational Land Imager (OLI) and Thermal Infrared (TIR) both provide re-ectance 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 reectances 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) [START_REF] Schmidt | Landsat ecosystem disturbance adaptive processing system (ledaps) algorithm description[END_REF] and Landsat-8 data with the Land Surface Reectance Code (LaSRC) algorithm (Version 1.5.0) [START_REF] Vermote | Preliminary analysis of the performance of the landsat 8/oli land surface reectance product[END_REF].

After 2017, Landsat-7's orbit drifted from its original schedule to an earlier overpass time. Changes in reectance and thermal data are limited until 2020 [START_REF] Qiu | Can landsat 7 preserve its science capability with a drifting orbit?[END_REF]. 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 data set.

f vg is estimated using the expression of [START_REF] Carlson | On the relation between ndvi, fractional vegetation cover, and leaf area index[END_REF] :

f vg =  NDVI -NDVI min NDVI max -NDVI min  2 (4)
with NDVI being the Normalized Difference Vegetation Index dened as the difference of near-infrared and red reectances 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 m and 100 m respectively. Both data sets were downsampled at 30 m resolution and provided over the same 30 m resolution grid of Landsat reectances. They were processed with the Landsat surface temperature algorithm (Version 1.3.0) to procude Ts maps from raw thermal radiances [START_REF] Cook | Atmospheric compensation for a landsat land surface temperature product[END_REF], [START_REF] Cook | Development of an operational calibration methodology for the landsat thermal data archive and initial testing of the atmospheric compensation component of a land surface temperature (lst) product from the archive[END_REF].

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 [START_REF] Boulet | Monitoring water stress using time series of observed to unstressed surface temperature difference[END_REF], [START_REF] Karnieli | Use of ndvi and land surface temperature for drought assessment: Merits and limitations[END_REF]. Meteorological conditions must also be relatively homogeneous within the study area [START_REF] Chirouze | Intercomparison of four remote-sensing-based energy balance methods to retrieve surface evapotranspiration and water stress of irrigated elds in semi-arid climate[END_REF]. Thus the study area extent is dened 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 Figure 2. On its northern part, it is delimited by the narrow Herault river and a plateau mainly covered by closed deciduous forests, were the ux site is set. Note that according to the data provider denition (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 Figure 2 a andb). The minimum elevation of 40 m is reached in the southern part while the maximum elevation of 441 m is located in the northern part (see Figure 2 b). Landsat remote sensing data were thus extracted over this area for the whole time series.

C. WDI method and image selection

The methodology to set the dry and wet edges of the Tsf vg space (II-C1 and II-C2) and to estimate the WDI (II-C3) is rst presented. Then, a quality control of the input data is undertaken to automatically lter out the dates when meteorological conditions weaken the WDI application (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 dened as the linear regression of the 99% quantiles of Ts evaluated in equal-length f vg bins against f vg as in [START_REF] Virlet | Stress indicators based on airborne thermal imagery for eld phenotyping a heterogeneous tree population for response to water constraints[END_REF] (see Figure 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 [START_REF] Jiang | A satellite-based daily actual evapotranspiration estimation algorithm over south orida[END_REF]. According to [START_REF] Colaizzi | Twosource energy balance model estimates of evapotranspiration using component and composite surface temperatures[END_REF], Ts is approximatively equal to aerodynamic temperature, so H is as a rst guess proportional to the difference between Ts and air temperature (Tair) [START_REF] Moran | Estimating crop water decit using the relation between surface-air temperature and spectral vegetation index[END_REF]. 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 Figure 3). This assumption has been successfully assessed in other previous studies undertaken in semi-arid environments [START_REF] Merlin | An original interpretation of the wet edge of the surface temperature-albedo space to estimate crop evapotranspiration (seb-1s), and its validation over an irrigated area in northwestern mexico[END_REF], [START_REF] Zhang | Surface soil water content estimation from thermal remote sensing based on the temporal variation of land surface temperature[END_REF], [START_REF] Amazirh | Retrieving surface soil moisture at high spatio-temporal resolution from a synergy between sentinel-1 radar and landsat thermal data: A study case over bare soil[END_REF].

3) WDI calculation: Figure 3 illustrates how WDI is computed in the Ts-f vg space from dry edge and wet edge distances [START_REF] Moran | Estimating crop water decit using the relation between surface-air temperature and spectral vegetation index[END_REF] :

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 denition, WDI is the complement in 1 of EF, i.e., WDI = 1 -EF. 4) Date selection: First, the use of thermal remote sensing 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 [START_REF] Chirouze | Intercomparison of four remote-sensing-based energy balance methods to retrieve surface evapotranspiration and water stress of irrigated elds in semi-arid climate[END_REF], [START_REF] Karnieli | Use of ndvi and land surface temperature for drought assessment: Merits and limitations[END_REF]. Hence, the study period is dened to include the hottest months from May to September for each year from 2015 to 2021.

Second, the presence of clouds, cloud shadows or recent rainfalls is likely to decrease Ts under the wet edge, dened 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 [START_REF] Chirouze | Intercomparison of four remote-sensing-based energy balance methods to retrieve surface evapotranspiration and water stress of irrigated elds in semi-arid climate[END_REF]. Therefore, an approach is proposed to automatically remove such dates. Only the pixels qualied as clear in the quality masks provided with Landsat-7/8 reectance products are used during this ltering step. [START_REF] Zhang | Surface soil water content estimation from thermal remote sensing based on the temporal variation of land surface temperature[END_REF] suggests 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 (qualied 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 lter out cloudy conditions. This limit quanties the clear or low cloudy conditions that are commonly used in the reference previous studies [START_REF] Moran | Estimating crop water decit using the relation between surface-air temperature and spectral vegetation index[END_REF], [START_REF] Merlin | An original interpretation of the wet edge of the surface temperature-albedo space to estimate crop evapotranspiration (seb-1s), and its validation over an irrigated area in northwestern mexico[END_REF], [START_REF] Jiang | Estimation of surface evaporation map over southern great plains using remote sensing data[END_REF], [START_REF] Tang | A moisturebased triangle approach for estimating surface evaporative fraction with time-series of remotely sensed data[END_REF].

D. Correction method of WDI for solar zenith angle 1) General approach: The correction method proposed in this paper 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 Figure 4.

Once the WDI is classically estimated (as described in section II-C3, see Figure 4 upper part), a relation is built between WDI and θ S :

∆WDI = a(θ S -b) (6)
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 ux 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 Figure 4 left part); 2) if ux 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,M in . This correction strategy is called self-calibrated correction (see Figure 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 Figure 4 lower part) to evaluate the correction performance. Each component of the correction method is described below.

2) Identication 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 signicant 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 dene the Very Dry dataset, as the dates of our dataset whose 15-day cumulative rain (P 15d ) is under the rst quartile of the 15day cumulative rain of all days during the 2015-2021 study period (q 0.25 (P 15d ) = 7mm). 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 dene the Dry dataset (Wet dataset), as the dates of our dataset whose 15-day cumulative rain (P 15d ) is under (over) the median of the 15day cumulative rain of all days during the 2015-2021 study period (q 0.5 (P 15d ) = 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 ux data are available, the difference between WDI and 1-EF regressed over θ S can be used to correct the WDI (see Equation 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 site-calibrated correction equation is :

WDI c,site = WDI -a site (θ S -b site ) (7)
4) Self-calibrated correction: In practice, ux data are rarely available and in general, they cannot be used for calibration purposes of Equation 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. rst by calibrating the slope a of Equation 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 inuenced by cast shadow effects and to calibrate the slope a = a self of Equation 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 ltering 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,M in ) 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 [START_REF] Guillevic | Directional viewing effects on satellite land surface temperature products over sparse vegetation canopies-a multisensor analysis[END_REF], [START_REF] Guillevic | Thermal infrared radiative transfer within three-dimensional vegetation covers[END_REF], 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,M in and that it can only decrease WDI when θ S increases. Therefore, the selfcalibrated correction equation is :

WDI c,self = WDI -a self (θ S -θ S,M in ) (8)

5) Evaluation of calibration and correction strategies:

To assess the quality and the improvements brought by both (site-and 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), coefcient 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 selfcalibrated correction method is proposed. A. Site-calibrated correction of WDI for solar zenith angle 1) Site calibration: As mentioned previously, Equation 1 is calibrated with the Very Dry data set (including the 25% driest dates) with the objective of identifying the dates when θ S signicantly drives the error ∆WDI = (WDI -EF + 1).

To assess the usefulness of this assumption, Figures 5 a andb plot the error ∆WDI versus θ S for All dataset and Very Dry dataset respectively. Table II summarizes the quality statistics of site-calibrated WDI when All dataset and Very Dry dataset is used for calibration separately. Figures 6 a andb plot the non-corrected 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 non-corrected WDI (see Table II and Figures 6 a andb). The slope and correlation coefcient 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 signicantly 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 [START_REF] Tang | A moisturebased triangle approach for estimating surface evaporative fraction with time-series of remotely sensed data[END_REF]- [START_REF] Nishida | An operational remote sensing algorithm of land surface evaporation[END_REF].

Even if the calibration regression using All dataset is significant (determination coefcient 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 Figures 5 a andb). These results suggest that θ S signicantly 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 Figure 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 Figure 5 b) that quanties the underestimation of 1-EF by WDI with the increase of θ S across the seasons. It also provides an intercept b site = 25.6 • (see Figure 5 b) that is near the minimum solar zenith angle (θ S,M in = 25.7 • ) reached in late June at the Puechabon site. This value implies that ∆WDI is 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 [START_REF] Guillevic | Directional viewing effects on satellite land surface temperature products over sparse vegetation canopies-a multisensor analysis[END_REF], [START_REF] Guillevic | Thermal infrared radiative transfer within three-dimensional vegetation covers[END_REF].

In brief, the site calibration provides an equation that efciently 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 non-corrected 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 coefcient) from 0.07 to 0.84 (from 0.08 to 0.59) for the non-corrected and site-calibrated corrected WDI, respectively. In parallel the bias decreases from -0.15 to 0.00 for non-corrected 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 signicantly 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 [START_REF] Anderson | Mapping daily evapotranspiration at eld to continental scales using geostationary and polar orbiting satellite imagery[END_REF] 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: Figures 7 (a andb) plots non-corrected 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 Figures 7 a andb). The slope decreases from 0.30 to 0.23 and the correlation coefcient from 0.53 to 0.28 for non-corrected WDI and site-calibrated corrected WDI respectively. Particularly, the correction for the Wet dataset tends to over-correct WDI. This year is the rainiest of the 7 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, [START_REF] Limousin | Drought acclimation of quercus ilex leaves improves tolerance to moderate drought but not resistance to severe water stress[END_REF] actually measured in 2015 the maximum of minimum pre-dawn summer water potentials -a proxy of plant and soil water limitation [START_REF] Rambal | Drought controls over conductance and assimilation of a mediterranean evergreen ecosystem: scaling from leaf to canopy[END_REF] -and the rainiest summer (June-August) over the 2003-2017 summers.

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 inter-annual 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 ux data. The calibration of the slope a self of Equation 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 efcient 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 non-corrected and self-calibrated corrected WDI, respectively. In contrast, using the Very Dry dataset the slope (and R) signicantly increases from 0.21 (0.32) to 0.50 (0.57) for the non-corrected 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 efcient correction on WDI (see Figures 5 c andd), 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 Figures 5 b andd). Even if a self is lower (in absolute value) than a site , the previous discussion in section III-A1 remains valid.

Only 7 years of data have been used for the self-calibration, 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 signicant results. In fact, the upper edge of the WDI-θ S space is strongly inuenced 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 dene the nal calibration Equation 8, the intercept is set to the minimum θ S between May and Septemberb self = θ S,M in = 25.7 • -which is close to the sitecalibrated 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,M in 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 Figure 6 and Table III indicate that the selfcalibrated corrected WDI outperforms the non-corrected WDI when evaluated on the Dry dataset. In particular, the selfcalibrated correction improves the slope and correlation co-efcient (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 non-corrected and self-calibrated corrected WDI respectively. However, the selfcalibrated 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 behaviour 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 Figures 7 a andc. 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 conrms 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 Equation 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 behaviour of the corrected WDI and the correction parameters (a self ) at the pixel scale within the study area. The study area (see Figure 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 5x5 pixels split window to take into account the actual spatial resolution of thermal Landsat data (60 m 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 selfcalibrated correction equation can nally be applied to correct WDI for cast shadow effects over the entire study area.

Figure 8 shows maps and density histograms of WDI and WDI c,self respectively, for data on the 22 th , 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 ux site is 0.88, site-calibrated corrected WDI c,self is 0.81 while non-corrected WDI is 0.67. It corresponds to an error reduction of 77% by correcting the WDI for trees casting shadow effects.

In Figure 8 a, the non-corrected WDI ranges from 0.2 to 1. In fact, the difference between the minimum Ts (299.7K) and Tair at 10:30 a.m. (297.7K, wet edge value) is 2 K, which implies strictly positive WDI values. The spatial distribution of WDI generally reects 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 efciently deal with low levels of soil water availability, a lower value of WDI lower than in vineyard/orchards would be expected [START_REF] Limousin | Drought acclimation of quercus ilex leaves improves tolerance to moderate drought but not resistance to severe water stress[END_REF]. Further north, the slope along the Herault river witnesses particularly low values of WDI.

Figure 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 nearzero 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 Figure 9 a and density histograms of a self values per dominant land use/cover class are provided in Figure 9 b. Values of a self are strongly related to the land cover map (see Figure 2 a) and to the elevation map (Figure 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 Figures 9 a andb). 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 ux tower is located, the Quercus Illex dominated forest is less open and witnesses more complex landscapes than orchards/vineyards, which implies a signicant decrease in the mean a self equal to -0.008 (see Figure 9 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 satellitederived WDI for trees casting shadow effects over large areas. To ll 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 (i) to quantitatively assess the impact of θ S on WDI during dry periods (high values of hydric stress), (ii) to build a self-calibration method that corrects WDI on a pixel basis from θ S and satellite data solely for operational use and (iii) 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 Southerneast of France. As a rst 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 coefcient (R) increases from 0.32 to 0.62 for non-corrected 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 (non-corrected) WDI respectively. Therefore, the self-calibrated method provides an operational workow 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 selfand 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 behaviour of forests, especially when the hydric stress is relatively low (e.g. thermal inertia [START_REF] Gaudio | Impact of tree canopy on thermal and radiative microclimates in a mixed temperate forest: A new statistical method to analyse hourly temporal dynamics[END_REF], [START_REF] Haverd | Air and biomass heat storage uxes in a forest canopy: Calculation within a soil vegetation atmosphere transfer model[END_REF]). 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 3D radiative transfer model like DART [START_REF] Gastellu-Etchegorry | Discrete anisotropic radiative transfer (dart 5) for modeling airborne and satellite spectroradiometer and lidar acquisitions of natural and urban landscapes[END_REF]. However, this would require the combined used of energy balance models to precisely solve the temperature of soil and vegetation components [START_REF] Malbéteau | Normalizing land surface temperature data for elevation and illumination effects in mountainous areas: A case study using aster data over a steep-sided valley in morocco[END_REF], [START_REF] Mwangi | Assessment of an extended sparse model for estimating evapotranspiration from directional thermal infrared data[END_REF]. Last but not least, this would also require a precise and exhaustive 3D 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 signicant viewing angle directional effects [START_REF] Buffet | The TIR instrument on TRISHNA satellite: a precursor of high resolution observation missions in the thermal infrared domain[END_REF], [START_REF] Lagouarde | The Indian-French Trishna Mission: Earth Observation in the Thermal Infrared with High Spatio-Temporal Resolution[END_REF]. Thus correcting for directional effects is especially relevant for the future operational use of these valuable data.

Fig. 1 .

 1 Fig. 1. Monthly cumulative rain and monthly mean temperature between 2015 and 2021 at the site of Puechabon.

Fig. 2 .

 2 Fig. 2. Land use and land cover map a) and elevation map b) over the study area (Source : BDtopo and BDAlti of IGN).

Fig. 3 .

 3 Fig. 3. Ts-fvg space and dry and wet edges plotted for Landsat-7 data collected over the study area on the 20 th , August 2015.

Fig. 4 .

 4 Fig. 4. Schematic diagram of the overall methodology to correct WDI for solar zenith angle effects, to calibrate this correction from in situ (site calibration, left) or from remote sensing (self-calibration, right) data solely and to evaluate the correction efciency at the Puechabon site (bottom).

Fig. 5 .

 5 Fig. 5. Scatter plot from May to September during the years 2015-2021 of a) (θ S ,WDI-EF+1) for the All dataset and the linear regression b) (θ S ,WDI-EF+1) for the Very Dry dataset and the linear regression ; c) (θ S ,WDI-EF+1) for All dataset and the linear upper hull d) (θ S ,WDI) for the Very Dry dataset and the linear upper hull.

Fig. 6 .

 6 Fig. 6. Scatter plot of a) non-corrected WDI vs 1-EF ; b) WDI c,site vs 1-EF; c) WDI c,self vs 1-EF. The Very Dry dataset is used for both site and self-calibration cases.

  This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2023.3288360 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

Fig. 7 .

 7 Fig. 7. Per year scatter scatter plot of a) non-corrected WDI vs 1-EF ; b) WDI c,site vs 1-EF; c) WDI c,self vs 1-EF.

Fig. 8 .

 8 Fig. 8. Maps and density histograms per dominant land cover class of a) and c) non corrected WDI, and b) and d) self-calibrated corrected WDI c,self for Landsat-7 data on the 22 th July, 2016.

Fig. 9 .

 9 Fig. 9. a) Map of the slope a self of the self-calibrated correction equation over the study area and b) density histograms of a self per land use/cover classe.

  a and b). The northern valley along the Herault river is characterized by very low a self values, consistent with topographic effects. Note that This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2023.3288360 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

TABLE I SEASONAL

 I CUMULATIVE RAIN AND THERMAL ANOMALIES FOR EACH YEAR SEPARATELY.

	Year 2015 2016 2017 2018 2019 2020 2021	Cumulative rain February-May (mm) 177 363 290 588 169 256 291	Cumulative rain June-September (mm) 267 271 74 98 142 211 190	Cumulative rain July-August (mm) 166 17 13 62 63 79 83	Mean Tair June-September ( • c) 21.8 22 21.7 22.8 22.6 22 21.7	Cumulative number of days Tair> 30 • c June-September 32 34 34 38 51 41 22	Cumulative number of days Tair> 35 • c June-September 2 0 3 5 3 1 1
	(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 remote sensing.

TABLE II STATISTICAL

 II RESULTS IN TERMS OF MAE, RMSE, R, SLOPE AND BIAS OF THE LINEAR REGRESSION BETWEEN SATELLITE WDI AND in situ 1-EF FOR DIFFERENT CORRECTION AND CALIBRATION STRATEGIES. THE ALL DATASET IS USED FOR EVALUATION.

	Correction method None Site Site Self Self	Calibration dataset All All Very Dry All Very Dry	MAE 0.14 0.10 0.10 0.13 0.11	RMSE 0.17 0.12 0.12 0.15 0.12	R 0.32 0.57 0.62 0.44 0.57	Slope 0.21 0.50 0.63 0.31 0.50	Bias -0.10 0.00 0.02 -0.07 -0.02

TABLE III STATISTICAL

 III RESULTS IN TERMS OF MAE, RMSE, R, SLOPE AND BIAS OF THE LINEAR REGRESSION BETWEEN SATELLITE WDI AND in situ 1-EF FOR DIFFERENT CORRECTION AND EVALUATION STRATEGIES.

	Correction method None Site Self None Site Self	Evaluation dataset Wet Wet Wet Dry Dry Dry	MAE 0.09 0.10 0.10 0.17 0.10 0.11	RMSE 0.11 0.14 0.12 0.20 0.11 0.12	R 0.21 0.36 0.35 0.08 0.59 0.48	Slope 0.14 0.37 0.29 0.07 0.84 0.59	Bias -0.02 0.08 0.04 -0.15 0.00 -0.05

This article has been accepted for publication in IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. This is the author's version which has not been fully e content may change prior to final publication. Citation information: DOI 10.1109/JSTARS.2023.3288360 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

V. ACKNOWLEDGMENTS

This study was funded in part by the Horizon 2020 ACCWA project (grant agreement # 823965) in the context of Marie Sklodowska-Curie research and innovation staff exchange (RISE) program. The authors would like to thank the principal investigator of the FLUXNET Puechabon site, Jean-Marc Limousin.