Methods

Site selection and data

A total of 13 hydrologic reference stations across Australia, focussing on stations with catchments areas that are less than 250km2 were selected from the Bureau of Meteorology (http://www.bom.gov.au/water/hrs/) (Figure 1 and Table 1). Using a basic land use data interpretation, we focussed on catchments that are essentially fully forested to minimise effects of landuse variation. By focussing on relatively small catchments, we also aimed to control for catchment size which influences storage and routing.
Rainfall and temperature stations were selected for each streamflow station based on proximity, data length and data continuity. Data for the selected rainfall and temperature sites for the period 1970-2010 were obtained from the Bureau of Meteorology. However, since climate stations were not necessarily close to the catchments, we also extracted the ANUCLIM 1km gridded rainfall product via the eMast NCI server (http://dapds00.nci.org.au/thredds/catalog/rr9/ANUClimate/catalog.html) for each of the catchment outlets.

Analysis of annual trends and elasticities following Chiew (2006)

As a first step annual summaries of the data were analysed, following Chiew (2006). Both the non-parametric approach for calculating the rainfall elasticity, εP, (Fu et al., 2007) and the modelling approach (Chiew, 2006) were used. The non-parametric approach calculates the elasticity as (Fu et al., 2007):
\(\epsilon_{p}=median\left(\frac{Q_{t}-\overset{\overline{}}{Q}}{P_{t}-\overset{\overline{}}{P}}\frac{\overset{\overline{}}{P}}{\overset{\overline{}}{Q}}\right)\)(1)
For the modelling approach, both GR4J (Perrin, Michel, & Andréassian, 2003) and SimHyd including the Muskingum routing routine (Chiew et al., 2009) were used on the daily data. The models were all calibrated on the full 40 year daily dataset (1970 – 2010). The calibration used shuffled complex evolution optimisation with 20 complexes, and this was repeated 10 times to capture model uncertainty due to parameter equifinality. The residuals of the modelling were calculated as the difference between observed and predicted values.
Following the original work (Chiew, 2006), the daily rainfall series (P ) was scaled by -15%, -10%, 0%, +10%, while the daily maximum temperature (MaxT ) series, rather than the monthly Potential Evapotranspiration (PET) as in Chiew (2006), was scaled by 0%, +5% and +10%. The resulting predicted daily flow series (Q ) were summarised to annual series and εP and the elasticity due to ET , εET, were calculated using regression (Chiew, 2006):
\(\delta Q_{\text{annual}}=\varepsilon_{P}\delta P_{\text{annual}}+\varepsilon_{\text{ET}}\delta\text{MaxT}_{\text{annual}}\)(2)
Where, MaxTannual is the average annual maximum daily temperature rather than the maximum annual daily temperature. The elasticity is a measure of amplification of the rainfall change in the streamflow, with elasticities greater than unity indicating amplification of the rainfall and ET signal (Chiew, 2006).

Mann-Kendall non-parametric trend test and linear trends

The presence of monotonic trends in weekly mean streamflow, rainfall and maximum temperature were initially assessed using non-parametric Mann-Kendall tests. The advantages of Mann-Kendall based tests are that there is no assumption that the data are normally distributed, nor does it assume that the trend is linear. However, there are two possible issues with simply using a Mann-Kendall test on the data. The first is outlined by Westra et al. (2012), which highlights that even in random data there is a probability of identifying a trend. As a result, we followed Westra et al. (2012) and ran a bootstrap with 500 samples to identify this error in the basic Mann-Kendall analysis. The second issue is outlined by Koutsoyiannis and Montanari (2007): A locally observed trend is not necessarily a true trend in the data due to scaling properties. Under the assumption of scaling, the Mann-Kendall trend can be further tested to identify if the Hurst coefficient is significant and the Mann-Kendall test can be corrected for this (Hamed, 2008), we indicate this test as the LTPMK test (long time period Mann-Kendall).
The daily data were de-seasonalised and summarised to weekly before running the Mann-Kendall and LTPMK tests using the package “deseasonalize” (McLeod & Gweon, 2013). Calculations for the standard Mann-Kendall and LTPMK tests were performed using the “Kendall” package (Mcleod, 2011) and “HKProcess” (Tyralis, 2016) in R (R Core Team, 2013), respectively.
In addition to analysing the data series, the trend in the residuals of the modelled series from the previous rainfall-runoff modelling step was also analysed using the same LTPMK test. If the modelled series deviates from the observed series for the post calibration period with a consistent trend, then this would be another indication of a trend in the data series. For this analysis, the residuals were summarised to weekly data to match the observed data. All analyses were repeated with both the station rainfall data and the gridded rainfall data.

Generalized additive mixed model (GAMM)

The basic starting assumption in this analysis is that streamflow (Q ), combining both surface runoff and groundwater inflows, in relatively small homogeneous catchments can be modelled as a simple water balance of rainfall (P ) and evapotranspiration (ET ):
\(Q=P-ET\) (3)
Due to the relatively small size of the catchments, storage effects based on the weekly data can be considered minor and the majority of the storage delays can be captured (in other words, residence times of water moving towards the stream are considered to be smaller than a week for the majority of the flow).
Generalized additive mixed models (GAMMs) were defined to model the weekly streamflow at each station as a function of rainfall and weekly mean maximum temperature. GAMMs may be considered as an extension of generalized linear mixed models whereby the predictors are additive, but can also be non-linear (Chen, 2000). The mixed model addition to the standard Generalized Additive Model (GAM), allows modelling of the autocorrelation in the residuals. Details of the statistical basis of GAMMs can be found in Chen (2000) and Lin and Zhang (1999).
Mean weekly streamflow and rainfall are indicated byQw and Pw , respectively. A smoothing function, s(), was applied to the rainfall term to account for non-linear responses in the rainfall-runoff relationship. Weekly mean streamflow data was used to remove any effects of filled in missing data by the Bureau of Meteorology, and to remove most of the stream routing effects.
The models were applied in several steps (Table 2). Initially, the linear trends in both rainfall and temperature were analysed using generalised least squares using an autoregressive order 1 (AR1) correlation structure fitted to the errors (ɛt ). As an added benefit, the AR1 error structure accounts for further storage effects in the catchment, including the effect of past rainfall.
Subsequently a GAMM with a smooth function of Pw , a linear trend and autoregressive order 1 (AR1) correlation structure fitted to the errors (ɛt ) was evaluated.
\(\text{Log}\left(Q_{w}\right)\sim s\left(P_{w}\right)+trend+\varepsilon_{t}\)(4)
Streamflow data were log transformed to stabilise the variance of the residuals. The trend in equation 3 can be interpreted as a change in the residual log (Qw ) over time. In order, to identify \(\frac{\text{δQ}}{\text{δt}}\) as a fractional change the trend needs to be back transformed using:
\(\frac{\text{δQ}}{\text{δt}}=\left(\exp\left(\text{trend}\right)-1\right)\). (5)
The first model (equation 4) removes the direct effect of weekly rainfall (Pw ) on the weekly streamflow (Qw ). Amplification can be defined as\(\frac{\text{δQ}}{\text{δP}}>1\), using the annual data. However, the trend in equation 4 can include other phenomena changing the rainfall-runoff response over time. As a second step, a GAMM with a smooth function of Pw and a smooth function of the interaction between maxTw andPw was fitted:
\(\text{Log}\left(Q_{w}\right)\sim s\left(P_{w}\right)+s\left(\text{maxT}_{w},P_{w}\right)+trend+\varepsilon_{t}\)(6)
The interaction between maxTw andPw is a proxy for the actual evapotranspiration. Equation 6 identifies the residual trends in the streamflow data after removal of the effects of rainfall and actual evapotranspiration. The difference in the trend in the second model (equation 6) and the first model (equation 4) highlights the impact of actual evapotranspiration effects on streamflow. The final trend remaining in the model therefore includes any “other” phenomena in the rainfall-runoff relationship. Examples of these could be changes due to tree physiology, CO2 fertilization or geomorphology and is the degree of amplification, i.e. the change in streamflow after changes in rainfall and evapotranspiration have been considered.
Using equation 6, a GAMM which removes the linear trend was also fitted, to test whether the assumption of a linear trend affected the results:
\(\text{Log}\left(Q_{w}\right)\sim s\left(P_{w}\right)+s\left(\text{maxT}_{w},P_{w}\right)+\varepsilon_{t}\)(7)
The residuals (εt ) of this final GAMM model (equation 7) were subsequently analysed for a trend using a LTPMK test.
The GAMM modelling was conducted using the “mgcv” package in R (Wood, 2006).
Given the complexity of all the methodology a flow chart of the different approaches and comparisons is presented in Figure 2.