The effect of year-to-year variability of leaf area index on Variable Inﬁltration Capacity model performance and simulation of runoff

This study assessed the effect of using observed monthly leaf area index (LAI) on hydrological model perfor- mance and the simulation of runoff using the Variable Inﬁltration Capacity (VIC) hydrological model in the Goulburn–Broken catchment of Australia, which has heterogeneous vegetation, soil and climate zones. VIC was calibrated with both observed monthly LAI and long-term mean monthly LAI, which were derived from the Global Land Surface Satellite (GLASS) leaf area index dataset covering the period from 1982 to 2012. The model performance under wet and dry climates for the two different LAI inputs was assessed using three criteria, the classical Nash–Sutcliffe eﬃciency, the logarithm transformed ﬂow Nash–Sutcliffe eﬃciency and the percentage bias. Finally, the deviation of the simulated monthly runoff using the observed monthly LAI from simulated runoff using long-term mean monthly LAI was computed. The VIC model predicted monthly runoff in the selected sub-catchments with model eﬃciencies ranging from 61.5% to 95.9% during calibration (1982–1997) and 59% to 92.4% during validation (1998–2012). Our results suggest systematic improvements, from 4% to 25% in Nash–Sutcliffe eﬃciency, in sparsely forested sub-catchments when the VIC model was calibrated with observed monthly LAI instead of long-term mean monthly LAI. There was limited systematic improvement in tree dominated sub-catchments. The results also suggest that the model overestimation or underestimation of runoff during wet and dry periods can be reduced to 25 mm and 35 mm respectively by including the year-to-year variability of LAI in the model, thus reﬂecting the responses of vegetation to ﬂuctuations in climate and other factors. Hence, the year-to-year variability in LAI should not be neglected; rather it should be included in model calibration as well as simulation of monthly water balance.


Introduction
The challenge of making accurate runoff predictions using hydrological models under changing or 'non-stationary' conditions, due to either changing climate and/or human intervention, is a significant issue in hydrology [5,26,33]. Rainfall-runoff models that lack representations of biophysical processes, such as vegetation dynamics, have been found to perform poorly when calibrated in a wet climate period and validated in dry climate period [6,25,43]. To address this problem different studies have suggested approaches including calibrating model parameters on a portion of the record with conditions similar to those of the future period to simulate [43], using temporal clusters [9] and adjusting the parameters according to the aridity of the catchment [38]. Rather than calibrating parameters that vary with the condition in the system, understanding the catchment processes and effectively incorporating them into the model may help us * Corresponding author. Tel.: +61 383449799.
E-mail address: ywei@unimelb.edu.au (Y. Wei). to improve model performance by considering the various processes that are modified with changing climate.
A large amount of evidence shows that vegetation is an important component of the hydrological process [14,32,35,40]. Vegetation has a significant role in the partitioning of rainfall into runoff and evapotranspiration (ET) mainly through canopy transpiration and interception loss [44]. Transpiration varies according to physiological (stomatal conductance) and structural properties, mainly leaf area index (LAI) of the vegetation [16], while interception varies according to structural properties of the vegetation and precipitation characteristics [28]. Changes in LAI not only affect evapotranspiration but consequent changes in soil moisture also impact other catchment processes including baseflow, recharge, infiltration excess, saturation excess, subsurface storm flow and catchment wetness [47]. Hence, lack of representation of the year to year variability of the monthly LAI in hydrological models may lead to lower monthly model performance due to underestimation of flow in dry periods and overestimation of flow in wet periods.
Remote sensing provides spatially and temporally variable LAI datasets that help to capture the vegetation dynamics and can  be incorporated in land surface models that include LAI in most evapotranspiration processes. There have been some efforts to exploit remotely sensed vegetation information into hydrological models [40,49]. In their study of the North American monsoon in Western Mexico, Tang et al. [40] applied year to year variable monthly LAI and mean monthly LAI obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) into the Variable Infiltration Capacity (VIC) model to predict evapotranspiration for the years 2001-2008 and validated their results with observations at two eddy covariance tower sites. They found that using mean monthly LAI in VIC biased evapotranspiration estimates by 10-30% due to not representing the year to year differences in vegetation greening onset and dormancy periods. Similarly, Ford and Quiring [14] investigated the effects of using observed monthly LAI compared with mean monthly LAI on simulated soil moisture from the VIC model for the period 2000-2009 in eastern Oklahoma, USA. The authors also compared VIC-simulated moisture results with in-situ soil moisture at different depths and locations and concluded that the models that incorporated observed LAI could better capture the intensity and duration of droughts. To date no studies have addressed the influence of observed monthly LAI on runoff simulation in the VIC model.
Sensitivity studies of VIC showed that the primary vegetation characteristic that affects the hydrological simulations is LAI [23]. Within VIC, LAI controls the partitioning of the canopy storage into canopy evaporation and through fall in vegetated areas. The default canopy storage (mm) value in VIC is estimated as 20% of the LAI [10], where field observations data are unavailable for each land cover type, and the storage is depleted using a nonlinear power function as given by Deardorff [7]. LAI is also used in VIC in the partitioning of water that infiltrates into the ground into evapotranspiration from the root zone soil column and recharge. This partitioning is made through the vegetation where LAI is used to scale the leaf to canopy transpiration as given by [1,11]. Hence, neglecting the year to year variability of LAI in the VIC model might lead to lower performance in runoff prediction and would be expected to be characterized by underestimation of flow in dry periods and overestimation of flow in wet periods, given the specific roles of LAI in the model. This paper addresses the question of whether including vegetation variations that reflect both variations in climate and human activities into a hydrological model improves model performance in terms of runoff simulation under changing climate and catchment conditions. Specifically, we examine the effects of using observed monthly LAI compared with long-term mean monthly LAI on VIC model performance in terms of runoff prediction in the Goulburn-Broken catchment for the period 1982-2012. Our approach is to calibrate and validate VIC using both observed monthly LAI and long term mean monthly LAI and compare the model performance in terms of simulated runoff for fourteen sub-catchments in the Goulburn-Broken catchment of south-eastern Australia. The paper is presented as follows: the study area is described in Section 2, followed by Section 3: the dataset and model setup, the results are presented in Section 4, followed by discussion and conclusion sections in Section 5.

Description of the study area
The Goulburn-Broken catchment (35.8-37.7°S, 144.6-146.7°E) is located in Victoria, south-eastern Australia and is part of the Murray-Darling basin (MDB). Fourteen sub-catchments were chosen for this study (Fig. 1). The Goulburn-Broken catchment covers approximately 24,000 km 2 representing about 10.5% of the total area of the State of Victoria, and 18% of the water supply for Victoria (www.riverfoundation.org.au). The Goulburn-Broken catchment contributes about 11% of the water resources of the MDB. The maximum altitude is approximately 1790 m above mean sea level (AMSL) on the southern side of the catchment and the minimum altitude is 86 m AMSL on the northern side of the catchment. The mean catchment elevation of the selected sub-catchments ranges from 155.83 to 1001 m AMSL.
The climate of the Goulburn-Broken catchment is influenced by mountain ranges with high precipitation in the southern part and lower precipitation in the flat plains of the northern part (declining precipitation from south to north). The long-term (1982-2012) mean annual precipitation peaks at 1632 mm yr −1 in the southern mountainous area, and reaches a minimum of 373 mm yr −1 in the north. In the selected sub-catchments mean annual precipitation ranges from 526 to 1407 mm yr −1 . About 60% of the total annual precipitation occurs in winter and spring; with about 45% occurring in the four months from June to September. The spatial variation in potential evapotranspiration (PET), using the Penman-Monteith formulation of the Food and Agricultural Organization (FAO56) method, is opposite to precipitation and varies from 775 mm yr −1 in the south to 1238 mm yr −1 in the north of the catchment, and in the selected sub-catchments it ranges from 903 to 1132 mm yr −1 ( Table 1). The catchment covers three climate zones based on the Köppen-Geiger climate classification as shown in Fig. 1. The north lowland part of the catchment experiences low annual precipitation and high potential evaporation and is semi-arid (BSk, 9%). The middle section of the Table 1 Characteristics of selected sub-catchments for this study: mean annual precipitation (P), potential evapotranspiration (PET), mean annual leaf area index (LAI) and the percentage of the three land cover types (crop, pasture and tree).  catchment has a hot summer temperate, without a dry season climate (Cfa, 35%). The southern part of the catchment has a warm summer temperate, also without dry season, climate (Cfb, 55%) [34].
Most of the southern part of the catchment is covered by trees: mainly open Eucalyptus tall trees and Eucalyptus woodlands (Fig. 1). The central part and most of the northern part of the catchment are covered by cropland and pasture with irrigated areas mostly found in the north. The land cover type is grouped into three dominant land cover types (trees, pasture and crop) which comprise 47%, 38% and 12% of the entire catchment, respectively, with the rest occupied by water bodies and urban areas.
The seasonal and year to year variability of the areal average LAI of the three dominant vegetation types in the study catchment is shown in Fig. 2. Crop (Fig. 2a) and pasture (Fig. 2b) show much higher LAI seasonality than tree (Fig. 2c), which is predominantly an evergreen genus. The minimum LAI for crop and pasture areas occurs during summer when trees reach their maximum LAI. The deviation of the observed monthly LAI from the long-term mean monthly LAI is observed to be significant in all the three vegetation types (Fig. 2d). The annual LAI time series show declines in annual LAI of crop and pasture and to a lesser extent trees during the recent prolonged drought (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)) in the study area. Climatic and land cover type characteristics for the 14 sub-catchments used in this study are presented in Table 1. There is no irrigation in any of these study sub-catchments.

The VIC model
VIC is a spatially distributed physical-based macroscale hydrological model that balances both water and energy budgets over a grid cell. It has been successfully applied in many settings, from global to river basin scale [23,30,37]. It simulates soil moisture, evapotranspiration, snow pack, runoff, baseflow and other hydrological properties at daily or sub-daily time steps by solving both the governing water and energy balance equations [20]. VIC independently simulates all processes in each grid cell, which are equally spaced. The infiltration and runoff are estimated using the Variable Infiltration Capacity model curve, which uses the soil moisture content of the upper two soil layers to approximate the spatial variability of surface saturation. VIC uses the Penman-Monteith equation to estimate potential evapotranspiration, which requires inputs of maximum and minimum temperature, vapour pressure, wind speed and solar radiation data. Overall, the model represents the spatial variability in climate, vegetation and physical properties of soil [2][3][4]. VIC has been successfully applied in many settings, from global to river basin scale [23,30,37]. In Australia, Sivapalan and Woods [39] used the VIC model to assess the effects of spatial variability of rainfall on soil moisture and Kalma et al. [18] compared the temporal trend in soil moisture storage between one derived from spatially distributed field measurement at catchment and the one simulated with the two layers VIC model and reported that VIC could simulated soil moisture status at catchment scale. Similarly, Western et al. [47] assessed VIC simulated soil moisture in a small catchment in southern Victoria, Australia and reported the time-series of spatially averaged internal total moisture storage was consistent with observations despite different assumptions in the statistical distribution of soil moisture storage used by VIC and their observed soil moisture dataset. After Liang et al. [20] modified the VIC model for three layers and spatially varied vegetation, Zhao et al. [50] assessed the model for prediction of runoff in a catchment located in south-eastern Australia. The advantage of the VIC model over a simple conceptual rainfall-runoff model is that it uses a "mosaic" scheme that allows spatial representation of gridded topography, infiltration rate, soil properties, climate variables and land cover which are important attributes in modelling runoff under spatially heterogeneous conditions.

Dataset
The input data used for the hydrological modelling include the daily precipitation, maximum and minimum temperature, vapour pressure and solar radiation data. They were obtained from the Australian Water Availability Project (AWAP) of the Bureau of Meteorology [17]. The gridded daily wind run data were obtained from [24] that was generated from point measurements. All data have a spatial resolution of 0.05°× 0.05°(approximately 5 km × 5 km), and the period 1982-2012 was selected for this study based on the availability of LAI data. The elevation data were collected from the GEO-DATA 9 Second Digital Elevation Model (DEM-9S) Version 3 [15] with a spatial resolution of 9 s (approximately 250 m). The elevation data were resampled to a resolution of 0.05°× 0.05°using the spatial average. The daily runoff data at the outlet of the selected calibration sub-catchments were obtained from the Victorian water resources warehouse (http://data.water.vic.gov.au/monitoring.htm).
The land cover input data were derived from the National Dynamic Land Cover Dataset which provides a land cover map for the whole of Australia at a resolution of 0.00235°× 0.00235°(approximately 250 m × 250 m) and can be accessed at (http://www.ga. gov.au/metadata-gateway/metadata/record/gcat_71071). The dataset was developed using the MODIS satellite and validated using a fieldgenerated land cover map [22]. For this study the land cover class was regrouped into three dominant classes: trees or forest, grass or pasture and crop; and resampled to the AWAP resolution to be spatially consistent with the other input data. Then the fraction of each land cover type inside each VIC model grid cell was computed and provided as an input to the VIC model. LAI data were collected from the Global Land Surface Satellite (GLASS) product which is available for download from Beijing Normal University (http://www.bnudatacenter.com). The dataset was derived by combining the Moderate Resolution Imaging Spectroradiometer (MODIS) and the Advanced Very High Resolution Radiometer (AVHRR) satellite products at 0.05°r esolution for the globe [19]. The dataset has been compared with other remotely-sensed LAI products and found to be in good agreement [12]. The root distribution in three soil layers was derived from the global ecosystem root distribution dataset [36]. The soil parameters in the VIC model running resolution were derived from the 5 min resolution Food and Agriculture Organization dataset [13]. The first soil layer was set to 0.1 m following Liang et al. [20] and the other two layer depths were calibrated. The empirical Arno curve was used to generate the baseflow based on the soil moisture content in the bottom layer [4]. The total runoff at each grid cell was routed through a defined river network that was generated from the digital elevation model using the algorithm developed by Lohmann et al. [21].

VIC model calibration and validation
After all the necessary input data for the model were collected and prepared, the VIC (version VIC 4.1.2g) model was calibrated on selected unregulated sub-catchments in the Goulburn-Broken catchment (Fig. 1). The seven most sensitive VIC model parameters (b, Ds, Ws, Ds max , d 2 , d 3 and exp) according to [8] were calibrated for each sub-catchments separately but were considered uniform within a sub-catchment. The physical meaning and possible ranges of values of these parameters are listed in Table 2. These parameters ranges were used as a boundary to guide the calibration algorithm.
This study employed the Multi-Objective Complex Evolution (MOCOM-UA) algorithm to optimize parameter values by minimizing, or maximizing, the objective function specified by the user [48]. The MOCOM-UA algorithm uses a multi-objective, rather than a single objective, function and is an advance over the Shuffled Complex Evolution Metropolis (SCEM-UA) global optimization algorithm [45]. The user sets the initial population size and the number of samples to be taken from that initial population to evolve towards a set of solutions stemming from a stable distribution Pareto set based on the concept of Pareto dominance [48]. The MOCOM-UA algorithm was implemented on each of the selected sub-catchments separately to calibrate the model against the observed runoff. The model was first calibrated for the entire period , then using the calibrated parameters as initial guesses, the model was re-calibrated for the period 1982-1997 and validated for the period 1998-2012. During the calibration, VIC ran on a daily basis but the objective function was calculated on a monthly basis. Three criteria (objective functions) were used to evaluate the model's performance during calibration: the Nash-Sutcliffe efficiency (NSE) [29] between observed and simulated flow Eq. (1), the logarithm of Nash-Sutcliffe efficiency (log NSE) which penalizes errors at peak flow Eq. (2), and the percentage bias (PBIAS) from the observed mean flow Eq. (3).
where Q i obs is the ith observed flow, Q i sim is the respective ith simulated flow from the model, Q mean is the mean of the observed flow for the calibration period and n is the total number of observed flows. Here the MOCOM-UA algorithm was set to maximize the NSE and log NSE and minimize the PBIAS in search for the optimal parameter set.

Assessing the two forms of LAI on VIC model performance
To assess the effects of using observed monthly LAI (hereafter LAI) compared with long-term mean monthly LAI (hereafter LAI mean ) on VIC model performance a systematic test was performed in validation mode for the data from 1998 to 2012 at each of the selected subcatchments. Model calibration was undertaken twice, once using LAI and once using LAI mean for the period 1982 to 1997. Each set of calibrated parameters for a given LAI were used to simulate runoff in VIC in validation mode forced with the matching LAI (LAI or LAI mean ) data. Model performance criteria were calculated from observed and simulated runoff for each run and each sub-catchment separately and compared to assess differences in model performance (Fig. 3).

Assessing VIC model performance when model and input errors are removed
The impact on model performance of using observed monthly LAI compared with long-term mean monthly LAI might be influenced by errors in input data and/or model structure. In order to minimise these effects on our results we adopted a synthetic data methodology similar to Oudin et al. [31], where they assessed the effect of using long-term mean monthly (simple) or observed monthly potential evapotranspiration (complex) on rainfall-runoff model performance. The basis of the synthetic data methodology is to take runoff generated from a calibrated rainfall-runoff model and then to re-calibrate the model against the simulated runoff using the same input data.
In this way all model errors or input errors are removed and the re-calibrated model should reproduce the simulated runoff perfectly (NSE = 100%). In order to assess any degradation in model performance between simple and complex data inputs, the model should initially be calibrated using complex input data to generate the synthetic runoff against which the model will be re-calibrated twice; once using complex and once using simple data inputs. The complex re-calibration is expected to have NSE = 100%, while any degradation in model performance (NSE < 100%) for the simple re-calibration is due to reduced information in the simple input data.
Here we used simulated runoff from the calibrated VIC model fed with observed monthly LAI (Section 3.3, complex data input) to generate the 'synthetic' runoff, against which we calibrated VIC twice, once using LAI (complex) and once using LAI mean (simple) input data. For the case where we calibrated using LAI the model performance criteria are expected to approach 100% as the model is capable of reproducing the synthetic runoffs. However, for the case where we calibrated with LAI mean the model performance is expected to be reduced and the degree of performance reduction reflects the impact of using LAI mean to calibrate and run the model rather than LAI.

Assessing influences on the sensitivity to LAI variability
To investigate the effect of using the two LAI inputs on simulated monthly runoff over time, two modelling experiments were conducted using the VIC model calibrated with LAI on the study subcatchments. In the first experiment VIC was fed with LAI and then fed with LAI mean to produce monthly runoff from each of the respective LAI inputs. The change in runoff (CS, Eq. (4)), the monthly leaf area index elasticity of runoff (ɛ i , Eq. (5)) and the root mean square of monthly leaf area index elasticity of runoff (ɛ rms , Eq. (6)) were calculated between the two simulated runoff series on a monthly basis as follows: where Q LAI is the simulated runoff from VIC using LAI, Q LAImean is the simulated runoff from VIC using LAI mean , i is the month and n is the number of all months. The root mean square of the leaf area index elasticity of runoff was then plotted in a series of scatter plots against the various subcatchment characteristics: mean annual precipitation, mean potential evapotranspiration, dryness index, percentage of tree cover, mean elevation, and catchment area.

Model calibration results
VIC was calibrated for the 14 selected sub-catchments with different climate and land cover composition that are representative of the main runoff generating regions of the Goulburn-Broken catchment. In Table 3 the calibrated model parameter values are listed for each sub-catchment together with the three model performance criteria during both the calibration (1982-1997) and validation (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) periods. Most of the calibrated sub-catchments have NSE of more than 70% during both calibration and validation periods.
In most of the selected sub-catchments the simulated runoff for both calibration and validation periods met the "satisfactory" criteria according to [27], with NSE > 50% and the percentage absolute bias is generally less than 25% during calibration and validation periods. The few sub-catchments that did not meet this criterion showed high biases in both calibration and validation. The temporal variability of runoff was well captured by the model in all calibrated subcatchments. The simulated and observed runoff for the three subcatchments where VIC performed best, average and worst are shown in Fig. 4. These selections did not include (Catchment 14) where VIC was less efficient. Although VIC captured the temporal variability of runoff well, there were some systematic biases in the runoff simulated. The model overestimates peak flow in a few cases and underestimates low flow in most of the sub-catchments. The sources of these biases need to be investigated in order to understand the performance of the model. To do this, the estimated monthly biases are plotted against the monthly climate inputs: precipitation, temperature and leaf area index (not shown here). The calibrated subcatchments showed no relationship between AWAP gridded climate data and simulated runoff biases. The biases are likely related to the model structure [18] rather than model inputs.
The spatial distribution of the calibrated parameters varies from sub-catchment to sub-catchment mainly due to large differences in annual precipitation. Sub-catchments with high annual precipitation (see Table 1) have higher values for the infiltration curve shape parameter (b), which indicates sub-catchments with lower infiltration and higher event runoff (see Table 1). This parameter varies from 0.0004 for catchment 13 to 0.244 for catchment 3 (Fig. 1) with an average value of 0.06 across all the calibrated catchments. The parameter for multiplying the exponent of the Brooks-Corey drainage equation (exp), which controls the vertical drainage between layers, varies from 3.2 (in catchment 12) to 1.15 (in catchment 3) with a spatial average of 2.4 times the value derived from the FAO soil map dataset. The smaller the value of the parameter (exp) the catchment has, the larger the drainage between the layers, for the same soil moisture, which then results in larger baseflow. The soil water storage capacity depends on the thickness the soil layers and the fixed threshold moisture contents (i.e. wilting point, field capacity and porosity). Plant access this soil water based on a root profile specified from the vegetation type. Hence, the thickness of the second layer (d 2 ) can influence the plant transpiration which contributes the largest component of evapotranspiration. The thickness of the second layer (d 2 ) varies from 0.52 m (in catchment 6) to 1.99 m (in Catchments 1 and 14). The third layer depth (d 3 ), which is not accessible to the plants, determines the baseflow and varies from 0.28 (in Catchment 3) to 2.0 (in Catchment 14). A weak inverse correlation was found between the thickness of the second layer and the third layer which might be due to a parameter identifiability issue.

Effect of the two forms of LAI on VIC model performance
In Section 4.2.1 the impact of the two forms of LAI input on VIC model performance in validation mode was investigated using observed runoff as a reference. Then in Section 4.2.2, the effect of LAI on VIC model performance in validation mode when model structure and input data errors are removed was investigated. Finally in Section 4.2.3, the sensitivity of runoff to LAI variability was investigated.

Effect of LAI on VIC model performance against observed runoff
The three assessment criteria (Nash-Sutcliffe, logarithm Nash-Sutcliffe and percentage bias) were calculated using observed and simulated runoff in validation mode with inputs of LAI and LAI mean ( Table 4). The three criteria have slightly different behaviours since the Nash-Sutcliffe is more influenced by peak flows, the logarithm Nash-Sutcliffe gives more emphasis to low flows, while the percentage bias provides the tendency of the model to overestimation or underestimation. The Nash-Sutcliffe model efficiency values show that calibrating VIC using LAI yields better or equivalent model performance than calibrating the model using LAI mean for all sub-catchments ( Table 4). The maximum improvement in Nash-Sutcliffe efficiency was found to be 25% in Catchment 11, which is covered predominantly by pasture, while no improvement was found in Fig. 4. Observed versus simulated runoff from using observed monthly LAI for selected calibration sub-catchments (1-13) with high, medium and low performance. Less biased a Differences: if the bias from observed monthly LAI is closest to zero = less biased, if the bias from observed monthly LAI is furthest from zero = biased.  Catchment 4, which is dominated by trees. In addition, some of the sub-catchments also show some improvement in the other two model efficiency criteria when LAI was used (Table 4). In interpreting these results it is important to recall that maximising Nash-Sutcliffe efficiency was the consistent calibration objective. The other two objectives (logarithm of Nash-Sutcliffe efficiency and percentage bias) were used when no difference in the Nash-Sutcliffe efficiency exist. The observed monthly LAI runoff produced less bias prediction in most of the sub-catchments, however highly forested sub-catchments showed some bias.

Effect of LAI on VIC model performance when model and input errors are removed
To assess the validity of the results in Section 4.2.1 we considered a separate analysis where VIC was re-calibrated against synthetic runoff data generated from the model calibrated with LAI, once with LAI and once with LAI mean . The reason for doing this is that it eliminates model input uncertainty and model structural uncertainty from the comparisons. The degradation in model performance when LAI mean was used instead of LAI is shown in Table 5. Sub-catchments located in the high annual precipitation zone covered pre-dominantly with trees showed less than 3% degradation in NSE model performance using LAI mean , whereas a tree dominated sub-catchment in a low precipitation area (catchment 14) showed 10.7% degradation in NSE. In contrast, catchments 11-13 (Tables 1 and 5) where the dominant land cover is pasture, showed more than 3% degradation in NSE model performance. The other performance criteria (log NSE and PBIAS) showed little change when the model was fed with either form of LAI data, which indicates that the low flows and runoff ratio are insensitive to the year to year changes in monthly LAI.

Influences on runoff sensitivity to LAI variability
The spatial pattern of LAI sensitivity is found to be related to the spatial patterns of precipitation and the distribution of land cover type. The Box-Whisker plots of change in monthly runoff simulated using LAI compared with using LAI mean are plotted in Fig. 5. The change in simulated runoff showed high seasonality. In a majority of the sub-catchments winter and spring runoff were most influenced by LAI input type. Overall, the grey range (25-75% quantiles) of most of the Box-Whisker plots in Fig. 5 are neither completely positive nor negative, suggesting that simulated runoff from using LAI mean was underestimated in years that have precipitation above mean annual precipitation and overestimated in years that have precipitation below mean annual precipitation. The spatial distribution of these changes in monthly runoff varied among study sub-catchments. Simulated monthly runoff using LAI mean was overestimated by up to 25 mm during wet periods and underestimated by up to 35 mm during dry periods when compared with simulated monthly runoff from using LAI. When sparsely forested sub-catchments (Catchments [11][12][13], are compared with tree dominated sub-catchments (Catchments [3][4][5][6][7][8], the later showed the smallest deviations in simulated monthly runoff whether LAI or LAI mean were used. The leaf area index elasticity of runoff is highly correlated to the dryness index (R 2 = 0.97), mean annual precipitation (R 2 = 0.95), percentage of forest cover (R 2 = 0.54) and mean annual potential evapotranspiration (R 2 = 0.89) (Fig. 6a-d). Sub-catchments located in the arid part of the study area showed the highest LAI elasticity of runoff. Mean catchment elevation also influences the elasticity of runoff to LAI (Fig. 6e), although these two variables are highly cross-correlated with mean annual precipitation (R 2 = 0.7). However the size of the catchment has been found not to have an influence on the leaf area index elasticity of runoff.
The results of sensitivity analyses for the simulated mean annual runoff of the calibrated sub-catchments to changes in ±10%, ±30% and ±50% of the mean monthly LAI are shown in Fig. 7. The sensitivity of simulated mean annual runoff to changes in mean monthly LAI for sub-catchments with a high proportion of area covered by trees (Fig. 7a) and a low proportion of area covered by trees (Fig.  7b) for low, medium and high mean annual precipitation are shown. When comparing Fig. 7a and b there is a significant difference in the sensitivity of mean annual runoff to mean monthly LAI between highly forested and sparsely forested sub-catchments. Highly forested sub-catchments (Catchments 1, 4 and 6) exhibited lower sensitivity than sparsely forested sub-catchments (Catchments 11, 13 and 12) for a small difference in their mean annual precipitations. Both land cover groups showed some increase in sensitivity to LAI at lower LAI values. A spatial difference in simulated mean annual runoff in response to the same change in LAI mean was observed in sub-catchments with similar vegetation cover, which is likely due to differences in dryness index (Fig. 6c) and difference in percentage of the forest cover (Fig. 6d).

Discussion
The performance of the VIC model was found to be good for the vast majority of sub-catchments in terms of Nash-Sutcliffe efficiency of flows and log flows (Table 3). In a few cases, bias in simulated runoff was significant which appears to be due to the model structure rather than the model inputs. Some sub-catchments only respond to precipitation events after the catchment becomes sufficiently wet and saturated areas develop [18] and become connected to the stream network [46,47]. The former can be addressed by modification of the relationship between soil moisture and runoff with addition of one parameter as suggested by [18]; however, the issue of connectivity is related to dynamic changes in soil moisture patterns, which implies that the soil moisture-runoff relationship changes over time and this would be harder to incorporate into VIC [47].
Previous studies have reported that rainfall-runoff models calibrated using a long period of record and tested for sub-periods with above long-term average precipitation perform well, but the performance of the rainfall-runoff model starts to deteriorate when tested for sub-periods with below long-term mean rainfall in the same region of this study [43]. This observation was not consistently evident for VIC in this study. In fact some sub-catchments (4 of 13) showed better performance in term of NSE when model parameters calibrated for a wet period (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997) were used in the predominantly dry validation period (1998-2012) (see Table 4). This might be due to using LAI, rather than LAI mean . It is clear that lower LAI during dry periods, due to drier moisture conditions, could further reduce actual evapotranspiration, while increased LAI during wet conditions might further increase actual evapotranspiration. Using LAI mean tends to underestimate (overestimate) runoff during winter in comparison with using the observed monthly LAI runoff in dry (wet) years (Fig. 5). The ability to allow monthly LAI to vary from year to year is lacking in most rainfall-runoff models but is possible in physicalbased hydrological models such as VIC. Previous findings [6,25,43] obtained with different catchment sets and models, emphasised the lack of robustness of conceptual rainfall-runoff models when calibrated during wet periods and validated in dry periods. In most rainfall-runoff models the inability to vary LAI of the vegetation might contribute to reduced model performance during drought due to not representing changes in LAI that we show improve model  performance in some sub-catchments when included. In support of this, Zhang et al. [49] showed an improvement in model performance when a rainfall-runoff model was coupled with actual evapotranspiration estimates from the Penman-Monteith equation and remotely sensed MODIS LAI.
Based on Nash-Sutcliffe efficiency the VIC hydrological models tested here revealed better performance when driven with LAI than LAI mean . Sub-catchments located in high annual precipitation zones that are covered predominantly with trees showed less than 3% degradation in model performance since those areas are more likely to be energy limited (Fig. 6c). In this case the actual evapotranspiration is triggered by energy availability rather than the leaf area index. Sub-catchments numbered 11-, where the dominant land cover type is pasture or sparsely forested sub-catchments, showed more than 3% degradation in NSE model performance criterion, which might be related to LAI mean not representing year to year variability in leaf area index due to fluctuations in climate [41] and change in phenological cycles or timing of planting or harvesting. The degradation in model efficiency from using LAI mean impacts soil moisture simulation in the VIC model as shown by Ford and Quiring [14].
A runoff sensitivity analysis was conducted by changing mean monthly LAI. The results indicated that decreasing mean monthly LAI has more effect on runoff than increasing mean monthly LAI consistently across all calibrated sub-catchments. The differences in the rate of response of mean annual runoff to increases and decreases in mean monthly LAI are related to the differences in above ground processes like throughfall, interception and canopy evaporation. Evapotranspiration and total soil moisture content simulated by VIC also showed sensitivity to the input LAI data type. In low precipitation pasture dominated catchments (Catchments 11 and 13) the model showed relatively higher sensitivity of evapotranspiration and total soil moisture content to the input LAI data type (a) Areal averaged sub-catchment monthly precipitation, (b) change in areal average sub-catchment monthly simulated evapotranspiration between the two forms of LAI data inputs (E LAI -E LAImean ), and (c) change in areal average sub-catchment simulated monthly total soil moisture content between the two forms of LAI data inputs (SM LAI -SM LAImean ) given separately for small, medium and larger impacted sub-catchments for the period 1982-2012.
( Fig. 8b and c). Whereas, Catchment 5 with relatively high annual precipitation and covered mainly with tree, showed less evapotranspiration and total soil moisture content sensitivity to the input LAI data type (Fig. 8b and c). For the same amount of precipitation more water reaches the ground under lower LAI than higher LAI since interception storage is directly related by LAI, which increases surface runoff. When climate-induced changes in LAI were represented in the hydrological modeling, reduction of LAI due to decline in precipitation decreases the evapotranspiration from vegetation that made the soil wet [42]. This effect of LAI on runoff is very important to model especially during prolonged drought when precipitation is low and unrealistic LAI input can result in unrealistic soil moisture status, with consequent impacts on runoff. The sensitivity of the model to mean monthly LAI was found to depend on the climatic conditions and vegetation type of the sub-catchment. The lower the mean annual precipitation the sub-catchment received, the higher the sensitivity of the mean annual runoff to change in the mean monthly LAI and vice versa. And also for a given mean annual precipitation, the lower the proportion of trees to pasture the higher the sensitivity of runoff to change in the mean monthly LAI and vice versa. Thus the effect of land cover change on mean annual runoff varies across sub-catchments with similar mean annual precipitation.

Summary and conclusion
The three layer Variable Infiltration Capacity (VIC) model was calibrated for 14 gauged selected sub-catchments located in the Goulburn-Broken catchment, south-eastern Australia. Two sets of experiments were conducted to assess the effect of using different LAI inputs on the VIC model performance and the simulated monthly runoff. In addition the impact of catchment characteristics including vegetation cover type and mean precipitation on the sensitivity of catchment runoff to changes in LAI was assessed. The most notable findings are: 1. The VIC model simulated runoff reasonably well with high Nash-Sutcliffe model efficiency in the Goulburn-Broken catchment with proper calibration of the seven sensitive model parameters. 2. For sub-catchments predominantly covered in pasture or crop, Nash-Sutcliffe model efficiency was improved in a range from 4% to 25% (Table 4) when VIC model was calibrated using observed monthly LAI instead of mean monthly LAI. Calibrating VIC model using observed monthly LAI showed less than 4% (Table 4) improvement in Nash-Sutcliffe model efficiency for subcatchments predominantly covered by trees. This implies that calibrating the model using observed monthly LAI is important if the catchment is dominated by pasture cover. However, in tree dominated sub-catchments using either mean monthly LAI or observed monthly LAI can give the same VIC model performance. 3. Applying the long term mean monthly LAI overestimated monthly runoff by up to 25 mm during wet periods and underestimated monthly runoff by up to 35 mm during dry periods when compared with the runoff simulated using the year to year variable monthly LAI. 4. The difference in spatial patterns of the effect of observed monthly LAI over the mean monthly LAI on runoff is most strongly related to differences in vegetation type but is also influenced by mean annual precipitation differences.
Making accurate runoff predictions using hydrological models under changing or 'non-stationary' conditions is a challenge in hydrology. It seems difficult to provide general guidelines for which model should be used for which purpose but our results do indicate the importance of accounting for potential vegetation changes in climate change impact studies. That is models that consider both first order forcing, precipitation and potential evaporation drivers and the second order forcing of vegetation, could provide better predictions than models that only consider the first order forcing for projecting runoff with projected climate inputs. This implies a need to predict vegetation changes in response to climate change in studies which assess the impact of future climate change on runoff.