Past and future rainfall change in sub-regions of Victoria, Australia

We examine rainfall variability and change in three sub-regions of the state of Victoria in Australia: the Murray Basin Victoria (MBVic), southeast Victoria (SEVic), and southwest Victoria (SWVic). These sub-regions represent three different hydrological super-catchments over Victoria and received average cool season rainfall for the 1997–2018 period, about 15%, 11%, and 8% less, respectively, than the 1900–1959 average. All three observed declines are shown to be very unusual in terms of historical variability. On analysing CMIP5 models under different forcing conditions (preindustrial, historical-all, historical-GHGs-only, historical-natural-only, RCP2.6, RCP4.5, and RCP8.5), we estimate that external forcing caused 30% of the observed drying in SWVic, 18% in MBVic, and 17% in SEVic. The external forcing contributions to the observed trend for the 1900–2018 period are estimated to be 56%, 17%, and 24% for SWVic, MBVic, and SEVic, respectively. Taken at face value, these figures suggest that only the 1900–2018 trend in SWVic was dominated by external forcing. Nearly all models underestimate the magnitude of the observed drying. This arises because models underestimate the magnitude of decadal variability, and because models might also underestimate externally forced drying, and/or the contribution of internal variability in the real world to the observed event was unusually large. By 2037, approximately 90% of the models simulate drying in SWVic, even under a low emissions scenario. Under a high emission scenario, the anthropogenically forced drying towards the late twenty-first century is so large in all three sub-regions that internal variability appears too small to offset it.


Introduction
Victoria, the second most populous state of Australia, is located on the southeast corner of mainland Australia, within the mid-latitudes of the global climate system with the subtropical ridge (STR) of high pressure to its north (~ 30° S) and the so-called Roaring Forties of westerly circulation to its south (40° S). Victoria's climate is influenced by a large range of climate drivers that operate at different time-scales. Collectively, they help to cause large variability in rainfall, and numerous flooding and several drought episodes over Victoria (Power et al. 1999a,b;DELWP 2020;Timbal et al. 2016;Hope et al. 2017). For example: the Australian Bureau of Meteorology gridded gauge based records which starts in 1900 shows three prolonged periods of below-average rainfall over Victoria: the socalled Federation Drought (1896)(1897)(1898)(1899)(1900)(1901)(1902)(1903)(1904)(1905), the World War II Drought (WWII;1936-1945, and the Millennium Drought (MD;1997-2009, each with different characteristics (Timbal and Fawcett 2013). However, the MD is regarded as the most severe, protracted drought of the instrumental period with a large decline in late autumn and early winter rainfall, whereas the previous droughts exhibited a larger reduction in winter-spring rainfall (Kiem and Verdon-Kidd 2010;CSIRO 2012;Grant et al. 2013;Timbal and Fawcett 2013;Cai et al. 2014;Dey et al. 2019).
In addition to natural climate drivers, Victorian climate is also responding to the increasing level of greenhouse gases (GHGs) in the atmosphere (CSIRO 2012;Cai et al. 2014;Delworth and Zeng 2014;Rauniyar et al. 2019;DELWP et al. 2020;Rauniyar and Power 2020). Previous studies showed that the magnitude of the Victoria-averaged cool season (April to October) drying since the late 1990s would not have been as large without anthropogenic forcing (Cai et al. 2014;Rauniyar and Power 2020). By analysing the global climate models from phase 5 of the Coupled Model Intercomparison Project (CMIP5; Taylor et al. 2012), Rauniyar and Power (2020) showed that there was an observed decline (~ 12%) in cool season over Victoria from an average of 449 mm during 1900-1959 to 394 mm during 1997-2018 period. They also estimated that about 20% of this decline (interquartile range (IQR): 40% to − 4%) was simply explainable by external forcing, suggesting that the recent decline, although largely driven by internal climate variability, was reinforced by external forcing to make it unprecedented (please see Rauniyar and Power (2020) for discussion on related studies focusing on attribution of regional changes in rainfall beyond Australia).
Looking into the future, various studies (e.g., CSIRO 2012; Timbal et al. 2010;Cai et al. 2014;Grose et al. 2015;Hope et al. 2015;Rauniyar and Power 2020) suggest that anthropogenic forcing is projected to strongly influence future rainfall in Victoria. For example, Rauniyar and Power (2020) concluded that there is an 88% chance that the next two decades will be drier than the pre-industrial average, irrespective of the future emission pathway. More recently, , exploiting the recent research of Rauniyar and Power (2020), developed a novel, model-based method to estimate future rainfall.  showed the multi-model median (anthropogenically forced) rainfall decrease monotonically over the remainder of the twenty-first century under a high greenhouse gas emissions scenario, called Representative Concentration Pathway 8.5 (RCP8.5;van Vuuren et al. 2011). They estimated that the probability of receiving very low rainfall below the observed 5 th percentile (using 1900-2018 as the reference period) becomes three times larger by the end of the century under RCP8.5. They also estimated that there is a 10% chance that All-Victoria rainfall could go beyond historical experience in any given year from 2025 onwards. These findings were used by the Victorian Government's Department of Environment, Land, Water and Planning (DELWP) in developing the current "Guidelines for Assessing the Impact of Climate Change on Water Supplies in Victoria" (i.e., DELWP 2020) to inform decisions it had to make for future water resources planning and management, including the baseline (historical reference climate) choices (DELWP et al. 2020).
Many of the previous studies on rainfall changes over Victoria were derived using all-Victoria area-averaged rainfall (e.g., Power 2020, 2022). This is not ideal because both long-term averages and externally forced signals are known to vary regionally Timbal et al. 2016). It would have been highly desirable to have the analysis done for the ten Catchment Management Authorities regions across Victoria as they are responsible for the integrated planning and coordination of land, water, and biodiversity management in each catchment and land protection regions. However, this is precluded due to coarse resolution of CMIP5 climate models, which means they may not have grid points in each region and will also have difficulty in capturing the physical features such as topography which differentiate rainfall at such scales. In addition, many decisionmakers would prefer information at a smaller scale. While this information can be estimated using downscaling tools, here we wish to learn as much as we can about the regional contrasts evident in the relatively coarse grid climate models. On the other hand, we know that confidence in models at small spatial scales is lower. So here we make a compromise and divide Victoria up into three sub-regions ( Fig. 1): the Murray Basin Victoria (hereafter, MBVic), southeast Victoria (hereafter, SEVic), and southwest Victoria (hereafter, Fig. 1 Spatial distribution of cool (April-October) season (a) rainfall climatology for the 1900-1959 period, (b) relative change, and (c) percentage change in rainfall for the 1997-2018 period relative to climatology shown in (a) and (d) the likelihood of occurring the 1997-2018 period observed change due to a random variability computed using bootstrapping method described in Sect. 2. The three sub-regions used in this study are: the Murray Darling Basin Victoria (MBVic), southeast Victoria (SEVic), and southwest Victoria (SWVic). The 500-m elevation contour is shown in red colour. Dotted grids represent the regions where the rainfall change is not statistically significant at 10% significance 92 Page 4 of 23 SWVic). These three sub-regions are meaningful from a climatic, vegetative, and hydrological perspective as they represent three different hydrological super catchments (Timbal et al. 2016) and also broadly align with the three climate regions identified in the study of Fiddes et al. (2021).
The rainfall across these sub-regions is influenced by fronts embedded in the largescale westerly flow, low pressure systems, particularly on the east coast, tropical influences modulated by ENSO and IOD, and interactions with the polar and sub-tropical jets (McKay et al. 2022;Murphy and Timbal 2008;Timbal et al. 2016). The relative importance of these weather features varies from sub-region to sub-region (Fiddes et al. 2021) and it is possible that the influence of climate change might also vary spatially. For instance, SWVic primarily receives much of its rainfall from the interaction between the frontal systems embedded in the large-scale westerly flow and the mountainous topography Hope et al. 2017). Rainfall from fronts has declined , and future rainfall changes will be influenced by changes in fronts and the rainfall from them (e.g., Blazquez and Solman 2019). SEVic and MBVic receive rainfall from thunderstorms that occur during the warm season, but the thunderstorm season is projected to lengthen (Allen et al. 2014;Singh et al. 2017). East coast lows are also particularly important in SEVic (Pepler et al. 2014, and it is as yet unclear how these will change in future. Given that each region's rainfall has different influences, and it is known that the influence of climate change varies spatially (e.g., Eyring et al. 2021). So, here we determine the extent to which forced changes vary from one sub-region to another.
The main objective of this study is to investigate past and future rainfall change in three sub-regions of Victoria due to climate change. Following the methods developed by Rauniyar and Power (2020), we will specifically address the following questions: • Are the characteristics (e.g., mean and distribution) of the cool season rainfall in each sub-region for the recent period  statistically different compared to any other multi-decadal period's characteristics? • What is the fractional contribution of externally forced climate change to the observed decline in cool season rainfall for the 1997-2018 period and does the response to external forcing vary from one sub-region to another? • When does the externally forced signal expel rainfall averaged over each sub-region beyond its pre-industrial and historical range (Power 2014) and what is the expected combined impact of both external forcing and internal variability over coming decades in each sub-region?
These questions have not been fully addressed previously for the selected sub-regions over Victoria, aside from question 1, that was partially addressed by Timbal et al. (2016), although they did not assess the recent cool-season rainfall decline in terms of its own stochastic variability and/or modelled internal climate variability. The 1997-2018 period is chosen in this study as there is a strong interest among stakeholders, both within and beyond the climate science community in Australia, in the role of climate change in this period as it encompasses the most severe, protracted drought (i.e., Millennium Drought) of the instrumental period. This also enables us to compare and contrast the results of this study with Rauniyar and Power (2020) who quantified the contribution of externally forced climate change in the decline of Victoria-averaged cool season rainfall using the same period. The paper is organized as follows: a brief description of data, climate models, and methods used is provided in Section. 2. Results are presented in Section. 3 including the results for the model evaluation. Section 4 summarizes the major findings, discusses possible explanations for the observed results, and suggests future work.

Data, climate models, and method
This study utilizes the same set of rainfall data used in Rauniyar and Power (2020) for both observations and CMIP5 models simulations. This facilitates a direct comparison between the results for the sub-regions and all-Victoria results. Specifically, we use the monthly resolution gridded rainfall observation for the period 1900-2018 from the Australian Bureau of Meteorology under the Australian Water Availability Project (AWAP; Jones et al. 2009) to analyse the past climate variability and changes over the sub-regions. Similarly, we use monthly rainfall simulations from the CMIP5 climate models under a variety of forcing conditions (Taylor et al. 2012). Unless otherwise stated, we only use those models that provide at least 500 years of rainfall simulations under preindustrial (piCTL) conditions to ensure a good estimate of the preanthropogenic climate (Delworth and Zeng 2014;Rauniyar and Power 2020) as the various forcing agents (e.g., concentrations of GHGs, aerosols, ozone, and solar irradiance) are prescribed (fixed) at the preindustrial level (year 1850). We use three different historical runs (i.e., all-forcing: histALL, GHG-only: histGHG, natural-only: histNAT) for the period 1900-2005 to examine how rainfall across Victoria responds to the different observed atmospheric forcings during the twentieth century. Furthermore, we use rainfall projections for the 2006-2100 period under three different future emissions scenarios (low, medium, and high), represented by three RCPs (i.e., RCP2.6, RCP4.5, and RCP8.5) to project changes in multi-decadal rainfall across the Victorian sub-regions during the remainder of the twenty-first century.
Each model shown in Supplementary Table 1 has different numbers of simulations, but we used only the first run (i.e., r1i1p1) in the first instance, to avoid weighting models with more than one simulation. In addition, not all models have rainfall simulation across all scenarios/ experiments. Therefore, most of our analysis are based on 24 models that have in common simulated rainfall across the piCTL, histALL, RCP8.5, and RCP4.5 scenarios to be consistent in our results across the different runs. For the low emission scenario (RCP2.6), we have only 19 models that provide simulations across the piCTL, histALL, and RCP2.6 and the number is further reduced to 13 for the piCTL, histALL, histGHG, and histNAT scenarios. However, we also examined the robustness of the results based on the first run only against the results of all available simulations and noted the similarities and differences where applicable. Prior to computing the area-averaged rainfall time-series for the three different sub-regions, all the data (observed and model) are interpolated from their native spatial resolutions to a common 1.5° × 1.5° grid resolution using a conservative interpolation method (Jones 1999;Rauniyar et al. 2017).
We applied a two-sided Kolmogorov-Smirnov (KS) test on the relative frequency distributions (RFDs) of cool season monthly mean rainfall for different 20-year blocks from 1900 to examine whether the RFD of 1997-2018 period is similar or different compared to RFD of any other multi-decadal periods or the long-term climatology of the twentieth century. Furthermore, to determine the occurrence probability of the recent decline in rainfall in terms of its own variability, we performed Bootstrap resampling (Efron and Tibshirani 1994). This is done by randomly selecting individual years (with replacement) from the full historical period of area-averaged cool season rainfall (119 years: 1900-2018) to create 119 samples. This is done 10,000 times. Then, the difference between the average rainfall of the last 22 years (i.e., 1997-2018 average) and the average of the first 60 years (i.e., 1900-1959 average) was computed each time. Here, the first six decades was chosen to minimize the influence of climate change, if any, on rainfall variability of reference period. These 10,000 samples of mean rainfall differences are used to estimate the probability distribution function (PDF). We repeated these steps to compute PDFs for the percentage differences and for difference scaled by the inter-annual standard deviation of cool season rainfall variability over the first 60 years. Finally, the actual observed change was compared against the distribution to estimate the likelihood of occurrence relative to its own statistical variability. This is justified because internal variability in regional rainfall does not exhibit statistically significant persistence in the sub-regionally averaged cool season rainfall (not shown).
Since CMIP5 model runs under historical simulations end at year 2005, we extended the histALL simulations by augmenting the corresponding model RCP8.5 data over the period 2006-2018 to evaluate the model's performance against the full observational period statistics (see Supplementary Material). We performed Fourier decomposition (harmonic analysis) on the mean annual cycle of observed and simulated rainfall to examine how well the models simulate the observed seasonal rainfall pattern. We also examined the consistency in models by comparing their trends against the observed trend for the 1900-2018 period (see Supplementary Material). We computed the multi-decadal rainfall changes relative to the statistics of 500 years long piCTL simulations for every 20-year block under piCTL runs and histALL simulations combined with the RCP scenarios. Estimates of 600 different multi-decadal rainfall change situations under piCTL condition from all 24 models (i.e., 24 models times 25 numbers of 20-year blocks in 500 years of piCTL) is then used to estimate the probability distribution that can arise due to internal variability only. Finally, the multi-decadal changes under historical and future experiments are compared against the piCTL PDF to approximate the timing of when the signal due to external forcing could fully expels the internal variability and also to estimate the relative proportion of the rainfall decline that was due to external forcing.

Observed characteristics
In general, mountainous areas in eastern Victoria receive the highest amount of rainfall during the cool season, followed by the southern coastal regions and substantially less rainfall occurs over MBVic (Fig. 1a). These distinct variations in rainfall arise from the fact that different sub-regions of Victoria are influenced by different combination of weather systems and large-scale climate drivers (Timbal et al. 2016;Pepler et al. 2020;Fiddes et al. 2021). For instance, the SWVic receives much of its rainfall from frontal systems embedded in the large-scale westerly flow Hope et al. 2017), whereas most of the rain over the SEVic comes from thunderstorm and east coast lows, particularly during autumn and winter seasons (Pook et al. 2014;Pepler et al. 2014Pepler et al. , 2021Dowdy and Catto 2017).
The spatial anomaly map (Fig. 1b) shows that the largest decline in cool season rainfall (> 90 mm) is concentrated over the highland areas where maximum rainfall also occurs (Fig. 1a). The decreases in rainfall over other sub-regions are comparatively small but show a similar magnitude (< 60 mm). However, the percentage decrease over MBVic is largest (Fig. 1c) and it also receives the lowest amount of total rainfall (Fig. 1a). Furthermore, the changes in rainfall are statistically significant at the 90% level over many catchments of Victoria, except over the central plain in SWVic and over the eastern parts in SEVic (dotted grids in Fig. 1b and c). Figure 1d shows that the probability of occurrence of such a low amount of rainfall in recent decades is less than 2% due to a random chance over the catchments where the declines are statistically significant.
Rainfall across Victoria also exhibits large temporal variations at inter-annual through to inter-decadal time-scales as shown in the left column of Fig. 2. However, it is notable that since the beginning of the MD in 1997, all three Victorian sub-regions have consistently received less rainfall during the cool season-i.e., the rainfall amounts have not "swung back". For the 1997-2018 period, the MBVic, SEVic, and SWVic received on average approximately about 15%, 11%, and 8% less rainfall than the reference period (i.e., 1900-1959) averages (Table 1). We find that the probability of occurrence of such low rainfall in recent decades relative to its own statistical variability is less than 1%, except for SWVic which shows a probability of less than 3%. Similarly, all three sub-regions exhibit  Table 1 Key statistics of area-averaged observed cool season (April-October) rainfall in the MBVic, SEVic, and SWVic regions for different historical periods used in this study. Means, standard deviations, and linear trends are in mm year −1 . Rainfall change (%) for the period 1997-2018 relative to reference period 1900-1959 is given in 4 th column and probability of occurrence of such change against its own statistical variability is given in 5 th column. Bold values are statistically significant at the 95% level. The likelihood of occruence of the observed magnitude of the recent decline by chance under all-forcings is computed by comparing the recent decline against the probability distribution of the 10,000 randomized differences which are generated by random resampling of the full historical records of cool season rainfall 10,000 times and each time the mean rainfall difference between the most recent 22 years period and the mean rainfall over the first six decades is computed. The occurrence probablility of observed change in rainfall without the contribution of externally forced drying (i.e., distribution of resample diffrences is shifted left by the multi-model mean value) is given in 6 th column.
The magnitude of externally forced drying for the period 1997-2018 is given in the last column of negative trends; however, they become statistically significant at the 95% level for timeseries segments beginning after 1950 but not over the full historical record. This is consistent with the previous studies that found the impact of climate change in Victoria is noticeable only after mid-twentieth century (Timbal et al. 2016;Hope et al. 2017;Jones and Ricketts 2017). Furthermore, the result shows both the median rainfall values (shown as vertical lines in varying colors in Fig. 2 right column), and the RFDs of different 20-year periods during the twentieth century are not statistically different at the 95% level with each other including the distribution of the twentieth century in all three sub-regions based on a two-sample KStest. In contrast, the distribution of 1997-2018 period differs substantially (p-value < 0.05, using the same test) from the rest of the multi-decadal distributions, including the distribution of twentieth century in MBVic and SEVic. However, the post-1997 period distribution for SWVic is statistically different only from the distributions of the 1960-1979 period and the distribution of the twentieth century. Nevertheless, the distributions of cool season rainfall for 1997-2018 in all three sub-regions have shifted to the left (i.e., to a mean lower monthly rainfall value) compared to the rest of the distributions and has experienced fewer very wet months during the cool season since 1997 compared to any other periods discussed here. These findings show that recent climate has been quite different (significantly drier) from earlier in the Century, suggestive of an emerging human fingerprint on Victoria rainfall. This is being examined in following sub-sections using climate models simulations under different forcings.

Past and future multi-decadal changes in simulated rainfall
In this section, we analyse changes in 20-year average simulated rainfall relative to 500 year-long piCTL simulations during the twentieth and twenty-first centuries (i.e., from 1900-1919 through to 2100) in all three sub-regions under historical and different emission pathways against the distribution of possible range of multi-decadal rainfall changes that can arise due to internal variability alone (Fig. 3). Multi-model median (MMM) values are selected instead of multi-model mean to reduce the influence of extreme outliers. The emergence of climate change signal is identified as the first occurrence of the 20-year period that shows the MMM change greater than one standard deviation of the preindustrial distribution and has at least 75% or more models agree with the sign of MMM change. However, before computing any statistics, we evaluated the CMIP5 models' ability to simulate the observed mean annual cycle of rainfall for the period 1900-2018 over the three sub-regions of Victoria to check if they are fit for purpose. The result showed that the majority of models were able to simulate the amplitude and timing of seasonal rainfall pattern reasonably well, but a few of them (CCSM4, CESM1-BGC, and NorESM1-M) exhibited a reversed annual cycle over MBVic and SEVic (see Supplementary Material). However, even after excluding the models with reversed annual cycle, we did not find any substantial differences in the magnitude of changes shown in Table 2 and Supplementary  Table 2. Furthermore, the trends in these models for the 1900-2018 period are found to be consistent with the observed trends in all three sub-regions (see Supplementary Material). Therefore, we included all the models in our analysis to increase the sample size. Figure 3 shows that the climate change signal is small during the twentieth century in all three sub-regions as MMM values are either close to zero or within one standard deviation of piCTL distribution's mean. Furthermore, there is very little agreement among models on the sign of the MMM change, specifically in MBVic and SEVic ( Fig. 3a and b). The externally forced drying becomes evident from 2010-2029 period in all three sub-regions but is very clearly visible in SWVic where more than 90% of models agree with the MMM sign change (Fig. 3c). Comparison of the boxplot against the distribution due to internal variability shows that substantial drying appears very strong after mid of twenty-first century under RCP8.5 in all three sub-regions. For example, the MMM change in all three sub-regions under RCP8.5 is more than four times the standard deviation of piCTL variability by the end of twenty-first century. Furthermore, it is also clear that the impact of emission reduction (RCP8.5 vs RCP2.6) is only visible from the mid-twenty-first century in terms of median change, since the MMM changes are similar under all the three RCP scenarios prior to 2060. This implies that substantial drying across Victoria may be inevitable after 2060 under high emission scenario because internal variability is, by comparison, too weak to offset the projected drying in decadal periods. Each boxplot represents the distribution of changes based on 24 models that are common between the preindustrial, historical and RCP8.5 simulations. The vertical line in the box indicates the median, the shaded box represents the 25th and 75th percentiles, and the whiskers indicate the minimum and the maximum values based on 24 CMIP5 models. The median values for RCP4.5 and RCP2.6 scenarios are overlaid on the boxplots as blue crosses and green stars. Proportions of models that exhibit the same sign as the median are shown in the square brackets in the order for RCP2.6, RCP4.5, and RCP8.5 scenarios. The color shading represents nothing more than time progression. The possible range of changes in mean rainfall in 20-year blocks in 500 years of piCTL runs is shown as the relative frequency distribution in gray shading. The recent decades  is highlighted using the horizontal dashed gray lines Throughout the remainder of twenty-first century (i.e., from now to 2100), the projected percentage decline under high and medium emission scenarios (RCP8.5 and RCP4.5) are largest in SWVic followed by MBVic and SEVic (Supplementary  Table 2). However, the percentage decline under low emission scenario (RCP2.6) is largest for MBVic towards the end of twenty-first century and is smallest for SWVic. The results show that SWVic, under RCP8.5, could experience a rainfall decline up to − 20% (IQR: − 27% and − 9%) by the end of twenty-first century relative to pre-industrial average. For the same period and scenario, the projected decline could be − 18% (IQR: − 30% and − 5%) and − 14% (IQR: − 24% and − 5%) for MBVic and SEVic. Furthermore, SWVic also exhibits the greatest agreement amongst the models (over 90%) to the MMM sign change while the same is slightly lower (83%) in MBVic and the lowest (79%) in SEVic. In addition, some of the models project increased rainfall in SEVic (i.e., beyond the range of internal variability) which is not the case for other two sub-regions (see Fig. 3b). Under low emission scenario, the projected decline could be almost half of high emission scenario for MBVic (MMM: − 9%; IQR: − 10% and − 4%) and SEVic (MMM: − 7%; IQR: − 10% and − 2%) regions while the same could be much less (only quarter) compared to high emission scenario in SWVic (MMM: − 5%; IQR: − 8% and − 2%).
The MMM value for near future (2018-2037) shows a change in rainfall of about − 5% with IQR between − 7% and − 2% across Victoria. However, our estimates show that there is only 7%, 11%, and 14% chance that externally forced drying averaged over this period could be completely masked by internal variability over SWVic, SEVic, Table 2 Contribution of anthropogenic forcing to observed drying for the period 1997-2018 relative to 1900-1959. In each standardization method, the first, second, and third rows are for the MBVic, SEVic, and SWVic regions respectively. Change (Δ), % change (100*Δ/µ), and standardized changes relative to internal standard deviation (Δ/σ) and decadal standard deviation (Δ/σ d ) are considered. Here Δ is the difference in area-averaged mean rainfall between 1997-2018 and 1900-1959, µ is the area-averaged rainfall over 1900-1959 period, and σ and σ d are the inter-annual and inter-decadal standard deviations of area-averaged rainfall over the same period. Observed change (ΔOBS) is shown in the 2nd column. The "Weighted mean" change (ΔMOD) at the 50 th percentile is the average of the multi-model median change (MMM) for all three RCPs, weighted by the number of models (i.e., 40, 38, 28, respectively). Similarly, the 25 th and 75 th percentile values are used to compute the weighted mean changes at the 25 th and 75 th percentiles. Best estimate of the relative contribution of anthropogenic forcing to the observed drying is shown in the last column in bold type as a percentage, together with IQR values and MBVic. This implies, according to models, that there exist a higher chance (93%) that the 2018-2037 period will be drier than the preindustrial average over SWVic compared to other two sub-regions. For the more recent period , the MMM changes are approximately -2% for SWVic, and are less than − 1% for other two subregions which suggests that externally forcing only makes a small contribution to the observed decline. However, our confidence in these numbers is lowered as the models tend to underestimate the magnitude of internal variability (see Section. 3.5 and discussion for details).

Role of greenhouse gases (GHGs) and natural processes.
In this section, we examine the influence of increasing levels of GHGs in the atmosphere and other natural forcings on rainfall across Victoria by analysing three different historical runs: histALL, histGHG, and histNAT. The smoothed (20-year running mean) version of percentage changes in rainfall relative to the piCTL long-term mean is shown in Fig. 4. As all historical simulations in CMIP5 models end in year 2005, the results shown are more representative until 1995 only (i.e., average of 1986-2005). Nevertheless, 10 out of 13 models simulate decreases in rainfall post-1970 under histGHG simulations with MMM value in all three sub-regions being significantly below zero at the 95% significance level towards the end of twentieth century (third row in Fig. 4). This implies that the increased concentration of GHG emission acting alone would have likely reduced rainfall across all sub-regions of Victoria. By contrast, no robust trend exists in rainfall under the histNAT simulation (i.e., including only volcanoes and solar forcing). It is likely that the drying seen towards the late twentieth century under GHGs simulation could have been masked out by natural processes, model internal variability, and aerosols. We examined this by removing (linearly) the impact from both histGHG and histNAT forcings from the rainfall response under histALL simulation. The result for histOTH (historical-other) forcings indeed shows an increase in MMM value towards the end of the twentieth century (last row of Fig. 4). As a result, the histALL forcing simulations do not show any significant decrease in rainfall towards the end of the twentieth century (first row of Fig. 4). However, we note that this might merely reflect sampling error arising from internal variability and further research is needed to clarify the relative importance and impact of the various individual forcings (e.g., sulphate-aerosols, ozone) to the observed changes. Table 2 shows that the observed declines in rainfall for 1997-2018 in all three sub-regions are significantly large compared to the inter-decadal variability of 1900-2018, which indicates that the recent decline is very large compared to internal variability. To estimate the probability of occurrence of the observed change due to internal variability alone, we used piCTL models with at least 200 years of simulations. For each piCTL model, we randomly resampled individual years from the 200-years of data 1000 times and computed the percentage difference between the average of the last 22 years and the average of the first six decades. This results into a sample of 35,000 values (35 models times 1000) of percentage rainfall changes which is used to estimate the PDF of multi-decadal rainfall changes due to internal climate variability alone as shown in Fig. 5 (left column). The results show that the magnitude of the recent observed drying over all three sub-regions are unusually large compared to internally generated multi-decadal rainfall variability in the models as the chance of such drying from internal variability alone is less than 0.5% for MBVic and SEVic and less than 2% for the SWVic. Figure 5 (left column) also shows the percentage changes in mean rainfall of the 1997-2018 period relative to the 1900-1959 period average for all the models under historical forcing and the three RCP scenarios. Many models forced under all three RCP scenarios struggle to capture the observed decline, with some exceptions for SWVic where only a few models able to simulate the magnitude of observed change. However, when the standardized data (i.e., difference scaled either by inter-annual or inter-decadal standard deviations) is considered, a larger number of models (depending on scenario) simulate dry conditions similar to the observed difference (not shown). This suggests that the models underestimate variability: this will be examined below.

Contribution of external forcing in recent rainfall decline
Following the methods described in Rauniyar and Power (2020), the magnitude of the externally forced drying in recent decades in models is estimated by averaging the MMM values of all three emission scenarios but each RCP is weighted by the number of models  Rainfall was randomized 1000 times to compute the change in mean rainfall of the recent 22 years period from the first six decades. The curves on the right represent the RFDs of 119-year trends arising from internal variability only, evident in models that have at least 200 years of rainfall simulations under pre-industrial conditions. The dashed bell curve is the same as the shaded curve, except that it is shifted left by the MMM value (i.e., externally forced response, see Table 2 and main text). The vertical dot-dash lines in the panels on the left indicate no change in average rainfall between the most recent 22 years and the first 60 years, while the dot-dash lines in panels on the right indicate zero 119-year trends. The observed % change and long-term trend are indicated using thick, black vertical dashed lines used under that scenario. Finally, the model "mean" estimate of contribution of externally forced drying in recent decades is computed by determining the proportional contribution of the averaged-MMM value to the observed change (Table 2). We did the same using the IQR values of three RCPs to capture the uncertainness around the best estimate. The results show that the percentage contribution of external forcing, according to models, is largest (MMM: 30% and IQR: 66 to 3%) in the SWVic followed by MBVic (MMM: 18% and IQR: 37% to − 9%) and smaller again (MMM: 16% and IQR: 39% to − 7%) in SEVic. No statistically significant differences are seen in the best estimates (i.e., MMM) of the contribution of external forcing in any of the three sub-regions even after considering all the available ensemble members (Supplementary Table 3 and Fig. 3). However, the IQR is found to be smaller in all sub-regions when the ensemble averages of models are being used.
Furthermore, the fractional contribution of anthropogenic forcing slightly increases when the standardized rainfall decline is considered which indicates that the model simulations underestimate the observed rainfall variability (Table 2, Supplementary Fig. 4). Nevertheless, irrespective of the standardization methods used, the magnitude of proportional contribution of external forcing suggests that the observed decline in rainfall, according to models, is largely dominated by internal climate variability in all three sub-regions. To examine whether the differences in the response to external forcing are statistically robust between the three sub-regions, we computed the RFD for each region separately using all the simulated differences (RCP8.5, RCP4.5, and RCP2.6). We found that the distribution for SWVic is statistically significant at the 95% level (shifted left) compared to other two sub-regions based on a two-sided KS test (not shown). This suggests that SWVic will dry more than the other two sub-regions, which is also evident in the projected change (Supplementary Table 2). Furthermore, we also computed the likelihood of occurrence of the observed decline, but with externally forced drying component removed (i.e., distribution of resampled observed differences is shifted left by MMM value) and found that the probability of occurrence is around 2% for MBVic and SEVic and about 9% for SWVic (6 th column in Table 1) under the assumption that no externally forced change has occurred. This again points to natural variability being an unlikely driver for the "residual" downturn in the three sub-regions, especially in MBVic and SEVic.
Up to now, we have considered multi-decadal changes. It is also of interest to know the degree to which observed and simulated trends over the period 1900-2018 ( Fig. 5 right column) agree. We find that the observed and simulated trends are relatively more consistent than the multi-decadal changes (Fig. 5 left column). The probability of occurrence of the observed trend due to the modelled internal variability alone for MBVic, SEVic, and SWVic is found to be 1.8%, 1.5%, and 15%, respectively. The likelihood is estimated from the distribution of all 119-year trends arising from internal variability evident in models with at least 200 years of pre-industrial simulation. We estimate the contribution of external forcing to the observed trend for MBVic, SEVic, and SWVic to be 14%, 19%, and 51%, respectively. These values further increase to 17%, 24%, and 56% for MBVic, SEVic, and SWVic, respectively, when all available ensemble members shown in Supp. Table 1 are considered (Supplementary Section 1.2 and Supplementary Fig. 2).

Projected multi-decadal changes relative to observed droughts
To better plan for future water availability across Victoria, it is prudent to examine how the projected multi-decadal rainfall changes for the twenty-first century under all three RCPs stand against the observed major historical rainfall changes, i.e., the World War II and Millennium Droughts and during the recent 22 years (i.e., 1997-2018). This cannot be done using preindustrial climate as a reference as historical records do not go that far back. Hence, we used the 1900-1959 period as a reference to compute the percentage change and results are shown in Fig. 6. MMM values during all periods from 1900 to 2100 under all emission pathways are smaller in magnitude than the drying that was observed during the MD, except in SWVic where MMM drying is slightly higher in magnitude from the 2060-2079 period under RCP8.5. This further illustrates-according to the models-that the internal variability was the dominant driver for the MD drying across Victoria.
Nevertheless, the projected MMM dryings in the late twenty-first century under a high emission pathway in all three sub-regions are larger than that observed for the WWII drought and approximately of the same magnitude of the drying observed during the last 22 years. Furthermore, the late twenty-first century projected MMM drying  Fig. 4, the distribution of changes in 24 models under RCP8.5 are represented as boxplots for each 20-year period. The horizontal line in the box indicates the median, the shaded box represents the inter-quartile range (IQR: 25 th and 75 th percentiles), and the whiskers indicate the minimum and the maximum values based on 24 CMIP5 models. The median values for RCP4.5 and RCP2.6 scenarios are overlaid on the boxplots as blue and green circles, with corresponding IQRs represented by the blue and green vertical lines, respectively. The color shading represents nothing more than time progression is smaller under RCP4.5 relative to RCP8.5, but it is still comparable to the WWII drought. Note that these MMM values are estimates of the impact of external forcing only. Changes are so large under RCP8.5 that any significant exacerbation by internal variability would result in levels that rival or exceed the MD drying. On the other hand, even significant offsetting by variability still results in conditions experienced during the WWII drought. Offsets larger than this appear unlikely. This suggests a significant increase in risk of experiencing events like the WWII and recent droughts under global warming.
However, we have only limited confidence in the modelled projections across Victoria because overwhelming majority of models underestimate the amplitude of the annual cycle of rainfall (see Supplementary Material) and the magnitude of the observed decline in rainfall in all three sub-regions (see Fig. 5). The reasons for these last underestimation-which are not mutually exclusive-are similar to those outlined by Rauniyar and Power (2020) and Power et al. (2017). More specifically, they are: (i) the observed drying is dominated by an unusual and large natural, internally generated event; (ii) drying in models in response to greenhouse gas forcing is underestimated; (iii) models underestimate the variability; and/or (iv) error in instrumental records. Figure 7 shows that the models do appear to underestimate the inter-annual and interdecadal variability in all three sub-regions. However, compared to other two subregions, the number of models underestimating the observed variability is less in SWVic. Nevertheless, the majority of models still tend to underestimate the observed change even if changes standardized with decadal variability are examined (Table 2 and Supplementary Fig. 4). So, while the weaker than observed decadal variability in models contributes to the apparent underestimate, (i) and/or (ii) might also play a role.

Summary and discussion
Many previous analyses used all-Victoria area-averaged rainfall to study the characteristics of rainfall changes, including detection and attribution of a climate change signal in recent cool season (e.g., April to October) rainfall decline over Victoria (e.g., Rauniyar and Power 2020). But rainfall across Victoria displays a strong climatological gradient (Fig. 1a) due to a complex interactions of various large-scale climate drivers and external forcing which modulates the weather systems differently over these sub-regions in the presence of contrasting orography. Therefore, the findings may not be directly applicable to catchments or sub-regional levels, which is very crucial for planning and management of water resources.
In this study, we analysed both observations and a suite of CMIP5 climate models under different forcings to investigate past, and future changes in cool season rainfall due to climate change over three sub-regions of Victoria: the Murray Basin Victoria (MBVic), southeast Victoria (SEVic), and southwest Victoria (SWVic). The division to these three sub-regions is based upon a contiguous response to weather systems, rather than simply by a political boundary. These sub-regions also represent three different super hydrological catchments and are meaningful from a climatic, vegetative, and hydrological perspective. The key characteristics of the observed cool season rainfall across the three sub-regions, how they have been influenced by the external forcing in recent decades, and how rainfall may change in future under different greenhouse gas emission pathways are summarized as below. Fig. 7 Comparison between the the observed and model-simulated cool-season rainfall standard deviations (mm month −1 ) for the period 1900-2018 for (a) MBVic, (b) SEVic, and (c) SWVic regions. Inter-annual and decadal standard deviations are shown along the X-axis and Y-axis. Crossing of dashed dark gray lines with the black filled circle represents the observations and the light gray lines represent the 5% and 95% of confidence intervals of the standard deviation using bootstrapping. The histALL simulations extended with the RCP8.5, RCP4.5, and RCP2.6 scenarios are respectively shown in red stars, blue crosses, and green plus symbols

Observations
In general, Victoria receives the highest amount of rainfall over its mountainous areas in eastern Victoria (i.e., SEVic), followed by the southern coastal regions (i.e., SWVic) and substantially less rainfall occurs over the north-west of Victoria (i.e., MBVic). However, since the beginning of the Millennium Drought (MD) in 1997, cool season rainfall is declined significantly by historical standards (p-value < 0.05) over most Victoria (Fig. 1d). The rainfall averages over the MBVic, SEVic, and SWVic for the 1997-2018 period are respectively, about 15%, 11%, and 8% less than the reference period (i.e., 1900-1959) average. The likelihood that such large declines could arise by chance alone are estimated as less than 1% for MBVic and SEVic while the same is less than 3% for SWVic. Furthermore, the distributions of cool season monthly mean rainfall for 1997-2018 period in all three sub-regions have found to be narrower and shifted to the left compared to the distributions of the twentieth century.

Climate model results
To understand the role of external forcing in the observed drying and to estimate the combined impact of external forcing and internal variability on future rainfall across Victoria, 24 CMIP5 global climate models runs under preindustrial, three different historical conditions (GHGs-only, natural-only, and all forcings), as well as three scenarios for the twentyfirst century: RCP2.6, RCP4.5, and RCP8.5 were analysed. We considered coarse resolution GCM results directly to avoid the uncertainties of downscaling.

Historical simulations
A majority of models simulated the timing of seasonal rainfall pattern reasonably well over SWVic, but a few of them (CCSM4, CESM1-BGC, and NorESM1-M) exhibited a reversed annual cycle over MBVic and SEVic. However, removing these models from analysis did not change the results substantially. On comparing models' historical rainfall simulations under different forcings, we found that the multi-model median (MMM) rainfall changes of GHGs-only runs are well below zero in all three sub-regions and a majority of models with GHG-only runs showed a clear decreasing trend from the mid-1960s. By contrast, the MMM rainfall changes of natural-only runs are close to zero with no clear trend in any direction towards the late twentieth century. These findings suggest that increasing levels of GHGs in the absence of other forcing changes generally lead to a decline in rainfall across Victoria. However, the signal is reduced or removed in models by natural variability as well as by aerosols in the historical all-forcings simulations. Further research is needed to clarify the relative importance and impact of the various individual forcings to the observed changes.
The multi-decadal rainfall changes in all three sub-regions are found to behave similarly, although the magnitudes vary. During the twentieth century, the MMM changes (relative to preindustrial average) are close to zero in all three sub-regions, indicating a very weak climate change rainfall signal across Victoria (at least compared with natural variability). For the 1997-2018 period, the majority of models simulated rainfall below the preindustrial average in all three sub-regions, regardless of emission pathways. However, the MMM drying is greater in SWVic compared to other sub-regions. By contrast, the observed rainfall decline in SWVic is smaller compared to other two sub-regions during this period (Table 2). According to the models, the externally forced contribution to the drying observed in the 1997-2018, relative to 1900-1959 is larger (MMM: 30% with IQR: 66 to 3%) in SWVic, followed by MBVic (MMM: 18% with IQR: 37 to − 9%) and then SEVic (MMM: 16% with IQR: 39 to − 7%). Noticeable differences are not observed in the MMM values of the first run of models vs all available simulations across all three subregions, except that the associated uncertainties (i.e., the IQRs) are found to be substantially reduced. Based on the multi-model ensemble means, the externally forced contributions to the observed long-term (1900-2018) drying trends for SWVic, MBVic, and SEVic are 56%, 17%, and 24%, respectively.
The changes from 1900-1959 to 1997-2018 are largely dominated by internal variability across Victoria and not by the external forcing. This is consistent with the conclusion of Rauniyar and Power (2020) for Victoria as a whole. We found that nearly all models underestimated the observed drying. We showed that part of the reason for this, perhaps the main reason, is that the models tend to underestimate the magnitude of decadal rainfall variability in all three sub-regions. In fact, we found that a larger number of models simulated drying as great or greater than observed if the observed change and simulated change were scaled by their respective levels of decadal variability ( Supplementary Fig. 4). However, most models still underestimated the drying. It therefore seems that the observed drying is dominated by a large natural, internally generated event and/or that the drying in models in response to greenhouse gas forcing is underestimated. Whatever the case, the above estimates of the external forcing contribution should be reassessed using the CMIP6 models. This conclusion is also true for long-term trends in MBVic and SEVic, but not in SWVic where the trend is dominated by external forcing.

Future rainfall
According to the model simulations analysed, the externally forced drying signal (estimated using the MMM change) is greater than one standard deviation of internal variability by 2029, and over 90% of models exhibit drying under all three scenarios. This indicates that the near-future (2018-2037) will be drier than the preindustrial average in all three sub-regions. We estimate that there is only a 7%, 11%, and 14% chance that internal rainfall variability will completely override externally forced drying averaged over SWVic, SEVic, and MBVic, respectively.
Under a high emission pathway, the externally forced change in all three sub-regions is so large by the late twenty-first century (relative to the variability) that drying equivalent to the WW2 drought or drier seems to be inevitable across Victoria. However, results here indicate that moving to a lower emission pathway would substantially reduce (by factor of 2) the magnitude of drying after 2060. Of course, internal variability over any particular decade may temporarily reinforce the externally forced drying to a level never seen before, as already seen in the case of recent 22-year drought. Under a low emission scenario (RCP2.6), the percentage decline is largest for MBVic towards the end of twenty-first century and is smallest for SWVic (Supplementary Table 2). In contrast, under high and medium emission scenarios (RCP8.5 and RCP4.5), the greatest drying occurs in SWVic, followed by MBVic and SEVic. This change in hierarchy may be related to the ozone hole repair, which is opposing GHGs related change and is most important in SWVic, so in RCP2.6 where GHG drying is weaker, the ozone-induced wetting is able to offset it more in SWVic. However, further research is needed to better understand the role of ozone on regional variation in rainfall across Victoria.

Summary
We have provided the best available estimates of the contribution of external forcing to the observed decline in rainfall in three Victorian sub-regions. However, confidence in our estimates is tempered by the fact that a majority of models tends to underestimate both the observed annual cycle during the cool season months, and the magnitude of the observed decline in rainfall. Furthermore, many models underestimate the observed inter-annual and inter-decadal variability during the cool season. As mentioned in Rauniyar and Power (2020), this may be related to the findings that CMIP5 models fails to reproduce some key features of atmospheric circulation and modes of climate variability and their relationship with rainfall across Victoria. This could be linked to the coarse resolution of the CMIP5 models (Chen and Dai 2019) or an underestimation of the Pacific Decadal variability (Power et al. 2017). Additional research is needed to clarify these issues using the next generation of climate models and downscaled simulations in future studies. Nevertheless, the methods described in this paper can be applied to any regions that are contiguous in their response to weather systems to examine the relative roles of external forcing and internal variability in the observed unusual change in rainfall by taking the advantage of GCM outputs directly.