Regional Equivalent Water Thickness Modeling from Remote Sensing across a Tree Cover / LAI Gradient in Mediterranean Forests of Northern Tunisia

The performance of vegetation indexes derived from moderate resolution imaging spectroradiometer (MODIS) sensors is explored for drought monitoring in the forests of Northern Tunisia; representing a transition zone between the Mediterranean Sea and the Sahara Desert. We investigated the suitability of biomass and moisture vegetation indexes for vegetation water content expressed by the equivalent water thickness (EWT) in a Mediterranean forest ecosystem with contrasted water budgets and desiccation rates. We proposed a revised EWT at canopy level (EWTCAN) based on weekly field measurements of fuel moisture in seven species during the 2010 dry period, considering the mixture of plant functional types for water use (trees, shrubs and herbaceous layers) and a varying vegetation cover. MODIS vegetation indexes computed and smoothed over the dry season were highly correlated with the EWTCAN. The performances of moisture indexes Normalized Difference Infrared Index (NDII6 and NDII7) and Global Moisture Vegetation Index (GVMI6 and GVMI7) were comparable, whereas for biomass vegetation indexes, Normalized Difference OPEN ACCESS Remote Sens. 2015, 7 1938 Vegetation Index (NDVI), Modified Soil Adjusted Vegetation Index (MSAVI) and Adjusted Normalized Difference Vegetation Index (ANDVI) performed better than Enhanced Vegetation Index (EVI) and Soil Adjusted Vegetation Index (SAVI). We also identified the effect of Leaf Area Index (LAI) on EWTCAN monitoring at the regional scale under the tree cover/LAI gradient of the region from relatively dense to open forest. Statistical analysis revealed a significant decreasing linear relationship; indicating that for LAI less than two, the greater the LAI, the less responsive are the vegetation indexes to changes in EWTCAN; whereas for higher LAI, its influence becomes less significant and was not considered in the inversion models based on vegetation indexes. The EWTCAN time-course from LAI-adapted inversion models based on significantly-related vegetation indexes to EWTCAN showed close profiles resulting from the inversion models using NDVI, ANDVI, MSAVI and NDII6 applied during the dry season. The developed EWTCAN model from MODIS vegetation indexes for the study region was finally tested for its ability to capture the topo-climatic effects on the seasonal and the spatial patterns of desiccation/rewetting for keystone periods of Mediterranean vegetation functioning. Implications for further use in scientific developments or management are discussed.


Introduction
Global to regional mapping and monitoring of drought is increasingly considered for food security assessment in agricultural lands [1] and for numerous ecosystem processes, such as carbon fluxes, plant productivity [2,3] or fire risk [4].Various indexes for monitoring drought have been developed on the basis of meteorological datasets (see [5] for a review).Meanwhile, drought detection from remote sensing has been increasingly used since the 1990s in natural and cultivated ecosystems, particularly with the increasing temporal frequency of satellite images assessing land cover status (SPOT-VEGETATION, AVHRR and MODIS).Thus, temporal tracking of water content with spectral measures has been widely investigated, especially in water-limited regions (e.g., [6][7][8][9]).In ecosystems belonging to transition zones, as is the case for the North African forests representing transition zones between sub-humid Mediterranean forests and the semi-arid region at the northern edge of the Sahara, the accurate assessment of drought by the use of remote sensing and field observations is essential to determine the major constraints of such ecosystems.Vegetation indexes, as the result of empirical combinations of reflectance in various wavelengths, had contributed in detecting changes in biomass status as changes in chlorophyll activity (biomass vegetation indexes, BVI) and vegetation water content (moisture vegetation indexes, MVI).Various studies have investigated the combination of BVI and MVI, especially to assess drought conditions as interannual precipitation deficit/surplus [10] or, more specifically, its impact on vegetation water stress (e.g., [11,12]) or fire risk [4].By combining the near-infrared (NIR) and shortwave infrared (SWIR) reflectances, variations induced by leaf internal structure and leaf dry matter content can be disentangled to improve the accuracy in retrieving the vegetation water content [13].In turn, simulated and remotely-sensed reflectances in the SWIR domain can capture leaf water content at the canopy level [11,14,15].
For biomass monitoring, various formulations of vegetation indexes were suggested, essentially to discriminate vegetation densities (combining reflectances in red and NIR bands), atmospheric effects (adding green and blue bands) and soil background adjustments.MVI computation is based on the combination of NIR and SWIR wavelengths (1200 to 2100 nm), where the leaf water content plays a decisive role on the leaf reflectance [14].We have summarized in Table 1 [15][16][17][18][19][20][25][26][27][28] the most used vegetation indexes, including their formulation and authors.BVI include, among others, NDVI (Normalized Difference Vegetation Index) and EVI (Enhanced Vegetation Index) for reducing atmospheric influence on the signal by the use of the blue spectral band, which does not saturate at high biomass compared to NDVI [16].The Soil Adjusted Vegetation Index (SAVI) was developed to consider the reflectance of the complex soil-vegetation interaction via an adjustment factor (L = 0.5), shown to reduce the soil signal to a non-significant noise for many vegetation density ranges [17].Qi [18] proposed a modified SAVI, named Modified Soil Adjusted Vegetation Index (MSAVI), where the L factor is adjusted automatically to vegetation density (through an iterative function) to decrease the soil background effects on the vegetation signal.The Adjusted Normalized Difference Vegetation Index (ANDVI) combines both the adjustment of soil and atmospheric effects by the blue and green bands [19].MVI indexes use multiple wavelengths in the middle infrared interval (B5 for the Normalized Difference Water Index (NDWI), B6 and B7 for the Normalized Difference Infrared Indexes (NDII6 and NDII7)) and for the Global Moisture Vegetation Indexes (GVMI6 and GVMI7), designed for reducing atmospheric effects [11,20,21].Vegetation water status may be related to various indicators, such as "equivalent water thickness" EWT [13,22,23], allowing the leaf amount of water computation or the fuel moisture content (FMC), the mass of water contained within vegetation in relation to the dry mass [4].Estimating EWT or FMC from remote sensing may be empirical (statistic-based) or physical, as indicated in [4].Empirical methods compare field measurements of vegetation water content to MVI and BVI computed from multi-temporal satellite images, whereas physical models, known as radiative transfer models (RTM), use rather simulated reflectances or hyper spectral data.Optical remote sensing was also combined in other studies with passive microwave observations known to better succeed in canopy water content, integrating the contribution of trunk and branch water content (e.g., [24]).Table 1.Biomass and moisture remote sensing vegetation indexes (bands correspond to the spectral interval of MODIS of Table 2).SAVI, Soil Adjusted Vegetation Index; ANDVI, Adjusted Normalized Difference Vegetation Index; NDII, Normalized Difference Infrared Index; GVMI, Normalized Difference Infrared Index; MSAVI, Modified Soil Adjusted Vegetation Index.

Biomass Vegetation Indexes
References NDVI = (B2 − B1)/(B2 + B1) [25,26] Table 1.Cont.The performance of these indexes (BVI and MVI or combined ones) has been largely investigated to link vegetation water content to remotely-sensed data in various climates and ecosystems.In Mediterranean-type ecosystems, vegetation water content was shown to be linearly related to a combination of BVI and MVI using either high or coarse spatial resolution images [29][30][31].In semi-arid ecosystems, Fensholt and Sandholt [14] found a strong relationship between the MVI and soil moisture (Sahelian zone in Africa); Ceccato et al. [32] and, more recently, Sow et al. [33] reported consistent relationships between field measurements of EWT and various BVI and MVI, allowing for a regional assessment of the seasonal dynamic of ecosystem dryness.Similar results were reported by [34] in the semi-arid region of Arizona.

Moisture Vegetation Indexes
However, the performance of VI as drought indicators was shown to be site-specific due to spatial differences in leaf and canopy characteristics, soil background, sensor characteristics and observation conditions [4].Methods were then proposed to overcome this fact by modifying the linear regression between VI and EWT variables, accounting for the seasonal and interannual vegetation variability [4,30,35], as well as the species mixture and their related biophysical properties [36].Considering that leaf area index (LAI) is an integrated proxy for canopy density and ecosystem water content capacity according to the ecohydrological equilibrium theory [37], canopy variables represented by LAI are also important parameters affecting the retrieval of water vegetation status from remote sensing vegetation indexes.Based on the Prospect simulation model [38] combined with canopy models, Ceccato et al. [20] showed that GVMI is suitable for retrieving vegetation water content when the LAI is equal to or greater than two due to the effect of soil background.Similar conclusions were reported in [33], who showed that the correlations between measured water vegetation content (EWT or FMC) and VI decreased moderately with decreasing aridity represented by increasing tree cover in the Sahel.Over areas with low LAI, there is a general consensus that red/NIR vegetation indexes saturate and, in turn, are less efficient than in foliar-water indexes for canopy water content estimation [21,39,40].
From these identified contrasted performances of VI for monitoring vegetation moisture across ecosystems with varying LAI, main objectives of our study are: (1) making the assessment of the performance of BVI and MVI in a Mediterranean natural ecosystem at the transition zone with semi-arid climate conditions in Northern Africa, characterized by a heavy rainfall gradient, where vegetation is characterized by a low LAI and a tree density gradient associated with a significant shrub or grassland understory layer; (2) comparing EWT computed from seasonal field measurements collected in the dry season of the year 2010 across the climate gradient with VI computed from the MODIS sensor and testing the hypothesis that EWT is differentially correlated with VI, according to LAI; and (3) investigating the ability of the inversed model to reproduce comprehensive regional scale EWT and seasonal drought across the region.

Study Area
The study area covers the Kroumirie forest territory (2500 km 2 ), belonging to the narrow Mediterranean forest band in North Africa and located at the extreme north of Tunisia (7°52′E-37°29′N to 9°33′E-36°22′N) between the Mediterranean Sea and the Sahara Desert (Figure 1).The region is mountainous with elevations ranging from 0 m on the northern coast to 1235 m.Regional climate is Mediterranean, where precipitation is mainly distributed in autumn, winter and spring, while summers are dry, with a strong north/south gradient.The northern bound of the region is influenced by its proximity to the Mediterranean Sea with wet (800 mm annually) and mild temperatures, the upper mountains are very wet (1250 mm annual rainfall) with low winter temperatures and mild summers, and the southern floodplain is dry (360 mm annually) with warm temperatures.Vegetation is mostly an evergreen Quercus suber forest, with deciduous Quercus canariensis on the upper coldest sites and evergreen needle leaf Pinus pinaster on the coast.Understory or burnt areas in the region are covered by a maquis shrubland composed of Cistus monspeliensis, Erica arborea, Arbutus unedo, Phillyrea latifolia, Pistacia lentiscus and Calicotome villosa, among others.

Field Measurements for Fuel Moisture Content (FMC)
We conducted a campaign of field measurements for FMC in 2010 during the dry summer period from 1 June to 28 September, according to the "reseau hydrique" method for fire risk and Mediterranean forest conservation in southern France [41].
Seven permanent plots were selected across the regional climate gradient (Table 3).For each plot, we selected four biological indicators of vegetation dryness: the surface leaf litter, the shallow-rooted evergreen shrubs (Cistus monspeliensis, Myrtus communis), the deep-rooted shrubs (Erica arborea, Phillyrea arborea, Pistacia lentiscus) and the deep-rooted evergreen cork oak tree (Quercus suber).We selected three individuals per plot and three apical twigs per individual (15/20 g fresh weight).Fresh vegetation samples were collected every week, around 12:00 P.M. at the hottest point during the day, and were inserted in a stainless metal box, hermetically sealed with a sticky and elastic band.The boxes were then carried to the laboratory located in Tabarka (Institut Sylvo Pastoral).Unsealed, but non-open boxes (box and vegetation humid BVH) were weighted with a precision balance.Once weighted, open boxes were inserted in an oven at 60 °C for 24 h and weighted again (box and vegetation dry BVD).Vegetation was extracted from the box, and the empty dry box was weighted on its own (BD).Fuel moisture content (% dry weight) was estimated according to Equation (1) (Figure 2).

Equivalent Water Thickness Calculation
A formulation of the equivalent water thickness (EWT) was suggested by [28] (Equation ( 2)) at the leaf level, to quantify the equivalent depth of water (in mm) contained in the plants, as relevant information for remote sensing studies using NIR and SWIR bands [4].At the canopy level, EWT is commonly computed by multiplying FMC by the amount of leaves per unit of ground area (LAI) and is defined as a hypothetical thickness of a single layer of water multiplied by the number of layers determined by the LAI [20].(2) EWT: equivalent water thickness of leaves (mm); FM: fresh matter (kg); DM: dry matter (kg); LA: leaf area (m 2 );ρw: density of pure water (1000 kg•m −3 ).Equation ( 2) is valid when the ecosystem is homogeneous in terms of vegetation functioning, the desiccation rate of plants and the seasonality of LAI through similar phenology.Our ecosystem is composed of evergreen trees with tree cover ranging from 25% to 100%, with an understory of evergreen shrubs and a ground layer covered by annual herbaceous species drying and dying off at the end of spring.We are confronted here with the peculiar aspect that: (1) the surface litter, the shrub layer and the tree layer have contrasted water budget functioning; and (2) the distribution of cover within each vegetation type varies spatially across the region.We reformulated Equation (2) into Equation (3) to account for EWT at canopy level EWTCAN in our mixed ecosystem composed of N species, admitting each a leaf mass area LMAi.

  
11 ρ As we have a unique LAI value for the whole MODIS pixel, which mixes tree/shrub LAI and herbaceous LAI, we decomposed LAI-MODIS based on previous observations that LAI-MODIS follows a seasonal pattern not matching the actual LAI variations measured on site for evergreen Mediterranean trees [42].We then stated that minimum LAI-MODIS is close to tree/shrub LAI when the herbaceous layer is dry, so that the seasonal variation between remotely-sensed LAI and LAImin characterizes the herbaceous LAI variations.We derived Equation (7) from Equation ( 6) by partitioning the herbaceous and the tree components of the pixel according to tree cover COV.Tree/shrub leaf area was estimated as [COV × LAImin] all along the year, and the litter LAI variations were estimated as the [(1 − COV) × (LAI − LAImin)].

MODIS Images Processing
We used MODIS time series of preprocessed LAI (freely available at [43]) from 1-km resolution of 16-day composite products, MOD15A2 (TERRA sensor), MYD15A2 (AQUA sensor), MCD15A2 (combined TERRA/AQUA), and the bidirectional reflectance (MOD09A1) 500-m resolution of 8-day composite product (Table 2).We compared the three MODIS LAI products to measured plots of LAI over the study region of Kroumirie in 2006 (refer to [44] for a detailed explanation of LAI measurements).The least RMSE was obtained with TERRA (RMSE = 2.36), compared to AQUA (RMSE = 3.24) and to TERRA/AQUA (RMSE = 3.49).However, the RMSE value remained high even for the TERRA product, especially for highly-covered plots (NDVI > 0.6).Thus, we suggested a calibration method for the TERRA-MODIS LAI product based on both ground measurements and SPOT images that were used as ground truth LAI maps (refer to [44] for the calibration process details).
From visualization and profile examination, we noticed important noise due to both cloud contamination (in the wet season) and sensor defaults that had to be removed from all reflectance bands of the MOD09A1 and the derived indexes; we also noticed a systematic noise (strips) in B5 (ρ1230-1250 nm).Therefore, it had been necessary to process time series data to suppress unusually high or low VI values.We used the Timesat software smoothing tool [45,46] to overcome MODIS product imperfections and contaminated images.This tool comprises two main processing stages.First, a preliminary processing is made to remove spikes and outliers due to clouds or imperfections of the sensor; then a filtering method is used.Two types of filtering methods are proposed in Timesat software: (1) Local model functions (the asymmetric Gaussians function or the double logistic function) that are being fit to data (VI at a given time "t") in intervals around maxima and minima in the time-series; (2) adaptive Savitzky-Golay method where each VI(t) value is replaced by a linear combination of nearby values for (2n + 1) points, with "n" being the number of time steps preceding or following the time "t" [47].The filter is adaptive since the window size "n" can be adjusted during the processing if there is a large variation in an interval around a given VI value; in this case, the filtering is redone with a smaller "n" size window.This latter filter has been used with success in minimizing noise in NDVI [48,49].Our time series data of various VI had been processed within the Timesat software tool; we also find a better performance of the adaptive Savitzky-Golay method over the local model functions.Thus, the smoothing significantly corrected the abrupt changes in VI, except for NDWI computed from B2 and B5 bands; this latter band represents a systematic strip noise over the smoothed images, and thus, NDWI was further excluded from the analysis.
Weekly reflectances MODIS images corresponding to a 2010 field campaign were processed as follows: we localized field plots on the referenced MODIS images and chose a relatively homogeneous area surrounding the plot (using Google Earth) to make comparisons between remotely-sensed data and EWTCAN based on field measurements.Then, within each region of interest (ROI), biomass and moisture indexes (BVI and MVI) were calculated in terms of basic statistics of the existing pixels (min, max, average, standard deviation).

Statistical Analysis
Relationships between EWTCAN and VI were evaluated by linear regressions with the "Lm" function in "R cran" software, from which we derived the intercept, slope, their corresponding p-values and R 2 coefficients.Regressions were performed for EWTCAN data assembled at the regional level (All Plots) and for each of the 7 study sites (Plot).We then evaluated the linear regression between plot-level regression slopes and LAI of these plots with the same statistical tool, to derive the intercept, slope, their p-values and R 2 for each remotely-sensed variable.
Statistical models at the All Plot and Plot level were then inversed to estimate EWTCAN at the plot level and at the landscape level using 8-day MODIS image time series VI.Root mean squared errors (RMSE) between field-derived and model-derived EWTCAN were determined for each model (All Plot and Plot with LAI-based relationships).

Spatial and Temporal Relationship between EWTCAN and Remotely Sensed vegetation indexes (VI)
The relationships between estimated EWTCAN from ground measurements and VI derived from MODIS spectral bands are shown in Figure 3, with their corresponding regression analysis statistics (Table 4).In general, all indexes performed satisfactorily well, except for EVI and SAVI, which provided the lowest correlations.All tested MVI show strong relationships (R 2 > 0.60), whereas NDVI, MSAVI and ANDVI show the strongest relation amongst the BVI (Table 4a).As a regional pattern (analysis of All Plots, Table 4a), we observe that, across the LAI gradient of the selected study plots with LAI varying between 1.10 and 2.90 (range of maximum values in Table 3), EWTCAN vary between 0 and 0.05 g•cm −2 , and all VI significantly follow a positive linear relationship with EWTCAN (R 2 varying from 0.43 for SAVI to 0.73 for NDII7, MSAVI and ANDVI), a result consistent with previously reported analysis for natural and cultivated vegetation [11,23,50].Based on these observations, we then tested if the plot-level seasonal course of EWTCAN follows the same linear relationship as the generic spatial regional relationship (Table 4b, regressions by plot), so that this latter one can be used to estimate the seasonal variation of EWTCAN.We found most values of R 2 to be much lower at the plot level, with the lowest value observed for the "Ain Draham" plot characterized by the highest LAI (Table 3).We noticed from Figure 3 that "Beni Mtir" is away from the general trend of the relationship for SAVI, which includes all of the sites.Actually, the SAVI general correlation (Table 4a) is the lowest, with R 2 = 0.43 quantifying this anomaly.This could be explained by the low LAI of "Beni Mtir" (Table 3) and the fact that the SAVI is a precisely sensitive index to soil effect.We also noticed that "Hammam Bourguiba" and "Khroufa" sites present somehow parallel scatterplots with the various indexes (Figure 3), but they admit different R 2 values (Table 4b), meaning that the dispersion of points around the trend line is different for these two sites while having the same slope.Globally, as we found most values of R 2 to be much lower at the plot level, our hypothesis that the regional regression can be applied to estimate the temporal course of EWTCAN across the season is then rejected.We then hypothesized that the slopes of the regression between EWTCAN and VI were correlated to LAI.Table 4c represents the slopes, intercepts, p-values and R 2 of the (LAI, Slope(A)) relationship for each study plot, with Slope (A) being the slope of the (EWTCAN, VI) linear relation given in Table 4a.Figure 4 represents the relationship between LAI study plots and the slope of the plot-level (EWTCAN, VI) relationship.First, we found that for LAI ≤ 2, the decreasing relation is significative at the 0.05 threshold (black dots in Figure 4) and becomes not significative for LAI greater than two.This suggests that for sparse vegetation (LAI ≤ 2), EWTCAN is driven by the LAI effect.For all BVI, we observed a significant decreasing linear relationship at the 0.05 threshold, indicating that the greater the LAI, the less responsive are VI to changes in EWTCAN.However, for MVI, the (LAI, Slope (A)) relationship is only significant for NDII6.This was also reported in [20], who mentioned that the GVMI6 (from simulated reflectances) is suitable for retrieving vegetation water content when the LAI is equal to or greater than two, which is not the case in all of the observed sites in the study region.In light of our experimental results, we can conclude that BVI are correlated with EWTCAN, and these relations are also driven by the LAI pattern for sparse vegetation (LAI ≤ 2) that ensures the ecosystem's water carrying capacity (or maximum EWTCAN); only NDII6 represents this pattern amongst the tested MVI.Our peculiar pattern of LAI-dependent (EWTCAN, VI) relationship points out the heterogeneous plant functioning in Mediterranean-type ecosystems.We measured contrasted species' desiccation rates, including litter, in accordance with previous studies on Mediterranean vegetation [51][52][53].These contrasted species-specific fuel moistures could partly explain the non-linearities of the (EWTCAN, VI) relationship as a response of LAI-dependent differential functional diversity.We might also wonder if the seasonal variations of LMA (see [54] for a review), as an increase in the canopy as a function of leaf aging [55] or the LMA decrease in the litter layer as a function of decomposition [56], could have substantially affected EWTCAN, so we calculated with a constant LMA, as already pointed out in [4,20], but which is hardly ever considered in regional-scale EWTCAN assessment from remote sensing.Indeed, an increase in LMA for aged leaves would enhance EWTCAN at the end of the dry season and, in turn, increase the (EWTCAN, VI) slope for high LAI sites.On the contrary, in lower LAI sites more dominated by grasslands, a decrease in litter LMA along the season could reduce EWTCAN and, in turn, decrease the (EWTCAN, VI) slope.These adjustments could rebalance the site-specific (EWTCAN, VI) relationships to the generic one to be used both spatially and temporally.Our approach relying on a plot-level (LAI, Slope(A)) relationship could then readjust EWTCAN, so we estimated with fixed LMA, and in turn, this would prevent time-consuming LMA measurements for each measured FMC.This hypothesis should be further evaluated with field campaigns.

Model Inversion for Regional EWTCAN
Based on the highest correlations (R 2 > 0.60), model retrieving of EWTCAN using VI is proposed.We called the MOD1 and MOD2 models, respectively, non-integrating and integrating LAI effects in the regression between EWTCAN and VI.Equation ( 8) makes use of the simple linear relation (with R1 2 ) between VI and EWTCAN, where "A" and "B" are, respectively, the slope and intercept (Table 4a).Equation ( 9) uses the slope varying with vegetation density expressed by a linear relation (with R2 2 ) with a decreasing slope when LAI increases (intercept and slope "α" and "β" from Table 4c and intercept "B" from Table 4a).As stated earlier (Section 3.1), the LAI pattern influences the retrieval of EWTCAN just in the range of LAI ≤ 2 (empty dots in Figure 4 for LAI > 2).Thus, MOD2 is appropriate for sparse vegetation, whereas MOD1 is proposed for LAI > 2.
Table 5 summarizes inversion models for EWTCAN estimation from BVI and MVI indexes and the conditions of use of these models to obtain positive values of EWTCAN (last column).The comparison of the two models' performance was carried out by the total RMSE computation showing the dispersion of both models in front of EWTCAN_FIELD and relative to the 1:1 line (Figure 5).Inversion models based on VI/LAI (MOD2, Equation ( 9)) having the strongest R1 2 and R2 2 at the same time produce slightly lower RMSE (underlined R1 2 , R2 2 and RMSE in Table 5) compared to MOD1.This could be explained by the fact that most of the plots admit LAI ≤ 2, and consequently, their EWTCAN is more appropriately simulated by MOD2 rather than MOD1.
The effect of LAI integration in the inversion model is highlighted in Figure 6 (comparison of inversion MOD1 and MOD2 for VI with both R1 2 and R2 2 > 0.60), showing that for a given value of VI, corresponding EWTCAN is decreasing as vegetation LAI is decreasing (for LAI ≤ 2).Consequently, at this range of LAI, EWTCAN is overestimated if it is determined from the VI-based inversion model (MOD1).Our results suggest that for a canopy characterized by reduced tree cover, LAI highly affects the EWTCAN value and had to be considered for a more realistic estimation in the Mediterranean natural vegetation areas to avoid overestimated EWTCAN.

Regional Analysis of EWTCAN
From the previously-established models between remotely-sensed VI and estimated EWTCAN from ground measurements, we built a multi-temporal set of weekly EWTCAN for the year 2010 based on inversion models using indexes for which both R1 2 and R2 2 are high, these indexes are NDVI, MSAVI, ANDVI and NDII6.We explored the effect of LAI integration (MOD2 compared to MOD1) at the regional scale, where EWTCAN were averaged according to two LAI classes (LAI ≤ 1 and LAI between one and two).Figure 7 shows the varying profiles of EWTCAN along the dry season of the year 2010.For sparse vegetation (LAI ≤ 1), EWTCAN is almost null because of the high desiccation of the herbaceous layer (Figure 7a) (cf.litter FMC in Figure 2); for LAI between one and two, all inversion models present coherent results with agreement of overestimated values of EWTCAN simulated by MOD1, which does not consider the LAI effect.The EWTCAN time-course from LAI-adapted inversion models (based on significantly related VI to EWTCAN ) is illustrated in Figure 7b (for both LAI between one and two and LAI > 2), showing close profiles resulting from the inversion models based on the various VI (NDVI, ANDVI, MSAVI and NDII6) applied during the dry season (June to September 2010).In summary, the inversion models for retrieving EWTCAN from VI have to consider the LAI effect when applied at the regional scale for LAI ≤ 2. On the other hand, for LAI > 2, inversion models based on only VI are rather recommended.We finally used EWTCAN from MSAVI (presenting the highest R1 2 and R2 2 and the least RMSE at the plot level) to extrapolate the weekly regional mapping of EWTCAN by applying MOD1 for LAI > 2 and MOD2 for LAI ≤ 2, as stated earlier.We presented the temporal difference in regional EWTCAN for spring, midsummer, end of summer and early fall, four keystone periods for Mediterranean vegetation functioning, representing, respectively, the start and the middle of the dry period, the early end of the dry period, when post-summer storms sporadically happen, and the end of the dry period, when large scale rainfall events refill soil water content up to field capacity.We observed contrasted regional patterns according to the season (Figure 8a-d) and to the topo-climate patterns represented by the digital elevation model (ASTER-GDEM) (Figure 8e).
In spring (Figure 8a), drier conditions appear in the valley bottom, while higher altitudes keep getting wetter.Regression analysis conducted between the spring deficit maps confirms the altitude-dependent effect on the desiccation process with a regression coefficient between spring deficit and altitude maps of 0.70.This contrasted pattern might be related to the differential evapotranspiration (ET0) rates across topographical features [57,58], so that, in the valley bottom, with higher temperatures, ET0 might be higher than rainfall amounts overall, leading to increasing drought, while lower ET0 at higher altitudes do not compensate for water inputs from rainfall.For this period, the regional pattern of drought is then driven by differential ET0 according to topographical effects, rather than the spatial pattern of rainfall.In the middle of summer (Figure 8b), the entire region, whatever its topographical situation, is submitted to increasing drought in the absence of any rainfall event.At this time, when drought intensity is at its maximum and soil water content is low, the desiccation pattern is homogenized, as pointed out for Pinus halepensis forest under the Mediterranean climate [59,60].
At the end of summer (end of August) (Figure 8c), sporadic rainfalls start to happen, mostly dominated by a west/north incoming wet air mass from the Mediterranean Sea warmed up by the summer season high temperatures.The altitude-dependent effect on rewetting is also present, with a regression coefficient between end of summer deficit and altitude maps of 0.67.The Atlas Mountains act as a barrier to cloud movements, leading to rainfall pouring only on the northern part of the region, while the most southern hills remain dry [61].We indeed observe a rewetting of the northern part of the region (green) from remote sensing for this period, while the most southern part of the region keeps drying off.During this period of the year, low intensity rainfalls (<10 mm), more sporadic than large-scale heavy events, can lead to short-term pulses with a significant impact on ecosystem water budgets [62].We observed here a contrasted pattern with a significant boundary between the rainfall-affected and non-affected parts of the region.
Finally, early fall changes in EWTCAN (Figure 8d) illustrate a homogeneous rewetting over the whole region, except high altitudes potentially already filled up to field capacity at that time and not submitted to rewetting.
This seasonal pattern of desiccation/rewetting with a high spatial resolution brings useful information for drought monitoring and forest management trough impact assessments on local-scale forest functioning [63][64][65], tree mortality [59,66] or fire risk [67].More particularly, this information brings fine-scale information on the spatial pattern of drought, hardly precisely assessed from rough interpolations of meteorological variables [68], particularly when rain gauges are sparse, as in Tunisia [69] and more generally in Africa [70], and still uncertain from global-scale remote sensing [71][72][73].

Conclusions
This study investigated the suitability of biomass and moisture vegetation indexes for vegetation water content expressed by the equivalent water thickness for Mediterranean vegetation across a tree cover LAI gradient with contrasted water budgets and desiccation rates.We proposed peculiar adjustments for low LAI and tree cover ecosystems.First, the comparison of remote sensing (MODIS) to vegetation water content required a formulation of EWT at the canopy level that had to consider the mixture of plant functional types for water use (trees, shrubs and herbaceous layers) and a varying vegetation cover.The performance of moisture indexes (computed from NIR and SWIR) was comparable (R 2 from 0.65 to 0.76), whereas for biomass vegetation indexes, MSAVI and ANDVI performed better than NDVI, EVI and SAVI (R 2 from 0.43 to 0.73).We particularly demonstrated and quantified the effect of LAI on retrieving EWTCAN mentioned in [16,24,27] for regional-scale drought assessment under a tree cover/LAI gradient from dense to open forest.LAI adaptive inversion models for retrieving EWTCAN are then proposed for Mediterranean forests drought monitoring by the VI/LAI inversion model for moderately-covered areas and the VI inversion model for denser forests.The regional scale investigation of EWTCAN in 2010 derived from the MSAVI/LAI inversion model (for which we found the best linear fits) revealed the temporal changes in regional patterns of ecosystem desiccation and rewetting along the summer season, including late summer rainfall pulses with scattered patterns.This integrated approach offered a fine-scale drought assessment for forest management when rainfall gauges are sparse and fine-scale topo-climates are hardly quantified.

Figure 3 .
Figure 3. Regional (All Plots) relationship between equivalent water thickness at the canopy level (EWTCAN; g•cm −2 ) and MODIS-VI indexes with the resulting fitted linear model (red line).Colored dots represent each of the seven study plots.

Figure 4 .
Figure 4. Relationship between LAI plot (X-axis) level and Slope(A) of (EWTCAN, VI) (Y-axis).(Red line, significant relationship p-value < 0.05; grey line, no significant relationship).Black dots represent significant plot-level relationships, and empty dots represent non-significant relationships at the 0.05 threshold for p-values.

Figure 8 .
Figure 8.Time course of regional maps of EWTCAN change (g• cm −2 ) in spring ((a); between DOY 137 and 169), mid-summer ((b); between DOY 209 and 217), late summer ((c); between DOY 241 and 257) and early fall ((d); between DOY 281 and 297).Reddish colors indicate drier conditions along the period, and greenish colors indicate wetter conditions along the period.(e) Map of the altitude of the Kroumirie region (processed from ASTER-GDEM).

Table 2 .
MODIS products and corresponding spectral bands.

Table 4 .
(a) Regression of linear relation (EWTCAN, VI) by plot; (b) regression of linear relation (EWTCAN, VI) for all plots; (c) regression between LAI and Slope (A) for all plots.p_inter and p_slope are p-values of intercept and slope.

Table 5 .
Inversion models for EWTCAN (g•cm −2 ) retrieval from biomass vegetation indexes (BVI), moisture vegetation indexes (MVI) and LAI.R 2 > 0.60 and the corresponding RMSE are underlined.Conditions imposed on VI are to avoid negative EWT values.