A simple ecohydrological model captures essentials of seasonal leaf dynamics in semi-arid tropical grasslands

A simple ecohydrological model captures essentials of seasonal leaf dynamics in semi-arid tropical grasslands P. Choler, W. Sea, P. Briggs, M. Raupach, and R. Leuning CSIRO Marine and Atmospheric Research, P.O. Box 3023, Canberra, ACT, 2601, Australia Laboratoire d’Ecologie Alpine, UMR CNRS-UJF 5553, Université J. Fourier, Grenoble I, BP53, 38041 Grenoble, France Received: 27 July 2009 – Accepted: 14 August 2009 – Published: 2 September 2009 Correspondence to: P. Choler (philippe.choler@ujf-grenoble.fr) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
Most recent advances in empirical models of leaf phenology, i.e. the time dependence of Leaf Area Index (LAI), have been made for temperate deciduous forests where temperature is the main controlling factor (Chuine, 2000).By contrast, leaf phenology in water-limited ecosystems remains poorly captured by current broad empirical approaches (Botta et al., 2000).Yet, water is the main controlling factor of ecosystem functioning for more than 50% of the land mass (Churkina and Running, 1998;Huxman et al., 2004).Models of leaf phenology are an essential component of land surface models (review in Pitman, 2003).Leaf phenology not only determines the favourable period for carbon uptake, but also modifies important surface parameters such as albedo and roughness, which feedback on exchanges of mass and energy between land and atmosphere (Pielke et al., 1998;Dickinson et al., 1998;Bounoua et al., 2006).
Two main reasons explain why capturing leaf dynamics in water-limited ecosystems.remains a challenge.First, soil water balance, which has long been recognized as a key driver of plant growth in water-limited ecosystems (Walker and Langridge, 1996;Farrar et al., 1994), is highly variable at the landscape scale and hence difficult to predict in global models.Spatial variability in precipitation and large-scale modelling of runoff and drainage are among difficulties encountered (Teuling and Troch, 2005).By contrast, mean air temperature exhibits smoother variations along broad latitudinal or elevation gradients.Second, there are immediate and important feedbacks of plant growth on soil moisture content through water extraction by roots during transpiration.Contrary to temperature or heat, water is a depletable resource exploited by plants at a rate dependent on both resource (soil water availability) and leaf biomass (Raupach, P. Choler et al.: Seasonal leaf dynamics in tropical grasslands 2007).The coupled dynamics between soil moisture and vegetation growth exhibit strong nonlinearities as highlighted in several conceptual papers in the emerging field of ecohydrology (Rodriguez-Iturbe et al., 2001;Daly and Porporato, 2005).So far, however, there has been no attempt to evaluate the performance of these nonlinear approaches using long data records of climate and remotely-sensed vegetation greenness.
Long time series of remotely-sensed data of phenology provide opportunities to relate vegetation cover or greeness with climatic controls (Pettorelli et al., 2005;Daly and Porporato, 2005).A large number of regional (Ji and Peters, 2004;Farrar et al., 1994;Richard and Poccard, 1998) and global (Potter and Brooks, 1998;Lotsch et al., 2003) studies have examined the response of greenness to rainfall or estimates of soil moisture.Most of these studies used linear modelling, with model performance usually maximised by smoothing or lagging the predictors, usually rainfall time series.Though these modelling approaches have proven to be efficient, they do not make explicit the processes linking precipitation, soil moisture dynamics, and water use by plants to sustain growth.
The motivation of the present study is to examine whether low-dimensional, nonlinear models outperform their linear counterparts for predicting leaf phenology in water-limited ecosystems.Our philosophy for nonlinear modelling is to capture plant-soil moisture coupled dynamics, while keeping the model tractable and comparable in complexity with a linear model that typically uses 3-5 parameters.This constraint in the number of parameters obviously leads to simplifications for soil water balance and plant growth.However, our purpose is to seek a simple, transparent and portable model that can be used with readily available large-scale climatic and remotely-sensed data and that will enable its inclusion in global land surface models.
Numerous studies have shown that the Normalized Difference Vegetation Index (NDVI) can be used as a proxy for green vegetation fraction (Myneni et al., 1995;Carlson and Ripley, 1997), especially for grassland ecosystems where the maximum Leaf Area Index (LAI) does not exceed 3-4 (Choudhury et al., 1994;Carlson et al., 1995;Carlson and Ripley, 1997).In this paper, NDVI calculated from reflectances in the red and near infrared wavebands from the MODerate Imaging Spectroradiometer (MODIS) was used to characterize the phenology of the Australian semi-arid Mitchell grasslands.The use of the Mitchell grasslands as a case study was motivated by three factors.First, these grasslands are known to exhibit a rapid response to precipitation (Christie, 1981).Second, phenology is driven by growth of one life-form, i.e. perennial grasses, which removes the problem of differentiating tree and grass contributions to the NDVI signature (Archibald and Scholes, 2007).Third, these grasslands are distributed along a large precipitation gradient in northern Australia, offering the opportunity to examine the ability of a single set of parameters to predict leaf dynamics under contrasting precipitation regimes.Our study has three steps: (i) calibrating linear and nonlinear models using time series of climate and NDVI-derived vegetation cover for 400 randomly chosen sites within the Mitchell grasslands; (ii) examining the ability of these models to capture the dynamics of test data sets of vegetation cover; and finally (iii) comparing the performance of these models to predict yearly integrated values of vegetation cover.

Study area
The study area of approximately 400 000 km 2 corresponds to the Mitchell Grass Downs ecoregion (Bailey, 2004), extending from central Northern Territory to central southern Queensland, Australia (Fig. 1).The lack of trees has been attributed to cracking clay soils and fire prior to European settlement (Orr, 1975).The grasslands are dominated by C 4 perennial grasses, among which Mitchell grass (Astrebla spp.), feather-top wire grass (Aristida spp.) and Blue grass (Dichanthium spp.) are the commonest.These grasslands support an extensive pastoral industry, and shifts in species composition and vegetation cover in response to sheep or cattle grazing have been documented (Orr, 1980a;Foran and Bastin, 1984).The climate is predominantly tropical semi-arid and exhibits a large north-south rainfall gradient.Annual rainfall is associated with pronounced wet-dry seasonality.Most of the Mitchell grasslands sites receive between 300 and 400 mm of precipitation during the wet period from November to April (Fig. 2a).Precipitation diminishes in amount and seasonality towards the southeast (Fig. 1).Maximum temperatures remain high throughout the year, though winter (July) minima are significantly lower in southern Queensland compared with the Northern Territory (Fig. 1).

Datasets
Four hundred sites were randomly chosen within the area of Mitchell grasslands, as classified in the digital map of Australia's Native Major Vegetation Subgroups (Australian Government, Department of the Environment, Water, Heritage and the Arts, (http://www.environment.gov.au/erin/nvis/mvg/).To avoid potential edge effects, we excluded from our dataset all 0.05 degree MODIS cells that were directly adjacent to other vegetation subgroups, usually Eucalypt open woodlands.Vegetation physiognomy for all sampled points was further checked using Ikonos images available on Google Earth version 4. 3 (http://earth.google.com/).
Annual vegetation often dominates in heavily grazed areas surrounding cattle watering points, causing the peak NDVI signal to occur more rapidly after the first rainfalls and to be narrower than for perennial-dominated communities (Pickup Biogeosciences, 7, 907-920, 2010 www.biogeosciences.net/7/907/2010/, 1998).The location of permanent and semi permanent sources of water for stock (bores and dams) was obtained from Geoscience Australia's Geodata Topo 250K vector product (http://www.ga.gov.au/mapspecs/250k100k/).These data were used a posteriori to test for a potential effect of the distance to the nearest watering point on model performance.
To calculate NDVI, we used reflectances in the red and near infrared wavebands from MODIS (MOD09A2 Collection 5) obtained from the Oak Ridge National Laboratory Distributed Active Archive Center (http://daac.ornl.gov/MODIS/modis.html).The period analysed was from the beginning of July 2000 to end of June 2008 for a total of 368 periods of 8 days each (except for the last period of each year where the period was shorter).Reflectances were averaged over a 2×2 km area around the target point.We tested different areas (4, 16 and 25 km 2 ) and found no significant effect on the analysis described below.Pixel values that did not have the highest quality flag value were discarded and missing values were interpolated using a two order polynomial fitting method.Because less than 2% of the data was missing, it is unlikely that the details of the gap-filling method will significantly affect the final results.
Following Carlson et al. (1995), we normalized the NDVI values to estimate the green vegetation fraction V , as follows: where NDVI obs is the observed pixel value, NDVI 0 is that for bare soil and NDVI ∞ is that for full vegetation cover.Different values for γ have been used in previous studies: γ =1 in a linear mixing model (Gutman and Ignatov, 1998) or γ =2 as in (Choudhury et al., 1994;Seaquist et al., 2003).
There was no qualitative change in the relative performance of the models when γ was varied between these two values and results below are given for γ =2.Averaged NDVI values during the dry season (July) were used to estimated NDVI 0 for each soil type as soil properties may affect reflectances (Montandon and Small, 2008).A value of 0.75 was assigned for NDVI ∞ , corresponding to the averaged maximum NDVI observed during the wet season in the most productive North-Western Mitchell grasslands near the Gulf of Carpentaria.Note that a lower value of 0.5 was reported for African tropical grasslands in Seaquist et al. (2003) using NOAA 8 km Pathfinder Land data archive.
Before modelling, time series of V were filtered to suppress unusually high or low values.We applied a weighted moving average filter in which the original data value V at t i is replaced by the value of a quadratic polynomial fitted for 2n + 1 points, with n being the number of time steps preceding or following t i (Savitzky and Golay, 1964).The Savitzky-Golay filter has been shown to perform well for minimizing noise in NDVI time series (Hird and McDermid, 2009).We used a moving window of length n = 2, which is narrow enough to follow rapid changes in NDVI while still able to reduce noise efficiently .
Daily time-series of precipitation (P ) and potential evapotranspiration (E) were retrieved from the 0.05 degree resolution SILO database (Australian Government, Bureau of Meteorology http://www.bom.gov.au/silo/).E refers to a Priestley-Taylor estimate of potential evapotranspiration given by 1.26s R n /(s + γ ), where R n is the net radiation absorbed by vegetation and soil, s is the slope of water vapour saturation versus temperature curve and γ is the psychrometric coefficient (Priestley and Taylor, 1972).Cumulative values of P and E matching the 8-day MODIS time periods were prepared for each site.Gridded data were produced by spatial interpolation of ground-based observations as de-scribed in Jeffrey et al. (2001).This dataset is useful for capturing the regional-scale and broad inter-annual patterns of climate variability, but it is of more limited use at finer scales, e.g.rainfall patterns of individual events.For the sampled area, the nearest weather station or rain gauge did not exceed 30 km for 80% of the sites.
Soil data were obtained from the Atlas of Australian Soils (McKenzie et al., 2000;McKenzie and Hook, 1992).The dominant soil type in the Mitchell grassland ecoregion is cracking clay soil.Most of the eastern and northern sites are deep cracking grey soils, coded Ug5.2 and hereafter named "grey soils", while the south-western sites in Queensland exhibit brown-red clay soils, coded Ug5.32 and hereafter named "brown-red soils".Maps for saturated soil water content (W sat ) and soil depth (Z) given in McKenzie et al. (2000) were used to calculate a common mean value of 450 mm for soil available water capacity (W cap ).We did not have enough data to discriminate soil types on the basis of W cap .

Modelling approach
We compared two classes of low-dimensional models, hereafter M1 and M2, with the aim of predicting 8-year long time series of soil water content (W in mm) and vegetation cover (V dimensionless).Both models are biophysical in origin, i.e. based on conservation principles, rather than being statistical models based on purely empirical curve fitting.We attempted to incorporate the processes and their feedbacks that govern the mass conservation balance of soil water (W ) and leaf carbon store represented by the proxy variable V .Models M1 and M2 differed in the way these processes are represented and coupled as detailed below.
In model M1, changes in V are linearly related to changes in W and there is no feedback between V and W .Because we assume a linear response of V toW , M1 is referred to as a linear model hereafter.This model is given by: Soil moisture content is calculated using a one-layer bucket model with precipitation (P ) as input, and evaporation and run-off as outputs.Total evaporation is represented as a function of potential evapotranspiration (E), relative soil water content (W/W cap ) and an exponential decay parameter α 1 .
Note that there is no distinction between bare soil evaporation and plant transpiration in M1.Initial values of V and W at t = 0 are specified by condition (i).Soil moisture content W is allowed to vary in the interval [0, W cap ].All precipitation is assumed to run off when W reaches W cap (condition ii).
The parameter α 2 determines the linear sensitivity of V to changes in We, i.e. the soil water that can be extracted by plant roots.In M1A, there is no threshold for determining We, i.e. we set α 3 = 0 and hence We = W , whereas in M1B, the threshold α 3 is an optimized parameter and We = W −α 3 ≤ W (condition iii).
V is allowed to vary in the interval [V min ,V max ] (condition iv).V max is a site-specific value that corresponds to the maximum value of V for the time series.Finally, M1 incorporates a lag in the number of time steps (L) to account for a delayed response of leaf growth to seasonal water availability.L has a considerable effect on model performance and all analyses were performed with L = 2 (16 days) since this consistently gave the best results.Note that including L adds one more parameter to M1.
In model M2, changes in V are non-linearly related to changes in W and there is a feedback between V and W through explicit incorporation of plant transpiration.Because we assume a non linear response of V toW , M2 is referred to as a nonlinear model hereafter.This model is given by: Soil moisture content is calculated using a one-layer bucket model with precipitation (P ) as input, bare soil evaporation, plant transpiration and run-off as outputs.Bare soil evaporation is represented as a function of the bare soil fraction (1−V ), soil relative water content (W/W cap ), potential evapotranpiration (E) and a parameter β 1 .Plant transpiration is modelled as a water-limited process dependent upon the fraction of plant cover (V ), the amount of extractable water by roots (We) and a parameter β 2 that accounts for plant water extraction ability.As for M1, we distinguish two cases for the calculation of We: in M2A, We = W , i.e. we set β 5 = 0, whereas in M2B, the threshold parameter β 5 is optimized and it is possible that We < W , i.e. the amount of water extractable by roots can be less than the total amount of water in the soil.
The dynamics of V is governed by a growth and a mortality term.Growth is represented by a logistic equation (Verhulst, 1838) with a growth rate, β 3 , dependent on the relative amount of extractable soil water content, i.e. the term We / (W cap − β 5 ) and a carrying capacity V max .Leaf mortality is a function of V and a decay rate parameter, β 4 .A non-zero minimum value of V , V min , is required to ensure initiation of leaf growth at the start of the wet season in M2 (see Dickinson et al., 2008 for further discussion).This small amount of persisting aboveground biomass may correspond to nearground leaf and shoot primordia of perennial grasses.Choosing contrasting V min values did affect parameter estimates in M2, but did not change the goodness-of-fit statistics.Here, we set V min = 0.001.For consistency, the same lower limit value for V was used in M1 and M2 although V min had no impact on M1 parameter estimates.
There is no need to include a lag parameter in M2 because the logistic growth model is able to simulate low responsiveness of V to initial increase of W .A similar coupled nonlinear relationship between V and W has been recently proposed by De Michele et al. (2008).

Model calibration and validation
Parameter estimation was performed using calibrating datasets comprising 100 randomly selected sites.Parameter estimation only relied on the available remotely-sensed data of greenness since we had no soil water content data to add further constraints on parameter estimates.Parameters were optimized using the R package rgenoud (Mebane and Sekhon, 2009).This optimization method combines a derivative-based quasi-Newton algorithm developed by Byrd et al. (1995) and a genetic algorithm.Derivative-based methods are efficient to find local optima.However, for difficult convergence problems where multiple local minima exist, these methods are notoriously poor at finding a global solution.Therefore, a genetic algorithm is necessary for a thorough search of the parameter space.The quasi-Newton algorithm allows setting of the lower and upper bounds for a given parameter.This flexibility was required to constrain all parameters to be positive.In addition, β 4 the fraction of existing V that disappears at each iteration (see Eq. 3) was constrained to be less or equal than one.For each of the 100 calibrated sites, a Coefficient of Variation of the Mean Absolute Error, hereafter CVMAE, was calculated as the Mean Absolute Error (MAE) normalised to the mean of observations (Eq.4).Because calibrating sites differ by their mean vegetation cover, MAE naturally tend to increase with increasing mean vegetation cover.We used CVMAE to ensure that each calibrating site had a similar weight on the final estimate of goodness-of-fit.To remove the undesirable influence of outliers, we minimized the median CVMAE value of the 100 training sites, instead of its mean or its sum.Therefore, the objective function, F , was: where t is time in periods of 8 days, n is the number of time periods, N=100 is the number of calibrating sites, V i,p and V i,o are the predicted and observed values of vegetation cover for site i, respectively, and Vi,o is the mean vegetation cover of site i.We followed recent recommendations on model performance evaluation (Willmott andMatsuura, 2005, 2006) and used MAE as the best measure of average error magnitude of a model.This is because goodness-of-fit statistics based on the squaring of an error term, e.g.Mean Square Error (MSE) or Root Mean Square Error (RMSE), are more sensitive to the distribution of errors and often lead to overweighting outliers.However, we have tested alternative cost functions with squared error terms and this did not change the conclusions of our study.We performed 30 different calibration runs, each conducted with a different subset of 100 sites.
The mean values of each parameter resulting from these calibration runs are given in Table 1.
Parameter estimates from each calibration run were used to predict time series of V and W for the remaining 300 testing sites.Several goodness-of-fit statistics were then calculated including R-squared, MSE, CVMSE and their rootmean-square RMSE and CVRMSE, but again MAE and CV-MAE were primarily used to assess model performance.The advantage of using the sum of squared errors as goodnessof-fit statistics is that they can be additively partitioned into a systematic and an unsystematic component as described by Willmott (1982).For example, in the case of CVMSE: with where a and b are the intercept and slope of the least square regression between the predicted and observed values of V for the calibrating site i, respectively.From Eq. 5, the percentage of systematic error was calculated as 100*(CVMSE s,i /CVMSE i ) and that of unsystematic error as 100*(CVMSE u,i /CVMSE i ).A systematic error approaching zero is indicative of a model structure that adequately captures the system dynamics.Using model II regression, or Standardized Major Axis (SMA) regression (Warton et al., 2006), we also tested whether the slope b was significantly different from 1, and the intercept a significantly different from 0.
To test whether the change from M1A to M1B, and from M2A to M2B, provides a better fit, we estimated the likelihood ratio: where MSE A (MSE B ) is the Mean Square Error of model M1A or M2A (M1B or M2B).The likelihood ratio significance test was based on the assumption that LR follows aχ 2 distribution with one df (i.e. the difference in number of parameters between models A and B).The Akaike's Information Criteria (AIC) was also evaluated for each model as: AIC = 2p + log (L), where p is the number of parameters and L is the maximum likelihood of the model.Finally, we conducted sensitivity analyses by varying one model parameter while keeping the other parameters constant and calculating the relative effect on CVMAE.We did not have any firm uncertainty to provide for NDVI input data, except one from subjective expert judgment (Raupach et al., 2005).Both models exhibited the same sensitivity to this source of uncertainty (data not shown).We also examined the relationships between model residual and distance to the nearest weather station as a way to assess sensitivity of model performance to accuracy of climate data.
Numerical simulations, statistical analyses and all graphics were performed within the R software environment (R Development Core Team, 2007).The source code is available upon request.All computations were performed on the cluster HealthPhy (CIMENT, Université J. Fourier -Grenoble I).

Results
Figure 2a shows that there is a positive relationship between the maximum vegetation cover V max and the mean precipitation during the wet summer period (R 2 = 0.31, P < 10 −4 ), though there was large scatter around 300-400 mm.No such relationship is observed for the dry winter period (Fig. 2b) and we did not find any correlation between V max and minimum, maximum or mean temperature (data not shown).
Table 1 summarizes the calibrated values of the parameters for the four models.Estimates from different calibration runs showed similar variation, i.e. coefficients of variation ranging from 10% to 20%, with the noticeable exception of β 5 for M2B and to a lesser extent of β 2 for M2A (Table 1).There was a strong discrepancy between M1B and M2B in the estimate of the threshold value for W , with α 3 = 72 mm and β 5 = 10 mm.With such a low value for β 5 , it is not surprising that the estimates for the other four parameters, from β 1 to β 4 , were very similar between M2A and M2B (Table 1).
Whatever the goodness-of-fit statistic, the performance of M1A was markedly lower than M1B (Table 2).The loglikelihood ratio test showed that there was a significant improvement of M1B over M1A by adding the extra parameter  where CV is the coefficient of variation resulting from the calibration runs (see Table 1).
α 3 in the optimization process (Table 2).While there was no marked contrast in the overall accuracy between M1B, M2A and M2B as measured by MAE, CVMAE or by CVRMSE, there are significant differences in error partitioning, with an inflated systematic error of 55% for M1B compared to M2A (31%) or M2B (29%).The higher bias in M1B was due to a consistent underestimation of V (SMA slope of 0.77).Ex-cept for M1A, no model had a SMA intercept that differed significantly from 0 (Table 2).The addition of one parameter, β 5 , in M2B did not lead to significant improvement of nonlinear model (  Sensitivity analysis indicated that both M1B and M2A were fairly robust, because a ±10% change in each parameter led to less than a 3% change in CVMAE (Fig. 3a and b).When ranges of parameter values obtained from the calibration runs were taken into account (Table 1), the sensitivity to parameter uncertainty became higher for M2A than for M1B (Fig. 3c and d).For M2A, sensitivity to β 2 was much lower than for the other three parameters even if β 2 exhibited a higher CV than the other parameters (Fig. 3d and Table 1).There was significant sensitivity of state output V to β 3 and β 4 i.e. the parameters describing the nonlinear response of plant growth to soil water content in M2.This tends to indicate that M1 and M2 represent two different solutions in the state space and that optimized parameters of M2 yield a solution that differs from a linear model.Two examples of predicted time series of V using M1B and M2A are shown in Fig. 4. The first example (Fig. 4a) corresponds to a relatively dry site with weak seasonality, while the second experiences heavier rainfalls and a more pronounced seasonality in V (Fig. 4b).Soil moisture dynamics predicted with M2A exhibited more rapid changes and was "peakier" than with M1B.This was particularly the case at the end of the wet season, for example in 2006 at the dryer site and in 2004 and 2006 at the wetter site.This affected the predicted V time series which also showed higher contrasts between dry and wet periods in M2A compared to M1B.Yearly maximum values of V tended to be underestimated by M1B, for example in 2001 and 2006 at the dryer site and in 2003 and 2004 at the wetter site.While the timing for leaf onset and leaf offset were correctly predicted by both models, they were less successful in reproducing the multiannual variations in peak amplitude, for example 2003 at the dryer site and 2002 at the wetter site.
The magnitude of error (CVMAE) tended to increase with decreasing summer rainfall for M1B (Fig. 5a), whereas CV-MAE of M2A was much more consistent throughout the rainfall gradient (Fig. 5b).Both models exhibited a significantly higher CVMAE for south-western sites with "brownred" soils than for northern sites with "grey soils" (data not shown), but because soil type and summer rainfalls covaried in space, we were unable to distinguish the relative importance of these two factors.There were no relationships between model performance and distance to the nearest rain gauge (Fig. 5c and d) and distance to the nearest stock watering point (Fig. 5e and f).
Finally, we explored the ability of both models to predict the integrated value of vegetation cover ( V ) from July to June.Consistent with model performance reported above (see Table 2), Fig. 6 shows that V predicted using model M1B were strongly underestimated in the most productive years (slope of SMA = 0.66), whereas the magnitude of error for M2A did not change with soil moisture availability (slope of SMA = 0.95).On the other hand, there were many cases for which M2A underestimated V compared to observations when V < 2 (Fig. 6b).As a result, the SMA slope of M2A shifted to the right and the SMA intercept was significantly different from zero (Fig. 6b).

Discussion
From our exploration of the ability of linear and nonlinear models to predict NDVI-derived vegetation cover in semiarid perennial grasslands, we conclude that that both types of models had a similar Mean Absolute Error.However, the nonlinear model had the following advantages: (1) it exhibited far less systematic error (Table 2); (2) it had a better ability to capture the sharp transitions in leaf cover, especially under high seasonality of precipitation (Fig. 4); (3) its performance did not deteriorate in the driest sites (Fig. 5a and  b); and (4) its parameters are more meaningful because the model captures the fundamental feedbacks between soil and plant growth through a process-based approach.The main caveat of the nonlinear model was its slightly greater sensitivity to parameter uncertainty (Fig. 3c and d).These points are discussed below.
Most previous studies that have examined relationships between NDVI and climate have done so with linear correlation or regression analyses.Model performance has usually been evaluated using Pearson correlation coefficients (du Plessis, 1999;Wang et al., 2001;Paruelo and Lauenroth, 1998;Chamaille-Jammes et al., 2006;Richard and Poccard, 1998;Wang et al., 2003) although this metric is a poor indicator of goodness-of-fit (Willmott, 1982;Willmott and Matsuura, 2006).In contrast, this paper compares two differ-P.Choler et al.: Seasonal leaf   ent process-oriented models of phenology and we did not include purely statistical models of NDVI-rainfall relationships.Data not shown indicate that, in our case, purely statistical models give very similar results to M1 but with a larger systematic error and strong deterioration in model performance for the driest part of the rainfall gradient.More than two thirds of the NDVI-vegetation cover variability was explained by the M1 and M2 models (see R-squared values in Table 2) indicating a high responsiveness of these grasslands to soil water balance.Similar values have been reported for other tropical grasslands found on clay soils (Nicholson and Farrar, 1994;Farrar et al., 1994).Our study expands on pre-vious literature because (1) we examined the appropriateness of different process-based model structure to address NDVIprecipitation relationships and (2) we provided a more comprehensive analysis of model performance.In particular, error partitioning indicated that in the case of linear models, no parameter set was able to give unbiased estimates at each time step and along the rainfall gradient (Table 2 and Figs. 4,  6).The nonlinear model greatly reduced the systematic component of the error but not the overall uncertainty in model predictions (Table 2).There is no single criterion for a model to be acceptable (Rykiel, 1996), but we advocate that for a given error magnitude it is preferable to retain the model that exhibits the minimum bias.
Several studies have focused on the relationships between climate and yearly or seasonally integrated values of NDVI, used as a proxy of productivity (examples in Paruelo and Lauenroth, 1998;Propastin and Kappas, 2008).Such an approach cannot identify fine timing of leaf onset, growth and offset events which are crucial for land surface modelling (Pielke et al., 1998;Dickinson et al., 1998;Bounoua et al., 2006).In the semi-arid tropical grasslands investigated here, the dynamics of leaf cover are characterized by abrupt transitions associated with seasonal water availability (Fig. 4).The amplitude of change in leaf cover in response to soil water content at the start of the growing season was better captured by a nonlinear logistic growth than by the linear model.The same holds for the senescence phase for which the feedback of plant cover on soil moisture was critical to simulate the decline in V .Sharp peaks of leaf biomass in response to resource pulses have long been viewed as archetypal examples of nonequilibrium dynamics (Seastedt and Knapp, 1993;Blair, 1997).The structure of nonlinear models provides much more flexibility to capture these transient maxima.Conversely, finding a common threshold value of W for peak initiation remains a delicate task for nonlinear models, as evidenced by the high coefficient of variation reported for β 5 in M2B (Table 2).Indeed, nonlinearities in M2 naturally tend to exacerbate errors associated with uncertainty in parameters (compare Fig. 3c and d).The nonlinear models satisfactorily predicted peak amplitude and integrated values of plant cover when soil available water was large (Figs. 4  and 6), but an unresolved issue is their tendency to predict lower values when seasonal soil water availability was more limited (Fig. 6b).
The concept of rain-use efficiency has arisen from a number of studies that examined how precipitation drives phenology in the tropics (Nicholson and Farrar, 1994).Empirical estimates of rain-use efficiency typically result from applying a low-pass filter to precipitation time-series and then correlating vegetation indices to these smoothed and lagged time-series of predictors.By contrast, a nonlinear modelling framework decomposes rain-use efficiency into different components.For example, it makes explicit the relationship between water extraction by roots (β 2 ) and plant growth rate (β 3 ), i.e. the increment of leaf area for a given amount www.biogeosciences.net/7/907/2010/Biogeosciences, 7, 907-920, 2010 of water extracted from the soil.This opens avenues for linking nonlinear phenological models to those used to describe plant physiology or biomass allocation in land surface models.
There are several reasons why both models fail to explain more NDVI-vegetation cover variability, including (1) observation errors, (2) neglect of other forcing variables and (3) failure of both model structures.As an example of possible observational biases, we found abnormally low NDVI values during the build up phase at the start of the wet season, possibly related to insufficient removal of cloud contamination of the NDVI data.It is likely that predictive models of NDVI time series are more sensitive to this source of errors than models using integrated values during the growing season.Second, there are a number of potential site-specific factors not included in the models that might affect plant phenological response such as contrasting soil water holding capacity, soil fertility, species composition, plant functional diversity, and grazing intensity.It is obviously impossible to translate all these effects into a four parameter model.Our preliminary studies did not provide any evidence for an effect of the distance to the nearest watering point, which can be viewed as a surrogate of grazing intensity.Neither did we find that temperature or photoperiod having any significant impact on the leaf dynamics of grasslands in the area covered.Third, the models did not take into account multiannual dynamics.For example, several drought years combined with overgrazing trigger a sharp decline in perennial grasses (Orr, 1980b).Alternatively, two or three consecutive wet years would favour recruitment of new grasses and hence increase the carrying capacity of a given site.Further refinements of the models should be able to combine short-term and long-term effects of the soil water balance.An option would be to add a third state variable representing root storage capacity of carbohydrates.
Both models were ecohydrological in essence because they included equations for soil water content and leaf dynamics.Regrettably, no data are available for soil moisture dynamics and so calibration and validation of models were only based on remotely-sensed estimates of NDVI.To match the scale and the spirit of our study, remotely-sensed estimates of soil moisture would be the most appropriate to provide an independent validation of the model.The available products (RADAR, estimates from NDVI and surface temperature, gravimetric methods) all have problems which prevent them from being used in the short term.These problems include data processing and validation, non-independence from NDVI and too coarse a spatial resolution.At this stage, our models have generated hypotheses on soil water dynamics that need to be further tested as remotely-sensed soil moisture data improve over time.

Conclusions
There have been very few attempts to calibrate and validate coupled dynamical systems model of soil moisture and plant cover by using time series of remotely sensed data (Seaquist et al., 2003;Hess et al., 1996).Indeed, there is still a gap between empirical, remote-sensing oriented modelling of NDVI-precipitation relationships and more conceptual and theoretical efforts towards process-oriented ecohydrological models (Porporato et al., 2002;Rodriguez-Iturbe et al., 2001).Our study contributes to bridging this gap by showing that simple nonlinear models of phenology can provide both elegant and mechanistic understanding on how precipitation variability affects vegetation growth in semi-arid grasslands.Further work should examine whether this mechanistic framework is also appropriate to modelling phenology at a larger scale, especially for other water-limited ecosystems of the world.

P
Fig. 1.Geographical distribution of Mitchell grasslands (grey) and location of the 400 random sites investigated.Left and right panels show monthly means (+ 1 sd) of precipitation (mm) and mean temperature ( • C) over the period 1961-1990 as recorded in 8 weather stations within or nearby the Mitchell grasslands area.Abbreviations are as follows: Bru: Brunette Downs; Bur: Burketown; Ric: Richmond; Ura: Urandangi; Vic: Victoria Rivers.Map uses Albers Equal Area Conic projection.

Fig. 2 .
Fig. 2.Relationship between the maximum vegetation cover (V max ) over the period 2000-2008 and the mean "summer", i.e. from November to April (A) and the mean "winter", i.e. from May to October (B) precipitation.Each point represents one of the 400 random sites investigated.Dotted lines represent median values.

P.
Choler et al.: Seasonal leaf dynamics in tropical grasslands with

Fig. 3 .
Fig. 3. Robustness of M1B (A, C) and M2A (B, D) to parameter uncertainty.Histograms show the effect of an increase (grey) and a decrease (black) of a parameter on the relative change of CVMAE.Magnitude of the change in parameter value is ± 10% (A, B) and ± CV (C, D), where CV is the coefficient of variation resulting from the calibration runs (see Table1).

Fig. 6 .
Fig. 6.Observed vs. predicted values of yearly integrated vegetation cover with model M1B (A) and model M2A (B).Each point represents one site for one year.Dotted lines show the 1:1 relationship and plain lines show the Standardized Major Axes (SMA).

Table 1 .
Calibrated values of parameters for linear (M1) and nonlinear (M2) models (see Eqs . 2 and 3).Mean (± se) and coefficient of variation (CV) were calculated from 30 calibration runs each using 100 randomly sampled calibrating sites.The mean values of the objective function, i.e. the median of CVMAE (see Eq. 4), are also indicated.

Table 2 .
Compared performance of linear (M1) and nonlinear (M2) models (see Eqs. 2 and 3).After each of the 30 calibration runs, goodness-of-fit statistics were calculated for the 300 remaining testing sites.Therefore, reported values are the means of 30 estimates.n is the number of time steps.See Material and Methods for details on statistics.
dynamics in tropical grasslands on model performance estimated by CVMAE.For each of the 400 sites, predicted values were calculated with model M1B (A, C and E) and model M2A (B, D and F) using parameter values of Table1.
Fig. 5. Effect of summer precipitation (A, B), distance to the nearest rain gauge (C, D) and distance to the nearest stock waterpoint (E, F)