Evaluation of root water uptake in the ISBA-Ags land surface model 1 using agricultural yield statistics over France 2 3

The simulation of root water uptake in land surface models is affected by large uncertainties. The difficulty in mapping soil depth and in describing the capacity of plants to develop a rooting system is a major obstacle to the simulation of the terrestrial water cycle and to the representation of the impacts of drought. In this study, long time series of agricultural statistics are used to evaluate and constrain root water uptake models. The inter-annual variability of cereal grain yield and permanent grassland dry matter yield is simulated over France by the Interactions between Soil, Biosphere and Atmosphere, CO 2 -reactive (ISBA-A-gs) generic land surface model (LSM). The two soil profile schemes available in the model are used to simulate the above-ground biomass ( B ag ) of cereals and grasslands: a two-layer force–restore (FR-2L) bulk reservoir model and a multi-layer diffusion (DIF) model. The DIF model is implemented with or without deep soil layers below the root zone. The evaluation of the various root water uptake models is achieved by using the French agricultural statistics of Agreste over the 1994–2010 period at 45 cropland and 48 grassland departements, for a range of rooting depths. The number of departements where the simulated annual maximum B ag presents a significant correlation with the yield observations is used as a metric to benchmark the root water uptake models. Significant correlations ( p value < 0.01) are found for up to 29 and 77% of the departements for cereals and grasslands, respectively. A rather neutral impact of the most refined versions of the model is found with respect to the simplified soil hydrology scheme. This shows that efforts should be made in future studies to reduce other sources of uncertainty, e.g. by using a more detailed soil and root density profile description together with satellite vegetation products. It is found that modelling additional subroot-zone base flow soil layers does not improve (and may even degrade) the representation of the inter-annual variability of the vegetation above-ground biomass. These results are particularly robust for grasslands, as calibrated simulations are able to represent the extreme 2003 and 2007 years corresponding to unfavourable and favourable fodder production, respectively.


Introduction
Modelling the land surface processes and the surface energy, water and carbon fluxes is an important field of research in the climate community, as soil moisture and vegetation play an essential role in the climatic Earth system (Seneviratne et al., 2010).Regular improvement and assessment of generic land surface models (LSMs) are also required.In particular, the seasonal and inter-annual variability of the vegetation interacts with hydrological processes, and must be represented well (Szczypta et al., 2012).Modern LSMs such as Interactions between Soil, Biosphere and Atmosphere, CO 2 -reactive (ISBA-A-gs) (Calvet et al., 1998;Gibelin et al., 2006) or ORganizing Carbon and Hydrology In Dynamic EcosystEms (ORCHIDEE) (Krinner et al., 2005) are able to simulate the diurnal cycle of water and carbon fluxes and, on a daily basis, plant growth and key vegetation variables such as the above-ground biomass (B ag ) and the leaf area index (LAI).In areas affected by droughts, soil moisture has a marked impact on plant growth, and the way root water uptake is represented in such LSMs may influence the simulated B ag and LAI values, in particular the maximum values reached every year.Therefore, long time series of observations related to the latter quantities, such as agricultural yields, have potential in the evaluation of the simulation of the available soil water content (AWC) and root water uptake in LSMs, provided their inter-annual variability is governed by climate and not by trends or changes in agricultural practices.
In Europe, a marked positive trend in crop yields has been observed in the last 45 years, due to the agricultural intensification and evolution of farming practices (Smith et al., 2010a, b).However, Brisson et al. (2010) and Gate et al. (2010) have shown that yields have been stagnating in Europe since the beginning of the 1990s, and particularly since 1996 in France.Therefore, it can be assumed that in the last two decades, the year-to-year change in the large-scale yield of a given rain-fed crop type is mainly driven by the climate variability.In Europe, Smith et al. (2010a, b) showed that the agricultural statistics can be used to assess crop simulations on a country level.On a finer spatial scale over France, Calvet et al. (2012), hereafter referred to as Ca12, have used agricultural statistics (Agreste, 2014) to benchmark several configurations of the ISBA-A-gs LSM through the correlation between yield time series and B ag simulations for the 1994-2008 period.The Agreste data are provided for administrative units (hereafter referred to as "départements").In ISBA-A-gs, the plant phenology is driven by photosynthesis: on a daily basis, plant growth is governed by the accumulation of the hourly net assimilation of CO 2 through the photosynthesis process, and plant mortality is related to a deficit in photosynthesis.The simulated annual maximum B ag and maximum LAI may differ from one year to another in relation to the impact of the weather and climate variability on photosynthesis.In regions where a deficit in precipitation may occur, soil moisture is a key driver of the photosynthesis and plant growth of rain-fed crops and grasslands.Although ISBA-A-gs is not a crop model and agricultural practices are not explicitly represented, Ca12 achieved a good representation of the inter-annual variability of the dry matter yield (DMY) for grasslands over many départements in France.On the other hand, representing the year-to-year variability of the grain yield (GY) of winter/spring cereals was more difficult.By performing a sensitivity study on different parameters of the model, they concluded that the maximum available soil water content (MaxAWC) and the mesophyll conductance under well-watered conditions (g m ) were the two key parameters driving the inter-annual variability of the simulated B ag .In particular, they showed that the model was markedly sensitive to MaxAWC (especially at low MaxAWC values).
In Ca12, an effort was made to benchmark two options of the vegetation model (drought-avoiding vs. droughttolerant).In this study, an effort is made to benchmark several options of the soil hydrology model.The main objective of this study is to assess to what extent using more refined representations of the soil hydrology and the root water uptake can improve the representation of the inter-annual variability of GY (and, possibly, DMY).The ISBA-A-gs model and the method proposed by Ca12 are used to evaluate a new option of the ISBA-A-gs model using a multi-layer soil model permitting a more detailed representation of soil moisture and soil temperature profiles, and of root water uptake.Since several options can be envisaged to implement the multi-layer soil hydrology simulations, a side objective of this study is to benchmark these options and learn about the representation of root water uptake.
The various versions of ISBA-A-gs are presented in Sect.2, together with the annual yield statistics of Agreste.The symbols used in this work are listed and defined in Table A1.The results obtained with the different set of simulations are shown in Sect. 3 and the differences in the interannual variability of the various simulations of B ag are presented, together with the hydrological variables.The results are analyzed and discussed in Sect. 4 and the conclusions of this study are summed up in Sect. 5.

Agricultural statistics in France
Agreste is an annually updated set of agricultural data over France (Agreste, 2014).An inventory of the land use in agriculture and of the crop, forage and livestock production is made on a yearly basis.The data are provided for département administrative units.For crops and grasslands, annual grain yields and dry matter yields (GY and DMY, respectively) are supplied.A new version of Agreste with recalculation since 1989 has recently been published.In this study, the new Agreste data set is used over the 1994-2010 period to examine the inter-annual variability of winter/spring cereal crop GY in 45 départements and of natural grassland DMY in 48 départements (Fig. 1).For cereals, we consider the following six crops: winter wheat, rye, winter barley, spring barley, oat and triticale.For grasslands, the DMY values of permanent grasslands are used.They correspond to natural grasslands or grasslands planted at least 6 years before.Figure 2 shows the inter-annual variability of the average GY and DMY time series derived from Agreste over the considered départements.Over the 1994-2010 period, no significant (p value < 0.01) trend is observed for any of the time series.A few anomalous years affected by particular climate events can be noticed.For example, Fig. 2 shows that the severe summer drought of 2003 impacted both crop and grassland yields.In 2007, the grassland production was the highest of the whole period.Conversely, it was one of the worst in terms of crop yield.The 2007 year was marked by a warm spring (favourable to permanent grasslands), followed by a slightly cold summer (detrimental to cereals).Furthermore, the rains were abundant over the grassland regions considered in this study, and have also contributed to the higher production (Agreste Bilans, 2007;Agreste Conjoncture, 2007;Agreste Infos Rapides, 2007).

The ISBA-A-gs land surface model
The Interactions between Soil, Biosphere, and Atmosphere (ISBA) model (Noilhan and Planton, 1989;Noilhan and Mahfouf, 1996) was designed to describe the daily course of land surface state variables into global and regional climate models, weather forecast models, and hydrological models.In the original version of ISBA, a single root-zone soil layer is considered.A thin top soil layer is represented using the Deardorff (1977Deardorff ( , 1978) ) force-restore approach.Soil characteristics, such as soil water and heat coefficients, the wilting point and the field capacity, depend on soil texture (sand and clay fractions).The stomatal conductance calculation is based on the Jarvis (1976) approach, and accounts for photosynthetically active radiation (PAR), soil water stress, vapour pressure deficit and air temperature.
The representation of the soil physics of the initial version of ISBA was upgraded gradually.A multi-layer soil model including soil freezing processes was developed by Boone et al. (2000) and Decharme et al. (2011).The multi-layer soil model explicitly solves the one-dimensional Fourier law and the mixed form of the Richards equation.The multi-layer representation is used to discretise the total soil profile.In each layer, the temperature and the moisture are computed according to the hydrological and texture layer characteristics.The heat and water transfers are decoupled: heat trans-fer is solely along the thermal gradient, while water transfer is induced by gradients in the total hydraulic potential.Hereafter, the two-layer force-restore model and the diffusion model are referred to as "FR-2L" and "DIF", respectively.
In addition to the simple Jarvis parameterisation of stomatal conductance, Calvet et al. (1998) and Gibelin et al. (2006) have developed ISBA-A-gs.ISBA-A-gs ("A" stands for net assimilation of CO 2 , and "gs" for stomatal conductance) is a CO 2 -responsive version of ISBA able to simulate photosynthesis and its coupling to stomatal conductance.This option was used in studies on the impact of climate change (Calvet et al., 2008;Queguiner et al., 2011) and on the impact of drought on the vegetation in the Mediterranean basin (Szczypta, 2012).
Under well-watered conditions, the A-gs formulation is based on the model proposed by Jacobs et al. (1996) (Calvet et al., 1998(Calvet et al., , 2004;;Gibelin et al., 2006).In this approach, the main parameter driving photosynthesis is g m .Under waterlimited conditions, a soil moisture stress function (F S ) is applied to key parameters of the photosynthesis model.For herbaceous vegetation, two parameters are assumed to respond to soil moisture stress (Calvet, 2000): the mesophyll conductance and the maximum leaf-to-air saturation deficit (D max ).Low (high) values of the latter correspond to high (low) sensitivity of the stomatal aperture to air humidity.These photosynthesis parameters are dependent on F S .Two contrasting responses of the model parameters to soil moisture are represented: drought-avoiding and drought-tolerant (see Fig. S1 in the Supplement).When F S is higher than the critical soil water stress F SC (F SC = 0.3 in our simulations), a drop in F S triggers an increase (decrease) in g m and a decrease (increase) in D max for the drought-avoiding (droughttolerant) parameterisation.The drought-avoiding parameterisation is used for cereal crops, and the drought-tolerant parameterisation is used for grasslands.This assumption was validated by Ca12.The drought response model is illustrated by Fig. S1.These parameters are then used to calculate the hourly leaf-level net assimilation of CO 2 and the stomatal conductance, in relation to sub-daily meteorological inputs such as the incoming solar radiation.A radiative transfer scheme is then used to upscale net assimilation of CO 2 and transpiration at the vegetation level.The plant transpiration flux is used to calculate the soil water budget through the root water uptake.The net assimilation of CO 2 serves as an input to the plant growth model, and LAI and B ag are updated on a daily basis.Figure 3 illustrates these mechanisms.For moderate soil water stress, the drought-avoiding response results in an increase in the water use efficiency (WUE).In the drought-tolerant response, WUE does not change, or decreases.It must be noted that another representation of the response to drought is used for forests (Calvet et al., 2004).
ISBA-A-gs contains a photosynthesis-driven plant growth model able to simulate LAI and the vegetation biomass on a daily basis.For herbaceous vegetation, the model simulates  the above-ground biomass.The B ag variable has two components (active biomass and structural biomass) related by a nitrogen dilution parameterisation (Calvet and Soussana, 2001).The leaf nitrogen concentration N L is a parameter of the model affecting the specific leaf area (SLA), the ratio of LAI to leaf biomass (in m 2 kg −1 ).The SLA depends on N L and on plasticity parameters (Gibelin et al., 2006).This version of ISBA-A-gs, called "NIT", is used in this study.
An assessment of the quality of ISBA-A-gs output variables has been performed in previous local studies with in situ data over France (Rivalland et al., 2005;de Rosnay et al., 2006;Sabater et al., 2007;Brut et al., 2009;Lafont et al., 2012).Gibelin et al. (2006) have shown that the LAI simulated by ISBA-A-gs on a global scale is consistent with satellite-derived LAI products.
Furthermore, a radiative transfer model within the vegetation canopy describes the attenuation of the PAR through a self-shading approach, and photosynthesis is calculated at three levels of the canopy using a three-point Gauss quadrature method (Jacobs, 1994).A new radiative transfer (hereafter referred to as "NRT") scheme was recently implemented in ISBA-A-gs by Carrer et al. (2013).The NRT is more detailed than the original model, and a vertical profile of ten layers within the canopy is represented.Because of the heterogeneity of the different vegetation canopies, distinct bottom and top canopy layer parameterisations are considered.Also, NRT has distinct representations of sunlit and shaded leaves, with two PAR calculations at each layer.Carrer et al. (2013) showed that NRT represents better the gross primary production (GPP) on both local and global scales.

Root density and the soil water stress
In the DIF simulations, the root density profile (Y ) is expressed by the following equation derived from Jackson et al. (1996): where Y (d k ) is the cumulative root fraction (a proportion between 0 and 1) from the soil surface to the bottom of a soil layer within the root zone, at a depth d k (m), d R is the rootzone depth (m) and R e the root extinction coefficients equal to 0.961 and 0.943 for crops and for temperate grasslands, respectively (Jackson et al., 1996).For a given value of d R , the lower value of R e for temperate grasslands corresponds to a cumulative root fraction higher than for crops close to the top soil layer, 15 % higher at d k = 0.36 m and more than 40 % higher at d k < 0.05 m.The cumulative root density is equal to 1 at the bottom of the root-zone soil layer (d R ).
The soil wetness index (SWI) of a bulk top soil layer of thickness d k , where k is the index of the deepest considered individual soil layer, and of a soil layer at depth d i (SWI TOP (d k ) and SWI i , respectively), are defined as where θ i is the volumetric water content (in m 3 m −3 ) at depth d i , d i is the thickness of soil layer i, and the subscripts "FC" and "WILT" indicate soil moisture at field capacity and at wilting point, respectively.Equation ( 2) is used to assess the soil moisture stress in a single soil layer or in several soil layers forming a bulk layer from the surface to a depth d k .Equation ( 3) is used to assess the soil moisture stress of an individual soil layer at depth d i .Equations ( 2) and ( 3) are used to calculate the stress function in FR-2L and DIF simulations, respectively.In this study, the same soil type is used for all the simulations, and a homogeneous soil profile is assumed with sand and clay fractions of 32.0 and 22.8 %, respectively, and θ FC i = θ FC = 0.30 m 3 m −3 and θ WILT i = θ WILT = 0.17 m 3 m −3 .Since the agricultural statistics we use concern rather large administrative units, it would have been illusory to try to use local soil texture properties.
The value of MaxAWC is expressed in units of kg m −2 , and depends on soil and plant characteristics: soil moisture at field capacity, soil moisture at wilting point (θ FC and θ WILT , respectively, in m 3 m −3 ), and rooting depth (d R , in m): where ρ = 1000 kg m −3 is the water density.The θ FC and θ WILT values are common to all the simulations, and the different MaxAWC values are obtained by varying the rootzone depth (d R ).
In the ISBA-A-gs simulations, the dimensionless stress function F S is used to calculate photosynthesis and the plant transpiration flux (F T , in kg m −2 s −1 ).The F S function varies between 0 (at wilting point or below) and 1 (at field capacity or above).Between these two limits, F S equals SWI TOP (d R ) in FR-2L, and plant transpiration is driven by the total soil water content in the root zone.In the case of DIF simulations, F S is the sum of the stress functions of each soil layer in the root zone F S i , i.e.SWI i , balanced by the root fraction R d i at depth d i : , and where N is the number of soil layers in the root zone.Once the F S stress index has been determined, the photosynthesis parameters can be updated, and the leaf-level and vegetationlevel fluxes can be calculated (Fig. 3).The F S value is used to calculate the photosynthesis parameters g m and D max under water-limited conditions (Fig. S1).
The root water uptake in layer i, S T i (in kg m −2 s −1 ), is calculated as (6)

Design of the simulations
In this study, the ISBA-A-gs LSM is used within version 7.2 of the SURFEX (SURFace EXternalisée) Earth surface modelling platform of Météo-France (Masson et al., 2013).For the first time, the NIT biomass option of the model and the NRT light absorption scheme are used together with the DIF multi-layer soil configuration.Two representations of the soil hydrology (FR-2L and DIF options) are considered, for both C3 crops and grasslands.The model simulations are offline (not coupled to the atmosphere) and driven by a meteorological reanalysis.We consider the vegetation cover fraction to be equal to 1 across seasons.We use the ISBA-A-gs default avoiding (tolerant) response to the drought for C3 crops (grasslands).Standard values of the model parameters used in this study are summarised in Table 1.
Six experiments are performed: -FR-2L is based on the force-restore representation of the soil hydrology, and is similar to the model configuration used by Ca12.The root zone corresponds to the whole soil layer.-DIF1 uses the new DIF capability of SURFEX v7.2 (Fig. 4).As in FR-2L, the root zone corresponds to the whole soil layer.The root profile reaches the bottom of the soil layer, and the total soil depth corresponds to d R .
-DIF2 includes additional subroot-zone base flow soil layers with respect to DIF1, and the deep soil layers contribute to plant transpiration through capillary rises.
It is assumed that MaxAWC is governed by the limited capacity of the plants to develop a root system in a deep soil, and the number of subroot-zone layers decreases when the rooting depth increases.A constant total soil depth of 1.96 m is prescribed, and d R is varied between 0.36 and 1.76 m (Fig. 5).
-DIF3 is similar to DIF1, as soil depth is the main limitation to root water extraction.However, two additional base flow soil layers contribute to transpiration through capillary rises.The total soil depth and d R are varied simultaneously, and two adjacent 0.1 m thick deep soil layers are represented (Fig. 6).
-DIF1-NRT permits assessment of the impact of a refined representation of the CO 2 uptake by the vegetation on the B ag inter-annual variability, as the NRT light absorption option is used together with DIF1.
-DIF1-Uniform permits assessment of the sensitivity of the ISBA-A-gs simulations to the shape of the root density profile.It corresponds to DIF1 simulations using a uniform root density profile instead of Eq. (1).These simulations are performed over the 61-Orne département (see Sect. 4.1).

Atmospheric forcing
The atmospheric forcing data required for our simulations are provided by the SAFRAN ("Système d'Analyse Fournissant  des Renseignements Atmosphériques à la Neige") mesoscale atmospheric analysis system (Durand et al., 1993(Durand et al., , 1999)).Precipitation, air temperature, air humidity, wind speed, incoming solar radiation and incoming infra-red radiation are provided over France at an 8 km × 8 km spatial resolution on an hourly basis.The SAFRAN product was evaluated by Quintana-Seguí et al. ( 2008) using independent in situ observations.One-dimensional model simulations are performed at the 8 km × 8 km spatial resolution of SAFRAN in grid cells corresponding to cereal and natural grassland départements (Fig. 1).These grid cells correspond to plots located within a département and with at least 45 % of their surface covered by either grasslands or crops, according to the average plant functional type coverage given by the 1 km × 1 km ECOCLIMAP-II global database (Faroux et al., 2013).

Optimisation of two key parameters
In this study, the method proposed by Ca12 is used: the values of two key parameters of the ISBA-A-gs simulations, MaxAWC and g m , are explored, and the parameter pair providing the best correlation coefficient (r) of the maximum annual value of B ag (B ag X ) and GY (DMY) is selected, for C3 crops (grasslands).For the FR-2L experiment, the optimisation of both MaxAWC and g m is performed for all the départements of Fig. 1.For the DIF1, DIF2, and DIF3 experiments, only MaxAWC is optimised, and the g m values derived from the FR-2L optimisation are used.In the case of crops, simulated B ag values after 31 July are not considered, in order to be consistent with the theoretical averaged harvest dates in France.Attempts were made to use other dates in July (not shown) without affecting the results of the analysis.
On the other hand, new optimal g m values are obtained together with MaxAWC for the DIF1-NRT experiment, as the representation of photosynthesis at canopy level differs from that of the other experiments.Moreover, major differences with Ca12 are that (1) a longer period is considered (1994( -2010( instead of 1994( -2008 in Ca12) in Ca12), and that (2) a more detailed screening of MaxAWC values is performed (12 values are considered, against 8 values in Ca12).For all the experiments, MaxAWC ranges from 50 to 225 mm, with a lower increment between the small values (50, 62.5, 75, 87.5, 100, 112.5, 125, 137.5, 150, 175, 200 and 225 mm; 12 in total).
For the g m parameter, the same range of values as in Ca12 is used (from 0.50 to 1.75 mm s −1 , six in total).For the three simulations DIF1, DIF2 and DIF3, the same values of optimal g m obtained for each département and vegetation type with the FR-2L version are used.

Metrics used to quantify the inter-annual variability
In Sect.4, the following metrics are used: the annual coefficient of variation (ACV), computed as the ratio of the standard deviation (σ ) of the simulated B ag X to the long-term and the scaled anomaly (A S ) of B ag X of a given year (yr) This metric is also called the z score, and can be applied to the Agreste cereal GY, For C3 crops (Fig. 7), B ag X values for FR-2L tend to reach slightly higher values than for DIF1.The largest difference is observed in 1996.Furthermore, some differences occur in the senescence period, especially in 2001 and 2009.Conversely, the simulated AWC values are higher for DIF1, especially in winter.For both simulations, the wintertime AWC is often higher than MaxAWC (set to 200 mm), in relation to water accumulation above field capacity, under wet conditions.This phenomenon is more pronounced for DIF1 than for FR-2L.Crop re-growth is simulated by both FR-2L and DIF1 during years with a marked summer drought, in 1995, 1996, 1998, 2006and 2010. During wet years (i.e. in 1994, 2000and 2007)), the two experiments provide similar AWC values in summertime.
For grasslands (Fig. 8), the two B ag simulations are also very close.However, contrary to C3 crops, the B ag values of the FR-2L simulation tend to be slightly lower than the DIF1 ones (e.g., in 1997, 2002, 2007, and 2009).The other difference with C3 crops is the systematic occurrence of regrowths.

ISBA-A-gs simulations vs. Agreste observations
The départements where FR-2L B ag X simulations present significant (p value < 0.01) correlations with the Agreste GY    Non-significant, significant at the 1 % level, and significant at the 0.1 % level correlations are indicated in red squares, yellow dots and black dots, respectively.and DMY time series are presented in Fig. 9, and the retrieved g m and MaxAWC median values are presented in Table 2 for all the experiments, together with the number of départements presenting significant correlations with Agreste, for C3 crops and grasslands.With FR-2L, 12 (5) départements present significant positive correlations at the 1 % (0.1 %) level for C3 crops.For grasslands, 34 ( 22) départements present significant positive correlations at the 1 % (0.1 %) level.Although the considered period is longer than in Ca12 (17 years instead of 15 years), these results are similar to those presented in Ca12, even if slight differences can be noticed, such as the number of départements with a significant correlation.In DIF simulations for C3 crops, DIF1 and DIF3 perform nearly as well as FR-2L, and they outperform DIF2: 10 (3) départements present significant positive correlations at the 1 % (0.1 %) level for both DIF1 and DIF3, against 6 (2) for DIF2.For the grasslands, a larger proportion of départements (among 48) presents significant correlations, from 27 (10) départements for DIF2 to 36 (20) for DIF1.The addition of deep soil layers below the root zone tends to degrade the results, especially in DIF2.Finally, the DIF1-NRT simulations perform as well as FR-2L or better, with 13 (4) and 37 (19) départements presenting significant positive correlations at the 1 % (0.1 %) level for C3 crops and grasslands, respectively.
Selecting the départements where the optimisation is successful, i.e.where the correlation between B ag X and GY or DMY is significant (p value < 0.01), the time series of the mean B ag X and mean GY and of the mean B ag X and mean DMY are compared in Fig. 10 for both the FR-2L and DIF1-NRT experiments.The inter-annual variability of the grassland DMY is better represented by B ag X than for the cereal GY, with R 2 = 0.83 and R 2 = 0.45, respectively.The FR-2L experiment presents slightly better R 2 values than DIF1-NRT.For C3 crops, it appears that the two experiments are not able to represent the lower GY in 2007, or the higher GY in 2004.For grasslands, the two experiments are not able to represent the lower DMY in 1996.

Optimal MaxAWC values
Table 2 shows that, for C3 crops, the median MaxAWC value is higher for FR-2L than for DIF1 (125.0 and 112.5 mm, respectively).For DIF2 and DIF3, the median MaxAWC is even lower (81.3 and 93.8 mm, respectively).For grasslands, the median MaxAWC is less variable from one experiment to another (from 68.8 to 81.3 mm).In Table 2, the median Max-AWC values are calculated irrespective of which Agreste cereal GY values are used to derive MaxAWC.Among the 10 départements with DIF1 simulations presenting significant correlations at the 1 % level with Agreste, 8 départements share the same cereal Agreste yields as FR-2L.
These eight départements are listed in Table 3 together with squared correlation coefficient (R 2 ) values and the Max-AWC for FR-2L and DIF1.The FR-2L R 2 is higher than the DIF1 R 2 , except for 08-Ardennes and 63-Puy-de-Dôme.Again, the median MaxAWC is higher for FR-2L than for DIF1 (118.8 and 112.5 mm, respectively).The FR-2L Max-AWC value is lower than the DIF1 MaxAWC value only once, for the 61-Orne département.This indicates that the  DIF1 root density profile tends to increase the impact of drought on plant growth for this département.Also, the largest difference in R 2 between FR-2L and DIF1 is observed for this département.

Plant growth
Table 2 shows that in DIF2 simulations the number of départements with a significant correlation at the 1 % level is lower than in other experiments.The use of DIF2 has a detri-mental impact on the representation of the inter-annual variability by the plant growth model.Figure 11 shows the impact of the root water uptake model on the simulated C3 crop B ag and root-zone soil moisture for the 08-Ardennes département during the growing season, from April to July 1996.
In the FR-2L, DIF1, DIF2, and DIF3 simulations shown in Fig. 11, the same g m = 0.5 mm s −1 and MaxAWC = 75 mm values are used.The growth period is longer in the DIF2 simulation than in the other ones, with senescence only starting during the second half of July.At the same time, the DIF2 root-zone soil moisture presents the highest values.It appears that in the DIF2 simulation, the additional water supplied by capillary rises from the subroot-zone soil layers has a marked impact on the phenology, with the date of maximum B ag shifted to the end of July and a much higher B ag X value than in the other experiments (1.02 kg m −2 for DIF2, against 0.62, 0.58, and 0.72 kg m −2 for FR-2L, DIF1, and respectively).The same phenomenon happens in the DIF3 simulation to a lesser extent.In particular, the DIF3 B ag X is not very different from the FR-2L one.The DIF1 simulation is closer to FR-2L.When the root-zone soil moisture reaches the wilting point (equal to 0.17 m 3 m −3 , as indicated in Fig. 11 by the dashed line), the senescence starts.A marked water stress occurs and impacts photosynthesis and biomass production.Since water is supplied by the subroot-zone soil layers of DIF2 and DIF3, the wilting point is reached later than for FR-2L and DIF1, and the senescence starts later.
In FR-2L, the growth of B ag is faster than in the other simulations.This leads to a slightly higher value of B ag X than for DIF1.This is related to the lower FR-2L root-zone soil moisture in May.In the drought-avoiding C3 crop parameterisation of ISBA-A-gs, a moderate soil moisture stress triggers an increase in water use efficiency (Calvet, 2000) and enhances plant growth.

Discussion
4.1 Are the Jackson root profile model (Eq. 1) and the resulting water availability (Eq.5) applicable on a regional scale?
In the DIF simulations, the stress function depends on the distribution of root density through Eqs. ( 5)-( 6).This allows the lower layers to sustain the transpiration rate to some extent when the upper soil layers dry out.However, one may emphasise that the approach used in this study to simulate the root water uptake is relatively simple, and may not be relevant to representing what really happens on a regional scale.Higher-level models are able to simulate the root network architecture and the three-dimensional soil water flow (Schneider et al., 2010;Jarvis, 2011).Also, the hydraulic redistribution of water from wetter to drier soil layers by the root system (hydraulic lift) is not simulated in this study.Siquiera et al. (2008) have investigated the impact of hydraulic lift using a detailed numerical model, and showed that this effect could be significant.
Another difficulty in the implementation of DIF simulations is that the proposed R e values in Eq. ( 1) are the result of a meta-analysis.A single R e value is proposed for a given vegetation type, while large variability of R e can be observed.This is particularly true for crops, and Fig. 1 in Jackson et al. (1996) shows that Y (d k ) and R e present much higher variability for crops than for temperate grasslands.This difficulty may explain the shortcomings of DIF1 simulations for the 61-Orne département described in Sect.3.2.1 (Table 3).In particular, the root density in the top soil layers has a large impact on the water stress modelling.This is demonstrated by performing an additional DIF1 simulation (DIF1-Uniform) using a uniform root density profile instead of Eq. (1). Figure 12 shows the evolution of B ag , SWI TOP (d R ) and SWI TOP (0.46 m) for the FR-2L, DIF1 and DIF1-Uniform simulations for the 61-Orne département over the period from April to July 1999.For all the simulations, g m = 1.75 mm s −1 and MaxAWC = 225 mm.The B ag evolution during the first 3 months is similar in the three simulations, with slightly faster growth for FR-2L.However, while senescence occurs in mid-July for DIF1, it occurs only at the end of July for FR-2L and DIF1-Uniform.Using the Jackson root density profile in Eq. ( 5) rather than a uniform profile has a marked impact on the simulated water balance.In situations where the top soil layers are drier (wetter) than deep soil layers (i.e.present lower (higher) F S i values), the total F S value is lower (higher) in DIF1 simulations than in FR-2L or DIF1-Uniform simulations.This tends to trigger an earlier senescence in DIF1 simulations.The early senescence for DIF1 is related to values of SWI TOP getting close to zero at the top fraction of the root zone: while SWI TOP (0.46 m) decreases below the 0.3 critical soil water stress value (Table 1) at the beginning of July, for DIF1, it never gets below 0.3 in July for DIF1-Uniform.It must be noted that Fig. 12 shows that root water uptake is reduced earlier with FR-2L than with DIF1, in relation to faster plant growth in the FR-2L simula- tion.For C3 crops, a drought-avoiding response to soil water stress is simulated, triggering an increase in WUE (and in the plant growth rate) as soon as θ is less than θ FC .the DIF1 simulations tend to accumulate water above the field capacity (i.e.θ remains longer above θ FC than for FR-2L), the increase in WUE tends to occur later than for FR-2L.Finally, the B ag X value for FR-2L and DIF1-Uniform is higher than for DIF1.This root profile effect also has an impact on the inter-annual variability and partly explains the lower R 2 value for DIF1 in Table 3 for this département.
Figure 12 shows that situations in which the top soil layers are drier than deep soil layers tend to be more frequent in DIF1 simulations than in DIF1-Uniform simulations, in relation to the enhanced root water uptake close to the soil surface.Therefore, for given MaxAWC and soil wetness conditions, the total F S values tend to be lower in DIF1 simulations than in DIF1-Uniform (and FR-2L) simulations.This results in less evapotranspiration and less GPP.The lower GPP in DIF simulations results in lower B ag X values, especially for cereals as illustrated in Fig. 10.As noted by Feddes et al. (2001), the limitation of transpiration in DIF simulations when a great deal of water is still available at depth is probably too severe.In the real world, plants are able to transfer water uptake to compensate for water stress in the top layers, and DIF simulations cannot adequately account for it.This fact probably explains part of why this model is not able to outperform the FR-2L simulations.

Do changes in the representation of photosynthesis have an impact on the model performance?
In this section, the impact of the revised vegetation radiative transfer scheme and the refreshed g m parameter (DIF1-NRT experiment) is discussed.Table 2 shows that while the DIF1-NRT results are close to those of DIF1 for grasslands, DIF1-NRT tends to outperform DIF1 for C3 crops.Figure 13 presents the simulated B ag of C3 crops and grasslands for the DIF1 and DIF1-NRT simulations in the 61-Orne département over the 1994-2010 period.The two grassland simulations are very similar.On the other hand, the two C3 crop simulations differ in B ag X values.The mean simulated B ag X values for C3 crops are 1.61 and 1.32 kg m −2 for DIF1 and DIF1-NRT, respectively.The lower B ag X values simulated by DIF1-NRT are related to the lowest gross primary production simulated by this version of the ISBA-A-gs model (Carrer et al., 2013).Also, DIF1-NRT simulates shorter growing periods and a slightly enhanced inter-annual variability: the ACV (see Sect.   2002, 2008, 20091997, 201020042001, 2007DIF1-NRT 2002, 20092008200119971998, 20042003, 2007Grasslands FR-2L 2007, 200820002003, 20101996DIF1-NRT 2000, 2007, 20082003, 20101996 tively, and ACV values for DIF1 and DIF1-NRT are both equal to 30 %.

Can the ISBA-A-gs model predict the relative gain or loss of agricultural production during extreme years?
ISBA-A-gs is not a crop model, and does not predict yield per se.The background assumption of this work is that the regional-scale above-ground biomass simulated by a generic LSM can be used as a proxy for GY or DMY in terms of inter-annual variability.The quantitative consistency between the simulated biomass and the agricultural statistics was discussed extensively by Ca12 (Sect. 3.3 and Figs. 12 and 13 in Ca12).For cereals, they considered the ratio of crop yield to the maximum above-ground biomass, called the harvest index.The later ranged from 20 to 50 %, and this was consistent with typical harvest index values given by Bondeau et al. (2007) for temperate cereals.The same result is obtained in this study (not shown).For grasslands, Ca12 simulated both managed and unmanaged grasslands.For managed grasslands, DMY was explicitly simulated, and ranged from 0.1 to 0.8 kg m −2 .The scatter of the simulated DMY was relatively small, with a standard deviation of differences with the Agreste DMY of 0.20 kg m −2 .ISBA-A-gs tended slightly to underestimate DMY values, with a mean bias of −0.08 kg m −2 .For unmanaged grasslands, the simulated B ag was 0.17 kg m −2 higher than the Agreste DMY values, on average.In this study, unmanaged grasslands were considered, only, and results similar to those of Ca12 were found (not shown).
The ISBA-A-gs model is optimised to maximise the correlation coefficient between Agreste GY (or DMY) and modelled B ag X .The resulting scores are used to assess the capability of a given model configuration to represent the interannual variability of B ag X over the 1994-2010 period.In studies where the objective of the model calibration is to improve the model prediction for operational applications, the model quality needs to be confirmed in an independent run with data not used during the calibration.An example of a rigorous calibration and validation procedure in hydrology can be found in Refsgaard (1997).In this study, a validation run was not performed, as the considered period was too short to apply a split-sample procedure and separate calibration and validation sub-periods.Moreover, the objective of this study is to benchmark DIF options, not to predict the agricultural yields.Therefore, using an independent data set to assess yield prediction is not needed.
While the main objective of this work is to evaluate contrasting root water uptake models using agricultural statistics, one can investigate how the resulting B ag X values react to extreme years (either favourable or unfavourable to agricultural production).The best simulations result from the optimisation of the MaxAWC parameter.Table 4 summarises the true and false detection of favourable and unfavourable years.The latter are defined as A S,B ag X or A S,DMY values higher (lower) than 1.0 (−1.0).The A S,B ag X or A S,DMY values are based on the mean time series of Fig. 10.The undetected favourable and unfavourable years are also listed in Table 4.The best detection performance is obtained by DIF1-NRT for grasslands, with only 1996 not detected as unfavourable.The worst detection performance is obtained by DIF1-NRT for C3 crops, with 2003 and 2007 not detected as unfavourable, 1998 and 2004 not detected as favourable, 1997 wrongly detected as unfavourable, and 2008 wrongly detected as favourable.For grasslands, the extreme years, defined as A S,DMY values higher (lower) than 1.5 (−1.5), are 2007 (favourable) and 2003 (unfavourable).These two cases are correctly identified in the two experiments.For C3 crops, the most favourable years are 2002 and 2009, and the most unfavourable year is 2007.While 2002 and 2009 are correctly identified in the two experiments, 2007 is not detected.The higher performance in the representation of extreme years for grasslands than for C3 crops is consistent with the results of Table 2 showing that significant correlations between B ag X and DMY are obtained more often than between B ag X and GY.This can be explained by the more pronounced inter-annual variability of the grassland DMY, with ACV = 30 % against ACV values less than 10 % for the cereal GY.The highest sensitivity of grasslands to climatic conditions is related to their growing cycle covering a longer period than cereals, and to their MaxAWC values, generally lower than for cereals (Table 2).Finally, ISBA-A-gs has no direct representation of agricultural practices and the cereal GY, and the consistency between B ag X and GY relies on the hypothesis that the harvest index (the ratio of GY to B ag X ) does not vary much from one year to another on the considered spatial scale.This issue is discussed in Ca12.For grasslands, the simulated B ag X is more directly representative of DMY.This explains why a better agreement of the simulations is found with the grassland DMY than with the cereal GY (Tables 2 and 4).

Prospects for better constraining MaxAWC
Ca12 have shown that MaxAWC is the main driver of the inter-annual variability of B ag in the ISBA-A-gs model.Representing the year-to-year B ag variability in a dynamic vegetation model is a prerequisite for correctly representing surface fluxes on all temporal scales (from hourly to decadal).Table 2 shows that significant differences in the representation of the B ag inter-annual variability are triggered by switching from one model option to another.Also, for a given model option, the median g m and MaxAWC values obtained for cereals contrast from those obtained for grasslands.This is very valuable information for guiding the mapping of the model parameters in future studies.It must be noted that using the inter-annual variability of plant growth to assess LSM parameters is a rather new idea.For example, Rosero et al. (2010) and Gayler et al. ( 2014) performed an assessment of key parameters of the Noah LSM, including a version with a dynamic vegetation module, using a set of experimental stations.However, they did not address the inter-annual variability of plant growth, as their simulations covered one vegetation cycle only.Such a short simulation period is not sufficient for constraining those model parameters which affect the inter-annual variability of plant growth (Kuppel et al., 2012).
In addition to the intrinsic limitations related to the use of a generic LSM unable to represent agricultural practices (see above), uncertainties are generated by the data sets used to force the LSM simulations.For example, the incoming radiation in the SAFRAN atmospheric analysis can be affected by seasonal biases (Szczypta et al., 2011;Carrer et al., 2012).Since phenology in ISBA-A-gs is driven by photosynthesis, biases in the incoming radiation can impact the date of the leaf onset.The impact of errors in the forcing data is probably more acute for cereals than for grasslands in relation to a shorter growing period.More research is needed to assess the impact of using enhanced atmospheric reanalyses (Weedon et al., 2011;Oubeidillah et al., 2014) and proxies for annual agricultural statistics such as gridded maximum LAI values at a spatial resolution of 1 km × 1 km derived from satellite products (Baret et al., 2013).
Another difficulty is that the coarse spatial resolution of agricultural statistics prevents the use of local soil properties (Sect.2.3).Models need to be tested on a local scale using data from instrumented sites.For example, the DIF version of ISBA was tested on a local scale by Decharme et al. (2011) over a grassland site in south-western France.However, the soil and vegetation characteristics at a given site may differ sharply from those at neighbouring sites.It is important to explore new ways of assessing and benchmarking model simulations on a regional scale.Remote-sensing products can be used to monitor terrestrial variables over large areas and to benchmark land surface models (Szczypta et al., 2014).At the same time, using in situ observations as much as possible is key, as remote-sensing products are affected by uncertainties.So far, the French annual agricultural yield data have been publicly available on a département scale only.In order to take advantage of the existing information on soil properties, an option could be to use satellite-derived LAI products at a spatial resolution of 1 km × 1 km in conjunction with soil maps at the same spatial resolution (e.g.derived from the Harmonized World Soil Database, Nachtergaele et al., 2012).Since these products are now available on a global scale, the methodology explored in this study over metropolitan France could be extended to other regions.
The ISBA-A-gs model is intended to bridge the gap between the terrestrial carbon cycle and the hydrological simulations (e.g.river discharge).In previous works, the ISBA-A-gs model was coupled to hydrological models able to simulate river discharge (e.g.Queguiner et al., 2011;Szczypta et al., 2012).While simulating vegetation requires a good description of the soil water stress, hydrological simulations are sensitive to changes in the representation of the surface water and energy fluxes.The latter are controlled to a large extent by vegetation.As suggested by Feddes et al. (2001) and Decharme et al. (2013), the obtained "effective root distribution function" could be validated using river discharge observations by coupling the LSM to a hydrological model.We will investigate this possibility in future work.Note however that the river discharge is often impacted by anthropogenic effects such as dams and irrigation.Such effects are not completely represented in large-scale hydrological models (Hanasaki et al., 2006).

Conclusions
The observed cereal GY and permanent grassland DMY production in France from 1994 to 2010 was used in this study to evaluate four contrasting representations of the root water uptake in the ISBA-A-gs land surface model within SUR-FEX.A simple representation of the root-zone soil moisture based on a single bulk reservoir (FR-2L) was compared with multi-layer diffusion models describing the soil water uptake profile.The latter used the Jackson root vertical distribution equation, with and without additional subroot-zone base flow soil layers.In order to limit the uncertainty related to the lack of knowledge of local rooting depth conditions, the Max-AWC quantity was retrieved by matching the simulated B ag X Hydrol.Earth Syst.Sci., 18, 4979-4999, 2014 www.hydrol-earth-syst-sci.net/18/4979/2014/ with the Agreste agricultural statistics, for given vegetation and photosynthesis parameters.The impact on the results of the representation of the vegetation was assessed using another representation of the light absorption by the canopy and using refreshed values of the g m photosynthesis parameter.The B ag X time series based on the multi-layer model without additional subroot-zone base flow soil layers presented correlations with the agricultural statistics similar to those obtained with FR-2L.On the other hand, adding subroot-zone base flow soil layers tended to degrade the correlations.Overall, a better agreement of the simulations was found with the grassland DMY than with the cereal GY in relation to several factors, such as (1) the more pronounced inter-annual variability of the grassland DMY, (2) the more direct correspondence between B ag X and DMY, and (3) less variability in the parameters of the Jackson model than for crops.More research is needed to map the MaxAWC parameter.In particular, long time series of satellite-derived vegetation products (e.g.GEOV1, Baret et al., 2013) could be used in conjunction with soil parameter maps to constrain MaxAWC.The next steps are (1) to verify that the new model parameters have a positive impact on the water and carbon fluxes derived from in situ flux-tower observations and satellite products, on a regional scale and on various timescales (hourly to decadal), and (2) to use a hydrology model coupled to SUR-FEX (Szczypta et al., 2012) to assess the impact of the new MawAWC maps on river discharge.
The Supplement related to this article is available online at doi:10.5194/hess-11-4979-2014-supplement.

Figure 1 .
Figure 1.Location of the 45 cropland and 48 grassland 8 km × 8 km grid cells (blue and green dots, respectively) and the corresponding département number.

Figure 2 .
Figure 2. Averaged annual statistics of Agreste over the 1994-2010 period of (top panel) grain yields of six cereals (winter wheat in black, rye in red, winter barley in blue, spring barley in green, oat in orange and triticale in purple) over the 45 départements of Fig. 1 and (bottom panel) dry matter yields of permanent grasslands over the 48 départements of Fig. 1.

Figure 3 .
Figure 3. Relation of biogeophysical variables to leaf-scale and vegetation-scale fluxes in the ISBA-A-gs simulations.

Figure 4 .
Figure 4. Soil profile of the DIF1 experiment.The soil depth within the root zone is in metres.Only two configurations are represented: for the minimum (left panel) and maximum (right panel) values of MaxAWC (50 and 225 mm, respectively).The cumulative root density profile for crops (Eq. 1 with R e = 0.961) is represented by a brown line.A top soil layer of 1 cm is represented.

Figure 5 .
Figure 5.As in Fig. 4, except for the DIF2 experiment.Subroot soil layers are added (blue lines) down to a constant soil depth of 1.96 m.

Figure 6 .
Figure 6.As in Fig. 4, except for the DIF3 experiment.Two subroot soil layers of 10 cm are added (blue lines).

Figures 7
Figures7 and 8show an example of the inter-annual variability of the simulated B ag and AWC (in kg m −2 ) as simulated by FR-2L and DIF1 for C3 crops and grasslands of the 61-Orne département.The optimal parameter values for C3 crops and grasslands are 1.75 and 0.5 mm s −1 for g m , and 200 and 50 mm for MaxAWC, respectively.For C3 crops (Fig.7), B ag X values for FR-2L tend to reach slightly higher values than for DIF1.The largest difference is observed in 1996.Furthermore, some differences occur in the senescence period, especially in 2001 and 2009.Conversely, the simulated AWC values are higher for DIF1, especially in winter.For both simulations, the wintertime AWC is often higher than MaxAWC (set to 200 mm), in relation to water accumulation above field capacity, under wet conditions.This phenomenon is more pronounced for DIF1 than for FR-2L.Crop re-growth is simulated by both FR-2L and DIF1 during years with a marked summer drought, in1995 , 1996, 1998 , 2006  and 2010 .During wet years (i.e. in 1994, 2000   and 2007)), the two experiments provide similar AWC values in summertime.For grasslands (Fig.8), the two B ag simulations are also very close.However, contrary to C3 crops, the B ag values of the FR-2L simulation tend to be slightly lower than the DIF1 ones(e.g., in 1997, 2002, 2007, and 2009).The other difference with C3 crops is the systematic occurrence of regrowths.

Figure 7 .
Figure 7. Simulations over the 1994-2010 period for C3 crops (g m = 1.75 mm s −1 , MaxAWC = 200 mm) in the 61-Orne département of (top panel) the above-ground biomass and (bottom panel) the available water content in the root zone, using the FR-2L and DIF1 configurations (black and red lines, respectively).

Figure 9 .
Figure 9. Best FR-2L simulation vs. Agreste statistical correlation levels obtained for (left panel) C3 crops and (right panel) grasslands.Non-significant, significant at the 1 % level, and significant at the 0.1 % level correlations are indicated in red squares, yellow dots and black dots, respectively.

Figure 10 .
Figure 10.Averaged simulated yearly B ag X values (ISBA-A-gs, solid lines) and averaged observed agricultural yields (Agreste, dashed lines) for départements with significant correlations (R 2 ) at the 1 % level, with both FR-2L (black solid line) and DIF1-NRT (red solid line) simulations for (top panel) C3 crop GY and (bottom panel) grassland DMY.

Figure 11 .
Figure 11.Simulations in 1996 for C3 crops (g m = 0.5 mm s −1 , MaxAWC = 75 mm) in the 08-Ardennes département of (top panel) aboveground biomass and (bottom panel) root-zone soil moisture in the DIF1, DIF2, DIF3 and FR-2L configurations (red solid, red dotted, red dashed, and black lines, respectively).The grey lines indicate the root-zone soil moisture values at field capacity and at wilting point.
2.7) is equal to 7.4 % for DIF1, and to 8.4 % for DIF1-NRT.For grasslands, the mean simulated B ag X values are 0.46 and 0.44 kg m −2 for DIF1 and DIF1-NRT, respec-

Table 2 .
Median MaxAWC value and median g m value under well-watered conditions, derived for each experiment ( 1 ) and number of départements where the simulated B ag X presents significant correlations ( 2 ) with the annual yields of Agreste statistics for six cereals (winter wheat, rye, winter barley, spring barley, oat and triticale) and for permanent grasslands in France over the 1994-2010 period.Median values are in bold.

Table 3 .
Optimal MaxAWC and squared correlation coefficient (R 2 ) between B ag X and Agreste for FR-2L and DIF1 simulations at départements where the same cereal Agreste data are used and where the correlation between B ag X values and the yields of Agreste statistics are significant at least at the 1 % level.The highest MaxAWC and R 2 values at a given département are in bold.

Table 4 .
Correspondence between simulated and observed extreme years for départements with significant correlations (R 2 ) at the 1 % level with both the FR-2L and DIF1-NRT simulations for C3 crops and grasslands as shown in Fig.10.Favourable (unfavourable) years are defined as z scores A S,B ag X or A S,DMY higher (lower) than 1.0 (−1.0).Years with A S,DMY higher (lower) than 1.5 (−1.5) are in bold.