How streamflow has changed across Australia since the 1950 s : evidence from the network of hydrologic reference stations

Streamflow variability and trends in Australia were investigated for 222 high-quality stream gauging stations having 30 years or more continuous unregulated streamflow records. Trend analysis identified seasonal, interannual and decadal variability, long-term monotonic trends and step changes in streamflow. Trends were determined for annual total flow, baseflow, seasonal flows, daily maximum flow and three quantiles of daily flow. A distinct pattern of spatial and temporal variation in streamflow was evident across different hydroclimatic regions in Australia. Most of the stations in southeastern Australia spread across New South Wales and Victoria showed a significant decreasing trend in annual streamflow, while increasing trends were retained within the northern part of the continent. No strong evidence of significant trend was observed for stations in the central region of Australia and northern Queensland. The findings from step change analysis demonstrated evidence of changes in hydrologic responses consistent with observed changes in climate over the past decades. For example, in the Murray–Darling Basin, 51 out of 75 stations were identified with step changes of significant reduction in annual streamflow during the middle to late 1990s, when relatively dry years were recorded across the area. Overall, the hydrologic reference stations (HRSs) serve as critically important gauges for streamflow monitoring and changes in long-term water availability inferred from observed datasets. A wealth of freely downloadable hydrologic data is provided at the HRS web portal including annual, seasonal, monthly and daily streamflow data, as well as trend analysis products and relevant site information.

stations will serve as 'living gauges' that record and detect changes in hydrologic responses to 93 long-term climate variability and other factors. 94 This paper presents a statistical analysis to detect changes or emerging trends across a range 95 of flow indicators, based on the daily flow data of 222 sites from the HRS network. The 96 objective of this study is to provide a nationwide assessment of the long-term trends in 97 observed streamflow data. Evaluation of past streamflow records and documenting recent 98 trends will be of benefit in anticipating potential changes in water availability and flood risks. 99 It is hoped that the findings from trend analysis presented in this paper will inform decision 100 makers on long-term water availability across different hydroclimatic regions, and be used for 101 water security planning within a risk assessment framework. 102 103 2 Site selection, data and methods 104

Hydrologic Reference Stations and data 105
The 222 Hydrologic Reference Stations (HRS) were selected from a preliminary list of 106 potential streamflow stations across Australia according to the HRS selection guideline (SKM 107 2010). These guidelines specified four criteria for identifying the high quality reference 108 stations, namely unregulated catchments with minimal land use change, a long period of 109 record (greater than 30 years) of high quality streamflow observations, spatial 110 representativeness of all hydro-climate regions, and the importance of site as assessed by 111 stakeholders. The station selection guidelines were then applied in four phases 112 (www.bom.gov.au/water/hrs/guidelines.shtml). The HRS network will be reviewed and 113 updated every two years to ensure that the high quality of the streamflow reference stations is 114 maintained. 115 Two features were considered in order to define the hydroclimatic regions in HRS: climatic 116 zones and Australia's drainage divisions. The climatic zones were defined according to 117 climate classification of Australia based on a modified Koeppen classification system (Stern 118 et al., 2000). Australia has a wide range of climate zones, from the tropical regions of the 119 north, through the arid expanses of the interior, to the temperate regions of the south (ABS 120 2012). The Australian Hydrological Geospatial Fabric (Geofabric) Surface Catchments (BOM 121 2012) were used to delineate 12 topographically defined drainage divisions approximating the 122 drainage basins from the Geoscience Australia (2004) definition. The selection of HRS 123 Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2015-464, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. The primary data used in this study were daily streamflow series of 222 gauging stations from 133 the HRS network. Table 1 lists the twelve drainage divisions and the number of stations in 134 each division. One third of the HRS stations are located within the Murray-Darling basin, half 135 of the rest are distributed along eastern coasts. This is the best compiled long-term quality 136 controlled data for Australia and the trends derived from this dataset constitute the first such 137 statement on long-term water availability across Australia. 138 The earliest record included in the data set is from 1950. Data prior to this has been excluded 139 due to the common existence of large gaps in the pre-1950 period. All stations included in the 140 HRS had a target of 5% or less missing data to meet the completeness criteria for high quality 141 streamflow records. Some stations were included with more than 5% missing data where they 142 excelled in other criteria such as stakeholder importance or spatial coverage. The periods of 143 data gaps were filled using a lumped rainfall-runoff model GR4J (Perrin et al. 2003), which 144 was found to perform well at most sites. The model was calibrated and forced with catchment 145 average rainfall and potential evapotranspiration from the Australian Water Availability 146 Project (AWAP) (Raupach et al., 2009). 147 The study examined sites with varying lengths of record depending on the data availability. 148 The daily flow data were aggregated into annual series based on a water year calculation. The 149 start month of the water year was defined as the month with the lowest monthly flow across 150 the available data period. In order to ensure the statistical validity of the trend analysis, all 151 stations had minimum 30 years of record, with an average time-series length of 45 years. The 152 longest record length was 62 years, 25% of the stations have 50 or more years of record. 153 Catchment sizes ranged from 4.5 to 232,846 km 2 with a mean size of 3108 km 2 . The majority 154 Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2015-464, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 18 January 2016 c Author(s) 2016. CC-BY 3.0 License.
(82%) of the stations had an upstream drainage area less than 1000 km 2 , and only three 155 stations had a drainage area larger than 50,000 km 2 . 156 The data and the long term series gathered in this study are the best compiled and quality 157 assured data for HRS catchments. The analysis and trends derived from the HRS datasets 158 constitute the first statement on long-term water availability across Australia. 159

Streamflow variables for trend analysis 160
Long-term climate variability can be reflected through trends in streamflow variables. To 161 understand the importance of the components of the hydrologic regimes and their potential 162 link to long-term climate variability, ten streamflow variables were chosen for statistical and 163 trend analysis. Two variables related to fluctuation of annual flows were annual total flow 164 (Q T ) and annual baseflow (Q BF ). Baseflow was separated from daily total streamflow using a 165 digital filter based on theory developed by Lyne and Hollick (1979) and applied by Nathan 166 and McMahon (1990). The trend analysis was applied to the ten hydrologic indicators of streamflow data at each 178 HRS station. 179

Trend and data statistical analysis 180
Changes in streamflow data can occur gradually or abruptly. Statistical significance testing is 181 commonly used to assess the changes in hydrological datasets (Helsel and Hirsch, 2002;182 Monk et al., 2011;Hannaford and Buys, 2012). The Mann-Kendall (MK) trend test (Mann,183 Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2015-464, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci.  1945;Kendall, 1975) was adopted in this study to identify statistically significant monotonic 184 increasing or decreasing trends (Petrone et al., 2010;Zhang et al., 2010;Miller and Piechota, 185 2008). In order to ensure the assumption of independence was met for the MK test, the non-186 parametric Median Crossing and Rank Difference tests (Kundzewicz and Robson, 2000) were 187 applied to entire datasets. When either of the randomness tests indicated that the time series 188 was not from a random process, the site was excluded from the MK trend assessment. As this 189 study attempted to examine patterns in historical streamflow records, no further adjustments 190 were made to account for the non-random structure of data. 191 The non-parametric MK trend test was used to detect the direction and significance of the 192 monotonic trend, and the trend line from a least squares regression (LSR) was used to 193 approximately represent the magnitude of the trend. The trend magnitude was standardised 194 [trend (mm/yr) / average annual flow (mm)] to make the change comparable across stations 195 by dividing the regression slope coefficient by the average annual flow over the data period. 196 All data were subject to step change analysis to detect any abrupt changes during the record 197 period. The distribution free CUSUM test (Chiew and Siriwardena, 2005) was applied to 198 identify the year of change in streamflow series. The significant difference between the 199 median of the streamflow series before and after the year of change was tested by Rank-Sum 200 method (Zhang et al., 2010;Miller and Piechota, 2008;Chiew and Siriwardena, 2005). More 201 information on the statistical tests used in this study can be found in Appendix A. 202 In addition to the trend analysis for the ten flow indicators, other statistical data analyses were 203 included to gain a broad understanding of hydrologic regimes. Aggregated monthly and 204 seasonal flow data were investigated for changes in flow patterns in different basins or 205 regions. Daily event frequency analyses were used to examine the variations in daily 206 streamflow magnitude, and daily flow duration curves were presented to examine changes in 207 daily flow among decades. 208 209

Development of the HRS web portal 210
A web portal has been developed to house the network of Hydrologic Reference Stations and 211 provide access to streamflow data, results of analysis, and associated site information. Figure  212 2 summarises the development process of the HRS network and website. Through a data 213 quality assurance process following the guidelines and stakeholder consultations, the final list 214 Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2015-464, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 18 January 2016 c Author(s) 2016. CC-BY 3.0 License. of 222 streamflow gauging stations was established. A suite of software tools, "the HRS 215 toolkit" was developed to undertake data aggregation, analysis, trend testing, visualisation and 216 manipulation. The toolkit is capable of automatically converting the flow variables to 217 monthly, seasonal and annual totals, and quantifying the step and/or linear changes in the 218 selected streamflow variables. The toolkit also generated and processed graphical products, 219 data, statistical summary tables and statistical metadata included in the web portal. 220 A snapshot of the HRS web portal is shown in Figure 3. The main page was designed with 221 three parts. A series of links on the top provide the project information. Below this is the 222 station selector, which facilitates searching for the site of interest by location. The third part is 223 the product selector containing the core information sections of the website. Several tabs are 224 offered for users to explore the web portal dependent on their needs and the level of 225 information they require. The daily streamflow data, graphical products, statistics and trend 226 analysis results are available for users to view and download. Information provided on the 227 HRS web portal will assist in detecting long-term streamflow variability and changes at the 228 222 sites, and therefore supports water planning and decision-making. More information can 229 be found at the website http://www.bom.gov.au/water/hrs. 230 This web portal provides public access to high quality data and information. It has more than 231 15,000 graphic products for display. It is carefully designed for the public to have synthesised 232 and easily understandable information on water availability trends across Australia. In order 233 to ensure currency of this web site, streamflow data are updated and reviewed every two 234 years. 235 236 4 Result and discussion 237 The study to detect long-term streamflow trends was performed on the 222 gauging stations 238 included in the HRS network. This section presents an overview bar-plot of the  test results for the selected ten hydrologic variables. Maps showing trend detection results and 240 step change analysis for the annual total flow are presented as well as a table listing the 241 stations with significant trends in annual total flow at 1% significance level ( < 0.01). In 242 addition, variations in trend among daily flow indicators and seasonal flows are examined. 243 Finally, regional patterns in long term trends, inter-annual and decadal variability are further 244 investigated for two feature stations. 245 Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2015-464, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 18 January 2016 c Author(s) 2016. CC-BY 3.0 License.

Overview 246
A stacked bar-plot is shown in Figure 4  The abrupt changes in a hydrologic time series could be due to hydrologic non-stationarity. Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2015-464, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. it was for annual total flow. However, some time-series of annual baseflow were non-random 283 and therefore not available for further trend testing. 284 Step change analysis was applied to all sites where the time series data was random to give 285 comparable results of gradual and abrupt changes in annual total flows. Table 2 gives the 286 Rank-Sum test results and lists the year of change for the 22 stations. Details of step changes 287 across Australia will be discussed in the following section. 288 Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2015-464, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci.

289
The Rank-Sum test was used to identify the presence of a step change in the median of two 290 periods, with the distribution free CUSUM method providing the year of change. Values were 291 reported for sites with Rank-Sum test at 0.1 significance levels or higher. Figure 6  in south-west West Australia had a key feature of 1975 step change, which might be partly 302 due to the observed rainfall decline since the mid-1970s. It was also noted that most stations 303 located on the south east coast of Queensland showed a significant step change in the 1980s. 304 The results from step change analyse imply that changes in streamflow and the consequent 305 hydrologic response are driven by changes in climatic forcing such as rainfall over the period 306 of record. Investigating this causative relationship and quantifying the relative impacts of 307 variations in climate on streamflow predictions is left for future work. 308

Spatial distribution of trends in daily flows and seasonal flows 309
Trend analysis maps shown in Figure 7 decompose trends of daily flow for Q Max , Q 90 , Q 50 and 310 Q 10 . In general, the identified trends were spatially consistent with the trend pattern in Q T : 311 with upward trends in the north-west and downward trends in the south-east, south-west and 312 Tasmania. The Q 50 and Q 10 series are notable for the number stations with non-random time-313 series and therefore an invalid MK test result, this can be seen most dramatically in Figure 7d, 314 and is due to the higher correlation of the time-series. This daily flow trend analysis indicated 315 similar results to previous studies (Tran and Ng, 2009;Durrant and Byleveld, 2009)  Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2015-464, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. Whilst it is a possible explanation, it is not explicit that climate change is the cause of 370 significant trends in streamflow. There are many other factors that may affect streamflow, for 371 example, natural catchment changes, climate variability, data artefacts and other influences. 372 Site specific comparison of rainfall, PE, and temperature may help to improve the 373 understanding of the underlying causes of trends in hydrological variables. Further 374 investigation would be required to discover the potential causes of detected trends, which was 375 beyond the scope of this study. Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2015-464, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 18 January 2016 c Author(s) 2016. CC-BY 3.0 License.

•
The analysis of step changes revealed definite regional patterns: stations in 413 southeast Australia were characterised with step changes in the mid-1990s, while a 414 key feature of a 1975 step change was identified for stations in south-west West 415 Australia. 416

•
The web portal (http://www.bom.gov.au/water/hrs) displays all the graphical 417 products, tables, and statistical test results of all 222 stations. It contains a 418 comprehensive unique set of graphical products for linear trends and step change. 419 The streamflow trends evident from the statistical data analysis showed some parallels with 420 climate variability patterns that the country experienced through recent decades. Long-term 421 trends in water availability across different hydroclimatic regions of Australia reported in this 422 study are derived purely from observations unlike other studies, they are not derived from 423 models which can invariably be influenced by biases. The high quality streamflow data of 424 HRS and the results from this analysis on streamflow variability provide critical information 425 for water security planning and for prioritising water infrastructure investments across 426 Australia. 427 428 Appendix A: Statistical tests 429

A1. Median Crossing Test 430
This method tests for randomness of a time series data. It is a non-parametric test. The n time 431 series values (X 1 , X 2 , X 3 ... X n ) are replaced by '0' if X i < X median and by '1' if X i > X median . If 432 the time series data come from a random process, then the count 'm', which is the number of 433 times 0 is followed by 1 or 1 is followed by 0, is approximately normally distributed with: 434 The z-statistic is therefore defined as: 437 Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2015-464, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 18 January 2016 c Author(s) 2016. CC-BY 3.0 License.

A2. Rank Difference Test 439
This method also tests for randomness of a time series data. It is a non-parametric test. The n 440 time series values (X 1 , X 2 , X 3 … X n ) are replaced by their relative ranks starting from the 441 lowest to the highest (R 1 , R 2 , R 3 … R n ). The statistic 'U' is the sum of the absolute rank 442 differences between successive ranks: 443 | | The z-statistic* is therefore defined as: 448

A3. Mann-Kendall Test 450
This method tests whether there is a trend in the time series. It is a non-parametric rank-based 451 test. The n time series values (X 1 , X 2 , X 3 … X n ) are replaced by their relative ranks starting 452 from the lowest to the highest (R 1 , R 2 , R 3 … R n where sgn(y) = 1 for y > 0 475 sgn(y) = 0 for y = 0 476 sgn(y) = -1 for y < 0 477 X median is the median value of the X i data set. 478 The time at which 'max|V k |' occurs is considered as the time of change. The distribution of V k 479 follows the Kolmogorov-Smirnov two-sample statistic (KS = (2/n) max|V k |). A negative value 480 of V k indicates that the latter part of the record has a higher mean than the earlier part and vice 481 versa. 482

A5. Rank-Sum Test 483
This method tests whether the medians in two different periods are different. It is a 484 nonparametric test. The time series data is ranked to compute the test statistic. In the case of Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2015-464, 2016 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 18 January 2016 c Author(s) 2016. CC-BY 3.0 License.