Manuscript: bg-2015-62 Title: Seasonal variation in vegetation water content estimated from proximal sensing and MODIS time series in a Mediterranean Fluxnet site

22 This study evaluates three different metrics of vegetation water content estimated in of an 23 herbaceous cover in a Mediterranean wooded grassland (dehesa) ecosystem. from proximal 24 sensing and MODIS satellite imagery: Fuel Moisture Content (FMC), Equivalent Water 25 Thickness (EWT) and Canopy Water Content (CWC) were estimated from proximal sensing 26 and MODIS satellite imagery. Dry matter (Dm) and Leaf area Index (LAI) connect the three 27 metrics and were also analyzed. in order to connect FMC with EWT and EWT with CWC, 28


Introduction
Water in leaves is a limiting factor for different physiological processes of vegetation and its deficit causes malfunctioning of different cellular processes.Water is involved in the thermal regulation of plant trough transpiration and also becomes crucial in the uptake of CO 2 for photosynthesis (Chaves et al., 2003).It is also fundamental to maintain turgor pressure, which controls different functional processes of plants like cell enlargement or gas exchange (Taiz and Zeiger, 2010).
Different metrics quantify vegetation water content.Fuel Moisture Content (FMC) (Desbois et al., 1997), defined as the mass of water per unit mass of vegetation, has been extensively used to estimate the fire risk occurrence and fire propagation (García et al., 2008;Yebra et al., 2008b).Equivalent Water Thickness (EWT) or Leaf Water Content (LWC), defined as the mass of water per leaf area, measures the thickness of the water layer with the same leaf area (Danson et al., 1992).Several studies showed that EWT can be retrieved from spectral information at leaf level as it is directly related to the water absorption depth of leaves (Ceccato et al., 2001;Datt, 1999).FMC and EWT are related each other since EWT can be expressed as FMC multiplied by the dry matter (Dm) and divided by 100.Dm is defined as the ratio of leaf dry weight and leaf area (Bowyer and Danson, 2004;Chuvieco et al., 2003).Finally, another metric is the Canopy Water Content (CWC), the mass of water in the canopy per ground area (Cheng et al., 2008;Trombetti et al., 2008).Therefore, CWC represent the product of EWT and Leaf Area Index (LAI), offering not just information at leaf level but at canopy level.
Field sampling of FMC, EWT or CWC relies usually on gravitational methods but this method is quite limited for estimates at regional to global spatial scales, since it requires interpolation to bridge the gaps in both time and space.Remote sensing is a powerful alternative data source to provide information on vegetation water content as it fills such temporal and spatial gaps.Monitoring vegetation water content with remote sensing benefits agriculture, to control crop production and prevent stress in plants (Peñuelas Introduction

Conclusions References
Tables Figures
To estimate plant water content with remote sensing, vegetation spectral reflectance has been primarily related to specific water absorption bands in the Short Wave Infrared region (SWIR, 1300-2500 nm) (Ceccato et al., 2001;Zarco-Tejada et al., 2003).Other studies related vegetation water content to spectral indices that do not include SWIR data.These indices monitor the vegetation water content by indirectly relating it to another biophysical parameter that is used as a proxy of water stress.This is the case of the Normalized Difference Vegetation Index (NDVI) (Tucker, 1979), with bands in the Visible (VIS) and Near Infrared (NIR) spectral region, that has shown a close relationship between vegetation biomass, chlorophyll and water content in grasslands (Chuvieco et al., 2003(Chuvieco et al., , 2004;;García et al., 2008;Yebra et al., 2008b).Least squares regression models have served to empirically relate observed measurements of water content to spectral indices.These models have their weakest point of being site dependent, requiring long datasets for calibration (Chuvieco et al., 2009) and showing different results when the models are extrapolated to other sites using different data sets, making difficult their applicability (Riaño et al., 2005;Yebra et al., 2008a).
Radiative Transfer Models (RTM) simulate vegetation spectra and are a sound alternative to empirical modeling.They can be applied to different locations to estimate different vegetation parameters, as long as the RTM is a true representation of the vegetation canopy.For example, Trombetti et al. (2008) predicted CWC for the continental US using RTM PROSAILH (Baret et al., 1992) simulations.Their model was calibrated with CWC from AVIRIS hyperspectral water absorption bands.Another example is Jurdao et al. (2013) who inverted the RTM GEOSAIL (Huemmrich, 2001) to estimate FMC that was validated with extensive field sampling in Spain.
This study compares the model performance using empirical relationships between three different metrics of vegetation water content (FMC, EWT, CWC) and nine spectral indices calculated from both proximal sensing and MODIS data.Dm and LAI were Introduction

Conclusions References
Tables Figures

Back Close
Full also analyzed in order to connect FMC with EWT and EWT with CWC, respectively.
Firstly, an analysis of the temporal and spatial variability of the different measurements collected in the field was conducted to evaluate which biophysical parameter offers more information.Secondly, the model performance was evaluated to select the most accurate models.Finally, these empirical models were compared to two RTM based models, one to estimate FMC (Jurdao et al., 2013) and the other for CWC (Trombetti et al., 2008).

Site description
The study site is located at Las Majadas del Tiétar (Spain) FLUXNET site (http: //fluxnet.ornl.gov/site/440,last accessed 05 June 2014) (Fig. 1).The area is a dehesa, a Mediterranean wooded grassland, which is an ecosystem that occupies about 4 % (2.5 Mha) of the Iberian Peninsula (Castro, 1997).Common tree species are different varieties of oaks, here mostly Quercus ilex ssp.ballota (L.), whose acorns and leaves are mainly used as forage for pigs and cows, respectively.The scattered oak trees have a 9 m mean height and 6 m mean crown diameter.Due to its deep and wide root system, this species is resistant to long drought periods (Camarero et al., 2012).Short grassland covers 86 % of the area that is managed for cow shepherding. of ecosystem in Mediterranean areas worldwide, the need to track the responses to water stress conditions, together with the presence of a FLUXNET eddy covariance flux tower justifies the selection of this site.

Vegetation sampling
Grasslands water status was destructively sampled every two weeks from April 2009 to April 2011 assuring no rain occurred in the two previous days to avoid sampling superficial water on the leaves.During the summer, grasslands become completely dry and samples were not collected in 2009.However, to ensure the time series continuity of at least one phenological cycle, sampling was restated throughout the summer of 2010.This strategy leaded to a total of 21 valid sampling days for the whole study period.
Six 25 m × 25 m plots were randomly located within the 500 m MODIS pixel that contained the eddy covariance flux tower (Fig. 1).Three grassland samples were collected from three 25 cm×25 cm quadrants randomly positioned within each plot.Rooted grasses were collected inside the quadrant using clippers.Additionally, a smaller sample was collected outside of the quadrant but nearby to estimate EWT -EWT Sample hereafter -in a cost-effective way, containing a representative proportion of surrounding species.All samples were placed in sealed plastic bags, weighed on a scale with 0.01 g precision and then transported in a cooler to the laboratory.Each EWT Sample and a sub-sample from each quadrant were scanned at 150 pixels per inch (ppi) in an Epson Perfection V30 color scanner (Epson American Inc., Long Beach, CA, USA).Leaf area was calculated automatically from the scanned images using the unsupervised classification algorithm ISOCLUS with 16 iterations in PCI Geomatica (PCI Geomatics, Richmond Hill, Ontario, Canada).All samples were then placed in an oven for 48 h at a constant temperature of 60 • C to obtain their dry weigh.Five biophysical variables were obtained from the collected vegetation samples: FMC, EWT, Dm, CWC and LAI.FMC was determined from the fresh and dry weights of both the whole quadrant sample (FMC Q ) and the EWT Sample (FMC E ) according to Eq. ( 1) where W Fresh is the fresh weight of the sample measured in the field and W Dry is the oven dried weight.
The EWT Sample permitted to calculate both, EWT and Dm, since fresh/dry weight and leaf area were all measured (Eq.2): where Area Leaf is the leaf area of the EWT Sample .
CWC was calculated from two different approaches.In the first one, quadrant and EWT Sample were combined using the following expression (Eq.4): where LAI is the leaf area index of the grass within the quadrant and EWT is obtained from Eq. ( 2).
The grass height was very short due to cow shepherding, so LAI, rather than optically estimated, was measured with gravitational methods (He et al., 2007).The biomass to leaf area ratio of a sub-sample inside the quadrant to the total quadrant's biomass provided LAI using the following expression (Eq.where W Dry is the total dry weight of the whole sample inside the quadrant, W Sub Dry is the dry weight of a sub-sample of W Dry , Area Sub Leaf is the sub-sample leaf area and Area is the total area of the quadrant. The second approach measured CWC from the fresh and dry weight difference inside the quadrant (CWC Q in Eq. 6):

Proximal sensing
Simultaneously to vegetation sampling, proximal sensing data were acquired using an ASD FieldSpec ® 3 spectroradiometer (http://www.asdi.com/)along NE-SW and NW-SE transects in each 25 m × 25 m plot.This instrument measures reflectance from 350 to 2500 nm.Before measuring along each transect, dark current was recorded, instrument settings were optimized and reference spectra were acquired using a Spectralon ® 99 % reflective reference panel (Labsphere Inc., North Sutton, NH, USA).All measurements were taken under clear sky within about ±2 h from local solar noon, to guarantee homogeneous illumination and maximum solar irradiance.Sky conditions were recorded in the field logs, and a quality control check removed the spectra where illumination changes may have occurred after calibration.The ASD was handled without fiber optics to reduce the directional effects on the spectroradiometer's fiber bundle field of view (FOV) (MacArthur et al., 2012).Spectra were acquired at approximately 1.2 m height, rendering a sensor footprint diameter of about 53 cm, since nominal FOV is 25 • .
An average of approximately 10 spectral measurements was calculated for each transect and then spectrally resampled to MODIS bands using ITT ENVI 4.7.(EX-ELIS, Boulder CO, USA).These bands were turned into spectral indices, according to the equations in Table 1.The indices selected to estimate the biophysical variables included bands in the water absorption SWIR region (Faurtyot and Baret, 1997)

Conclusions References
Tables Figures

Back Close
Full bands sensitive to vegetation greenness in the NIR region (Paltridge and Barber, 1988;Yebra et al., 2008b).Finally, two RTM based algorithms were run to compare with these spectral indices.GEOSAIL RTM model inversion estimated FMC testing two different look up tables that constrained the simulations to either a grassland or a mixed tree-grassland cover based on Jurdao et al. (2013).PROSAILH RTM model inversion predicted CWC assuming a pure grassland cover following Trombetti et al. (2008).This model applies the same look up table to all land covers but different calibration coefficients which will render the same R2 independently of the land cover considered.

MODIS data
MODIS Terra daily surface reflectance (MOD09GA) data from 1 April 2009 to 15 April 2011 were downloaded from NASA Land Processes Distributed Active Archive Center (LP DAAC) at the USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota, USA.This product includes the reflectance of bands 1 to 7, from 469 to 2130 nm at 500 m spatial resolution, as well as sensor and solar observation angles and quality flags at 1 km.A script programmed in Matlab (Mathworks, Batick, Massachusetts, USA) extracted the MODIS pixel value of our study site from all the images to build the time series.The impact of angular effects on reflectance was reduced by removing images with sensor zenith angles wider than 20 • , which also assures the accuracy of the geometrical location of the pixel (Wolfe et al., 2002).In addition, the quality flag layer eliminated images under clouds, cloud shadows and/or with high atmospheric aerosol content.The algorithm selected the closest valid MODIS image to the field sampling day within +5 day window, or the MODIS image acquired before the sampling day in case they were equal.Minimal time lag between sensor and field data reduces the chances of discrepancy, as grassland grazing could affect LAI in a short period of time.This leaded to a total of 14 days of MODIS data with coincident proximal sensing measurements and field data.Similarly to proximal sensing, MODIS estimated the biophysical variables from the spectral indices in Table 1 and from the RTM model inversions.

BGD Introduction Conclusions References
Tables Figures

Back Close
Full

Data analysis
Intra-group, inter-group and overall R 2 values between FMC Q , FMC E , EWT, CWC Q , CWC E or LAI, and each of the proximal sensing spectral indices were calculated to investigate their variability within the 500 m MODIS pixel.More specifically, the intragroup R 2 offers information about the spatial variability, due to the collection of samples at different plots within the MODIS pixel.A linear regression model was created for each sampling day where the biophysical variable and the spectral index were the dependent and the independent variable, respectively.The average R 2 of all the regression models for each day provided the intra-group R 2 .Instead, the inter-group R 2 explains the temporal variability due to the collection of the samples on different days.In this case, the biophysical variable and the spectral index for all plots were averaged for each day.The linear regression model of these averaged values determined the intergroup R 2 .To explain temporal and spatial variability together, the overall R 2 fitted in a single regression model including all plots and sampling days for each spectral index and biophysical variable.
Later, using the mean values of each biophysical variable and the proximal sensing spectral indices, a univariate linear regression model was applied.The same procedure was repeated for MODIS data.Bootstrapping techniques evaluated the empirical model robustness, which is a valid alternative to traditional leave-one-out methods to validate regression models predictability according to Richter et al. (2012).Root Mean Square Error (RMSE), Relative Root Mean Square Error (RRMSE), Nash-Sutcliffe Efficiency index (NSE), determination coefficient (R 2 ) and Taylor's diagrams evaluated the models' performance.As recommended in Steyerberg et al. (2001), two hundred bootstrap model simulations were run for each model and the median value of each statistics represented its performance.The RMSE measured the error in the estima-Figures

Back Close
Full tion of the biophysical variable by each model: where V i est is the estimated variable and V i obs is its observed field measurement.RMSE cannot compare the error of different variables with different units.To address this limitation in order to compare the model performances between different variables, RRMSE divides RMSE by the average of the observed values V obs (Richter et al., 2012): The NSE indicates the model predictive power which ranges from −∞ to the best predictive power value of 1 (Richter et al., 2012).It establishes if the model performs at least as accurate as the average of observed values through the following expression: The R 2 measures the proportion of variance explained by the model and is calculated as: where σ 2 r represents the residual variance and σ 2 is the variance of the dependent variable.Introduction

Conclusions References
Tables Figures

Back Close
Full Taylor diagrams allowed the comparison between the spectral indices and the RTM based algorithms of Jurdao et al. (2013) andTrombetti et al. (2008) for FMC and CWC, respectively.In these plots the observed variable and its SD are plotted in the x axes.RMSE is represented as semicircles centered at the observed data.The correlation coefficient (r) is displayed in the azimuthal position.Best models are closer in the plot to the observed measurement; therefore they will have a high r, a low RMSE and a SD similar to the observed values.

Results
All variables showed similar temporal evolution with a peak in spring and second minor peak in the fall except Dm (Fig. 2).Dm fluctuated throughout the year and exhibited its highest values in the summer.The 47 % Coefficient of Variation (CV) for Dm was less than for CWC Q (CV = 95 %), CWC E (CV = 0.95 %), FMC Q (CV = 60 %) and FMC E (CV = 56 %), but higher than for EWT (CV = 38 %).A higher precipitation in the spring of 2010 vs. the previous year translated into higher FMC, CWC and LAI values.FMC E and CWC E , calculated from the EWT Sample , presented similar trends but in some cases higher values than FMC Q and CWC Q , calculated from the quadrant sample.
A low intra-group R 2 for all the variables indicates a low spatial variability between plots (Fig. 3).Contrary, the high inter-group R 2 also for all variables points to the high temporal variability between sampling dates.The main differences between variables occurred for overall R 2 .Similar overall and inter-group R 2 values for CWC E and FMC E indicated that the combination of the temporal and spatial factors matched in importance each factor on its own.Instead, overall R 2 for CWC Q and FMC Q laid in between the inter-group and the intra-group R 2 underling the temporal factor as the main source of variation.GEMI offers the best R 2 for all variables while VARI had the weakest R 2 .The most accurate univariate empirical bootstrap model to retrieve each variable differed from proximal sensing (Fig. 4) to MODIS (Fig. 5).FMC E and FMC Q showed the best correlations with GEMI from proximal sensing data but EVI performed better 5515 Introduction

Conclusions References
Tables Figures

Back Close
Full from MODIS.EWT was more accurately estimated with GEMI from either sensor but presented the lowest R 2 and NSE of all the biophysical variables.NDII and GVMI were the most accurate predictors for LAI, CWC E and CWC Q with proximal sensing.When using MODIS, the most accurate results for LAI were achieved with NDII and GVMI, but EVI did so for CWC E and CWC Q .When the quadrant sample was used instead of the EWT Sample , both FMC and CWC showed more accurate results with lower RRMSE and higher NSE and R 2 .
Figures 6 and 7 compare spectral indices and RTM using the Taylor diagrams.In the case of FMC Q from proximal sensing (Fig. 6), RTM presented higher RMSE and lower r than the spectral indices whereas RTM SD was more similar to the observed values.
In the case of FMC Q estimated from MODIS (Fig. 6), RTM was closer to the empirical models.For CWC Q (Fig. 7), the RTM from proximal sensing overestimated the SD of the observed CWC Q .At MODIS scale (Fig. 7), RTM showed a very high overestimation of the CWC Q .
Temporal evolution of the biophysical variables estimated using the most accurate model for proximal sensing and MODIS in Figs. 4 and 5 are shown in Fig. 6.Both sensors predicted well EWT, FMC Q and FMC E but showed an overestimation, especially during the dry season.Contrary, the models for LAI, CWC E and CWC Q adjusted well even during the dry season.

Discussion
Results revealed that Dm varies significantly throughout the year (CV = 47 %) with high values in the summer.These changes could be related to the temporal variation in plant community structure, species composition and diversity in this ecosystem (Casado et al., 1986).Summer should be the best time of the year to invert However, our study suggests that, due to the high seasonal variation in Dm, a constant annual value would not be recommended here.
The high inter-group and low intra-group R 2 implies that the temporal trend is much more critical than the spatial variation within the MODIS pixel (Fig. 3).Therefore, the strategy to capture better the variability of vegetation water content in this ecosystem should be to sample more times but fewer plots.In addition, CWC Q and FMC Q , generated from larger sample sizes than CWC E and FMC E , presented higher inter-group R 2 values, which indicate a better characterization of the temporal variability.This suggests the need to standardize sampling protocols for the estimation of vegetation biophysical parameters to ensure data quality, repeatability and to facilitate accurate cross comparison from different studies.Some initiatives already exist to facilitate this standardization, as the Global Terrestrial Carbon System (GTOS) guidelines in support of carbon cycle science (Law et al., 2008).However, currently there is no international backbone that ensures this and an agreement in the protocols is needed in order to validate remote sensing products.
CWC was better predicted than the other two water content variables, FMC and EWT (Fig. 3).CWC depends on LAI which is even higher correlated than those two variables.It is possible to have the same FMC and EWT for different LAI and hence different CWC and amount of soil background, which will change its reflectance.Yebra et al. (2013) demonstrated through PROSAILH simulations how a very different CWC for the same EWT based on changes of LAI translates into a huge range of NDII values.Our results confirm this theoretical assumption described in Yebra et al. (2013).This issue is especially critical over areas like this one with large dynamic annual growth.
The empirical methods estimated FMC and CWC with slightly different results for proximal sensing and MODIS (Figs. 4 and 5).While GEMI and NDII were the most accurate for FMC and CWC respectively from proximal sensing in our study; EVI was the most accurate estimator of both variables from MODIS.In addition, it is remarkable that MODIS estimations were more accurate than proximal sensing.Roberts et al. (2006) also observed different correlations between indices and platforms and the discrep-

BGD Introduction
Full ancies here need further investigation.The MODIS pixel included not only grasslands but also trees, their shades, and other marginal covers like bare soil and a water pond (Fig. 1), and its view angle could be up to 20 • .Proximal sensing measures only two transects within each of the six plots and provides only nadiral measurements which could be more affected by the soil signal.
Similarly to this study, Casas et al. (2014) reliably predicted water content variables in California (USA) from GEMI, NDII and EVI using simulated MODIS spectral response from airborne hyperspectral AVIRIS instrument.In their case, it was actually VARI the most accurate for grasslands (FMC and CWC), chaparral (EWT, FMC and CWC) and a Mediterranean oak forest (EWT).Contrary, VARI showed very poor accuracies in our case to estimate FMC, EWT and CWC, but was still capable of capturing the variability in LAI (Fig. 3).This fact also contradicts other studies that predicted FMC from VARI on chaparral (Peterson et al., 2008;Roberts et al., 2006;Stow et al., 2005Stow et al., , 2006)).VARI was developed to detect vegetation fraction in homogenous wheat crops (Gitelson et al., 2002b), but Gitelson et al. (2002a) nor the above studies have tested this spectral index to detect vegetation water content on sites like ours, with strong seasonal changes in species composition and LAI.
The empirical methods calibrated for this specific site outperformed the physical RTM estimates for CWC and FMC (Figs. 4 and 5).This confirms the results in Casas et al. (2014) where the CWC algorithm based on RTM inversion developed by Trombetti et al. (2008) also failed to improve results from empirical estimates.RTM only overcome empirical approaches when structural information constrains the model inversion (Yebra et al., 2008b;Casas et al., 2014).Such ancillary information is key to successfully extrapolate a RTM inversion at broader scale.

Conclusions
This work showed a complete analysis of three metrics, EWT, FMC and CWC, to measure vegetation water content at two different spatial scales by using proximal sensing Figures from a field spectroradiometer and MODIS images.The temporal changes in these metrics are more critical than their spatial variation within the MODIS pixel.Results indicated that larger samples collected using quadrants are more representative than the small EWT Sample in order to follow the temporal trends in FMC and CWC.Protocol standardization in order to make different datasets comparable should be considered to make different dataset comparable both spatially and temporally.Due to the high seasonal Dm variability, a constant annual value should not be used to estimate EWT from FMC in this ecosystem.The dependence of CWC on LAI makes this vegetation water content variable easier to predict than FMC or EWT.GEMI, NDII and EVI reliably predicted vegetation water content.The best empirical estimator differed between sensors.Results were a bit better for MODIS than proximal sensing, probably due to differences in view angles, sampling strategy and canopy observed.These empirical methods still exceed RTM inversions developed for other sites to predict FMC (Jurdao et al., 2013) and CWC (Trombetti et al., 2008).Figures

Back Close
Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | It is mainly composed by Rumex acetosella L., Plantago carinata Schrad, Trifolium subterraneum L., Cynodon dactylon (L.) Pers., Taraxacum dens-leonis Desf.and Vulpia myuros (L.) C. C. Gmel.During the summer, grassland dries out rapidly and turns into dead matter.Summers are hot and dry, with 30 • C daily average temperature and only 67 mm total precipitation, which are not representative of mean annual 16.7 • C and 572 mm.The average altitude is 256 m above mean sea level.Soils are lixisols with an average thickness of 80 cm.Due to the presence of a clay layer in some of the areas, small water pools may appear in winter after rainy periods.The occurrence of this type Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 5): LAI(cm 2 cm −2 ) = Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | RTM and predict Dm, since leaves are drier and therefore EWT does mask the Dm spectral absorption signal (Riano et al., 2005).Casas et al. (2014) applied a constant annual Dm value from the literature to successfully predict seasonal variations in EWT and CWC.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 1 .Figure 4 .
Figure 1.Plots sampled near the FLUXNET tower within the 500 m MODIS pixel at Las Majadas del Tiétar (Spain) study site.