Comparing hydrological modelling , linear and multilevel 1 regression approaches for predicting baseflow index for 596 2 catchments across Australia 3

Junlong Zhang, Yongqiang Zhang, Jinxi Song, Lei Cheng, Rong Gan, 4 Xiaogang Shi, Zhongkui Luo, Panpan Zhao 5 1 CSIRO Land and Water, GPO Box 1700, ACTON 2601, Canberra, Australia 6 2 Shaanxi Key Laboratory of Earth Surface System and Environmental Carrying Capacity, 7 College of Urban and Environmental Sciences, Northwest University, Xi'an 710127, China 8 3 State Key Laboratory of Soil Erosion and Dryland Farming on the Loess Plateau, Institute 9 of Soil and Water Conservation, Chinese Academy of Sciences, Yangling 712100, China 10 4 Lancaster Environment Centre, Lancaster University, Lancaster, UK, LA1 4YQ 11 5 CSIRO Agriculture Flagship, GPO Box 1666, ACTON 2601, Canberra, Australia 12 6 State Key Laboratory of Hydrology-Water Resource and Hydraulic Engineering, College of 13 Hydrology and Water Resources, Hohai University, Nanjing 210098, China 14


Introduction
Baseflow, the outflow from the upstream aquifers when the recharge is ceased (e.g., precipitation or other artificial water supplies) (Brutsaert and Lopez, 1998;Brutsaert, 2005), is an important indicator of catchment hydrogeological characteristic (Knisel, 1963).
Various methods have been developed to separate baseflow from streamflow (Lyne and Hollick, 1979;Rice and Hornberger, 1998;Spongberg, 2000;Furey and Gupta, 2001;Eckhardt, 2005;Tularam and Ilahee, 2008;Lott and Stewart, 2016), which can be categorized to tracer based and non-tracer methods (Gonzales et al., 2009).However, tracer based method is only applied to experimental catchments due to expensive the high consumption of both experimental time and materials (Koskelo et al., 2012).The alternative is non-tracer methods (e.g., digital filter methods) (Zhang et al., 2017), which are widely used because of their high efficiency and repeatability in estimating BFI (Arnold et al., 1995).More importantly, they perform well when the digital-filtering parameters (e.g., recession constant and maximum baseflow index) are appropriately estimated (Zhang et al., 2017).The non-tracer methods can only be used for catchments with streamflow observations.For ungauged catchments, hydrological models and regression approaches can be used to separate baseflow form total streamflow.Their accuracy can be evaluated against ensemble estimates from the non-tracer methods at gauged catchments.
One group considers baseflow as a linear recession process for groundewater reservoir, including SIMHYD (simplified version of the HYDROLOG model) (Chiew and McMahon, 1994;Zhang et al., 2016), 1LBY (Abdulla et al., 1999;Stoelzle et al., 2015), HBV (Ferket et al., 2010) models; the another group takes baseflow as a non-linear recession process including Xinanajing (Zhang and Chiew, 2009), PDM (Ferket et al., 2010) and ARNO (Abdulla et al., 1999) models.It is expected that BFI obtained from the hydrological models is largely uncertain as a result of different model structures, model calibration and parameterisation schemes (Beven and Freer, 2001).There are few studies in the literatures to evaluate the accuracy of baseflow estimation from the hydrological models at a regional scale.This study evaluates two hydrological models (SIMHYD and Xinanjiang models) for predicting BFI against the ensemble BFI estimates from the non-tracer methods.
Linear regression approach is another commonly used method to predict hydrological signature indices, including baseflow index (Gallart et al., 2007;Longobardi and Villani, 2008;Bloomfield et al., 2009;van Dijk et al., 2013).This method uses catchment physical characteristics (i.e.descriptors) and BFI obtained from the gauged catchments to establish linear regressions that are then used to predict BFI in ungauged catchments (Bloomfield et al., 2009;Beck et al., 2013).Several studies show some catchment characteristics have important control on BFI.For instance, geological characteristics such as soil properties were found to be key for accurate BFI estimates (Brandes et al., 2005;van Dijk, 2010).Other studies also used climate-related indices, such as mean annual precipitation and mean annual potential evaporation, to simulate BFI (van Dijk, 2010;Beck et al., 2013).In similar studies, mean annual precipitation, slope and proportion of grassland are used for building the regressions for predicting BFI (Haberlandt et al., 2001;Brandes et al., 2005;Mazvimavi et al., Hydrol. Earth Syst. Sci. Discuss., https 2005;Gebert et al., 2007;Bloomfield et al., 2009;van Dijk, 2010).Beside BFI, linear regression is also an useful approach in estimating other hydrological signatures (e.g., runoff coefficient, runoff seasonality, zero flow ratio and concavity index (Zhang et al., 2014)) and understanding the catchment hydrology behaviour (Zhang et al., 2014;Su et al., 2016).One limitation of the linear regression approach is that it uses constant parameters to predict BFI, and cannot handle cross-interactions at different spatial scales (Qian et al., 2010), which could result in large errors for catchments located in a wide range of climate regimes.This limitation can be overcome by the multilevel regression approach that provides a robust tool to establish the relationships between BFI and catchment attributes.The basic idea of this approach is that higher level variables vary within a lower level (Berk and De Leeuw, 2006).This approach can handle the variables with various solutions using random effects (i.e., hierarchical structure) (Dudaniec et al., 2013).This approach has been extensively used to understand interplay of ecosystem dynamics (i.e., carbon cycle across different ecosystem (McMahon and Diez, 2007;Luo et al., 2015) and N2O emissions from agricultural soils (Carey, 2007)).However, no literatures have been reported to use this approach for hydrological signature (such as BFI) predictions.This study, for the first time, explores the possibility of using multilevel regression (Qian et al., 2010;Luo et al., 2015) to predict BFI across widely distributed Australian catchments.Catchment characteristics are used here as lower level (i.e., individual-level) predictors, and the effect of these predictors is assumed to vary across higher level predictors (i.e., climate zones) (Gelman and Hill, 2006).Details of the multilevel regression approach are elaborated in section 3.3.
The main aim of this study is to improve the large-scale BFI prediction.To achieve this, we compare the three BFI prediction methods (hydrological modelling, classic and multilevel regression approaches) against ensemble average estimates from four non-tracer baseflow separation methods.The objectives of this study are to i. Obtain "benchmark" BFI using the four non-tracer baseflow methods (Lyne-Hollick, UKIH (United Kingdom Institute of Hydrology), Chapman-Maxwell and Eckhardt) for 596 Australian catchments (Figure 1); ii.Introduce the multilevel regression approach for FBI predictions across large regions; iii.Assess relative merits of the three approaches for BFI predictions; and iv.
Investigate good BFI predictors for the multilevel regression approach.
Figure 1 is about here 2 Data sources

Streamflow
There are 596 catchments selected across Australia for assessing the three methods (hydrological modelling, linear regression and multilevel regression) used in this study to predict BFI.Streamflow measurements and related catchment attributes were collated by Zhang et al. (2013).Following criteria are used to filter the streamflow data for each catchment: i.It is a small catchment with catchment area 50 to 5000 km 2 ; ii. Streamflow was not subject to dam or reservoir regulations; iii.Except for the climate forcing data, the two models also require remote sensing leaf area index, land cover and albedo data that were used to calculate actual evapotranspiration (ETa) using the Penman-Monteith-Leuning model (Leuning et al., 2009;Zhang et al., 2010).The leaf area index data from 1981 to 2011, derived from the Advanced Very High Resolution Radiometer (AVHRR), were obtained from Boston University (Zhu et al., 2013).The temporal resolution is half-monthly and its spatial resolution is ~8 km.The land cover data required to estimate aerodynamic conductance came from the 2000-2001 MODIS land cover product, obtained from the Oak Ridge National Laboratory Distributed Active Archive Center (Friedl et al., 2010).The dataset has 17 vegetation classes, which are defined according to the International Geosphere-Biosphere Programme.The albedo data required to calculate net radiation were obtained from the 8-day MODIS MCD43B bidirectional reflectance distribution function product at 1 km resolution.All of the forcing data were reprojected and resampled using nearest neighbour approach to obtain 0.05 o gridded data.

Baseflow separation algorithm
The benchmark BFI data were estimated using four baseflow separation methods.They are Lyne-Hollick (Lyne and Hollick, 1979), UKIH (Gustard et al., 1992), Chapman-Maxwell (Chapman and Maxwell, 1996) and Eckhardt (Eckhardt, 2005) respectively.It is found that estimates of the recession constant and maximum baseflow index are the key to improve the performance of the digital-filtering methods (Zhang et al., 2017).This study used the Automatic Baseflow Identification Technique (ABIT) for the recession analysis, which was developed by Cheng et al. (2016) based on the recession theory provided by Brutsaert and Nieber (1977).Figure 2 demonstrates how the recession constant is estimated using the ABIT method.
In order to eliminate uncertainties raised from different algorithms, the ensemble mean from the four methods was taken as the benchmark (denoted as 'the observed BFI').The observed BFI was used either to evaluate the two hydrological models for BFI prediction, or to build the linear and multilevel regression approaches together with the catchment attributes.
Figure 2 is about here

Hydrological models
The SIMHYD and Xinanjiang model are two conceptual rainfall-runoff hydrological models.
Since developed by Chiew and McMahon (2002), SIMHYD has been widely applied in runoff simulation and regionalization studies (Chiew et al., 2009;Vaze and Teng, 2011;Li and Zhang, 2016;Zhang et al., 2016).Four water stores are used in this model to describe hydrological processes, namely the interception store, soil moisture store, groundwater store and channel store (Chiew and McMahon, 2002).Detailed model structure can be found in Chiew and McMahon (1994).The modified SIMHYD model by Zhang and Chiew (2009), which uses remote sensing data and contains nine model parameters, is used in this study.
The Xinanjiang model was developed by Zhao (1992) and has been widely used in humid and semi-humid regions (Li et al., 2009;Lü et al., 2013;Yao et al., 2014).This model reproduces runoff by describing three hydrological processes including ETa, runoff generation, and runoff routing.Details of Xinanjiang model are available from studies conducted by Zhao (1992) and Zhang and Chiew (2009).Here we use the modified Xinanjiang model proposed by Zhang and Chiew (2009), in which ETa was estimated using remote sensed LAI and the model parameters were reduced from 14 to 12.
The revised version of those two models is denoted as original models.The details of two hydrological models and regionalization approaches are described by Zhang and Chiew (2009).We used three types of BFI estimates from hydrological modelling: calibration, Efficiency of the daily square-root-transformed runoff data and minimise the model bias (Li and Zhang, 2017).
For the spatial cross-validations, two regionalisation approaches, spatial proximity and integrated similarity approaches (Zhang and Chiew, 2009) were used.The spatial proximity approach is where the geographically closest catchment is used as the donor basin to predict the ungauged catchments; integrated similarity approach is derived from combination of the spatial proximity and physical similarity approaches.

Linear regression and multilevel regression approaches
Traditionally, BFI was predicted using one set parameters for all catchments.The details are: where BFIi is the baseflow index for each catchment i=1,..., 596, N is normal distribution function, α is the intercept, β is slop, X is the variables (i.e., catchment attributes), and ε is 2006; Longobardi and Villani, 2008;Bloomfield et al., 2009).However, catchment attributes vary with hydrometerological conditions, therefore the constant parameters are not adequate to reflect the catchment characteristics.This approach ignored variability of catchment characteristics in various backgrounds.In order to reduce the uncertainties of prediction using one set of parameters, one level reflects hydrological background should be introduced.
In this study, we assumed that BFI associates with the climate variables (annual precipitation, potential evapotranspiration) and terrain attributes (area, elevation, slope, land cover and available soil water holding capacity in top soil) in each catchment (i.e., i = 1, 2, 3, …, 596).
We further assumed that the effects of those predictors on BFI vary with climate zones including arid, tropics, equiseasonal and winter rainfall (i.e., j = 1, 2, 3, 4).In this process, the catchments were divided into multiple datasets based on climate zones, then individual linear regression model were built for each subset.
where j is catchment in each climate zone, BFIji is the baseflow index for catchment in each subset j = 1,2,3,4.N is normal distribution function, α is the intercept, β is slop, X is the variables (i.e., catchment attributes), and ε is variance in each subset.However, hydrological processes in a catchment have close connections with other catchments, interactions crossing various group levels are primary drivers to influence baseflow processes.Therefore, an approach should be developed to consider cross level effects for predicting hydrological signatures.
Thus, we introduced the multilevel regression approach (Gelman and Hill, 2006;Qian et al., 2010;Luo et al., 2015) to improve the prediction of BFI and quantify the relative importance of predictors under different climate zones.Comparing the traditional linear regression approach, the hierarchical structure of the multilevel regression approach allows the  , 2, 3, 4).In this model, the intercept and slope vary with the group level (i.g., climate zone).The details of the approach is elaborated as follows: where Xi is the catchment attributes for each basin, and its intercepts and slopes can be decomposed into terms for climate zone, where µα and σα are the mean and standard deviation of variable intercept α, µβ and σβ are the mean and standard deviation of variable slope β, ρ is the correlation coefficients between the two variables αj and βj.The Eq. ( 3) can be rearranged as block matrix of the details of Eq. ( 5) ) , ( ~  N A can be described as: the Eq. ( 4) can be calculated individually by: The density function of the normal distribution N is (for example, α variable): This model considers variation in the αj's and the βj's and also a between-group correlation parameter ρ (Gelman and Hill, 2006;Qian et al., 2010).In essence, there is a separate regression model for each climate zone with the coefficients estimated by the weighted average of pooled (which do not consider groups) and unpooled (which consider each group separately) estimates, i.e. partial pooling.When fitting the model, all predictors are standardized using z-scores.
Where x' is the new catchment attributes using function z-scores.

Leave-one-out cross-validations
We apply leave-one-out cross-validation to assess the ability of the two regression approaches to predict BFI in 'ungauged' catchments where no streamflow data are available.
For each of the 596 catchments, the data from other 595 catchments are used to predict its BFI.This procedure is repeated over all 596 catchments.This cross-validation procedure explores the transferability of the two regression approaches from known catchments to the ungauged and particularly evaluates the value of the between-catchments information.
4 Model evaluation

Bias
The absolute percentage bias was used to evaluate model performance, which is calculated as: where BFIo is the observed BFI derived using the ensemble average from the four non-tracer baseflow separation approaches (i.e., Lyne-Hollick, UKIH, Chapman-Maxwell and Eckhardt), BFIs is the simulated BFI from the two hydrological models or the two regression approaches.And n is the total number of catchment.The unit of bias is a percentage (%), the larger of the absolute bias, the worse of the simulation.The bias is 0 indicates that simulated value is the same as the observed value.The Nash-Sutcliffe efficiency (NSE) is a normalized statistic that measures the relative magnitude of the residual variance ("noise") compared to the measured data variance ("information") (Nash and Sutcliffe, 1970).It is a classic statistical metrics used for evaluating the model performance.The closer NSE is to 1.0, the better the simulation is.

Spatiality of observed BFI
It can be seen from Figure 3 that BFI varies dramatically across Australia (location, i.e. coordination and distance away from ocean).Within the latitude ranges from 20º S to 30º S, which is smaller than that of the regions beyond this latitude range.Catchments located in latitudes higher than 30º S tend to have larger BFIs in general.Yet this is not the case for  Figure 4 is about here We further compared the observed and simulated in scatterplots (Figure 5). Figure 5 is about here

Comparison of traditional regression and multilevel regression approaches
Figure 6 compares the observed BFIs and simulated BFIs using traditional linear multivariate regression and multilevel regression approaches across four different climate zones.The result shows that the multilevel regression approach generally outperforms the traditional linear regression approach, evidenced by the NSE from multilevel regression approach being 0.31, 0.30, and 0.18 higher than that from linear regression in arid, tropics, and equiseasonal regimes respecitively, and the percentage bias from multilevel regression approach being 8, 7, and 8 lower than that from the linear regression.The two approaches show no significant difference in winter rainfall climate zone, indicated by same bias or NSE.
Figure 6 is about here We further check the leave-one-out cross-validation results obtained from the two approaches (Figure 7).It is clear that there exists noticeable degradation from calibration to cross validations for the traditional regression in the three climate zones: arid, tropics, and equiseasonal regimes.Compared to that, there is no noticeable degradation for the multilevel regression approach for the three climate zones.In the winter rainfall zone, the both approaches do not have degradation, and perform similarly.The leave-one-out cross- validation results further demonstrate the multilevel regression approach outperforms the traditional linear regression.
Figure 7 is about here Figure 8 summarises parameters of the multilevel regression approach.It can be seen that precipitation has the most positive impact on BFI, which does not greatly vary across climate zones.ETP has the most negative effect among all climate zones, and has significant large effect in equiseasonal zone.The H and Kst also have the noticeable positive effect on all the climate zones.Other three characteristics A, S and F have slope close to zero, suggesting small impacts on BFI.
Figure 8 is about here

Discussion
Our results suggest there are large biases to use hydrological models to simulate and predict BFI.It is understandable since hydrological models are not designed to simulate baseflow directly, but the baseflow component, in order to better simulate streamflow.It seems that model structure is more important than parameterisation since the three parameterisation schemes (calibration, spatial proximity and integrated similarity) obtain similar BFI, and SIMHYD has larger bias than Xinanjiang model as summarised in Figures 4 and 5.However, both hydrological models are calibrated against total streamflow, rather than its components, such as baseflow.This suggests that better estimate streamflow.This issue has been well recognised in other hydrological models as well (Fenicia et al., 2007;Lo et al., 2008).In fact, baseflow is designed as an integrated store combined with the river recharge (Chiew and McMahon, 2002).This structure feature of hydrological models tends to overestimate baseflow and therefore leads to unsatisfactory BFI prediction.Interactions of catchments crossing group level would influence the baseflow processes.BFI is affected by catchment attributes, and in relevance with terrain and climate factors (Gustard and Irving, 1994;Longobardi and Villani, 2008;van Dijk, 2010;Price, 2011).However, how to predict the effect of BFI response to such various environmental conditions remains challenging.In order to improve our understanding of BFI, interaction of catchment attributes within different climate zones should be considered (Berk and De Leeuw, 2006).
Climate influences the hydrological process and thus leads to changes in baseflow generation.
Implementation of multilevel regression approach in this study, P and ETP have the most significant effects on BFI, are the two essential elements controlling baseflow processes.The effect of these two factors varies across climate zones.As studies by Santhi et al. (2008) and Peña-Arancibia et al. (2010), they have shown that climate attributes can be used to best predictors for recession constant.The increase of the precipitation can cause the more saturation of the soil, and lead to the baseflow increase (Mwakalila et al., 2002;Abebe and Foerch, 2006).In addition, the ETP is related to the baseflow discharge over the extended period (Wittenberg and Sivapalan, 1999).ETP has the adverse effect on BFI for all climate zones.This result agrees well with the conclusion drawn by Mwakalila et al. (2002).The influence is relative smaller in arid zone than other climate zones.In general, ETP is related to the baseflow discharge over the extended period (Wittenberg and Sivapalan, 1999), catchment with low evapotranspiration will have higher BFIs (Mwakalila et al., 2002).
Comparing to climate attributes, F tends to have smaller effects and with various effects with climate zones (i.e., positive effect in arid and winter rainfall zones).F associates with quick flow generation and thus leads to the changes in the baseflow.The influence comes from vegetation regulation of water flux through moist conditions and ETP (Krakauer and Temimi, 2011).The plant on the ground can cover the land surface and influence the ETP and then increase the baseflow.Studies have shown that vegetation cover has a strong control on ET in catchments, and thus influences baseflow generation (Schilling and Libra, 2003).Wittenberg (2003) found that water consumption of deep-rooted vegetation has significant influence on baseflow generation where faster recession is usually found.Furthermore, baseflow is more closely related to the water storage of the saturated zone in plant root zone drainages (Milly, 1994).Studies have also shown that higher watershed forest cover usually corresponds well with lower BFI (Price, 2011).This is particularly significant during dry seasons, where the reduction of vegetation cover can lead to increase baseflow in dry seasons (Singh, 1968;Price, 2011).In the tropic zone, the proportion of the forest cover within a catchment has negative effect on BFI.This is because of the high water loss through ET in forests, and the vegetation draws heavily on the artesian leakage and contacts the spring flow (Meyboom, 1961;Knisel, 1963).Although a relatively close correlation between forest cover and BFI is found for most catchments, there are exceptions in some catchments.For instant, BFI was found to have a weak correlation with forest area in the Mediterranean region (Longobardi and Villani, 2008) and a case study in the Elbe River Basin (Haberlandt et al., 2001).
Our study demonstrates that those two topographic features are insignificant impact on the BFI cross Australia, and have different effects on various climate zones (i.e., slope has positive impact on arid but negative on other climate zones).However, some studies found that S and H have positive correlation with the recession timescales (Peña-Arancibia et al., 2010;Krakauer and Temimi, 2011).When interactions crossing level have been implemented, adding those two factors can greatly improve performance of multilevel regression approach.
Other studies show that the watershed area and slope are highly associated with the baseflow statistics (Vogel and Kroll, 1992).This can be a result of the catchments in their study are under the 150 km 2 .The effect of the slope will be induced when the catchment area are larger (Peña-Arancibia et al., 2010).However, the study conducted in southeaster Australia found that the topographic parameters have no significant relationship with the BFI (Lacey and  Grayson, 1998), this may be groundwater is relatively deep reducing connections between groundwater and streams (Mazvimavi et al., 2005).Besides, Kst is positively related with BFI for all catchment across climate zones.This may be explained by the strong interactions between soil water content and P as well as ETP (Milly, 1994).
Our result shows that multilevel regression approach, this approach can better understand the hydrological dynamics within different systems.To be specific, this method considers climate controls on catchment BFIs cross continental scale (Figure 8). Figure 9 shows the different coefficients in each climate zone.The hydrological processes are controlled by various climate conditions at large scale as has been proved by a number of studies (Lacey and Grayson, 1998;Abebe and Foerch, 2006;Merz and Blöschl, 2009;Ahiablame et al., 2013).
The baseflow processes will have the interactions at different climate zones (within and between group).The multilevel regression approach considers the cross-level interactions, and the prediction not only influenced by predictors at one scale (i.e., continental scale) but also different spatial scale (i.e., climate zones) (Qian et al., 2010), incorporates the group level information, and this approach takes the fixed and random effects into account one single model, the coefficients of the model for the whole data and the group has the variances.Prediction of BFI using group level information (i.e., climate zones) will help capturing the climate spatial variability at different regional scales.
According to the good performance as illustrated above, it is promising that this method can be used as a robust tool to estimate BFI across changing backgrounds (i.e., climate zones), and can promote improved understanding of hydrological processes.

Conclusion
This study estimated ensemble baseflow index from four well-parameterised baseflow including precipitation and potential evapotranspiration are proven to be most significant.
The multilevel regression approach can provide insights into the control factors of baseflow generation.This approach has the potential of being used to estimate baseflow index.We proposed the framework of using this approach to estimate hydrological signatures of under various backgrounds.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-737Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 21 December 2017 c Author(s) 2017.CC BY 4.0 License.regionalisation from spatial proximity, and regionalisation from integrated similarity.Herein, a short description of these three kinds of estimates is given below.For model calibration, a global optimisation method, the genetic algorithm from the global optimisation toolbox in MATLAB(MathWorks, 2006), was used to calibrate the model parameters for each catchment.This optimiser used 400 populations and the maximum generation of 100 for searching the optimum point, which converges at approximately 50 generations of searching.The model calibration was to maximise the Nash-Sutcliffe variance.This model ignores the potentially different effects of the same variable on BFI across different climatic zones.That is, α and β are constant irrespective of the climatic zone to which the BFI belongs.To be specific, many studies have conducted the baseflow prediction at large area, yet constant α and β are used in the model (Abebe and Foerch, Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-737Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 21 December 2017 c Author(s) 2017.CC BY 4.0 License.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-737Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 21 December 2017 c Author(s) 2017.CC BY 4.0 License.assessment of the variation in model coefficients across groups (e.g., climatic zones) and accounting for group-level variation in the uncertainty for individual level coefficients.The multilevel regression approach could be written as a data-level model (the predicted BFIi belonging to climate zone j), allowing the model coefficients (α and β) to vary by climate zone (j = 1 Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-737Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 21 December 2017 c Author(s) 2017.CC BY 4.0 License.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-737Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 21 December 2017 c Author(s) 2017.CC BY 4.0 License.Tasmania, where catchments with latitude higher than 40º S have smaller BFI values in the southeastern region within this island.This indicates that the BFI spatiality is distinct from the main continent to island.It is also interesting to notice that beyond the range of 20-30º S, observed BFI increases from inner land to coastal catchments, especially in southeast region within the main Australian continent.

Figure 3
Figure 3 is about here Figure 5(a) and 5(d) compares the observed and simulated BFIs from calibrated SIMHYD and Xinanjiang models, respectively.Figure 5(b)-(c) and 5(e)-(f) show the regionalisation results (i.e., spatial proximity and integrated similarity) of these two hydrological models.Notably, BFI estimated using SIMHYD model is much larger than the observed values (Figure 5(a), (b), and (c)), with the majority catchment BFIs dotted above the 1:1 line.SIMHYD model under calibration, spatial proximity, and integrated similarity gives NSE being -8.30, -8.42 and -8.44 respectively, and gives percentage bias being 146, 152 and 152 respectively, indicating similar poor model performance.In comparison, BFI estimated from Xinanjiang model tends to scatter a larger range around 1:1 line regardless of the parameterisation method (Figure Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-737Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 21 December 2017 c Author(s) 2017.CC BY 4.0 License.5(d), (e), and (f)), and is closer to the observed BFI.Xinanjiang model under calibration, spatial proximity, and integrated similarity give NSE being -2.75, -2.70 and -2.58 respectively, and gives bias being 84, 81 and 83 respectively, indicating still similar poor model performance in prediction of BFI.The results obtained from Figures 4 and 5 indicate that parameterisation has much smaller impact on BFI estimates, compared to model structure.

Figure 1 .
Figure 1.The location of 596 selected unregulated small catchments in this study and climate classification based on Köppen-Geiger (2006) classification schemes in Australia.

Figure 2 .
Figure 2. Estimation of the recession constant (Log (-dQ/dt) versus log (Q)) using automated baseflow identification technique (ABIT) for Endeavour catchment (station ID 107001).The black line is 5 % lower envelope line has a slope 0.983 and the estimate of the characteristic drainage time scale K = 57.1 days.

Figure 3 .
Figure 3. Spatial distribution of the observed baseflow index across Australia.

Figure 4 .
Figure 4. Baseflow index duration curves obtained from the observed, SIMHYD model and Xinanjiang model.Calibration and two regionalisation results are shown for each hydrological model, where R1 and R2 represent spatial proximity and integrated similarity approaches, respectively.SIMHYD is simplified version of the HYDROLOG model.

Figure 5 .
Figure 5. Scatterplots of observed baseflow index versus simulated baseflow index using SIMHYD and Xinanjiang models, where calibrated and regionalised model results are presented in (a) and (d) (calibration), (b) and (e) (spatial proximity regionalisation) and (c) and (f) (integrated similarity regionalisation), respectively.The blue ellipses represent the confidence level at 0.95.The full naming of SIMHYD is introduced in Figure 4.

Figure 6 .
Figure 6.Scatterplots of observed and simulated baseflow index using traditional linear regression ((a)-(d)) and multilevel regression ((e)-(f)) approaches that are built using the full catchment samples in four climate zones, with (a) and (e) for arid, (b) and (f) for tropics, (c) and (g) for equiseasonal and (d) and (h) for winter rainfall, respectively.The blue ellipse is drawn at 0.95 confidence level.The black line represents 1:1 line.

Figure 7 .
Figure7.As same as Figure6, but using the leave-one-out cross validation approach.

Figure 8 .
Figure 8. Parameter values using multilevel regression approach, fixed and random variables758

Figure 1 .
Figure 1.The location of 596 selected unregulated small catchments in this study and climate classification based on Köppen-Geiger (2006) classification schemes in Australia.

Figure 2 .
Figure 2. Estimation of the recession constant (Log (-dQ/dt) versus log (Q)) using automated baseflow identification technique (ABIT) for Endeavour catchment (station ID 107001).The black line is 5 % lower envelope line has a slope 0.983 and the estimate of the characteristic drainage time scale K = 57.1 days.

Figure 4 .
Figure 4. Baseflow index duration curves obtained from the observed, SIMHYD model and Xinanjiang model.Calibration and two regionalisation results are shown for each hydrological model, where R1 and R2 represent spatial proximity and integrated similarity approaches, respectively.SIMHYD is simplified version of the HYDROLOG model.

Figure 5 .
Figure 5. Scatterplots of observed baseflow index versus simulated baseflow index using SIMHYD and Xinanjiang models, where calibrated and regionalised model results are presented in (a) and (d) (calibration), (b) and (e) (spatial proximity regionalisation) and (c) and (f) (integrated similarity regionalisation), respectively.The blue ellipses represent the confidence level at 0.95.The full naming of SIMHYD is introduced in Figure 4.

Figure 6 .
Figure 6.Scatterplots of observed and simulated baseflow index using traditional linear regression ((a)-(d)) and multilevel regression ((e)-(f)) approaches that are built using the full catchment samples in four climate zones, with (a) and (e) for arid, (b) and (f) for tropics, (c) and (g) for equiseasonal and (d) and (h) for winter rainfall, respectively.The blue ellipse is drawn at 0.95 confidence level.The black line represents 1:1 line.

Figure 8 .
Figure 8. Parameter values using multilevel regression approach, fixed and random variables are represented.Error bar represents standard error of each parameter.The abbreviations of catchment attributes are introduced in Table1.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-737Manuscriptunder review for journal Hydrol.Earth Syst.Sci.Discussion started: 21 December 2017 c Author(s) 2017.CC BY 4.0 License.Chapman-Maxwell and Eckhardt), and found that the baseflow index varies significantly in corresponding to climate zones across Australian continent.Multilevel regression approach is introduced to improve BFI estimate for 596 catchments across Australia.BFI obtained from this new method is compared to that of traditional linear regression method and two hydrological models.When compared to observed BFIs, the multilevel regression approach outperforms both linear regression approach and hydrological models.Traditional linear regression approach fails to considerate the interactions across group levels.The two hydrological models have good performance for simulating runoff yet fail to separate baseflow.In contrast, the multilevel regression approach indicates that annual precipitation, potential evapotranspiration, elevation, land cover and available soil water holding capacity in top part of the soil all have strong control on catchment baseflow, where climate factor separation methods (Lyne-Hollick, UKIH (United Kingdom Institute of Hydrology),