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.