Soil Drought Anomalies in MODIS GPP of a Mediterranean Broadleaved Evergreen Forest

The Moderate Resolution Imaging Spectroradiometer (MODIS) yields global operational estimates of terrestrial gross primary production (GPP). In this study, we compared MOD17A2 GPP with tower eddy flux-based estimates of GPP from 2001 to 2010 over an evergreen broad-leaf Mediterranean forest in Southern France with a significant summer drought period. The MOD17A2 GPP shows seasonal variations that are inconsistent with the tower GPP, with close-to-accurate winter estimates and significant discrepancies for summer estimates which are the least accurate. The analysis indicated that the MOD17A2 GPP has high bias relative to tower GPP during severe summer drought which we hypothesized caused by soil water limitation. Our investigation showed that there was a significant correlation (R = 0.77, p < 0.0001) between the relative soil water content and the relative error of MOD17A2 GPP. Therefore, the relationship between the error and the measured relative soil water content could explain anomalies in MOD17A2 GPP. The results of this study indicate that careful consideration of the water conditions input to the MOD17A2 GPP algorithm on remote sensing is required in order to provide accurate predictions OPEN ACCESS Remote Sens. 2015, 7 1155 of GPP. Still, continued efforts are necessary to ascertain the most appropriate index, which characterizes soil water limitation in water-limited environments using remote sensing.


Introduction
Water availability has been shown to have dominant or co-dominant effects in the productivity of most terrestrial biomes on Earth [1,2].This is particularly true today of Mediterranean-type (MT) climate areas and will become crucial with global climate change [3,4].In these areas, natural vegetation has to cope with a strong seasonality of environmental conditions.The five relatively small, isolated regions of the world with Mediterranean-type climate, i.e., respectively parts of the Mediterranean Basin, California, South-Western and Southern Australia, Central Chile and Southern Africa, are approximately located between latitudes ranging from both 30° to 43° North and South with a total area of about 2.75 Mkm 2 .The Mediterranean-type climate can be coarsely considered as a transition between dry tropical and temperate climates.They are characterized by a distinctive annual climate sequence in which a hot dry summer alternates with a cool to cold, humid period lasting of 5-10 months from fall through winter to spring.In these five regions, vegetation is dominated by short to tall shrubs or trees with small, evergreen, sclerophyllous leaves [5].These shrubland and woodland formations are fairly similar in their physiognomy and provided the ideal testing vegetation type for studying low productive ecosystems with evergreen cover that undergo severe water limitations during the growing season.They also constitute a good model region for evaluating the generality of plant functions in the so-called evergreen broadleaf forests (EBF) of the Global International Geosphere-Biosphere Programme (IGBP) land classification that grow in some subtropical areas.
Terrestrial ecosystems affect carbon dioxide (CO2) fluxes to and from the atmosphere through an exchange between gross primary production (GPP) and whole respiration.The net balance determines whether an ecosystem is sequestrating carbon (C) or releasing it into the atmosphere.This balance regulates C sink strength or net ecosystem productivity and changes continuously with climate drivers.Forest ecosystems are a major sink in the global carbon cycle, sequestrating large amounts of carbon annually and thus slowing the rise of CO2 concentration in the atmosphere.Forests store ca.(standard abbreviation for latin allocution "circa", which means "around") 45% of terrestrial carbon and are responsible for half of terrestrial net primary production (quoted by [6]).Unfortunately, measuring and describing how trees whatever their functional type transfer C from the atmosphere into terrestrial biomass remains poorly quantified.We observe a lack of mechanistic understanding on this question and [7] suggest developing generalizable models of C allocation to biomass growth of plant parts, respiration, nonstructural C reserve, reproduction and defense, as a challenging endeavor.
Global modeling, remote sensing and ecosystem flux derived from eddy-flux towers have demonstrated that climate and climate variations influence biogeochemical processes that control net ecosystem productivity at multiple timescales.The C sink strength is affected by differential sensitivity of respiration components and storage in both short-and long-lived C pools with environmental drivers.
The MOD17A2 product is one of the primary sources of remote sensing-based GPP at global scale.It provides an 8-day mean GPP at 1 km spatial resolution for the entire land surface [8].However, several recent studies have highlighted limitations of this model [9][10][11].GPP MOD17A2 (further called GPP MODIS) is an extension of the light use efficiency (LUE) model of Monteith previously used for estimating net primary production [12,13], where minimum temperature Tmin and vapor pressure deficit VPD were retained as the two scalars directly modifying maximum LUE, LUEmax [14].The limitation arises from: (1) the uncertainties of meteorological reanalysis used as input data, (2) the realism of LUE attributed to a given pixel derived from a biome properties look-up table (BPLUT) and more crucially (3) the ability of both scalars to correctly represent the abiotic controls for any particular ecological conditions.In MODIS, estimates of LUE are obtained from lookup tables based on vegetation type, which may contain errors either in the assignment of vegetation type to a pixel or in the initial estimates of LUEmax.
We tested for the accuracy of this simple model on an evergreen forest largely dominated by Quercus ilex where flux data are available.Q. ilex is an emblematic Mediterranean tree species growing over more than 6 Mha around the Mediterranean Sea [15].For testing we will use long-term GPP derived from an eddy-flux tower (further called EC GPP) set up in our research site.This site underwent recurrent summer droughts and demonstrated the key-role of soil water deficits on its carbon exchanges [16].We suspect biased estimates of the active season GPP impacted by summer drought because VPD seems an insufficient surrogate of the soil dryness in such MT climate.The main objectives of this work are: (1) to compare the time course of EC GPP and to those derived from MODIS, (2) to assess the anomalies of MODIS GPP to between-year variations in drought severity; and (3) to quantify the limitations of such model.

Study Area
The study site is located 35km north-west of Montpellier (southern France), on a flat plateau in the Puéchabon State Forest (3°35′45″ E, 43°44′29″ N, 270 m a.s.l.).This forest has been managed as a coppice for centuries with the last clear cut performed in 1942.Vegetation is largely dominated by a dense overstory of the evergreen oak Quercus ilex.The top canopy height is approximately 5.5 m.In 2010, the density of the resprouted stems was 6070 stems•ha −1 .Stems with diameter at breast height (DBH) < 4 cm represented 6% of total stems, whereas stems with DBH > 10 cm represented 20.6%.Understory species Buxus sempervirens, Phyllirea latifolia, Pistacia terebinthus and Juniperus oxycedrus, compose a sparse shrubby layer with a percent cover less than 25% and height less than 2 m.
The area has a Mediterranean-type climate.Rainfall mainly occurs during autumn and winter, with almost 80% of the precipitation taking place between September and April.The mean annual precipitation (MAP) is 903 mm, with a range of 556-1549 mm recorded for the 1984-2010 period.Mean annual temperature (MAT) over the same period was 13.0 °C, with a minimum in January (5.5 °C) and a maximum in July (22.9 °C).The soil is extremely rocky consisting of hard Jurassic limestone origin; on average, the volumetric fractional content of stones and rocks is about 0.75 for the top 0-50 cm and 0.90 deeper soils.The stone-free fine fraction of the soil is a homogeneous silty clay loam (USDA texture triangle) within the top 0-50 cm layer (38.8% clay, 35.2% silt and 26% sand).The fine fraction fills up the space between stones and rocks and thus provides a source of stored water throughout the long dry summers for the deep-rooted Q. ilex.The highly permeable soil prevents any surface runoff to occur even during downpours.

Flux and Ancillary Data
As a CarboEurope and Fluxnet site, eddy covariance fluxes of CO2, sensible heat, latent heat and momentum have been measured continuously since 2001 from the top of a 12 m high tower that is approximately 6 m above the canopy.Our eddy covariance facility included a three-dimensional sonic anemometer (Solent R3, Gill Instruments, Lymington, England) and a closed path infrared gas analyzer (IRGA, model LI 6262, Li-Cor Inc., Lincoln, NE, USA), both sampling at a rate of 21 Hz.Flux data were processed with protocols in accordance with the CarboEurope network [17].Processing schemes of Fluxnet have been used for filling data gaps and partitioning net ecosystem productivity into GPP and ecosystem respiration [18,19].Half-hour GPP fluxes were summed at a daily time step for further comparison with MODIS GPP and called EC GPP.Half-hourly data of air temperature, VPD, incoming flux of photosynthetically active radiation PARtop were recorded at the top of the flux tower.The fraction of PAR absorbed by the canopy FPAR was derived from 14 PAR-sensors randomly set up in understory locations measuring PARbelow as:

MODIS Land Products Subsets
Collection 5 MODIS land product subsets for Fluxnet sites are readily available from the DAAC (Distributed Active Archive Center) at Oak Ridge National Laboratory [20].These subsets represent MODIS product cutouts of 7 km by 7 km around the eddy flux towers, greatly simplifying data availability.We used both collection 5 MOD15A2 (FPAR canopy product [21]) used for computation of the collection 5.1 MOD17A2 (GPP) for the period ranging from 2001 to 2010.For Puéchabon, we selected all pixels from the 7 km × 7 km subset to represent the eddy tower footprint (see quality analysis and control for our site in [22] and [23]) with its mean and standard-deviation values because our site is representative of vegetation cover that extends largely across our study region [24] considering the surrounding area reduce errors with gridding artifacts and with spatial data with nested source of variations [25].
where LUEmax is the biome-specific maximum light use efficiency (g•C•MJ −1 ), PAR is the photosynthetically active radiation (MJ), and the potential GPP is linearly scaled between 1 (no limitation) and 0 (whole limitation = no GPP) using biome-specific functions for minimum temperature, f(Tmin) and vapor pressure deficit, f(VPD).The fraction of PAR absorbed by the vegetation cover FPAR is from MOD15 while meteorological inputs are from large spatial-scale meteorological data, obtained from NCEP-DOE Reanalysis II dataset (see supporting online material in [26]).Due to acknowledged effects of input climate data on MODIS17A2 GPP [21,27,28], we also tested for f(VPD) sensitivity to bias in climate data source (NCEP vs. tower climate).

Modeling Soil Water Balance and Daily Course of Relative Water Content
Soil water storage integrated over the rooting depth (SWSobs), that is ca.4.5 m [29], was measured, during the vegetative periods, at approximately monthly interval, with a neutron moisture gauge in 1984-1986 [30] and since July 1998.Discrete measurements were interpolated at a daily time step with a soil water balance model proposed in [31] and further used in [32].Soil water storage and soil water potential were related by a Campbell-type retention curve [33] whose parameters are strongly dependent on soil texture [16].The drainage curve relating depth drainage to soil water storage depends on the stone content over the whole-soil profile [34].The model was driven by daily values of incoming solar radiation, minimal and maximal temperature and rainfall amount.Comparison of measured against simulated soil water storages, in mm, and predawn leaf water potentials, in MPa, displayed very good agreement.Leaf water potential values came from discrete measurements performed on the study site [35].Reduced major axis (RMA) regressions yielded SWSsim = RMA SWSobs + RMA with RMA ± standard-error (SE) = 0.94 ± 0.03, RMA ± SE = 6.0 ± 4.4, R 2 = 0.93, F = 1137, p < 0.0001 and n = 91.For the predawn potential, pdsim = RMA pdobs + RMA with RMA ± SE = 0.93 ± 0.05, RMA ± SE = −0.09± 0.09, R 2 = 0.840, F = 273.3,p < 0.0001 and n = 54.The continuous daily course of soil relative water content, RWC, was derived from SWSsim divided by the soil water storage at field capacity.We choose to fix it at 205 mm.It corresponds to the value we observed after 2 days of free drainage after a substantial rain event falling in cool wet period.For characterizing the whole-year water limitation, we calculated the water stress integral (WSI) as the yearly sum of pdsim.For days with RWC  1, pdsim is fixed to −0.03 MPa.The WSI are expressed in MPa day.

Annual and Seasonal Pattern of GPP: Comparing MOD17A2 with EC
Yearly correlations coefficients (R) between the 8-day means of MODIS and EC GPP at the Puéchabon study site from 2001 to 2010 are presented in Table 1.The yearly correlation coefficients (R) vary between 0.36 in 2006 and 0.87 in 2004, with a mean correlation of 0.67 for the whole period, this somehow is a weak correlation when compared to other non water-limited ecosystems [36].
Mean annual GPPs are 1240 g for respectively EC and MODIS.The yearly totals of MODIS GPP consistently exceeded the tower averages, with a relative error of 41%.When comparing cumulated annual GPP on a yearly basis, annual EC GPP varied between a minimum of 945 g ), highlighting nonlinear discrepancies between the two datasets (Figure 1).Intra-annual GPP pattern is highly seasonal with maxima at the end of spring and early summer, minima in winter as shown in Figure 1.In spring (DOY 65-152), we observed a constant increase of GPP in both datasets before reaching the highest values at the end the season.For most years between 2001 and 2010, MODIS GPP starts to overestimate EC GPP in the middle of spring.In summer (DOY 153-248), MODIS GPP reaches its highest value while EC GPP decreases almost linearly and usually reached its lowest point at the end of summer.MODIS GPP and EC GPP exhibited the maximum difference in 2006 (6.21 ).In fall (DOY 249-335), the first period is still influenced by summer discrepancies between the two datasets which progressively reach the same final GPP pattern with a smooth decrease finally reaching the low winter values (around 1 g•C•m −2 •d −1 ).Indeed, MODIS GPP was fairly low-correlated with EC GPP for the summer season, with correlation coefficient varying between 0.07 in 2007 and 0.12 in 2008.In winter (DOY 0-64 and 336-366), both MODIS GPP and EC GPP have similar low values with a higher correlation coefficients, except for the year 2006 (R = 0.25).We will note here that the MODIS GPP product consistently underestimated GPP during the winter season compared to the other seasons.As we observed significant differences in the seasonal pattern of MODIS GPP compared to measurements from the tower, we explored the potential candidates for this bias based on variables used in Equation ( 2).

Annual and Seasonal Pattern of FPAR: Comparing MOD15A2 with Tower
MOD15A2 FPAR provides a key and seasonally-varying input to the MODIS GPP algorithm, and is derived from remote sensing.As FPAR is also measured on site with no effect from atmospheric components, we tested for its accuracy in the MODIS GPP calculation.The FPAR estimates from MODIS ranged from 0.75 to 0.92, in fair agreement with the corresponding measurement at the tower with values ranging from 0.77 to 0.97.The seasonal pattern is moderately correlated between the two datasets (Figure 2 and Table A1), with an average relative error (RE = 100 × (MODIS − EC)/EC) of 6%.REs reached up to 16% in autumn or winter.This low seasonality represented by FPAR for this type of evergreen forest was also evidenced by [37].The tower FPAR pattern decreases in mid-summer and reaches a maximum in winter; in contrast, the MODIS FPAR exhibits a reversed seasonal pattern, with the maximum FPAR reached in summer.Indeed, MODIS FPAR displays good accuracy of prediction during the spring and summer periods with RE varying between 0 and 9% (Table A1).The annual maximums of RE are all concentrated in winter when GPP is low (due to low air temperatures during this period which leads to low photosynthetic activity), showing that errors in FPAR MODIS are a weak contributor to anomalies in GPP MODIS.

Climate and Drought Patterns
As the main driver of GPP, we analyzed the seasonal pattern and between-year variability of climate variables registered at the flux tower on a half-hour basis.We also used the subsequent soil water status from a water budget model as an index of drought which is identified as the main constrain for carbon assimilation in the dry season for water limited ecosystems [38].
Rainfall amounts and temperatures were used in the analysis for the 2001-2010 periods (Figure A1).Rainfall amounts present both high seasonal and inter-annual variations, with annual amounts varying between 1234 mm in 2003 and 570 mm in 2007.Mean rainfall in spring (184.3 mm) is higher than in summer (97.5 mm), with the most extreme dry summers occurring in 2006 (37 mm)  We also analyzed vapor pressure deficit (VPD), minimum temperature (Tmin) and photosynthetically active radiation (PAR), as drivers for LUE model calculation in MODIS GPP.The three variables show similar yearly patterns (Figure 3).The annual average value of VPD is 6.7 hPa, varying between 2.5 hPa in winter and 15.5 hPa in summer.Tmin usually drops below 5 °C in winter to reach 18.6 °C in summer, a seasonal pattern in phase with VPD variations.Table A2 presents the interannual variability of yearly mean and maximum PAR for the 2001-2010 periods, with a coefficient of variation of 10% showing little interannual variability.Regarding the seasonal pattern of PAR, we observed a time lag of 25 days in reaching the maximum values compared to Tmin.The maximum value of PAR is reached between DOY 150 and 200, while it is reached between days 175 and 225 for VPD and Tmin.
Furthermore, we characterized daily drought intensity susceptible to affect ecosystem functioning by the soil relative water content (RWC) and the yearly integral of water stress (WSI) of Quercus ilex ecosystems by capturing plant water deficit through daily estimates of predawn water potentials, and also integrating soil properties and the species functional traits controlling water extraction (leaf area index, root profile, stomatal conductance).When comparing the seasonal pattern of VPD controlling MODIS GPP algorithm with soil RWC (Figure A2), we observed a concomitant seasonal signal with highest VPD in summer when RWC is the lowest, but with a time lag > 20 days for the maximum values, except for 2003 and 2004 where the two peaks coincide in time.The VPD trend is, however, much smoother than RWC because heavy rainfalls at the end of the dry season lead to a rapid increase in the soil water content to values close its field capacity.The seasonal day of occurrence of maximum Tmin, maximum VPD and minimum RWC are presented in Table A3.We observe that the day of occurrence for the lowest value of tower GPP is generally correlated with the day of occurrence for the minimum RWC timing rather than both VPD and Tmin, suggesting that RWC could be a better control for GPP estimates.

Control of RWC on GPP MODIS Anomalies
Before testing our hypothesis that RWC would be a better GPP control than f(VPD), we checked if f(VPD) used in MOD17A2 GPP is affected by potentially biased climate inputs from CRU NCEP or a biome type misclassification as already pointed out by [21,27,28] (Figure A3A,B).We obtained enhanced f(VPD) effects with the tower climate data when compared to NCEP, leading to lower root mean square error (RMSE) between MODIS/EC GPP anomalies and f(VPD) (RMSE = 0.247 for NCEP and RMSE = 0.202 for tower climate).Using the mixed forest vegetation type (MF) parameters also increased the performance of f(VPD) for GPP estimates with RMSE = 0.214 with NCEP and RMSE = 0.179 for tower data (Figure A3B).However, the temporal matching between MODIS/EC GPP anomalies and f(VPD) remained low (year 2008 Figure A3A for example).We then tested for the relationship between soil RWC and GPP errors observed between MODIS and EC (Figure 4).We obtained a significant relationship following a non-linear decay (Equation (3)): with a = 36.7,b = 0.39, c = −70.2,R 2 = 0.77 and p < 0.0001.The curve exhibits a vertical asymptote for RWC  0.4, and the RE starts to deviate from 0 when RWC = 0.7.The larger anomalies of GPP were mainly concentrate in the 0.4-0.5 interval of RWC.When comparing RE with GPP anomalies between MODIS and EC, we obtained RMSE = 0.162 (Figure A3A), and RMSE = 0.149 when combined with f(VPD) calculated with the tower climate, both outperforming the f(VPD) correction alone.We finally used the yearly WSI, derived from daily simulation of soil RWC and ecosystem water budget, to characterize the dry season for each year in terms of both the length and intensity of drought.
We investigated this annual drought index as an empirical driver controlling the anomaly between MODIS GPP and the actual tower GPP cumulated for each year, by plotting these two variables in Figure A4.We obtained a significant linear relationship (RE = −26.65 − 0.4WSI, R 2 = 0.71 and p = 0.0021) between the relative errors of annual GPP (between MODIS and EC), and WSI.We will note that the years 2003 and 2006 were the most out-of-the-trend years, both years respectively corresponding to an extreme heat wave where high temperatures might have abnormally modified carbon assimilation rates, and the driest year of the period with very low RWC for which water potential calculations are very sensitive and potentially biased.

Discussion
Mediterranean-type ecosystems cover a small fraction of the terrestrial component of the Earth.They are often neglected in Global modeling exercises or included in vegetation classes that comprise various vegetation types growing under contrasted abiotic and biotic constraints.Nevertheless, Mediterranean-type ecosystems remain a good system model for studying broadleaved evergreen plant covers that have to tackle recurrent severe droughts during their vegetative season and under the threats of the on-going climate changes.Among the dominant key-species in this area, Quercus ilex trees are present across more than 6 Mha around the Mediterranean Sea.They occur in dense woodlands, such as the one we analyze here, as well as open woodlands or even savanna-like formations.
Concerning the climate constraints, all simulations performed with high resolution coupled atmosphere-ocean regional climate models or with coarser resolution GCMs under SRES A2 and A1B scenarios gave a collective picture of a substantial drying and warming of the regions surrounding the Mediterranean Sea, especially in the warm season [4,39,40].The comparative analysis of [41] placed this region among the most responsive regions to global climate change, and one of the main climate change "hotspots".This is also true for California.Global changes have the potential to deeply modify local climate patterns mostly by lowering spring soil water availability leading to reduced summer convection and large land-sea contrasts in warming.All this is leading to increased air dryness and decreased precipitation [42].As a result, we are expecting a warming from +3.2 °C in Winter (DJF) and exceeding 5 °C (+5.2 °C) in Summer (JJA).For our study site, precipitation projections yielded decreases of ca.0.1, 0.4, 0.5 and 0.35 mm•day −1 for DJF, MAM, JJA and SON, respectively [4].The consequence for the Spring MMA rains will be a significant decline of 18%.More drastically, the Summer JJA amount, currently 108 mm, will lose 45 mm with expected negative impacts on ecosystem functions or more frequent perturbations and particularly fires [43][44][45].Inter-annual variability is projected to increase as the occurrence of extreme heat and drought events [46,47].Return periods for drought durations lower than 4-6 months will be multiplied by 3 while return periods for drought durations longer than one year will be multiplied by 7 [3,48].
In such an already water-limited climate going to drier conditions, we have to focus on the mechanisms that facilitate or constrain the ability of Mediterranean terrestrial ecosystems to cope, adapt or recover from various disturbing forces and particularly to severe droughts.Factors related to climate change are likely to play a key role in predisposing ecosystems to tipping points [49][50][51], while vulnerable ecosystems are already being influenced by multiple environmental drivers.This reinforces a general view that synergisms among stressing elements can be extremely important, predisposing ecosystems to serious environmental changes.GPP is the first step of C-use cycle which has been shown to be drastically affected by heat waves [9] and drought [16,52,53].GPP may be useful for evaluating regional carbon budget while also contributing to detect early warning sign of dysfunctions.

Soil Drought as the Main Driver of GPP Anomalies
The comparisons of EC GPP and MODIS GPP have already been performed across a large range of climate, land use, and vegetation class and have been found to provide rather good agreements [10,54,55].Discrepancies have been explained by the validity, or not, of internal hypothesis concerning either the maximum LUE derived from the BPLUT or the response sensitivities to both scalars Tmin and VPD.In our case we observed good agreements between the two cool wet seasons.The underestimations of FPAR during these periods (Figure 2) may explain why EC is slightly greater than MODIS.The Tmin scalar seems effective.It increases linearly from zero at −8 °C to reach 1 at 9.1 °C.This scalar is rarely questioned as a potential source of error in the literature.The VPD scalar does not act during these periods because it begins to control GPP at only 11 hPa and across the cool periods' 8-day mean, VPD were always lower than this threshold (Figure A2).
At the peak of productivity, EC GPP was always lower than MODIS GPP.We suspected first an overestimation of LUEmax.These over or underestimations have been already discussed as the main source of errors in [55], [10] or in [56].Leuning, R. et al. [36] observed this overestimation by comparing multiyear EC measurements between two contrasted ecosystems in Australia, an Eucalyptus spp.dense forest and a tropical savanna both ecosystems dominated by broadleaved evergreen tree species.Kanniah, K.D. et al. [57] carried out an extensive analysis showing the strong control of rainfall on the inter-annual variation of GPP as estimated using a LUE-based GPP model in a water-limited region.As [38] did further, they suggested applying site-specific estimates.Kanniah, K.D. et al. [56] proposed a modifier based on evaporative fraction (EF) to replace VPD in water limited regions.They found good agreement between tower GPP and GPP estimated using site specific parameters and EF via a LUE model.In our case, for the two well-watered, highest WSI years, 2002 and 2004, the lowest discrepancy was observed in 2002.This year displayed a very limited effect of VPD and an insignificant control by soil water limitation (Figure 5).At a VPD of 15 hPa the scalar is affected by only 14%.At 25 hPa it declines by half.The more drastic effects on peak GPP was observed in 2006 the driest year of our study period with a WSI = −358.6MPa day.This year revealed how early severe drought that started in February is not taken into account by the VPD scalar.Interestingly, in 2003, a heat wave took place from June to August.It was associated with a period of high VPD in phase with the lower relative soil water storage.At this level of VPD, the scalar drastically declined by 80%, but insufficiently to closely mimic EC GPP.In a MT climate, there exists a substantial delay between peak air VPD and lower soil RWC (Table 4).As a consequence, this delay, even changing the threshold values in the VPD scalar, prevents the direct use of VPD as an efficient surrogate of soil water limitation.
Among others, there is an increasing debate about the use of MOD15A2 FPAR, both in terms of quality of this dataset when compared to other available products, and in term of effective radiation absorbed by chlorophyll for photosynthesis (see for a substantial account [58][59][60][61][62] and also [63][64][65]).The widely used FPAR is in turn a canopy level index, accounting for both the photosynthetic active vegetation and the non-photosynthetic parts including stems, branches or senescent leaves.This debate also applies to ground estimates of the leaf area index [66] and biases linked to remote sensing signals in general [67].Increasing efforts are being devoted to obtaining more accurate information concerning the vegetation photosynthetic capacity retrieved from space borne measurement [68][69][70][71].Now, when focusing on our study site, and quantifying the potential errors related to the FPAR canopy instead of fraction of photosynthetic active radiation absorbed by the photosynthetic elements, we are not alarmed by this bias as GPP MODIS and measured GPP are actually in close agreement during the non water-limited periods.In addition, as our model explained the error rate with R 2 = 0.71, we can also attribute to an additional error rate of 0.29 for the FPAR bias.In turn, we could also claim that refining drought is twice more important than the bias due to erroneous FPAR.The comparison of different GPP outputs when using FPAR or refined photosynthetic active interception in [64] illustrates large discrepancies in croplands/grassland test sites, but a lower bias for evergreen or mixed forests.
Leuning, R. et al. [36] were the first to report that MODIS GPP was significantly overestimated during the water-limited periods, but gave reasonable estimates during the wettest seasons in the savanna ecosystem where rainfall amount and timing exclusively controls plant productivity.In the Eucalyptus spp.forest, they found that rainfall was a key-control factor in explaining inter-annual variations of productivity, while evaporation fluxes were less affected by drought than carbon uptake.The lack of a soil water limitation term in MODIS GPP algorithm may result in significant overestimation in productivity particularly during extreme drought.Kanniah, K.D. et al. and Hwang, T. et al. [56,72] advocated for the development of such soil water modifier.Leuning, R. et al. [36] and Pan, Y. et al. [73] examined the possible benefits of modifying the MODIS algorithm by adding simple water balance scalar from the ratio of antecedent rainfall and potential evapotranspiration PET data.Leuning, R. et al. [36] proposed a modifier based on the ratio of the sum of rainfall on the sum of PET over a given time window.Their estimate of PET is the equilibrium evapotranspiration and the optimal time window is three months.Further, Coops, N.C.et al. [74] successfully applied this modifier with the same integration period for a Douglas-fir needle leaf forest.A similar time window of 100 days is also adopted by [75] in developing their new drought index for grasslands, the "local dryness".Alternatively, Pan, Y. et al. [73] proposed a soil water correction index at a yearly scale.They first calculated the monthly course of soil water availability with a soil bucket model receiving monthly rainfalls minus PET.For evergreen coniferous, their correction index is integrated over the whole-year because they assumed that photosynthesis occurs throughout the year.For deciduous forests only the growing season is concerned.All found that the modifiers significantly improved the predictions.
For our Q.ilex forest, we develop a modifier depending on daily values of the soil RWC integrated over the rooting depth that we further average over eight days to facilitate comparison with MODIS GPP.This modifier explains more than 75% of the variance (Figure 4).The use of soil water content as driver for a modifier has already been suggested by [74].They observed that [36]'s modifier paralleled the eight-day averaged measurements of the relative extractable soil water content REW over the 0-60 cm layer.This layer contains most of the roots of the dominant species of their study forest, Pseudotsuga menziesii.Interestingly, we observed two characteristic values in the course of GPP anomalies against RWC.First there exists a plateau for values of RWC < 0.7 during which we do not observe significant effect of soil water limitation.Most of the errors sources from the previously discussed parameters.Finally, we observed an asymptote for RWC = 0.4.For 0.7 < RWC < 0.4 large anomalies may be related to the soil water limitation.Both values, 0.4 and 0.7, are in agreement with [76], who found in a dry Pinus radiata plantation that the canopy conductance zeroed at RWC = 0.35 with a plateau till 0.6-0.65 and that C-flux seemed affected at higher RWC.Similar values have also been obtained in [16] for Q. ilex ecosystem which concords with our results.Other publications papers validated our results.Granier, A. et al. [77,78] used a generic soil water balance model for studying both water and C-fluxes from numerous forest ecosystems ranging from boreal to Mediterranean climates.They found that GPP were dramatically reduced during the drought when soil relative extractable water (REW) dropped below a threshold comprised between 0.35 and 0.4, a threshold independent of tree species and soil type.If we fix the zero GPP value at 0.4 as the limit for extractable water, we obtain close results to those derived in using REW.

How Does Circumvent the Problem of Soil Drought Limitation?
A central question still remains unsolved: How to include soil drought limitations in the MODIS GPP algorithm.As seen previously, FPAR is relatively easy to evaluate from remote sensing.On the other hand, LUE has been proved to be more difficult to estimate because: (1) it is based on large spatial-scale meteorological data available from the NCEP-DOE Reanalysis II dataset [21]; (2) LUEmax and parameters of the climate scalars are obtained from lookup tables on the basis of a limited number of vegetation types and (3) soil drought is not included in the calculation.Three axes have been followed to tackle the deficiency of MODIS GPP in soil water-limited ecosystems.The first one is to find remotely-sensed proxies of GPP, the second one is to find a remotely-sensed surrogates of LUE and finally the third one is to develop parallel algorithms that estimate evaporation fluxes and derive a soil water limitation scalar or directly measure the soil water storage.
It would be easier and more direct if we could base GPP of water-limited ecosystems only from remotely-sensed data and thus have continuous estimates at the spatial resolution of the satellite data.A lot of works concerning a large range of vegetation types have been interested in the extent to which GPP could be estimated directly from greenness indices such as the enhanced vegetation index (EVI) without direct estimation of LUE [79,80].This approach has been shown to be very useful in up scaling EC fluxes at continental scale [81].Unfortunately, poor correlations between EVI and EC GPP have been identified for water-limited sites [11,82].Moisture vegetation indices (MVI) combining NIR and SWIR in the middle infrared interval wavelengths (1200 to 2100 nm) have also been widely investigated to assess drought conditions as inter-annual precipitation surplus/deficit [83], or more specifically its impact on vegetation moisture content and water stress [84,85] and in turn fire risk [86], but with decreasing efficiency when tree cover increases [87].Differential plant responses to soil water deficit and desiccation rates however make it difficult to directly relate leaf moisture content to stomatal closure and subsequent controls on GPP.Some studies suggested that inclusion of MODIS land surface temperature (LST) or the land surface water index LSWI [11,82] can partly address this limitation.LST could be used to restrict the active vegetation period and provide a measure of water limitation through its correlation with vapor pressure deficit in climates where VPD pattern captures the seasonality of soil drought.
In EBF ecosystems including our site, the photochemical reflectance index (PRI) has been successfully applied [24,88] and constitutes a promising basis for further developments.PRI has been shown to correlate well with LUE.A strong limitation of PRI is that the originally proposed reference band for PRI is not available on MODIS.We tested the reference bands 1 (620-670 nm), 4 (545-565 nm), 12 (546-556 nm), 13 (662-672 nm), and 14 (673-683 nm) [24].MODIS spectral band 1 turned out to be the most suitable reference band, followed by the narrow red bands 13 and 14.The strongest correlation between LUE and PRI was also found when considering only a narrow range of viewing zenith angles at a time reducing drastically the potential availability of useful data.However, this indicates that, at site level, MODIS-based PRI is very competitive as a proxy for LUE.Despite the potential advantages of using PRI to estimate LUE at site-level, we could not establish a universally applicable LUE model based on MODIS PRI.Models that were optimized from a pool of data from several contrasted sites did not perform well [89].
Rather than researching remotely-sensed surrogates of LUE, the last option is to evaluate evaporation flux in parallel and derive a water limitation index of the vegetation cover or measure directly the soil water storage.Development of global evapotranspiration algorithms based on both remotely-sensed MODIS and global meteorology data may provide method to better express ecosystem water limitation without using precipitation amounts or soil moisture.Previous algorithms [54,90] have been further improved and yielded satisfying estimates of actual evapotranspiration [91].Yuan, W. et al. [92] proposed a water stress factor or evaporation fraction as the ratio of latent heat flux to the sum of latent and sensible heat fluxes.
Existing coarse resolution soil water content data from active and passive microwave sensors were found beneficial for climatology and hydrology at global to regional large scales (see for instance [93][94][95]).It can be anticipated that further applications would become feasible with medium (<1 km) scales.This downscaling procedure may be critical.These include forest and shrubland ecosystems and soil moisture measurements over the complex landscape mosaic we observed in MT areas.The benefit of a C-band Advanced Synthetic Aperture Radar has already been demonstrated for mapping top soil moisture generally over few centimeter depths [96].For instance, Sentinel-1 will carry onboard a C-band radar instrument [97] with a configuration and high temporal sampling rate providing great interest for the operational soil water storage estimation we need to calculate a drought scalar for the GPP algorithm.
However, land surface features such as dense forest cover, surface roughness and rock outcrops significantly limit its sensitivity [97].In parallel, considerable efforts are required in order to have a reasonable access to the whole water storage covering the rooting depth for deep-rooted tree or shrub species [98,99].Finally, in both cases, evapotranspiration or soil moisture estimates, we advocate the use of a multi-sensor approach to address the practical problem of introducing a scalar describing the soil drought limitations in the MODIS GPP algorithm [100].

Conclusions
We found that GPP as modeled from the MOD17A2 algorithm is biased during the summer period in a water-limited Mediterranean climate.We screened the seasonal or within-year and between-year patterns of contributive variables and identified that relative soil water content appears to be a useful modifier for MOD17A2 GPP in summer and for the Puéchabon Quercus ilex forest, despite the influence of FPAR, PAR, Tmin and VPD on daily GPP variation.Anomalies in MOD17A2 GPP compared to actual measurements are more pronounced during summer drought rather than the vapor pressure deficit due to high temperature leading to an overall yearly 41% overestimate.Large scale applications of this result require information on plant water status resulting from climate, soil properties and plant functioning or integrated estimates from remote sensing.A1).
2006 and a maximum of 1472 g•C•m −2 •y −1 in 2004.Extremes in cumulated annual GPP for MODIS follow a different timing and intensity, the lowest value being reached in 2003 (1637 g•C•m −2 •y −1 ) and the highest value in 2007 (1915 g•C•m −2 •y −1
and in 2003 (43.2 mm).Autumn and winter are the wettest seasons, with an extremely wet autumn in 2003 accounting for 82% of the total annual amount.Temperatures are also highly seasonal with hottest temperatures reached during the dry season (summer) and varying between 24.6 °C in 2003 and 20.6 °C in 2007.Absolute maximum daily mean temperature can reach up to 40 °C and even more, as in the extreme heat wave of summer 2003.Mean winter temperatures vary around 6 °C.Therefore the climate in the Puéchabon forest is typical Mediterranean characterized by cool, wet winters and hot, dry summers and experiences a strong inter-annual variability of most climate variables.During the 2001-2010 periods, the most extreme events were excessively high temperatures as in the summer of 2003 and sustained moisture deficits in spring and summer of 2006, two key years for GPP anomaly analysis.

Figure 4 .
Figure 4. Nonlinear regression between the relative soil water content and the weekly relative error of GPP between MODIS estimates and measurements at the tower flux.

Figure A2 .
Figure A2.Comparison between 8-day mean vapor pressure deficit (VPD, in blue) and 8-day mean relative soil water content (RWC, in red) from 2001 to 2010.Blue dashed lines represent the yearly maximum for VPD and red dashed lines represent the yearly minimum value for RWC.

Table 1 .
Correlation coefficients (R) between 8-day means of MOD17A2 GPP and EC GPP from 2001 to 2010.