Unprecedented High Northern Australian Streamflow Linked to an Intensification of the Indo‐Australian Monsoon

Streamflow in Australia’s northern rivers has been steadily increasing since the 1970s, most likely due to increased intensity in the Indo‐Australian monsoon. However, because of limited data availability, it is hard to assess this recent trend and therefore contextualize potential future climatic changes. In this study, we used a network of 63 precipitation‐sensitive tree‐ring chronologies from the Indo‐Australian and Asian monsoon regions to reconstruct streamflow in the Daly catchment in the Northern Territory of Australia from 1413 to 2005 CE. We used a novel wavelet‐based method to transform the variance structure of the tree‐ring chronologies to better match the hydroclimate prior to reconstruction with a hierarchical Bayesian regression model. Our streamflow reconstruction accounts for 72%–78% of the variance in the instrumental period and closely matches both historical flood events and independent proxy records, increasing confidence in its validity. We find that while streamflow has been increasing since the 1800s, the most recent 40‐year period is unprecedented in the last ∼600 years. Comparison to an independent coral‐based streamflow record shows regional coherency in this trend. Extreme high flows were found to be linked to La Niña events, but we found no significant relationship between streamflow and El Niño events, or streamflow and other regional climatic drivers. More work is therefore needed to understand the drivers of the recent streamflow increase, but, regardless of the cause, water managers should be aware of the paleoclimatic context before making decisions on water allocations.

peoples, are dependent on the wet/dry seasonality in flows (CSIRO, 2009;Russell-Smith et al., 1997). The IASM is closely tied to the El Niño-Southern Oscillation (ENSO). Both climatic phenomena, the IASM and ENSO, exhibit substantial inter-annual and decadal variability (Smith et al., 2008;Suppiah, 1992). This variability can cause both droughts and floods, with severe impacts on the lives and livelihoods of people living in the region.
Since 1950, monsoonal Australia has experienced an increase in summer rainfall (Taschetto & England, 2009), and annual streamflow (Zhang et al., 2016), thought to be a result of changes in the IASM. Historical meteorological records from ship logs show that this recent trend is part of a longer-term intensification of the IASM, occurring at least since the early 1800s (Gallego et al., 2017). Despite this intensification, rainfall variability remains high; consecutive years of monsoon failure occurred during the 2018-2019 and 2019-2020 wet seasons, resulting in very low river levels and water restrictions over much of Australia's north (BOM, 2019(BOM, , 2020. It is highly uncertain whether the recent trend in rainfall will continue. Projections of future monsoon rainfall are challenging because of the strong climate coupling between the land, atmosphere, and oceans. Climate models predict a future thermodynamic increase in mean monsoonal rainfall (Christensen et al., 2013), however, there are large inter-model uncertainties, most notably for the IASM (Brown et al., 2017). The ensemble mean of the Climate Model Intercomparison Project Phase 5 projects an increase in IASM rainfall variability at daily to decadal timescales (Brown et al., 2017). Increased variability has important implications for the management of water resources, as existing infrastructure, management plans, and entitlements may be inadequate for future changes (De Loe & Kreutzwiser, 2000), particularly in regions like monsoonal Australia where water storage potential is limited (CSIRO, 2009).
Palaeohydrological reconstructions are useful tools for supporting water resource management and planning under uncertainty Rice et al., 2009). There has been a recent interest in palaeohydrological reconstructions for monsoon rivers (e.g., Nguyen et al., 2020 and references therein) due to the inability of short instrumental records to fully characterize hydroclimatic variability, and uncertainty around future changes in streamflow and flood frequency (Rao et al., 2020). Long reconstructed records can increase our understanding of variability in monsoon streamflow and floodplain dynamics, as well as their relationship with major drivers of climate variability such as ENSO. However, very few streamflow reconstructions are available for the IASM region (see D'Arrigo et al., 2011;Verdon-Kidd et al., 2017) in part due to the limited availability of high-resolution (sub-annual to annual) precipitation and streamflow proxies across the region D'Arrigo et al., 2008). There is an urgent need for multi-centennial reconstructions of streamflow variability in monsoonal Australia to reduce the uncertainty over the impact of future climate change and increasing resource demands.
The Daly River, a monsoonal river system in the northwest of the Northern Territory, is one of the northern Australia's largest rivers and one of the few with perennial flows (Morrison, 1970). Annual streamflow is highly dominated by monsoon rainfall, with dry season flows, fed by groundwater discharge, accounting for less than 10% of total annual flow. Reliable stream gauging stations and monitoring bores for the Daly River system are sparse. Very limited data exists prior to the mid-1970s and there is low confidence in dry season flow records for many gauges. The paucity of data greatly inhibits the potential to assess the linkages between the hydrological regime, resource availability and dependent ecosystems (CSIRO, 2009). Declining rainfall trends in southern Australia and the perception of abundant land and water resources in the monsoon tropics have recently renewed interest in agricultural development in northern Australia and the Daly region (Petheram et al., 2008;Yeates et al., 2013). With limited storage potential in seasonally dynamic aquifers, increased dry season agricultural demand is likely to change surface and groundwater regimes with consequences for the environment (CSIRO, 2009). Extending short hydrological records is essential for improving the assessment and management of water resources in the Daly catchment and the IASM region more broadly.
Here we present a paleohydrologic reconstruction of Daly River streamflow that extends available instrumental observations by more than five centuries and discuss its linkages with both changes in the IASM and Pacific atmosphere-ocean climate variability. Our Daly River streamflow reconstruction was developed using a tree-ring network from across the IASM region and remote sites with strong teleconnections to the IASM. We demonstrate the benefits of a new method for paleohydrologic reconstruction, in which a wavelet-based variance transform maximizes the ability of the tree-ring predictors to characterize the hydroclimate, resulting in a streamflow reconstruction that better matches the instrumental data. We then compare our record to other observational and proxy records for the region, specifically low-resolution sediment cores from the Daly River and a coral-based river . The inset map shows the location of the Daly catchment (gray) and Burdekin catchment (beige) within the Australian tropical monsoon region (light blue), and the locations of Callitris intratropica and C. columellaris tree-ring chronologies within and proximal to the monsoon region (black diamonds). The approximate southern monsoon region boundary was adapted from the Australian Bureau of Meteorology climate classification.
the Oolloo Dolostone (CSIRO, 2009). Variations in dry season baseflow are due to higher or lower aquifer levels in the preceding wet season (CSIRO, 2009).

Hydrological Data
Daily streamflow data (ML/day) for gauges in the Daly River catchment were downloaded from the Northern Territory Government Water Data Portal at https://water.nt.gov.au/. Gauges on perennial rivers with at least 40 years of data before 2005 (the last year for many tree-ring chronologies) were selected. Suitable gauges were identified on the Douglas River downstream of the Old Douglas Homestead, the Katherine River at Railway Bridge, and the Daly River upstream of Dorisvale Crossing (Figure 1). For the Daly River at Mount Nancar, the gauge nearest to the catchment outlet, the time series was extended from 1961 to 1971 using data from nearby gauges and a regression model (see Supporting Information S1).
Years with more than 15% of daily values missing during the peak discharge period (January-March), were excluded from analysis (Table 1). For other years, gap-filled daily streamflow data was downloaded from the Australian Bureau of Meteorology (BOM) Hydrologic Reference Stations database (http://www.bom.gov.au/ water/hrs/index.shtml). Annual streamflow (km 3 ) was then calculated for the Northern Territory water year of September to August (CSIRO, 2009). Streamflow between the four gauges is highly correlated (Pearson r = 0.78-0.97; Figure S2 in Supporting Information S1). Including all four gauges in the reconstruction model meant that 1962 was the only year without any streamflow data (i.e., missing data from all four gauges).  55 19611962, 1966, 1977, 1978, 1985G8140067 Daly Upstream Dorisvale Crossing 33,227 5.79 19661968, 1975, 1980, 1989, 1991G8140063 Douglas Old Douglas Homestead 831 0.25 19591962, 1968, 1969, 1971, 1976, 1977, 1993, 2002 Adler et al., 2003) monthly gridded rainfall from 1979 to 2020, used to calculate the correlation between Daly catchment and regional rainfall, and the Extended Reconstructed Sea Surface Temperature (ERSST) v5 dataset Huang et al. (2017), were both obtained from the NOAA/OAR/ESRL PSL website at https://psl.noaa.gov/.

Historical Flood Records
The largest settlement in the Daly catchment is the town of Katherine on the Katherine River, also the location of gauge G8140001. Three significant flood events, in 1998, 2006, and 2011, have occurred since streamflow was reliably recorded at Katherine. However, only the 1998 event falls within the reconstruction period. Pre-instrumental data sources, predominantly newspaper articles, provide a record of major flood events since the permanent European occupation of the region began in 1872. Large floods at Katherine occurred in December 1897 (reconstruction year 1898), 1914, 1931, 1940, and 1957(Water Studies Pty Ltd., 2000.

Tree-Ring Proxy Records
Very few tree-ring chronologies exist for the tropics, including northern Australia, due to the difficulty in finding tropical tree species with annual growth rings (Rozendaal & Zuidema, 2011). However, recent efforts have identified the dendrochronological potential of both Callitris intratropica and C. columellaris (cypress pine) species (Allen et al., 2019(Allen et al., , 2020Baker et al., 2008;O'Donnell et al., 2015). While these species have great potential to inform understanding of the hydroclimate of monsoonal Australia (Allen et al., 2020;D'Arrigo et al., 2008), the chronologies developed to date are short, with the majority ∼100-250 years in length. Longer reconstructions can be developed by utilizing remote tree-ring proxies  from regions with strong teleconnections to monsoonal Australia. For the initial predictor pool, we downloaded all publicly available tree-ring chronologies from Australia, New Zealand and Monsoonal Asia located between 70°E, 40°N and 180°E, 47°S ( Figure S4 in Supporting Information S1) from the International Tree Ring Databank (ITRDB). The wood property series from each site were converted to chronologies by the process of standardization (i.e., detrended and transformed into dimensionless growth indices). The process aims to remove growth trends thought to be largely unrelated to climate using the "signal free" (Melvin & Briffa, 2008) and "Regional Curve Standardization" (Briffa et al., 1992) methods of tree-ring standardization.

Methods
Figure 3 provides a schematic overview of the methods used to reconstruct Daly catchment streamflow. The methods are described briefly below and in further detail in the Supporting Information S1.

Streamflow Reconstruction
Tree rings are useful predictors of streamflow because both tree growth and variations in streamflow are controlled by soil moisture, which is influenced by rainfall and evapotranspiration (Meko et al., 1995). However, tree rings provide only indirect information about climate variations. In addition, factors including local proxy noise, seasonal sensitivities, and standardization methods result in biases in the high-to low-frequency spectrum of proxy records compared to the instrumental datasets they are attempting to reconstruct (Franke et al., 2013). Spectral biases are likely to be compounded when using remote tree-ring chronologies in reconstructions as only a portion of the climatic information contained in the proxy-related to the teleconnection-is relevant to the reconstruction.
To address these issues, after selecting the tree-ring chronology predictor pool, we applied a novel method for paleohydrologic reconstruction, in which a unique wavelet-based variance transform (Jiang et al., 2020(Jiang et al., , 2021 was used to modify the spectral characteristics of each tree-ring chronology to better match Daly catchment average annual rainfall (see Supporting Information S1). A nested hierarchical Bayesian regression model with partial pooling (Devineni et al., 2013;Rao et al., 2018), which is suitable for our short data records with missing years, was then used to produce the streamflow reconstruction from the transformed chronologies. The partial pooling framework combines the regression strength of the model across the four gauges, which can result in lower uncertainty in estimated parameters and reconstructed discharge, as well as improved final model skill (see Supporting Information S1). We used a moving block calibration-verification scheme (Nguyen et al., 2020) in which 7-year, overlapping periods of the instrumental record were successively withheld to independently test the modeling (see Supporting Information S1).
To evaluate the reconstruction, we used standard verification tests in dendrochronology: the calibration period coefficient of multiple determination (CRSQ or R 2 ), the validation period reduction of error (VRE), and the validation period coefficient of efficiency (VCE; Cook & Kairiukstis, 1990). The VCE is equivalent to the Nash-Sutcliffe efficiency test (Nash & Sutcliffe, 1970) and is the most stringent verification criteria. The Bayesian R 2 , a data-based estimate of the proportion of the variance explained for new data Gelman et al. (2019), and the Sign Test, which calculates the number of years that the reconstruction correctly (+) or incorrectly (−) tracks the sign of observations during the calibration period (Cook & Kairiukstis, 1990), were also used in verifying the model.

Historical Flood Events
To further verify our results, the ability of the streamflow reconstruction to identify historical flood events recorded at the Katherine township was assessed via superposed epoch analysis (SEA), a compositing technique that tests the probability of an association between high flows and flood years occurring by chance (Haurwitz & Brier, 1981;Rao et al., 2019). To apply SEA to flood identification in our reconstruction, the six historical flood events which occurred during the reconstruction interval were used as key event years for SEA. A composite matrix was created by selecting reconstructed streamflow at Katherine (G8140001) from 10 years prior to 5 years following each event. Streamflow was normalized to the 10 years before an event to remove noise unrelated to the event signal, then averaged across all events to determine the composite signal (see Supporting Information S1). The significance of the flood event-streamflow relationship was tested by comparing the key event composite to that generated from 10,000 draws of 6 years at random without replacement ("pseudo-flood years") from the reconstruction between 1888 (first event year-10) and 2005 CE.

Extreme Event Analysis and Climate Forcing
Extreme dry/wet events were classified as those exceeding below the 5th and above the 95th percentiles of flows, respectively, during the reconstruction period (1413( -2005. Changes in the frequency of extreme events were estimated using a nonparametric Gaussian kernel function (Mudelsee et al., 2004) with bandwidth selected via the method of Sheather and Jones (1991) and 95% confidence intervals estimated based on 1,000 bootstrap simulations. Trends in the recurrence rate of extreme events were tested against the null hypothesis of constant recurrence rate using the Cox-Lewis statistic (Cox & Lewis, 1966): which is a standard normal distribution under the null hypothesis (Mudelsee et al., 2003), where T(1), T(n) are the first and last years of the observation interval, and m is the number of extreme events.
Bootstrapping was used to assess whether extreme events are significantly associated with different ENSO phases (Allen et al., 2020). The occurrences of El Niño/La Niña events since 1875 were determined based on December-February (DJF) Nino3.4 sea surface temperature (SST) anomalies above/below 0.5°C and referenced against the events listed on the Australian Bureau of Meteorology website (www.bom.gov.au). Bootstrap replicates were drawn with replacement from the SST data to match the same number of years as the wet and dry extremes. The co-occurrence of extremes and ENSO phases was considered significant if the number of actual events lay outside the 95% bias-corrected and accelerated (BCA) bootstrap confidence interval (Efron, 1987) of the distribution based on 30,000 bootstrap replicates. The number of bootstrap replicates was increased from 1,000 until the confidence interval remained unchanged. Figure 4 shows the correlations between Daly catchment average rainfall and rainfall from regions with strong teleconnections to the IASM, during both the dry and wet (monsoon) seasons. Daly catchment rainfall is significantly correlated to rainfall in southeastern Australia during the May-September dry season when Indian Ocean sea surface temperatures affect Australian rainfall (McBride & Nicholls, 1983;Risbey et al., 2009). Conversely, Daly catchment rainfall is significantly anti-correlated to the Asia Summer Monsoon during the November-April wet season, reflecting coherence between the regional monsoons (Wang et al., 2014). Correlations between Indonesian and northern Australian rainfall are significant throughout the year but are stronger during the dry season. Figure 4 also shows the locations of the 63 tree-ring chronologies with significant correlation to Daly catchment streamflow retained for the final model (see Table S1 in Supporting Information S1 for details). Nearly 70% of the retained predictors are from regions impacted by monsoon rainfall. Five chronologies are within the IASM region, of which three are proximal to the Daly catchment (within 20 km, see Figure 1). A further 36 are within the Asian Summer Monsoon region, and the remaining 22 from subtropical to temperate Australia, and New Zealand. 8 of 17

Tree-Ring Predictor Selection
Absolute Pearson correlations between streamflow and significant tree-ring predictors vary from r = 0.28 to 0.61 with a median value of r = 0.42. This is a substantial increase in correlation compared to the pre-transform predictors, which varies from r = 0.0 to 0.51 with a median value of r = 0.21. Transforming the spectral variance of the tree-ring chronologies is therefore shown to strengthen their correlation to Daly catchment streamflow.

Streamflow Reconstruction
The final nested reconstruction based on the variance transformed tree-ring chronologies spans 1413 to 2018 CE, with instrumental data appended after 2005. Figure 5a compares observations with reconstructed streamflow at the most downstream gauge, G8140040 at Mount Nancar, using the tree-ring data for the best replicated (1898-2005) nest. Results for the other three gauges are shown in Figure S9 in Supporting Information S1. This reconstruction accounts for 72%-78% of the variance in streamflow observations at each gauge over the calibration period from 1959 to 2005. Table 2 shows the calibration-validation statistics used to evaluate the validity and stability of the nested models for the best-replicated nest, and median statistics of all nests between 1413 and 2005 CE. Timeseries results for all nests are provided in Figure S10 in Supporting Information S1.
The calibration and verification statistics demonstrate that the reconstruction is reliable, despite the relatively short calibration period. Despite the decreasing number of tree-ring predictors moving back in time, the statistical results remain strong for the earliest part of the record, with the least replicated nest accounting for 59%-62% of the calibration period variance. Uncertainty intervals ( Figure S9 in Supporting Information S1) show that the fifth percentile VRE and VCE values are above zero throughout the reconstruction period, indicating that the reconstruction contains meaningful information over its entire length (Cook & Kairiukstis, 1990). In addition, the significant (p < 0.01) Sign Test results are a good indication that the tree-ring chronologies are accurately tracking the direction of year-on-year variability in streamflow during the calibration period. The full reconstruction for G8140040 is shown in Figure 5b; results are similar for the other gauges (not shown). The reconstruction shows that streamflow during the instrumental period represents a period of unusually high flow that is unprecedented in the preceding five centuries. The 57-year gauge instrumental period from 1961 to 2018 has mean streamflow of 7.79 km 3 /year, which is significantly higher (p < 0.001, two-sided Student's t-test) than the pre-instrumental mean of 4.60 km 3 /year, and any period of equivalent length in the reconstruction. Similarly, the gauge calibration period from 1961 to 2005, which has mean streamflow of 6.99 km 3 /year, exceeds any previous 44-year period, but the difference is not significant for periods between 1873 and 1922 CE. These results are consistent with an intensification of IASM rainfall starting in the early to mid-1800s.

Historical Comparisons
The largest recorded flood event at Katherine occurred in 1998 after ex-Tropical Cyclone Les brought 400-500 mm of rain within three days (Skertchly & Skertchly, 1999). The rainfall characteristics of other major historical flood events are similar, with several hundred millimeters of rain falling within a few days (Water Studies Pty Ltd., 2000). There is a moderately strong and significant positive correlation (r = 0.60; p < 0.001) between the annual maximum 5-day rainfall, which represents the approximate length of rainfall events resulting in historical floods, and total water year rainfall at Katherine (Figure 6). In most cases, major floods occurred during very wet years, with total rainfall in the top ∼11% of the record. The exceptions are 1914 and 1931, which were average rainfall years. These results are indicative only as extreme rainfall events, especially those related n is less than the number of calibration years due to data gaps. to cyclonic rainfall, are not evenly distributed over the catchment, and there is a paucity of rainfall gauges in the early part of the record.

Table 2 Calibration and Verification Statistics for the Best-Replicated Nest (1898-2005) and Median Values for the Whole Reconstruction Interval (1413-2005) for the Four Daly Catchment Gauges
The significant positive relationship between high-intensity events which may cause flooding and total annual rainfall suggests flood years should be identifiable by higher-than-average reconstructed Katherine River streamflow. We tested the probability of random association between flood years and high streamflow at gauge G814001 using superposed epoch analysis. We found that mean reconstructed flows across these 6 years are significantly higher (p < 0.01) than would be expected by chance ( Figure 6). The ability of the reconstruction to identify pre-calibration historical flood events increases confidence in the reconstruction model outside of the instrumental period.

Extreme Events
Surface warming of the equatorial Pacific associated with El Niño events drives an intensification of the Walker circulation, resulting in decreased average deep convection over Australia (Holland, 1988;McBride & Nicholls, 1983). The reverse occurs during La Niña events, with weakening of the Walker Circulation and increased convection. However, this relationship is not symmetric, with La Niña events generally more highly correlated with rainfall over Australia than El Niño events (Power et al., 2006;Risbey et al., 2009). Teleconnections between ENSO phases and Australian rainfall also vary with the phases of the Interdecadal Pacific Oscillation, with warm (positive) phases decreasing the strength and spatial coherence of ENSO-rainfall relationships (Power et al., 1999). Recent research (Allen et al., 2020) has shown that precipitation at the end of the wet season (March-May) in monsoonal Australia is asymmetrically linked to ENSO, with increasingly extreme wet events associated with increasingly cooler SSTs, but less clear relationships between dry extremes and El Niño phases.
We tested the relationship between annual Daly River streamflow extremes and Pacific SSTs by comparing the number of high and low streamflow events that co-occurred with a known ENSO event. Figure 7 shows that extreme high streamflow events, defined as exceeding the 95th percentile of reconstructed flow, are significantly linked to Pacific SST anomalies, as demonstrated by the significant co-occurrence with La Niña events. Conversely, low flow events, defined as below the 5th percentile of reconstructed flow, were found to significantly not co-occur with La Niña events. El Niño events were not found to be significantly related to either high or low flow extremes, with the number of co-occurring events very similar to the mean of the bootstrap distribution. Not all El Niño events are associated with weak IASM rainfall. While lower-than-average monsoon rainfall occurs during canonical El Niño events, there is little change in total rainfall during El Niño Modoki events, which result in a shorter but more intense monsoon (Taschetto et al., 2010). We tested the relationships between El Niño events and Daly River low flow extremes for canonical El Niño and El Niño Modoki events separately, but again the results were not significant ( Figure S11 in Supporting Information S1).
These findings are consistent with the relationship between cool Pacific SSTs and rainfall over Australia in observations, and are similar to the findings of Allen et al. (2020), despite the differences in the study season. Figure 7, panels (e) and (f), shows the composite December-January SSTs over the 9 extreme low and 15 extreme high flow years. High streamflow years are associated with warmer-than-average SSTs off the coast of Australia which lead to greater convection and monsoon rainfall, whereas low flow extremes are associated with cool SST anomalies and suppressed convection (BOM, 2012). The mechanisms behind the asymmetrical relationship between ENSO and northern Australian rainfall and streamflow is an area of active research, with the interaction between ENSO and other modes of variability a likely cause (Cai & van Rensch, 2013;Heidemann et al., 2021;Power et al., 2006).
While high total annual streamflow is more likely to occur during a La Niña than El Niño phase, the rainfall events that cause significant flooding at Katherine cannot be linked to ENSO phasing based on the available data. While the flooding at Katherine during 2011 was associated with a strong La Niña that also caused extensive flooding in Eastern Australia, other major historical floods at Katherine have occurred during years with average Pacific SSTs and El Niño years. As most of the region's rain falls as intense, intermittent, tropical showers, monsoon troughs, or ex-tropical cyclone lows, rain events during average years are sufficient to cause local flooding. Tropical cyclones in the instrumental record are ∼20% more likely to occur in northern Australia during La Niña years (Kuleshov et al., 2009), but ENSO phases do not significantly impact the occurrence of monsoon bursts which are more likely to occur during active phases of the Madden-Julian Oscillation (Berry & Reeder, 2016). Rainfall event timing linked to antecedent conditions is also a factor, with heavy rainfalls occurring later in the wet season when the subsurface is saturated more likely to result in extreme discharge than equivalent events earlier in the season (Chappell & Bardsley, 1985).
Indian Ocean SSTs also affect northern Australian rainfall, with positive phases of the Indian Ocean Dipole (IOD) linked to decreased winter-spring rainfall. The influence of the IOD begins in May-June, and peaks in September-October, diminishing rapidly at the start of the monsoon. Therefore, the IOD has very little impact on the peak of the monsoon (Jourdain et al., 2013;Taschetto et al., 2010), and, as expected, we found no relationship between the IOD and Daly River annual average or extreme streamflow. An IASM-region streamflow reconstruction for the Citarum Basin in Java, Indonesia, which specifically targets spring (Sept-Nov) streamflow, has shown an increase in the frequency of low streamflow events accompanying the increase in the frequency and magnitude of positive IOD events since the 1960s (D' Arrigo et al., 2011). This illustrates the importance of the seasonality of reconstructions when considering the impacts of different climate drivers on extreme events in the IASM region.
Changes in the frequency of extreme low and high streamflow events in the Daly catchment were investigated using a Gaussian kernel technique. Results for the Katherine River and Daly River gauge G8140040 are shown in Figures 8a and 8b and remaining gauges in Figure S12 in Supporting Information S1. Persistent weak monsoons have been proposed as the cause of the sixteenth-century megadrought (1560-1587 CE) in monsoon Asia (Cook et al., 2010). We also see a peak in the occurrence of low flow events in the mid-sixteenth century, noting that the Daly River reconstruction is not independent of the drought reconstruction of Cook et al. (2010).
The mid-1600s to early 1800s shows low reconstruction variance and a low frequency of both low and high streamflow events. Paleoclimate studies suggest two potential mechanisms that could explain this result. The low variance period broadly coincides with the Little Ice Age (LIA), a period of lower-than-average temperatures in the Northern Hemisphere, believed to result from both low solar activity (Maunder and Dalton grand solar minima) and high volcanic activity (Grove, 1988 ton et al., 2016;Yan et al., 2015). Proxy reconstructions of ENSO by D' Arrigo et al. (2005) and McGregor et al. (2010) have also identified a period of low ENSO amplitude over the 17th and 18th centuries, although their findings disagree with other studies (Cobb et al., 2003). Either smaller north-south movements of the ITCZ or lower ENSO amplitude could plausibly result in lower monsoon streamflow variability during the LIA. Figure 8 shows that the occurrence of high flow events has increased markedly over the last ∼100 years. A significant (p < 0.001) increasing trend in 95th and 90th percentile flows at Daly River gauge G8140040 was identified over the period 1800 to 2018 CE using the Cox-Lewis test against the null hypothesis of a constant occurrence rate. Noting that the Daly catchment gauge reconstructions are not independent, trends in both 95th and 90th percentile flows at Katherine River gauge G8140001 were also significant over this interval (p < 0.01; p = 0.05, respectively), although the increase in high flows began slightly earlier at ∼1770 CE. The results did not change if either the reconstruction end date of 2005 or the observation end date of 2018 were used, and increasing trends calculated over the entire reconstruction interval from 1413 CE were also significant for both gauges.
A sand splay deposit taken at the Nancar Hideout on the lower Daly River, close to the location of streamflow gauge G8140040, shows that low rates of deposition occurred from ∼1420 to ∼1760 CE, indicating a dry period, followed by an increase in deposition and streamflow in the period from ∼1760 to 2005 CE (Wasson et al., 2010). Flood frequency derived from the sediment cores has risen steadily over the last ∼160 years, with a doubling in the last ∼60 years (Wasson et al., 2010). Noting the dating uncertainties (±10-60 years), this record is highly consistent with our tree ring-based reconstruction, providing independent verification of the recent trend in extreme streamflow.

Regional Coherency
To test whether recent increases in monsoon streamflow are regionally coherent or unique to the Daly catchment, we compared our terrestrial proxy reconstruction to a coral luminescence-based reconstruction of the Burdekin River from 1648 to 2011 (Lough et al., 2015). The Burdekin River, located in the dry tropics of Northeast Queensland (Figure 1), also experiences wet season rainfall related to the IASM. Lough et al. (2015) found an increase in the magnitude and frequency of high flow events from the mid-1900s compared to the preceding century in their Burdekin streamflow reconstruction.
We assessed changes in Burdekin River streamflow using the same methods described above for the Daly catchment (Figure 8c), extending the reconstruction from 2011 to 2018 using instrumental data. As for the Daly catchment, trends in both 95th and 90th percentile Burdekin River flows increased significantly over the interval from 1800 to 2018 CE (p < 0.05), with trends also significant since the beginning of the reconstruction period at 1648 CE (p < 0.05).
Both catchment and rainfall changes may have contributed to the increasing streamflow trends. Agricultural development in the Burdekin Catchment since the settlement period (1851-1900) has likely contributed to higher discharge through vegetation cover and soil compaction changes to rainfall-runoff ratios (Lough et al., 2015). However, there are currently low levels of development and intensive agriculture in the Daly, with less than 10% of land under intensive use (Álvarez-Romero et al., 2016). Land surface changes, therefore, are unlikely to have had a significant contribution to the observed increase in streamflow in the Daly River or the coherent trend between the catchments.
Recent increasing trends in northern Australia summer rainfall are well documented. The increasing trend has been observed using rainfall data from 1900 onward, but may have begun in the early 1800s (Gallego et al., 2017). The trend has intensified since the 1950s (Nicholls, 2004;Rotstayn et al., 2007;Suppiah, 1992;Taschetto & England, 2009), which matches the observed changes in the frequency of extreme events in the Daly and Burdekin catchments. Most of the increasing rainfall trend can be explained by an increase in the frequency of multi-day rainfall events during active monsoon phases, rather than changes in the intensity of individual rainfall events (Clark et al., 2018;Dey et al., 2020). This implies that the weather systems causing heavy rainfall patterns are developing more often (Clark et al., 2018).
The mechanisms behind increasing monsoon rainfall are still unclear. Recent changes in ENSO are not a contributing factor as the increased frequency of El Niño events in the late twentieth Century should result in fewer wet years and fewer high annual streamflow events (Shi et al., 2008). Nor does the trend follow the thermodynamic increase in rainfall event intensity expected with climate change (Berry et al., 2011;Clark et al., 2018;Smith et al., 2008). Rotstayn et al. (2007) showed that increased concentrations of anthropogenic aerosols can drive temperature and pressure-induced changes in monsoonal winds and increased precipitation in climate modeling, but these findings were later disputed (Shi et al., 2008). Changes in regional sea surface temperatures (Shi et al., 2008), land-ocean temperature differences (Wardle & Smith, 2004), and timing of the monsoon onset (Taschetto & England, 2009), are also proposed mechanisms. It follows that the causes of the recent increasing trend in Daly and Burdekin annual streamflow totals are currently unknown, and further research effort is required.

Conclusions
A long-term perspective on past northern Australian streamflow variability is provided by our tree-ring reconstruction for the Daly catchment, one of the few high-resolution streamflow proxy reconstructions for the IASM region. Despite the relatively short length of the instrumental records used for calibration and verification, little decrease in model statistics was observed moving back in time, indicating that the reconstruction provides useful information throughout the reconstruction period. Our Daly catchment reconstruction extends the record by more than five centuries. The length and robustness of the reconstruction, ability to identify historical flood events and coherence with other proxy reconstructions demonstrate the utility of the variance transform methodology to isolate streamflow signals from noisy tree-ring proxies, particularly when proxies far from the reconstruction target location are used. This method can be applied to reconstructions from other catchments in regions with few local proxies.
The reconstruction shows that high annual streamflow in the Daly catchment is associated with La Niña events, but low streamflow is not associated with El Niño events, confirming previous findings that IASM rainfall responds asymmetrically to ENSO. More generally, warm (cool) SST anomalies near the Australian coastline are associated with high (low) annual streamflow due to enhanced (suppressed) convection.
The recent magnitude and frequency of high streamflow events is unmatched over the past five centuries, regionally coherent, and closely follows observed trends in summer monsoon rainfall. The mechanisms behind the increasing trend in monsoon rainfall, and thus streamflow, are currently unknown. Increasing annual streamflow cannot be directly interpreted as a trend in flood hazard, because the increased frequency of rainfall events rather than increased event intensity lies behind this trend. The short duration high-intensity rainfall events that cause flooding at Katherine can occur in otherwise dry years, and antecedent catchment conditions and tropical cyclone events also contribute to flood occurrence.
Although it is unclear how average IASM rainfall will change with continued global warming, rainfall variability will likely increase. The consecutive failure of two monsoon seasons in 2019 and 2020 demonstrate that despite robust and regionally coherent trends in high streamflow, multi-year dry events still occur. Increased variability could have serious implications for water resources in the Daly catchment, as ecological and social functions rely on dry season baseflow from aquifers recharged during the previous summer's monsoon. Our reconstruction shows that current resource allocations in the Daly have been set during a period of unprecedented high river and aquifer levels which should be carefully considered by water managers when deciding on sustainable future allocations have contributed tree-ring chronologies to the ITRDB, allowing for studies such as this one to be undertaken. Our thanks also to Ze Jiang for his helpful discussions on methodology and to Hung Nguyen and two anonymous reviewers, whose comments have improved this manuscript. PAH is supported by an Australian Government Research Training Scholarship and the UNSW Scientia PhD Scholarship Scheme. MPR is supported by a NOAA Climate and Global Change Fellowship under UCAR CPAESS award #NA18NWS4620043B. FJ is supported by the UNSW Scientia Program. Further support provided by the ARC Centre of Excellence in Australian Biodiversity and Heritage (CE170100015). Open access publishing facilitated by University of New South Wales, as part of the Wiley -University of New South Wales agreement via the Council of Australian University Librarians.