Trends in the seasonal cycle of modelled streamflow across Australia, 1980–2018

Communities rely on a degree of certainty in the seasonality and magnitude of streamflow for consumptive use and cultural and spiritual practices (Jackson et al. 2022). Concurrently, engineers and agriculturalists determine hydraulic infrastructure (such as reservoirs, culverts, and river crossings) and water allocation based on streamflow variability as evidenced in historical records, and statistical measures thereof (for example, the empirical probability distribution and derived return periods). Extreme changes and unexpected variations in streamflow cycles impact human use and ecosystem health and sustainability. Australian Indigenous seasonal calendars highlight these interlinkages between communities, infrastructure, ecosystems, and streamflows through relationships with the country that encompass acute localised awareness and caring for land (which includes both land and water resources) (Prober et al. 2011; Woodward & McTaggart 2019). The knowledge of interlinkages across colonially stratified domains demonstrates a long and deep understanding of the history of ranges of hydrological seasonal variation and the impacts of these variations (Prober et al. 2011).

Large-scale climate patterns and climate change have been implicated in recent hydrologic changes in streamflow across regions of Australia (Zhang et al. 2016), threatening ecosystems, infrastructure, and human cultural experiences. Historical annual streamflows across south eastern Australia, namely in the states of New South Wales (NSW) and Victoria, have been significantly decreasing, while there have been significant annual increases in northern Australia (Zhang et al. 2016), and significant declines in both streamflow and rainfall in south western Australia (WA) (Petrone et al. 2010; Liu et al. 2019). Streamflows across the Australian continent are affected by coupled ocean-atmospheric climatic phenomena. The South Pacific Convergence Zone is associated with having a stronger effect on annual variability of streamflow in the northern regions of Australia than the streamflows in more southerly regions (Shams et al. 2018); while variations in rainfall and resulting streamflow more broadly across Australia have been indicated to have a strong relationship with the El Niño–Southern Oscillation, typically months in duration and the Interdecadal Pacific Ocean Oscillation of inter-annual duration (Verdon et al. 2004). The inter-annual and longer-term variations in the monsoon, presenting as increasing frequency and intensity of rainfall, strongly influence streamflow variability in the northern region of Australia (Higgins et al. 2022).

Altered streamflow can have a diverse array of impacts on the human and environmental systems linked to riverine flow. Streamflow data can be used as an indicator of episodic to seasonal and annual variability in the occurrence of these impacts and trends in streamflow changes. The above-noted studies offer a glimpse of the climate influences and trends on mean streamflow across Australia; however, the manifestation of trends over the entire distribution of streamflow is not known. Reliance on the mean to inform streamflow trends, most often through linear trends, can be easily influenced by extreme low and high flow values while concurrently losing clarity in streamflow behaviour at the low and high flows. Mann–Kendall trend tests are another popular method for assessing streamflow trends; these, however, do not delineate the low and high flows, where changes can have the greatest impact on human and environmental systems. It is therefore of interest and value to examine whether changes are occurring at low and high flows, which may be obfuscated by observing only the mean value or using other common methods. To this end, we seek to quantify trends in modelled daily streamflow across Australia by carefully considering the entire distribution of streamflow based on quantile regression and clustering. This will provide a detailed view of modelled streamflow behaviour at the low and high flows, which are obfuscated when considering streamflow through mean trends. Trends across streamflow quintiles and resolved at daily scales are especially important for a side-by-side assessment of the potential impacts on human and environmental systems.

Specifically, the analysis will be conducted on modelled streamflow and observed rainfall for daily and annual variability and trends. Analyses will include comparisons of the low, median, and high streamflows to provide a more complete picture of identified variances, followed by spatial analysis (clustering) to delineate groups of streamflow stations with similar trends. The annual median seasonal cycle for each station will demonstrate the annual periods of low and high streamflow, while quantile regressions will identify significant and near-significant changes in streamflow at daily scales. Significant changes of streamflow trend within the cluster-average high streamflow period will be identified, while the same will be performed on the identified high rainfall period, as well as the low rainfall period. The last portions of the analyses will focus on the observed rainfall data, identifying the annual mean low and high rainfall periods, and trends in substantial rainfall during these periods.

Daily observed rainfall data from the Bureau of Meteorology (BOM) from 1980 to 2018 (Bureau of Meteorology 2022) was also analysed for comparison with the modelled streamflow results. Rainfall stations were selected based on their proximity to modelled streamflow station coordinates and the completeness of daily records for the time period (Table 1). Using these matching criteria, two sets of duplicates resulted from the close proximity of streamflow stations, both of the Borroloola Crossing and MIM. Pump streamflow stations were matched with the McArthur River Mine Airport rainfall station, and both the Glenore Weir and Walkers Bend streamflow stations were matched with the McAllister Station for rainfall observations. The resulting observed rainfall dataset consisted of 35 stations, equivalent to the modelled streamflow dataset, including the two repeated rainfall stations. Rainfall leap day received the same treatment as was performed in the streamflow dataset prior to analysis.

The primary analysis methods used in this study were quantile regression (Koenker & Hallock 2001), clustering (Lattin et al. 2003), and Poisson regression (Maindonald & Braun 2006). Cluster analysis enabled the identification of data with similar features into groups, while quantile regression revealed the direction and significance of trends in low, median, and high flows. Poisson regressions were used to determine changes in rainfall during low and high rainfall seasons, in particular analysis of count data, such as the number of rain days.

The data used in this study come from diverse physiographic and climate regions, so it was of interest to identify if the streamflow stations could be grouped by common characteristics, such as the seasonality of streamflow. Euclidean clustering was chosen since it would classify the streamflow stations with the most similar patterns of the annual cycle of streamflow. Data from each streamflow station were delineated into three quintiles to further analyse the difference in streamflows across the probability distribution. Low streamflow was represented by quintile 0.2, median streamflow by quintile 0.5, and high streamflow by quintile 0.8.

Rolling pentad averages were a common starting point throughout the analyses due to the smoothing of outliers provided by using rolling averages. The average daily streamflow and rainfall values by station and by median cluster were used to compare patterns between the two datasets and median clusters. Proportions of station streamflow throughout the year were examined using the percent of station median streamflow rolling pentad averages.

Data preparation began with getting the rolling pentad average of daily streamflow values across 39 years. A basic quantile function was then applied at τ = 0.2, 0.5, and 0.8, resulting in 365 values for each station. Plots of the stations at this stage revealed a wide range in streamflow magnitudes between stations, which had the potential to skew clustering outcomes. The large difference in magnitudes between stations was standardised by rank ordering each station's values from 1 to 365, dividing by n + 1, with n being the number of days in a year, (366) to estimate the empirical quantile and then estimating the z-variates using qnorm in R (R Core Team 2022).

These values were then clustered using the K-means function from the ‘stats’ R package with the parameter nstart = 55 to stabilise cluster membership. The number of clusters, represented by a variable k, for each quintile was selected using silhouette plots, with two clusters being optimal for all three quintile levels. When k > 2, the overall average silhouette coefficient (represented by variable S) for the clusters decreases rapidly (k = 2, S = 0.31; k = 4, S = 0.16); the additional individual clusters have very low or negative silhouette coefficients, indicating that k > 2 is less effective at grouping clusters. The average daily streamflow and rainfall values by station and by median cluster were used to compare patterns between the two datasets and median clusters. The daily average across years was used to get 365 values for each location from the streamflow data and rainfall data. The locations were then grouped by median streamflow cluster membership and the daily values were averaged again for one overall cluster-average for visual comparison.

The 5-day rolling mean for each day across the 39 years was used to compute the representative daily median values, resulting in 365 daily median values for each streamflow location. These daily medians were then summed to provide the total annual median streamflow, from which the daily percent of annual median streamflow was found.

Annual seasonal variation in streamflow across each station was assessed using the daily percent of annual streamflow. The 5-day rolling mean for each day across the 39 years was used to get daily median values resulting in 365 daily median values for each streamflow location. These daily medians were then summed to provide the total annual median streamflow, from which the daily percent of annual median streamflow was found.

Analysis of the streamflow data consisted of quantile regressions to assess long-term trends in daily streamflow, the direction of these changes, and analysis of seasonal variation through the median daily percent of annual streamflow. A quantile regression affords the flexibility to assess changes at select quantiles and as such, enhances the interpretation of trend assessment. In the context of assessing hydrologic change, one benefit of this method over linear regression is that changes in direction and magnitude in the quintile extremities can be revealed, whereas they would otherwise be obscured through linear regression methods (Koenker 2005). In the application of quantile regression, a no-cross provision was employed to prevent invalid results from intersecting lines for different quintiles (Bondell et al. 2010).

Analysis of the streamflow data used quantile regressions to assess daily streamflow change significance and the direction of these changes. The following steps were applied to each of the 35 stations. The modelled daily streamflow values were averaged over 5-day to provide rolling means across the 39 years of data. For each day, the 5-day rolling average values across the 39 years were then sampled with replacement 1,000 times to generate a comparison distribution for tau levels of 0.2, 0.5 (median), and 0.8, resulting in 365 sample distributions per year. The actual 5-day rolling average regression values for each day were then located on the sample probability distribution. The slope and significance values were retained for subsequent analyses. The tail ranges of τ = 0.2 and τ = 0.8 were chosen to allow examination of changes in the extreme values over 39 years.

Significance values from the quantile regressions were used to identify daily increases and decreases in streamflow that were significant and near-significant at each τ = 0.2, 0.5, and 0.8. Increases in streamflow are represented by significance values above the 0.90 thresholds, while decreases in streamflow are represented by significance values below the 0.10 threshold. For example, a span of days with p < 0.05 at τ = 0.2 is indicative of significant decreases in streamflow when streamflow is already low and could be interpreted as representing an increase in hydrologically low flow periods.

Seasonal windows for low and high rainfall and low and high streamflow using the original modelled streamflow values were found using the following method. A 75-day rolling sum of the daily annual values were calculated for each location and the first and last years of data were removed due to the subsequent 37 days of missing data. The daily median across years for each location was then found, resulting in 365 days of data per location. Locations were then separated into their 0.5 cluster groupings and the minimum/maximum median day for each location found. These overall cluster minimum/maximum medians were then identified as the low/peak value days and 37 days were added on either side to get the full 75-day low/high season window for each cluster. Quantile regressions using τ = 0.5 were then performed on the 75-day low/high streamflow volume values against the number of years of data. The slope and significance from the regressions were used to identify significant increases and decreases in streamflow for each station within their clusters’ 75-day low/high-volume period of streamflow.

Rainfall trends across stations were analysed using Poisson regressions to identify increases and decreases in rain during the 75-day low and high rainfall windows over the years of available data. Poisson regression on the counts of no rain days and substantial rain days compared with the number of years of rainfall data will indicate whether the number of no rain days or substantial rainfall days is significantly increasing or decreasing during the identified 75-day low and high rainfall windows at each location.

Within each cluster, each station's original rainfall values within the identified 75-day high rainfall season were then dichotomised according to whether they met the threshold for substantial rainfall, that is, 0.1 inch (2.54 mm) or greater in a day. The number of days with rain during a high season was used against the number of years of available rainfall data in Poisson regression for each station. The analysis was then repeated for the low season, except the number of days without substantial rainfall was counted instead and used with the number of available years of rainfall data for each station. The slopes and significances from the Poisson regressions of both the low and high rain seasons were then used to identify if the day of lowest or peak rainfall had experienced significant increases or decreases over recent years.

The parallel coordinates plot (Figure 2(b)) demonstrates the changes to group membership across the quintile clusters with size and membership staying fairly stable in the largest two clusters across quintiles. Two stations changed cluster membership. Station 24 (Penrith, NSW) was in cluster one for quintile 0.2, changed to cluster two for quintile 0.5, and then returned to cluster one for quintile 0.8. Station 28 (Singleton, NSW) was also in cluster one for quintile 0.2 and changed to cluster two for quintile 0.5, but then remained in cluster two for quintile 0.8. These two stations are located close to each other on the southeastern coast, near the border between the two clusters (Figure 3(b)), possibly indicating that the distinction between the two clusters is less clear in this region and that these streamflows may exhibit characteristics present in both clusters or alternate between having characteristics of each cluster at different flow levels. This may be due to the stations’ proximity to the coast in this region, causing different streamflow behaviour compared with nearby and inland stations.

The accuracy of the modelled data was evaluated using Kling–Gupta Efficiency (KGE) validation values, which ranged from −1.97 to 0.84. These KGE values were included in the dataset provided by Alfieri et al. (2020). For the purposes of this analysis, KGE values greater than −0.41 are considered as improving upon mean streamflow (Knoben Freer & Woods 2019; Ayzel et al. 2021; Nickles & Beighley 2021). The KGE values were split into four ranges with less than −0.41 being very low, −0.41 to zero being low, zero to 0.5 being moderate, and good being greater than 0.5. These are represented in the left side column of Figure 5 with the corresponding streamflow station on the right. Ten locations had KGE values greater than 0.5, 14 locations had moderate KGE values, and 7 locations had low KGE values. The remaining four locations had very low KGE values (KGE < −0.41).

The southern cluster displayed very different significant daily trend patterns, with most days having zero stations with a significant positive trend. The lower quintile 0.2 had some days early and in the middle of the year with one station per day having a significant positive streamflow trend; the median had a similar pattern but with greater sparseness, and the upper quintile 0.8 had mostly days at the end of the year with one station per day experiencing a significant positive streamflow trend. The significant negative streamflow trends were far more prevalent than the significant positive trends throughout the year in the southern cluster, both in daily frequency and station frequency. Across τ = 0.2, 0.5, and 0.8, there was a lower frequency of stations through approximately March and April, with up to five–eight stations per day recording significant negative streamflow trends throughout the rest of the year, primarily in January through February and May through November.

This trend changes across quintiles throughout the year demonstrate the changing behaviour of streamflow across the year. For example, increases in high streamflow trends indicate an increased risk of floods, while decreasing trends in low flow could restrict the water supply available for human activities and impact biota processes. It bears mentioning that DS Burrendong is a dam with the streamflow heavily regulated, so the decreasing trends are likely related to conserving water in the dam while increases at low flow may be related to releasing stored water during drier periods in the middle of the year.

Poisson regression results of rain event frequency, wherein events are defined as those exceeding the 2.54 mm trace amounts, within the high and low rain seasons were grouped by streamflow cluster (Table 2). Out of the 24 stations in the median northern cluster during the high rain seasons, 4 had significant positive slopes and 3 were significantly negative. The rainfall stations with significant negative values were located in the northern areas, while significant positive Poisson regression values were more spread out across the northern and eastern coasts (Table 2). These indicate that some of the locations in the northern areas are experiencing fewer days of rain during their respective high rain seasons, while increases in days of rain during high rain seasons are increasing inconsistently across the region. The Poisson regression results in the southern cluster are another reflection of the decreasing streamflow trends in the region, with negative trends in days of rain during the high rain seasons across all stations except one (Penrith) during the high rain season. However, Drik Drik is the only station with a statistically significant negative trend.

The change in the number of days of rain during the respective low rain seasons was mixed positive and negative but non-significant for both the northern and southern cluster locations. The prevalence of within-cluster increases and decreases in days of rain in the northern and southern locations indicates the potential for impacts on ecological and environmental seasonal processes, while the absence of coherent, widespread change in days of rain during the low and high seasons points to the highly variable nature of the seasonal monsoon extent, intensity, and episodic nature of moisture incursions. Changes in the days of rain during the low and high rain seasons can function as a stressor on environmental and ecological processes, as insufficient and excesses of rain days or changes to interstorm intervals can parch environments or inundate them beyond infiltration capabilities. Meanwhile, ecological processes that have developed to particular rhythms of rainfall may be prematurely initiated, delayed, or interrupted, which has carry-on effects throughout ecological systems.

Spatial patterns of streamflow variability were delineated robustly by examining coherent scales of variability across three representative quintiles of streamflow resolved to daily time scales. The identified groups of streamflows belonging to the two clusters show similar variability and thus represent regional scales appropriate for analyses of human-environmental system linkages, future changes induced by climate change and scales over, which may be drivers of hydrologic variability, for example, monsoonal weather may evolve and affect regional resources.

The northern and southern regions exhibit markedly different streamflow and rainfall behaviour. The northern group demonstrated strong seasonality with high rainfall and streamflow from November through May followed by much-reduced rainfall and streamflow in the later months. The streamflow trends were mixed increases and decreases but increases were more frequent early in the year across the low, median, and high flows. In comparison, the southern region experienced a slight increase in rainfall and streamflow in the middle to later portion of the year; however, streamflow trends were almost exclusively negative with the greatest frequency in the early and later months of the year.

Communities develop in tandem with their environment, particularly around their water sources, which are so vital for health and life. Excesses or shortages in water can be devastating to human health, life, and infrastructure. The northern streamflow stations demonstrate strong seasonality with a dramatic spike in streamflow volume followed by a much-decreased volume; however, communities in these northern regions are not devastated by the considerable changes in streamflow volume and associated meteorological events. Rather, their lives are adapted to these seasons. This is particularly evident among Indigenous peoples within the northern region and across Australia, whose relationships with the country have produced detailed knowledge of local seasonal activity, which some Indigenous nations have made available in seasonal calendars (Green Billy & Tapim 2010; Prober et al. 2011; Woodward & McTaggart 2019). Many of these Indigenous seasonal calendars have observations for the timing, duration, and severity of hydroclimatic events. Recent climate variability and change are driving an increase in variations from what has historically been consistent seasonal cycle within an expected range of variation; however, the seasonal calendars, through respectful and equitable collaboration with Indigenous knowledge holders, can assist in identifying and anticipating variation in hydroclimate events at local scales. These fine-scale observations of estimated hydroclimate deviations can be further quantified through the analysis methods used in this study.

The much-reduced scale of seasonality in the southern cluster has meant that communities in these regions have had less experience with regular water excesses and shortages, resulting in a historically reduced need for communities to develop adaptation strategies to such extremes. This prior lack of need to manage water deficits could result in greater vulnerability to increased occurrences of such events, which is of particular concern since the significant negative streamflow trends across the southern cluster indicate reduced water availability throughout the year. Adaptation strategies used in the northern cluster locales may be applicable to locations in the southern cluster for improving resiliency to the increasing water uncertainty coinciding with climate change.

Alternatively, knowledge of how to live and manage resources sustainably is likely already present among the Indigenous peoples in the region of the southern cluster through their greater holistic knowledge of their local region (Woodward 2013; Woodward & McTaggart 2019; Yang et al. 2019). Knowledge available within Indigenous seasonal calendars is indicative of rich and detailed knowledge of signs of dry periods and management methods to reduce the severity of these periods (McKemey & Rangers 2020). There is a great opportunity for collaboration efforts prioritising respect and equitable outcomes through relationship building between Indigenous peoples and colonial benefactors. Such collaboration could result in real changes to resource management, policies, and attitudes towards relationships with country, which could enhance adaptation methods for climate-driven declining water trends in the southern cluster.

Streamflow and rainfall trends are changing across Australia in accordance with the hydro-meteorological response to climate change, as demonstrated in this study through (a) Characterisation of the full range of streamflow variations and trends resolved to daily timescales; and (b) Cluster analysis applied to results from three different quintiles, revealing the relative coherence in groupings as well as streamflow magnitude-related regional variations. Clustering and quantile regression methods delineated locations into groups of streamflow seasonality and presented at daily resolution how seasonal trends are changing at each location at low and high extremes and not just the average, while Poisson regressions identified changes in substantial rainfall during low and high rainfall seasons. The increases and decreases in rainfall and streamflow trends are representatives of advancing hydrological uncertainty in the region, which can culminate in floods and droughts with devastating immediate and long-term effects on communities and regions.

The analyses presented in this study sought to present quantitative estimates of long-term changes in streamflow to further better understanding and adaptation to the impacts of climate change. Similar approaches could collaborate with local knowledge holders, such as Indigenous peoples with their detailed knowledge of the expected range of variations in local seasonal hydroclimate events, to enable estimates of unusual variations in hydroclimate at fine scales. For example, when examined simultaneously within a region, the daily-resolved streamflow trends and Indigenous seasonal calendars (Coleman 2022) have a higher likelihood of uptake of knowledge for pre-empting, mitigating, and adapting to the localised impacts of climate change. In addition, collaboration with community members in regions that have historically extensive experiences of similar seasonality presentations to the current and future seasonality trends of other regions may be able to offer insight into resiliency methods.

All relevant data are available from an online repository or repositories. (Alfieri streamflow reanalysis data: https://data.jrc.ec.europa.eu/collection/id-00288 and BOM rainfall data: http://www.bom.gov.au/climate/data/).

The authors declare there is no conflict.