Continental-scale bias-corrected climate and hydrological projections for Australia

The Australian Bureau of Meteorology has developed a national hydrological projections (NHP) service for Australia. The NHP aimed to provide nationally consistent hydrological projections across jurisdictional boundaries to support planning of water-dependent industries. NHP is complementary to those previously produced by federal and state governments, universities, and other organisations for limited geographical domains. The projections comprise an ensemble of application-ready bias-corrected climate model data, derived hydrological projections at daily temporal and 0.05° × 0.05° spatial resolution for the period 1960–2099, and two emission scenarios (Representative Concentration Pathway (RCP) 4.5 and RCP8.5). The spatial resolution of the projections matches that of gridded historical reference data used to perform the bias correction and the Bureau of Meteorology's operational gridded hydrological model. Three bias correction techniques were applied to four CMIP5 global climate models (GCMs), and one method was applied to a regional climate model (RCM) forced by the same four GCMs, resulting in a 16-member ensemble of bias-corrected GCM data for each emission scenario. The bias correction was applied to fields of precipitation, minimum and maximum temperature, downwelling shortwave radiation, and surface winds. These variables are required inputs to the Bureau of Meteorology's landscape water balance hydrological model (AWRA-L), which was forced using the bias-corrected GCM and RCM data to produce a 16-member ensemble of hydrological output. The hydrological output variables include root zone soil moisture (moisture in the top 1 m soil layer), potential evapotranspiration, and runoff. Here we present an overview of the production of the hydrological projections, including GCM selection, bias correction methods and their evaluation, technical aspects of their implementation, and examples of analysis performed to construct the NHP service. The data are publicly available on the National Computing Infrastructure (https://doi.org/10.25914/6130680dc5a51, Bureau of Meteorology, 2021), and a user interface is accessible at https://awo.bom.gov.au/products/projection/ (last access: 24 November 2023).

1 Introduction

Australia's climate has large natural variability encompassing many different climate zones and multiple climate drivers, leading to frequent floods and extended periods of droughts. In addition, Australia's climate is changing with increasing temperatures and shifts to rainfall patterns impacting all components of the hydrologic cycle. Seasonal precipitation has changed across the country, broadly increasing in the north and decreasing in the south (Dey et al., 2019), with similar patterns evident in streamflow (Zhang et al., 2016). The State of the Climate 2022 (CSIRO and Bureau of Meteorology, 2022) reports that cool-season (April–October) precipitation has declined by around 15 % since 1970 in Australia's southwest. Across the same region, May–July rainfall has seen the largest decrease, i.e. around 19 % since 1970. In the southeast of Australia, precipitation started to decline around 1990, and the average cool-season precipitation from 2000 to 2019 was 10 % less than the previous century. Along with this observed decline in precipitation, streamflow has declined substantially in both the southwest and southeast, with changes in streamflow typically disproportionally larger than changes in precipitation (Chiew, 2006; Wasko et al., 2021; Zhang et al., 2016).

To deal with changes in water availability, Australia's water policy and infrastructure investment decisions require high-resolution hydroclimate information representative of past, current, and future variability. Of prime interest are precipitation changes and how those changes will impact runoff, soil moisture, and evapotranspiration across Australia. At present, many water utilities, planners, and operators use in-house hydrological models to manage their resources at strategic time frames of 20–50 years (or more), whereby the simulations of streamflow feed into decision-making, planning, and implementation requiring significant investment. Traditionally, the modelling approach for these studies has been to use the historic climate record of a variable (e.g. rainfall, runoff) for a region or catchment of interest and extend it using stochastic methods to increase the sample of extreme events (Crosbie et al., 2009; CSIRO, 2009). However, this approach is insufficient to adequately represent the range of plausible climate futures. Several Australian states have commissioned and produced downscaled climate projections of rainfall (and in some cases runoff). However, a key issue is that the methodologies, ensemble sizes, emissions scenarios, and sets of host global climate models (GCMs) used vary considerably, meaning the data are too heterogenous for use in applications that intersect jurisdictional boundaries.

Detailed climate projections for Australia have been produced by the Bureau of Meteorology and the Commonwealth Scientific and Industrial Research Organisation (CSIRO), resulting in the Climate Change In Australia service (CSIRO and Bureau of Meteorology 2015; hereafter referred to as CCiA). The CCiA service provides application-ready data and information on climate projections for Australia using a subset of eight GCMs selected from all CMIP5 GCMs for their ability to reproduce climate features that are important for the Australian climate and suitable for use in impact models. The application-ready data were produced using a mean scaling for all variables other than precipitation for which quantile–quantile scaling was applied. However, the application-ready data are not continuous¹, and explicit hydrological output was not provided. Data currently available for runoff projections that are sufficient to cover multiple catchments from one data source have been developed by Chiew et al. (2018) (see also Zheng et al., 2019b). This is an ensemble of national projections of runoff using a delta change methodology that applies the mid-century change signal from 42 GCMs to historical time series, which are then run through GR4J rainfall–runoff models (Perrin et al., 2003) for each catchment. As with any scaled dataset, applications are limited to uses for which continuous series are not required and assume the variability and persistence in the observed period (for example, the daily distribution and sequencing of precipitation) is representative of a future climate undergoing warming. Ensembles of daily hydrological projections produced from bias-corrected climate model output are available globally, including for Australia, via ISIMIP (https://esg.pik-potsdam.de/search/isimip/, last access: 24 November 2023) at a spatial resolution of 0.5° (approx. 50 km). The application of this data is limited by the relatively coarse spatial resolution, and there is no information describing the suitability of the bias correction reference data for Australian conditions.

To address these deficits in hydrological projections, the Australian Bureau of Meteorology (the Bureau) has produced an ensemble of daily projections of hydrological variables at 0.05° × 0.05° (approximately 5 km) spatial resolution spanning the period 1960–2099. In addition, application-ready datasets of the variables required to produce the hydrological output at the same temporal and spatial resolution are also available. These were produced as a central component of the Bureau's National Hydrological Projections (NHP) project. The overarching goal of the NHP was to provide a nationally consistent set of projections, including publicly accessible datasets and communication tools that will be useful to water-sensitive sectors of the Australian community to prepare for and adapt to the many possible impacts of climate change, especially those sectors for which the hydrological cycle is of importance. By developing these datasets we aim to address these limitations and provide nationally consistent hydrological projections for the Australian community.

Development of the application-ready climate data and the derived hydrological projections required several considerations:

• GCM selection and the associated choice of emission scenarios,

• selection of a regional climate model (RCM) to dynamically downscale the GCM output;

• selection of techniques to correct the biases of the GCM and RCM output;

• creation of an observational dataset, which is required for both implementation and evaluation of bias correction;

• a choice of a hydrological model to run using the downscaled climate data;

• development of communication tools for users of the projections data.

Each of these decisions was guided by both pragmatic and scientific considerations. The various components of these decisions and the resultant modelling chain are shown in Fig. 1. The objective of this paper is to detail the decisions made and methodologies employed to develop a national continental-scale hydrological projections service.

Figure 1 Modelling flow diagram for the National Hydrological Projections. Individual components are detailed in the main text (Sects. 2, 3, 4 and 5).

The structure of the paper broadly follows points listed above. Section 2 details the GCM selection process, including how the selected GCMs compare to the larger collection of CMIP5 models, and the model used to dynamically downscale the GCMs. The reference observational data are described in Sect. 3. The bias correction methods implemented are reviewed in Sect. 4. The AWRA-L hydrological model is described in Sect. 5. A brief evaluation of the bias correction methods against the reference observational data is detailed in Sect. 6. The GCM and hydrological projections and some example analyses are detailed in Sect. 7, with a focus on extreme rainfall and runoff events in Sect. 7.4. We provide details of how the datasets were used to provide communication tools for the projections to end users in Sect. 8. Finally, conclusions, considerations around the limitations of the projections, and details of code availability are presented.

2 GCM selection

The GCM selection was guided by both pragmatic and scientific considerations. The AWRA-L hydrological model (see Sect. 5) requires daily values of precipitation (pr), maximum and minimum temperature (tasmax, tasmin), surface wind speed (sfcWind), and downwelling shortwave radiation (rsds) as inputs. Furthermore, the Australian research community uses CMIP5 GCM data stored on the National Computing Infrastructure (NCI) supercomputer, meaning that the ready availability of these variables (at daily temporal resolution) was of paramount practical importance. Ideally, an exhaustive set of GCM data, for instance, the complete set of 47 CMIP5 simulations assessed in CCiA, would have been utilised; however, that was not feasible, requiring us to choose a subset. The CCiA report guided this decision-making, wherein it recommended the use of eight GCMs to prepare application-ready data for Australia. They were chosen primarily based on their ability to represent the climate drivers of significant importance for Australia (El Niño–Southern Oscillation (ENSO), the Southern Annular Mode (SAM), and monsoonal rainfall patterns), as well as to provide a reasonable representation of the wet–dry and cool–warm variation between CMIP5 models (CSIRO and Bureau of Meteorology, 2015). Using this guidance and the practical considerations around data availability, the final suite of models chosen was ACCESS1-0, CNRM-CM5, GFDL-ESM2M, and MIROC5 (see Table 1).

Table 1 List of GCMs used in this study.

A choice of emission scenarios was also required, and simulations produced using the Representative Concentration Pathway (RCP) 4.5 and RCP8.5 emission pathways were chosen. These emissions pathways were chosen to provide a high (RCP8.5) and moderate (RCP4.5) set of temperature projections, noting that the set of modelled greenhouse gas emission pathways provided in CMIP5 have relatively minimal deviation before 2050 and that the observed climate change trends for CO₂ emissions and temperature in recent decades indicate that the high-emission pathway (RCP8.5) has been followed more closely than other emissions pathways (e.g. RCP2.6) (Schwalm et al., 2020; Stocker et al., 2013). However, we also note that RCP8.5 may be at the high end of emission scenarios for future risk assessment (Peters and Hausfather, 2020).

2.1 Regional climate model (CCAM)

To extend the GCM data, the Conformal Cubic Atmospheric Model (CCAM; McGregor, 2005; McGregor and Dix, 2008) was used to dynamically downscale the four GCMs. CCAM outputs contribute to the CORDEX (Coordinated Regional Downscaling Experiment)² international intercomparison of downscaling models. It is a global atmosphere-only climate model that has the ability for its grid to be stretched such that it has a finer resolution over the domain of interest (Australia) and coarser resolution in more distant regions (McGregor, 2005). CCAM can be forced via two methods, one using bias-corrected sea surface temperatures (SSTs) and another using spectral nudging to update the circulation to that of the host GCM. It was the former configuration used for NHP, whereby the mean and variance SSTs of the host GCM are bias corrected to provide boundary conditions for CCAM. This produced 50 km resolution projections of the atmospheric state over the Australian continent (Clarke et al., 2019; Hoffmann et al., 2016). As for the host GCMs, the simulations were run for the period 1960 to 2100 for both RCP4.5 and RCP8.5 scenarios. The resulting 50 km simulations represent a reconstruction of the atmosphere for which the GCM SST biases have been removed and the atmosphere is allowed to evolve independently of the host GCM.

The only information incorporated into CCAM from the host GCM is the bias-corrected SSTs, so the resultant atmospheric evolution may depart quite markedly from the host GCM. This has the effect that the range of the projections, and by implication the uncertainty, will differ from those that only include the host GCM. There are systematic differences between the GCM and CCAM projections, which we will explore in Sect. 7; however, the general effect of the inclusion of the CCAM simulations in the NHP ensemble is to increase the ensemble range. The approach to forcing the downscaling model is the same as that used for the Victorian Climate Projections – 2019 (Clarke et al., 2019) but differs from that used in some other projects, e.g. the Energy Sector Climate Initiative (ESCI), where a spectral nudging approach was used.

2.2 Comparison to the CMIP5 ensemble

It is informative to understand how our selection of GCMs fits within the larger ensemble of CMIP5 model simulations. Figures 2 and 3 show the ranking of 40 GCMs used in CCiA, for surface temperature (tas) and precipitation (pr). The ranking was evaluated by calculating the areal (in this case all of Australia) and temporally averaged difference in the projection period and the same quantity for the historical period; we will refer to this as a “change factor”. To examine the change factor, we define four future 30-year time slices: 2016–2045, 2036–2065, 2056–2085, and 2070–2099, which are centred on 2030, 2050, 2070, and 2085, respectively, while the historical period was defined as 1976–2005. Due to Australia's large size and geographical location, the climate of Australia varies markedly from the tropical north to the temperate south. As such, nationally averaged precipitation may not provide meaningful insight from a climatological perspective; nevertheless, it does impart interpretation of how our choice of GCMs occupies the phase space spanned by the CMIP5 models. For surface temperature, our selection falls approximately midway in the range of the CCiA ensemble; however, it still spans a relatively large range of change factors. The spread amongst the change factor values generally increases as a function of time and is also slightly larger for the RCP8.5 scenario. The ACCESS1-0 and GFDL-ESM2M models have an increased warming signal compared to MIROC5 and CNRM-CM5. The change factor for precipitation consists of both wetting and drying signals across each time slice and scenario, although generally two models (ACCESS1-0 and GFDL-ESM2M) show a drying signal and two show increased precipitation (CNRM-CM5 and MIROC5). GFDL-ESM2M has a consistent and large drying signal, while the precipitation change signal from the ACCESS1-0 model varies between time slices and scenarios; CNRM-CM5 and MIROC5 both indicate enhanced precipitation.

Figure 2 The ranking of surface temperature (tas) for the GCMs used in this study (darker shading) with regard to the CCiA ensemble (lighter shading). The horizontal bars indicate the change signal (the difference in the nationally averaged quantity from the climatology for the period 1976–2005). Four 30-year periods are shown centred on 2030, 2050, 2070, and 2085 for RCP8.5.

Figure 3 The same as Fig. 2 but for precipitation. Note that percentage change in precipitation is shown.

Figure 4 Change plot of surface temperature (a, c) and precipitation (b, d) for all of Australia. The NHP models are indicated by the larger blue symbols, while the CCiA models are shown as smaller grey-filled circles. The thick vertical lines represent the ensemble median, and the left and right extent of the rectangles show the 10th and 90th percentile of the model ensemble, respectively. Time is displayed on the ordinate, and the magnitude of change (from the historical period) is shown on the abscissa; RCP4.5 is shown in blue, and RCP8.5 is shown in red. The larger rectangles and vertical grey lines, and dots show the spread of the CCiA ensemble, while the smaller hatched rectangles, vertical blue lines, and shapes display the NHP ensemble.

Although focussing on a particular model is useful, when using the output from several models the range of uncertainty (provided by the differing simulations) can be examined by constructing the model ensemble and calculating ensemble statistics. As shown in Fig. 2, the change factor is calculated for each model separately, which is then aggregated to calculate ensemble statistics, including the median and the 10th and 90th percentiles. Here, we introduce “change factor plots”, an example of which is shown in Fig. 4. It displays the model ensemble change factor spread for surface temperature (left) and precipitation (right) for the NHP ensemble and the CCiA ensemble, averaged over Australia, for both RCPs and for the four future time slices. The NHP ensemble median is cooler than that of the CCiA ensemble and has a decreased ensemble range. The precipitation signals between the NHP ensemble and CCiA ensemble are often of opposite sign, especially for RCP4.5. Although the ensemble range is generally smaller for the NHP models, it is comparable to that of CCiA, and in the case of the 2030 and 2050 time slices for RCP4.5 it is larger.

Figure 4 indicates that the NHP ensemble provides a subsample of the CCiA ensemble. To focus on how the climate in a particular region is projected to change in the future it is useful to divide a large geographical region, such as Australia, into sub-regions. In keeping with the recommendations of CCiA, the Natural Resource Management (NRM) areas are used in this report, and we will refer to changes across NRM clusters (see Fig. 5). These regions were chosen to reflect locally similar climatic conditions, biophysical factors, and broad patterns of climate change (CSIRO and Bureau of Meteorology, 2015). CCiA also included further sub-divisions that may be better suited for impact and assessment studies. However, our intention here is not to provide such assessments but rather a broad overview of the data available that can be used for such studies in the future.

Figure 5 The eight NRM regions used in the NHP to present results relevant to the differing climate regimes across Australia. Image sourced from CSIRO and Bureau of Meteorology (2015).

Change plots were also constructed for the NRM clusters (not shown) and the discrepancy between the NHP and CCiA ensembles changes depending on the geographical region. For instance, the NHP warming signal in the Southern and Southwestern Flatlands (SSWF) NRM is not as large as the CCiA ensemble; however, the precipitation change signal for both the CCiA and NHP ensembles indicates a consistent drying trend despite the NHP ensemble drying signal being lower. As for all of Australia, the range of temperature and precipitation projections of the NHP ensemble was evaluated against the 40-member CCiA ensemble, and it was found that the NHP ensemble spanned the complete ensemble range remarkably well. This means that the four GCMs used in NHP include both cooler and warmer and dryer and wetter projections, giving us confidence that the mean changes in the NHP ensemble reflect that of the complete ensemble. Nevertheless, the spread of the NHP ensemble is less than that of CCiA, suggesting that the NHP GCM ensemble spread may be less than the CCiA GCM recommendations. However, we note that this analysis only considers the four host GCMs and does not include CCAM simulations. Since CCAM develops a climate independent of the host GCM, it has the effect of increasing the ensemble spread compared to that of the host GCMs. While the inclusion of the CCAM ensemble increases the uncertainty in the projections, rather than considering this a limitation, it is instead advantageous as it enables sampling more of the phase space represented by the full suite of CMIP5 models (see Figs. 2, 3, and 4), thereby decreasing the reliance on the interpretation of the results for GCM selection.

3 Reference data

The outputs from GCMs and RCMs have systematic biases that need to be corrected before using their data for climate change impact and risk assessment studies. This is achieved using bias correction methods that are designed to preserve certain moments between the observations and the simulations. Implementation of a bias correction method requires three main components: (1) GCM and/or RCM simulations that are to be corrected for bias; (2) reference data, which can be direct observations or a derived product (e.g. gridded observations or reanalyses); and (3) a bias correction algorithm. The reference datasets required are those necessary to run the AWAR-L hydrological model (see Sect. 5): precipitation (pr), maximum and minimum temperature (tasmax$/$tasmin), downwelling shortwave radiation (rsds), and surface wind speed (sfcWind). Here we describe the reference datasets used to bias correct the GCM output over the historical period (defined to be 1976–2005).

3.1 Precipitation and temperature

Daily gridded precipitation and temperature data were obtained from the Australian Water Availability Project (AWAP)³ climate dataset, which consists of surface air temperature (daily minimum and maximum) and daily precipitation from 1 January 1911 to the present (Jones et al., 2009). The rainfall and temperature data are interpolated from daily station records using a Barnes successive-correction scheme that applies weighted averaging to the station data and then provided on a 0.05° (approximately 5 km) spatial grid across Australia. Topographical information is included by using anomalies from long-term (monthly) averages in the analysis process. For temperature, the daily anomaly is represented by the difference from the monthly mean, while for rainfall the ratio of the daily observed precipitation to the monthly average precipitation is used.

Precipitation data represent the total recorded over the 24 h beginning (and ending) at 09:00 LT (local time) and is assigned to the day from which the recording began. The highest temperature over the 24 h before the observation at 09:00 LT is recorded as the maximum temperature (T_max) for the previous day. The lowest temperature (T_min) for the 24 h before 09:00 LT is recorded for the day on which the observation was made. This convention usually guarantees that T_min and T_max will have occurred on the same calendar day; however, in some instances it is possible that $T_{min}>T_{max}$. For the ISIMIP2b method (see Sect. 4.1), the surface temperature (tas) is corrected (not T_min and T_max) and the minimum and maximum are derived (see Eq. 26 in Hempel et al., 2013). As such, for the ISIMIP2b method, tas was derived from the AWAP data using their average (i.e. $\mathrm{tas}=(T_{min}+T_{max})/\mathrm{2}$). In some circumstances tasmin can exceed tasmax ($T_{min}>T_{max}$), for which we set $T_{min}=T_{max}$. The other bias correction methods (multivariate recursive nested bias correction and quantile matching for extremes) correct T_min and T_max directly; therefore, no calculation of tas was required for these methods. Note that while observations are made from 09:00 to 09:00 LT the next day, the bias corrections are calculated by calendar day.

Since AWAP is a gridded dataset it has errors compared to measurements from the original station data. This is particularly true in data-sparse areas for rainfall since it is highly variable both spatially and temporally. For this reason, some regions, particularly in central Western Australia, have a mask applied and are excluded from the analysis. Furthermore, in areas of steep topography and a sparse gauge network, particularly along the Great Dividing Range, AWAP has been shown to underestimate rainfall (Chubb et al., 2016). It has also been shown to underestimate extreme rainfall percentiles when compared to station data (King et al., 2013).

3.2 Solar radiation and wind speed

Daily surface shortwave downwelling solar radiation flux (rsds) is obtained from geostationary satellites (Grant et al., 2008) and aggregated to the same 0.05° AWAP grid. The solar radiation record is from 1990 to the end of 2005 (the calibration period). From 1 January 1990, the GMS-4, GMS-5, GOES-9, and MTSAT-1R satellites were used, although these were replaced by the Himawari-8 satellite from 23 March 2016 onward. Daily climatological averages (for each day of the year) are used for solar radiation before 1990.

Site-based daily near-surface wind speed observations (sfcWind) collated by the Bureau of Meteorology and interpolated nationally (McVicar et al., 2008) were used for daily average wind speed from 1975 onwards (when sufficient site observations were being collated by the Bureau) for the GCM bias correction. The near-surface wind dataset is for the wind speed recorded at a height of 2 m; however, since the GCM data are reported at a height of 10 m, the observed 2 m winds were transformed to 10 m using the following logarithmic transform (Garratt, 1992),

$$\begin{array}{}\text{(1)}& u\left(z_{\mathrm{2}}\right)=u\left(z_{\mathrm{1}}\right){\displaystyle \frac{\ln \left[\right(z_{\mathrm{2}}-d)/z_{\mathrm{0}}]}{\ln \left[\right(z_{\mathrm{1}}-d)/z_{\mathrm{0}}]}},\end{array}$$

where u(z₂) is the observed wind speed at height z₂, u(z₁) is the wind speed at height z₁, z₀ is the roughness length to account for the effect of the roughness of a surface on wind flow, and d is the zero-plane displacement to account for flow around obstacles. The bias correction was applied to the derived 10 m wind speeds; however, AWRA-L (Sect. 5) requires wind speed at 2 m height, so Eq. (1) was again used to transform the near-surface wind speed back to 2 m height as input for the AWRA-L simulations.

4 The bias correction methods

Many bias correction (BC) techniques have been developed, and several authors have provided detailed summaries and critiques of them (e.g. Fowler et al., 2007; Hewitson et al., 2014; Maraun, 2013, 2016; Maraun and Widmann, 2018; Teutschbein and Seibert, 2012); however, they fall into three general categories: (1) empirical scaling techniques, (2) non-parametric quantile matching (QM), and (3) parametric QM (Potter et al., 2020). Empirical scaling techniques include linear methods that preserve only the mean and apply a constant correction factor (Lenderink et al., 2007), while non-linear scaling, such as power transformations, preserves the mean and variance (Leander and Buishand, 2007). Non-parametric QM techniques map the simulated quantiles of the cumulative distribution function (CDF) to the observed CDF quantiles without any underlying assumptions that the variable can be modelled by a mathematical distribution, whereas for parametric QM distributions are fitted to the variables before the application of QM (Heo et al., 2019). One of the advantages of using a parametric method is that it can provide a useful estimate based on the distribution when matching simulated fields that fall outside the calibration range of the historical observations (while noting that the extrapolation of a parametric distribution to values outside the range used to produce that distribution involves some uncertainties). For instance, the simulated field of projected temperature may not have any historical precedent, in which case assumptions have to be made as to how it will be associated with a quantile based on observations (Piani et al., 2010a, b).

It is important to note that bias correction cannot add climate change information to that of the GCM to which it is applied. It instead modifies the coarse-scale GCM projections at a finer resolution to produce localised projections (Ekström et al., 2015). Furthermore, the application of bias correction may correct the bias at a local scale while ignoring the interdependency of variables, which has the unintended consequence of interfering with the underlying physics of the GCM (e.g. conservation of energy and mass) such that the resulting temporal and spatial representation of fields is unrealistic (Maraun et al., 2017).

Three statistical bias correction methods were applied to the GCM output, and one (ISIMIP2b)⁴ was applied to the CCAM output (see Fig. 1). In addition, the original GCM data spatially interpolated to the resolution of the AWAP were archived to provide data before the application of bias correction; we refer to these data as “NOBC” and “NOBC-CCAM” for the GCM and CCAM data, respectively. At a fundamental level, they are all variants of quantile matching (e.g. Maraun and Widmann, 2018); however, they have important distinctions concerning their implementation. Despite the various pros and cons discussed in the literature of various BC methods, the choice was also (as for GCM selection) guided by pragmatic considerations, the foremost being ready availability to source code, thorough documentation in the scientific literature, and access to technical support. Here we present a mostly qualitative overview of the individual methods; a complete mathematical formulation of the methods can be found in the respective citations.

4.1 ISIMIP2b

The Inter-sectoral Impact Model Intercomparison (ISIMIP) Project (https://www.isimip.org/, last access: 24 November 2023) was designed to offer a consistent framework for projecting climate change impacts at different global warming levels across different sectors and spatial scales. Some of the sectors catered for in ISIMIP include regional water, fisheries, agriculture, biomes, and terrestrial biodiversity. To maintain consistency across so many diverse impact assessment models a common framework for bias correction and downscaling of GCM data was implemented. There have been several iterations of the bias correction method adopted for ISIMIP. At the commencement of the NHP, the most recent version of the bias correction available for download was the one developed for the ISIMIP2b protocol. The complete description of the bias correction used in ISIMIP2b is a modification of the method described by Hempel et al. (2013), which extends the method described in Piani et al. (2010a, b). Although ISIMIP is a set of protocols and output data across different impact models, we will refer to the bias correction as ISIMIP2b to be explicit. The ISIMIP2b method is classified as a parametric quantile matching (QM) method as it assumes variables of interest can be described by a particular parametric distribution. It uses Gaussian distributions to model daily temperature and gamma distributions to model daily precipitation probability distributions. For temperature, which can take both negative and positive values, the correction applied is additive, while in the case of fields such as precipitation, solar radiation, and surface wind speed, which have positivity constraints, the correction applied is multiplicative.

The ISIMIP2b method is comprised of the following steps: (1) correction of the monthly mean, (2) correction of the daily variability, and for precipitation (3) correction of the frequency of dry days and (4) correction for the intensity of wet days. Steps (1) and (2) have different implementations depending on whether the correction to be applied is additive or multiplicative. Furthermore, steps (1) and (2) are designed to preserve the absolute trend of simulated temperature and the relative trend of the other variables. The preservation of the absolute long-term trend in the simulated temperature field preserves the climate sensitivity of the GCM and further ensures that there is consistency between the projected temperature and the bias-corrected climate change signal. For variables that have positivity constraints, the ISIMIP2b method conserves the relative trend. The use of an additive correction for temperature and a multiplicative correction for precipitation ensures that the ratio of the relative change in precipitation to an absolute change in temperature is conserved and thus does not alter the hydrological sensitivity of the GCM (Hempel et al., 2013).

4.2 MRNBC

The multivariate recursive nested bias correction method (MRNBC) is an extension of QM to include inter-variable dependencies and in addition corrects across multiple timescales. It was progressively developed from the nested bias correction (NBC; Johnson and Sharma, 2012) and recursive nested bias correction techniques (RNBC; Mehrotra and Sharma, 2012). The NBC corrects the distribution (mean and standard deviation) and persistence (lag 1 autocorrelation coefficient) at monthly, seasonal, and annual timescales using a standard autoregressive lag 1 model (Srikanthan and Pegram, 2009). The RNBC method is an extension of the NBC method, whereby the method is applied three to five times repeatedly so that it significantly reduces the biases in mean, variability, and persistence-related attributes in GCM and RCM simulations. The MRNBC method is a multivariate version of the above RNBC method. It simultaneously corrects many GCM and RCM variables, using a multivariate first-order autoregressive model at daily, monthly, quarterly, and annual timescales to impart observed distributional and persistence properties of the input fields (Mehrotra and Sharma, 2015).

The implementation of a multivariate autoregressive model means that the MRNBC may better capture the joint dependence among input variables and hence be more effective in capturing the behaviour of natural processes that contribute to the variability in a field, particularly precipitation. Furthermore, as the MRNBC also corrects for the biases at seasonal and annual timescales, the bias-corrected atmospheric fields may provide a better representation of long-term variability in the downscaled simulations. This is of particular importance for hydrological applications where the representation of variability and persistence of precipitation is important for the simulation of the occurrence and intensity of extreme events such as floods and droughts.

4.3 QME

The quantile matching for extremes method is described in Dowdy (2020, 2023). As its name implies, it has a focus on matching the extremes of the CDF of a variable. It was originally developed to examine future projections of fire weather conditions for Australia based on the Forest Fire Danger Index (FFDI) (McArthur, 1967). In general, QM methods divide the quantiles of the observed and simulated fields into discrete bins (for example 100 percentile bins) so that events in the top percentile of the simulated field are associated with those in the same quantile bin from the observational period. However, having uniform resolution quantiles does not capture extreme events well, especially those that fall outside the historical record. Rather than using uniform percentiles, the QME method applies a logarithmic transformation so that there is higher resolution in the percentile bins at the upper and lower ranges. This has the effect of being able to assign percentiles in the upper (and lower) tails (e.g. the 99, 99.5, and 99.9 percentiles) without the use of interpolation.

The QME method allows for simulated values that fall outside the historical reference values to be bias corrected while still avoiding the influence of extreme outliers that may have occurred in the historical period. As an example, the temperature recorded near Melbourne on Black Saturday (7 February 2009) had no historical precedent in the prior 109 years of records (Trewin et al., 2020). An exceptionally rare event that happens to occur in the training period may not be ideal to use in isolation for calibrating the model data. In the QME method, an additional step is included to help represent extremes, where the top and bottom five most extreme values are used to calculate the mean bias between the model and the reference data. By using the mean of the five extreme values (rather than fewer than that) the QME method is intended to provide a reasonable estimate of extremes up to about the 10–20-year average recurrence interval (ARI), while not being overly influenced by individual exceptional events that may occur. This is done individually for the top five values, as well as for the bottom five values, where the mean biases are subtracted from the extreme values for the model data. The QME method is trained using the 43 years from 1975 to 2017, such that estimates of extremes should be reasonably well calibrated up to values of around the 10-year return period (i.e. ARI near 10 years) based on a sample size of about four to five events on average for estimating those extreme values. Additionally, the mean bias for the five most extreme values is also subtracted from simulated values that lie outside the historical range of occurrence, noting the relevance of this in our changing climate.

4.4 Effect of bias correction on GCM output

Insight as to how the bias correction and the dynamical downscaling have modified the original GCM biases are shown in Figs. 6 and 7. Areal averages of the input and output data were produced for each ensemble member for all of Australia and each NRM region. Time series of maximum temperature for each GCM, averaged over all of Australia and smoothed with a 10-year running mean, are shown in Fig. 6. The NOBC data (shown in magenta) are significantly cooler than the AWAP reference data (shown in black), and the bias correction has brought the GCM data in closer agreement with the AWAP data. Consequently, the bias-corrected projections are warmer than the NOBC data for all GCMs. In contrast, the NOBC-CCAM data are warmer than the AWAP reference data, and the application of ISIMIP2b has generally decreased the warming bias; however, CCAM-ISIMIP2b is still significantly warmer than both the original GCM (NOBC) and bias-corrected data. The three bias correction methods follow very similar trajectories, and only slight differences are discernible for the GFDL-ESM2M model. Since temperature is a smooth field (both spatially and temporally) and all the bias correction methods are variants of quantile matching, this indicates that (at least for the mean and large geographical aggregation) the bias correction imparts less uncertainty than GCM selection in the projections.

Figure 6 Time series of maximum temperature averaged over Australia for each GCM in the ensemble. Different colours correspond to the bias correction methods. The labels NOBC and NOBC-CCAM indicate the GCM and CCAM data before the application of bias correction. Time series have been smoothed with a 10-year running mean. The AWAP reference data are shown in black.

Figure 7 The same as Fig. 6 but for precipitation.

A contrasting bias, shown in Fig. 7, is present for precipitation. All the original GCMs except ACCESS1-0 are significantly wetter than the AWAP reference data and that bias continues in the projections. The bias correction methods have lowered the average daily Australia-average precipitation. The effect of the bias correction methods on precipitation shows more contrast among the methods than for maximum temperature. The NOBC-CCAM simulation is wetter than the reference data; however, the application of ISIMIP2b has produced a generally drier bias than when applied to the GCM data. The bias-corrected data correlate with both the NOBC and NOBC-CCAM data; however, the data have translated the time series to lower values. Finally, the application of a 10-year smoothing window demonstrates how the CCAM model produces quite distinct projections from its host GCM.

5 The AWRA-L hydrological model

The Bureau's operational Australian Water Resources Assessment Landscape (AWRA-L) model was used to produce the National Hydrological Projections outputs. AWRA-L is a daily semi-distributed water balance model (Frost et al., 2018) that is run daily using climate input data (pr, tasmax$/$tasmin, rsds, sfcWind) to produce hydrological outputs (runoff, evapotranspiration, soil moisture) at the 5 km × 5 km grid resolution. Land use and vegetation are represented within each grid cell by the fractional coverage of two hydrologic response units (HRUs) representing shallow-rooted vegetation (grass) and deep-rooted vegetation. The terrestrial water storage component is represented by three soil layers, namely, top (0–10 cm), shallow (10–100 cm), and deep (100–600 cm) layers. Shallow-rooted vegetation is assumed to access soil moisture up to the shallow layer (up to 100 cm), while deep-rooted vegetation is assumed to access soil moisture up to and including the deep layer (up to 600 cm). The soil moisture content modelled AWRA-L in the three soil layers is denoted by s0, ss, and sd, respectively. For many hydrological applications, the root zone soil moisture (sm), which is the sum of the soil moisture in the top two soil layers (0–100 cm), is of most interest (i.e. sm = s0 + ss). Many fluxes and stores are output by the AWRA-L model; however, for the NHP, we focus on runoff (qtot), potential evapotranspiration (e0), and root zone soil moisture (sm).

The AWRA-L model was calibrated by using 295 unimpaired catchments and validated against 291 evaluation catchments (Frost and Wright, 2018). The median values of the relative bias and daily and monthly Nash–Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970) were −0.01 and 0.49 and 0.67, respectively, for runoff. Wasko et al. (2020) evaluated the performance of AWRA-L to simulate trends in catchment-aggregated surface runoff and found good correspondence with observed trends in the southeast of Australia, particularly in the Murray Basin and Southern Slopes NRM regions. They also found reasonable correspondence between modelled and observed trends in streamflow in the northern parts of Australia for the warmer months.

For the NHP, the AWRA-L model is forced using the historical daily gridded reference dataset to output a gridded observation-based dataset of the three hydrological output variables (qtot, e0, and sm); this dataset constitutes the reference hydrological data and, as for the GCM data, the historical period is defined as 1976–2005. To produce the historical AWRA-L hydrological data for the reference period (1976–2005), it was necessary to use historical GCM simulations (bias corrected to the reference period) beginning in 1960. This is because AWRA-L requires a “spin-up” time to equilibrate its water cycle since soil moisture is initialised at zero. The modelling of the deep-layer soil moisture requires several years of simulation for surface water to drain into the deeper levels. Starting the AWRA-L simulations in 1960 allows 16 years for the deep-layer soil content to reach a pseudo-equilibrium before the beginning of the historical period (1976–2005).

6 Evaluation of the bias correction methods

Each of the bias correction methods has been published previously and comprehensively evaluated (see references in Sect. 4). A further detailed evaluation of both the GCM and AWRA-L output was performed for the historical simulations (1976–2005) against the gridded reference datasets (Vogel et al., 2023). They examined the spatial and temporally aggregated (annual, seasonal, and monthly) bias for each NRM region and multiple individual grid locations around Australia. A complete evaluation is outside the scope of this paper; however, we provide an overview of the spatial biases for the hydrologically relevant variables (pr, qtot, sm and e0) and annual and seasonal biases aggregated over Australia for all variables. The bias is calculated as the difference between the bias-corrected historical simulations and the reference gridded data. Figure 8 shows the area-averaged annual biases (seasonal biases are shown for water-related variables for ACCESS1-0 bias corrected with the MRNBC method in Figs. S3–S6 in the Supplement). Absolute and relative biases are shown to aid the interpretation of the maps. For instance, an absolute bias of about 20 mm in the north of Australia during austral summer (DJF), which is also during the northern Australian “wet season”, is on the order of a −2.5 % relative difference. Conversely, in southern Australia, absolute deficits of about 5–10 mm translate to relative deficits of 5 %–25 %. The mostly dry relative biases in precipitation manifest as corresponding dry biases in the runoff and root zone soil moisture. Potential evapotranspiration most likely increases due to the increases in minimum and maximum temperature. Similar maps (not shown) were produced for all 16 ensemble members.

Figure 8 Bias maps of annual precipitation (a), runoff (b), root zone soil moisture (c), and potential evapotranspiration (d). The maps are for the ACCESS1-0 GCM with the MRNBC bias correction method applied. The bias is the difference between the historical (1976–2005) GCM simulations from the historical (1976–2005) gridded reference data. The mean bias is shown at the top of each map in the panel. Note that the scale for the absolute biases is dynamic, while for the relative bias it is fixed in the range −50 %–50 %. Data-sparse regions have been masked (see Sect. 3.1).

Figure 9 Bias histograms for all of Australia for precipitation (a), runoff (b), root zone soil moisture (c), and potential evapotranspiration (d). Annual and seasonal (four seasons) biases are shown. The legend indicates the various bias correction and GCM combinations.

Figure 10 The same as for Fig. 9 but displaying the biases for maximum temperature (a), minimum temperature (b), downwelling solar radiation (c), and near-surface wind speed (d).

The relative biases averaged over Australia are shown in Fig. 9 (for water-related variables) and Fig. 10 (for the remaining GCM output variables). Similar plots are shown for the SSWF NRM in Figs. S7 and S8 in the Supplement. Each subpanel corresponds to one of the bias correction methods implemented on each GCM. These summarise the information that was shown in the spatial bias maps (Fig. 8). The bar plots presented here are for Australia; however, similar analyses were conducted for each of the NRM regions to quantify their spatially aggregated bias and assess the contribution from each GCM and bias correction method. For instance, the application of MRNBC results in negative biases in rainfall, runoff, root zone soil moisture, and downwelling solar radiation for most GCMs on annual and seasonal timescales. The various model and bias correction combinations have differing bias magnitudes and signs between NRM regions. The reasons for this are due to the various GCMs' ability to simulate major climate drivers and also due to the AWRA-L model's ability to simulate the diverse hydroclimate of Australia. For example, each NRM will be dominated by differing land use, agricultural, and natural reserve regimes, which in turn impact the hydrological response in an NRM (Frost and Wright, 2018). However, it was found that the application of the bias correction methods resulted in satisfactory confidence in the ability to simulate the hydroclimate across each NRM on annual, seasonal, and monthly timescales. Note that the numbers displayed in Fig. 8 (and Figs. S3–S6 in the Supplement) differ slightly from those in the bar plots displayed in Figs. 9 and 10 due to slightly different methods in their calculation. For the former (bar plots), the Australia-wide mean was calculated first and then the relative bias (as a percentage) was calculated from the two regionally averaged values, while for the spatial maps, area-averaged biases were calculated by first calculating the relative bias (as a percentage) for each grid cell and then averaging it over Australia.

7 National Climate and Hydrological Projections

The core datasets produced comprise gridded daily time series at 5 km × 5 km spatial resolution spanning the period 1960–2099. They are composed of bias-corrected climate datasets of precipitation (pr), maximum and minimum temperature (tasmax and tasmin), surface winds (sfcWind), and surface downwelling solar radiation (rsds), which were used to produce National Hydrological Projections of runoff (qtot), actual and potential evapotranspiration (etot and e0), and three layers of soil moisture (s0, ss and sd). For communication purposes, only projections of potential evapotranspiration are displayed, and soil moisture is referenced by root zone soil moisture (sm), which is the linear addition of the top two layers (sm = s0 + ss) and represents the soil moisture in the top 1 m soil layer. For each variable, 16 ensemble members exist for each RCP, corresponding to the four GCMs with the three bias correction methods applied (12 ensemble members) and the CCAM RCM with ISIMIP2b applied (4 ensemble members). We will refer to the bias-corrected GCM data as input data and the AWRA-L simulations as the output data. The bias-corrected climate input data complement the existing climate projections for Australia and can be used as input to other impact and assessment models.

To guide the interpretation of nationally consistent hydrological projections for Australia (and the NRM regions), several lines of evidence were explored: area-averaged time series, spatial change maps, and area-averaged change plots for each RCP scenario and 30-year time slice. Here, we show examples of each of these, choosing the SSWF as an exemplar as its westernmost portion (southwestern Western Australia) has a well-documented historical rainfall decrease signal (CSIRO and Bureau of Meteorology, 2015). We also note that the rainfall declines in southwestern Western Australia have a robust attribution to an anthropogenic influence (Timbal et al., 2006), and the scientific consensus is for those trends to continue (Dey et al., 2019). However, our intention is not to provide a complete hydrological overview of a particular region as these are provided elsewhere (see Sect. 8) but rather to demonstrate the products and lines of enquiry that were developed.

7.1 Area-averaged time series

To examine the projections of the hydroclimate of a particular region, consider Fig. 11, which shows projections for precipitation, runoff, soil moisture, and potential evapotranspiration for the SSWF NRM region. Total area-averaged precipitation is expected to decline with concomitant declines in runoff and soil moisture, while potential evapotranspiration projections show increases. The ensemble statistics have been obtained by concatenating the historical and projected data and applying a 30-year running mean. This results in 15 years of data being excluded from the beginning and end of the time series, such that they cover the period 1975–2085. This emphasises the long-term decadal-scale signal in the projections. Interannual variability, shown as yearly averaged data using the ACCESS1-0-ISIMIP2b ensemble member, is shown to demonstrate the large departures from the long-term signal for all variables.

Figure 11 Projected ensemble change in precipitation, runoff, root zone soil moisture, and potential evapotranspiration for RCP4.5 (blue) and RCP8.5 (red) scenarios in the SSWF NRM region. The shaded area represents the 10th to 90th percentile range for all ensemble members. The time series for ACCESS1-0 (rcp8.5) with the ISIMIP2b method applied is included after 2020 to show the variability projected for an individual ensemble member. The AWAP reference gridded observations are shown prior to 2020. The ensemble has been smoothed with a 30-year running mean.

7.2 Spatio-temporal average change signal

Further information can be obtained by examining the change signal of future periods from the historical period (cf. Fig. 4). Spatial and temporal averages of the change signal were calculated for Australia and individual NRM cluster regions for each of the input and output variables. An example for the SSWF NRM at annual timescale is shown in Fig. 12. A general decrease is projected for all variables except potential evapotranspiration. The following general observations can be made: (1) there is a large spread among the ensemble members; however, the spread among the individual BC methods is less than that of the individual GCM to which they were applied; (2) the change signal may be of opposite sign between the various GCM/bias correction combinations; and (3) the median projected change is not necessarily monotonic, especially for the RCP4.5 scenario. We note that the median projected change signal for temperature (not shown) exhibits monotonic increases for all NRMs and both RCPs. In addition, we note that the (ensemble median) change signal for RCP8.5 for all variables is monotonic; however, that is not the case for all NRMs and is likely affected by the strong change signal for precipitation in southwestern WA (which is a subset of the SSWF NRM).

Figure 12 Relative change in annual (a) precipitation, (b) runoff, (c) root zone soil moisture, and (d) potential evapotranspiration for the SSWF NRM. Changes are indicated for each GCM and BC method combination for four 30-year time periods centred on 2030, 2050, 2070, and 2085. The change is the difference in the future from the historical period (1976–2005) expressed as a percentage. The left and right extents of the shaded boxes represent the 10th and 90th percentile of the GCM and BC ensemble, while the vertical solid line shows the median (50th percentile). RCP4.5 is shown in blue, and RCP8.5 is shown in red. The ensemble spread for each RCP has been offset from the centre of the 30-year time slices for legibility. Note that time is descending on the ordinate and that the change (in percent) is shown on the abscissa.

The first observation implies that, in general, the choice of GCM contributes more to the ensemble spread than the application of the BC techniques and that the inclusion of the CCAM model has added further diversity to the ensemble. The differing sign of individual members is a consequence of the model selection as indicated in Fig. 3, whereby some GCMs differ in their wetting or drying signal. The non-monotonicity of the change signal (more apparent for other NRMs) could be due to differing multi-year timescale processes (e.g. ENSO) represented by each GCM, especially in those areas where these processes dominate rainfall variability. Combined, these uncertainties make communication of the projections for end users a challenging problem; however, we address some of these in Sect. 8.

The diversity of climate in Australia makes the use of annually averaged changes unsuitable for most hydrological applications of the projections. As such, change plots were also constructed considering both four seasons (DJF, MAM, JJA, SON) and two seasons (a wet season of November–April and a dry season of May–October) to gain insight into the projected changes in seasonal hydrological projections. An example is given for SSWF in Fig. 13 for four seasons, which is appropriate for this region of Australia. It can be seen that the major signal of the decrease in precipitation, runoff, and soil moisture is due to their decreases in winter (JJA) and spring (SON), while the annual increase observed in potential evapotranspiration has no clear contribution from any particular season. The wetter signal is generally due to the MIROC5 and CNRM-CM5 models; however, the signal of CCAM forced with these models tends toward a drying of these signals. Note that in Fig. 13 the annually averaged absolute (rather than relative) values are shown, indicating that small absolute changes may express as large percentage changes, particularly in low-rainfall regimes. For example, the annual change in runoff ranges from an approximately 20 % to 60 % decrease; however, most of this is expressed as winter and spring runoff decreases of approximately 2 to 8 mm yr⁻¹.

Figure 13 The same as Fig. 12 but displaying seasonally averaged absolute change values for four seasons (DJF, MAM, JJA, SON) for each variable: precipitation (a), runoff (b), root zone soil moisture (c), and potential evapotranspiration (d).

7.3 Time-averaged spatial plots

The final line of evidence used was that of time-averaged spatial plots, where the temporal average of the daily data for the historical and future time slices was calculated and the change (future minus historical) plotted. An example of the winter (JJA) ensemble median of variables is shown in Fig. 14. They present the same temporal information displayed in Fig. 13; however, they include the spatial variability in the change signal. For instance, the decreases in precipitation, runoff, and soil moisture for SSWF, displayed in the winter (JJA) sub-panels of Fig. 13, are evident along the southern and western regions of Australia. The spatial variability indicates that the change signal of decreases in precipitation, runoff and soil moisture is more prominent in the Western Australian portion of SSWF than in the South Australian sector. The projections for increased potential evapotranspiration apply to all time slices and RCPs. Similar plots were produced for all time slices and RCPs for annual and seasonal (both two and four) temporal aggregations for both relative and absolute changes.

Figure 14 Ensemble median relative change in winter (JJA) precipitation (a), runoff (b), root zone soil moisture (c), and potential evapotranspiration (d). The change signal is only shown for the 2030 and 2070 time slices for RCP8.5.

7.4 Extreme wet events

The previous analyses have focussed on changes in the mean, whether it be the temporal and spatial mean (e.g. Figs. 12 and 13) or a temporal mean (Fig. 14). Mean changes are useful quantities; however, for many impact studies it is the quantification of events with a low probability and high impact (or extreme events) that are most crucial. As the climate warms, heavy-rainfall events are expected to increase due to the increased moisture content of the atmosphere (Sherwood et al., 2010). Model projections have indicated that the frequency of extreme rainfall events in Australia is also expected to increase (Alexander and Arblaster, 2009; Perkins et al., 2014; Watterson et al., 2017). Extreme rainfall may be expected to manifest itself as an increased probability of flooding events; however, the risk of flooding is also determined by antecedent soil moisture conditions. Wasko and Nathan (2019) suggested that in Australia (and other parts of the world) that following heavy rainfall events soil moisture deficits are first augmented such that despite increases in heavy rainfall events there has been a decrease in flood magnitudes. However, it was also noted that soil moisture plays a decreasing role in flood severity as catchment size decreases (Wasko and Sharma, 2017). Here we examine some future flood scenarios based on changes in low-probability, heavy-rainfall and runoff events, which we denote as extreme events. Characterising changes in flood frequency and intensity at a large spatial and temporal scale is challenging; flood risk is often dependent on local topography, sub-daily rainfall intensity, and antecedent conditions. A set of threshold-based indicators using precipitation and runoff have been calculated here to capture changes in flood risk on a broad scale. We have analysed these extreme rainfall and runoff events using the generalised extreme value (GEV) distribution (Kharin et al., 2007; Perkins et al., 2014), where the changes are determined using the projected annual maximum daily rainfall and runoff and the estimate of the 20-year return period of the annual maximum.

Figure 15 shows the results of the application of the GEV analysis in the SSWF NRM to the change in the daily mean, maximum, and 1-in-20 return period rainfall and runoff event to the reference period (1976–2005). Each statistic was evaluated at each grid cell, and then the mean of each statistic was evaluated over the NRM. Even though the mean daily rainfall shows little change or a slight decrease, the maximum daily and 20-year-return-period rainfall indicates a substantial increase under both RCPs, and the increases are projected throughout the ensemble. The pattern of a decrease in mean rainfall and an increase in rainfall extremes is found in almost all other NRM clusters and is supported by results from other studies (Alexander and Arblaster, 2009; CSIRO and Bureau of Meteorology, 2015; Wasko and Sharma, 2017). The projections for the extremes in runoff are less emphatic than those of rainfall; however, like rainfall, mean daily runoff shows decreases (relative to the reference period), while maximum daily runoff and 20-year-return-period runoff indicate a substantial increase under both RCPs (except for 2085 under RCP8.5). Similar analyses were conducted for each of the individual NRM regions.

Figure 15 Relative changes in mean daily, maximum daily, and 20-year-return-period daily rainfall (a) and runoff (b).

8 Communication for end users

To understand user requirements for the National Hydrological Projections service, a user-centred design process was undertaken using an expert consultant approach (Wilson et al., 2022). Over 8 months, interviews were held with 56 potential users of NHP information from 20 organisations across Australia. The organisations represented included water utilities, government departments responsible for infrastructure planning, agriculture and water resources, hydro-electricity generators, and emergency management services. They comprised existing customers of the Bureau of Meteorology and those anticipated as being interested in using a hydrological projections service. Interviews were recorded and transcribed for documentation purposes. Following this process, key user requirements were grouped into three different types: (1) general information about the impacts of future climate change on Australian water resources in the context of historical hydro-climatic trends (general audience); (2) local- and/or regional-scale impact information, including key messages of regional impact water resource assessment and guidance on the use of hydrological projections information (technical and policy audience); and (3) access to application-ready hydrological projections datasets to be able to run their hydrological models and localised assessment analysis (planners and technical users).

To be able to attend to all three customer groups, we developed the following strategy: (1) the provision of application-ready datasets, (2) the creation of a user interface to provide a national picture of changes to key water variables and, (3) the provision of guidance material and regional hydroclimatic assessments.

8.1 Application-ready NHP dataset

The complete foundational NHP dataset is hosted on the Australian National Computational Infrastructure (NCI) via the NCI Data Collection and its Thredds server. The data are made publicly available free of charge under the following link (under the Creative Commons – Attribution 4.0 International – CC BY 4.0): https://doi.org/10.25914/6130680dc5a51 (NCI Data Catalogue, 2024). The dataset consists of 16 datasets per greenhouse gas emission scenario ordered by climate input variable and hydrological output variable (Table 2). Each variable is stored in individual netCDF files following a structured approach that is ordered by the host GCM, historical or RCP scenario, and bias correction method (including data with no bias correction applied). Metadata information is made compliant with CMIP5 metadata standards. The data are accessible either as a OpenDAP, netCDF, HTTP, WMC, or WMS data provision service.

Table 2 Overview of NHP foundational dataset variables available in the NCI Data Collection.

8.2 Australian Water Outlook user interface

The publicly facing portal created by the Bureau of Meteorology called the Australian Water Outlook (AWO) service (https://awo.bom.gov.au, last access: 9 April 2024) connects people with interpreted key information on conditions and changes to water balance components. The AWO brings together information about the current and historical state of water in the Australian landscape with near-term and seasonal forecasts and information about the impacts of climate change on water resources (Wilson et al., 2022). The portal enables users to explore projections information in a meaningful way that is selectable for their region of interest. Hydroclimatic changes are expressed for rainfall, evapotranspiration, soil moisture, and runoff as absolute and relative change signals from the reference period of 1976–2005. The underlying NHP data are visualised as maps of mean changes for rainfall, runoff, soil moisture, and evapotranspiration for multiple future periods (30-year periods centred at 2030, 2050, 2070 and 2085) and for the two greenhouse gas emission scenarios. The portal also provides the ensemble statistics showing the spread across the NHP member ensemble. The data in the interface is supported by regional reports (see Sect. 8.3) that contextualise the information against the current climate and the development of case studies to demonstrate the potential application of the information. In addition, non-technical users are provided with digestible information on scientific and technical background information about the hydrological projections through guidance material (e.g. FAQs and “About” sections). Feedback from a range of key users during multiple stages of the portal development ensured that the information is fit for purpose. This co-design process also guided the user-friendly accessibility of the information presented. Maps, charts, information, and data can be freely downloaded from the AWO interface under the Creative Commons – Attribution 4.0 International – CC BY 4.0 license.

8.3 National hydrological projections assessment reports

To understand future impacts on Australia's water resources, we prepared eight tailored hydrological change assessment reports on plausible future changes in rainfall, potential evapotranspiration, soil moisture, and runoff, including an analysis of wet and dry extreme conditions. These assessment reports are based on the eight clustered NRM regions (see Fig. 5) and use the information and analysis presented in Sect. 7. These NRM regions broadly represent groups of similar climate and biophysical settings in Australia and corresponding natural resources. The reports build on and are consistent with similar reports of CCiA, which provided comprehensive, robust information on future climate changes for each of the NRM regions (https://www.climatechangeinaustralia.gov.au/en/overview/methodology/nrm-regions/, last access: 9 April 2024). Our work builds a complementary picture in the context of changes to the regional hydrologic cycle and its future impacts.

The assessment reports provide an overview of the key findings for each of the NRM regions. They describe each NRM region, including background information about physiographic and hydroclimatic characteristics, recent conditions, and long-term hydroclimatic trends. The ability to simulate such changes is also examined by investigation of several questions: (1) are the climate models chosen able to represent the region's climate, (2) how well does the hydrologic AWRA-L model perform in the region, and (3) how does the evaluation of the bias correction methods affect the interpretation of the projections? Together, addressing these questions provides important context when assessing each NRM's hydroclimate. Results from the 16-member ensemble (for each RCP) are presented as plausible future representations of change in magnitude and sign for each of the hydrological output variables.

We use a storyline approach to demonstrate the applicability of the NHP data for future water resource impact analysis across Australia and to address the uncertainty in hydrological projections (Shepherd et al., 2018). The storyline approach describes an internally consistent evolution of plausible future events, allowing for a way to focus on assessing only uncertainties that relate to a specific impact of interest. By establishing a set of storylines that represent plausible changes in risk, the impacts of a particular risk can be explored (as opposed to a forecast or likelihood of specific outcomes). They are effective for assessments where the change to risk is triggered by the interaction of multiple variables, such as water supply and demand. For example, applying the storyline concept to the Wet Tropics region in northern Australia, we used the wet-season runoff as an indicator for water storage (and hence supply) since almost all infilling occurs in the wet season. Changes to soil moisture in the dry season were used as an indicator of changes in demand due to increased agricultural irrigation and domestic and consumptive use. Readers are referred to the Australian Water Outlook (https://awo.bom.gov.au/about/overview/assessment-reports, last access: 9 April 2024) to explore the regional assessment reports in further detail.

9 Limitations of the hydrological projections

A factor not considered in the NHP projections was vegetation responses to increased CO₂ and the resultant effects on projections of runoff. Offline models generally indicate an increase in terrestrial drying and drought metrics, especially those based on potential evapotranspiration (Naumann et al., 2018; Scheff and Frierson, 2015). They also predict decreases in runoff that exceed those projected by GCMs (Milly and Dunne, 2016). This has been attributed to potential evapotranspiration being overestimated in offline models since they do not capture stomatal conductance reductions in response to elevated CO₂ (Roderick et al., 2015; Yang et al., 2019). AWRA-L was chosen as the hydrological model based on the evaluation and benchmarking of the available national models (Frost and Wright, 2018). It was also selected for its ability to model projected hydrological variables when considering the larger uncertainties inherent to the climate models and the static soil and vegetation inputs (Azarnivand et al., 2022). Importantly, those evaluations considered runoff, soil moisture, and actual evapotranspiration in the assessment of the models. However, AWRA-L is not a coupled model and was run independently using the bias-corrected GCM climate data as input.

The lack of feedback between the GCM's and AWRA-L means that the potential roles of increased CO₂ levels on vegetation growth and evapotranspiration rates are not considered (Yang et al., 2019). It is also uncertain how important vegetation feedbacks are in a mostly water-limited environment such as Australia, as the sensitivity of evapotranspiration parameterisation is more important in humid regions (Zheng et al., 2019a). In addition, recent observations indicate a robust increase in the partitioning of rainfall to evapotranspiration rather than runoff in response to the warming over land (Pascolini-Campbell et al., 2021). Finally, future land use change and vegetation change due to future temperature and water availability changes are also not considered and contribute to the uncertainties in the hydrological projections (Prestele et al., 2017). Despite these caveats, we have developed these datasets and complementary analysis for the wider community. In particular, they have been developed for those who wish to use hydrological projections for climate adaptation and for those in the scientific community who wish to further interrogate these datasets for research purposes.

10 Conclusions

This article has documented the National Hydrological Projections project (NHP), which resulted in the production of datasets and communication tools for end users. We developed a 16-member ensemble of bias-corrected GCM output which was used to force an offline hydrological model (AWRA-L). The ensemble was constructed by applying three bias correction techniques to the output of four GCMs and one of the bias correction methods to the output of the CCAM RCM forced by the same set of four GCMs. The four GCMs were CMIP5 model simulations (Taylor et al., 2012) using mid-range (RCP4.5) and high-emission (RCP8.5) scenarios (IPCC, 2000). The bias-corrected GCM variables were those required as input to force the AWRA-L hydrological model, namely precipitation, minimum and maximum temperature, surface wind speed, and downwelling solar radiation. The subsequent hydrological output consisted of runoff, soil moisture, and potential evapotranspiration. The bias correction methods were applied to historical (1960–2005) simulations using 0.05×0.05 (approximately 5 km × 5 km) gridded reference data at a daily timescale. The resulting calibrations were then applied to the GCM output at a daily timescale for the period 2005–2099, and the subsequent AWRA-L simulations were produced on the same spatial and temporal scales. This is the first time that a nationally consistent set of hydrological projections have been produced for the entire Australian domain.

The bias-corrected GCM data and AWRA-L outputs were evaluated against historical reference datasets and shown to produce simulations that matched the historical record (i.e. the bias correction methods had allowed the GCM output to be suitable to force the AWRA-L hydrological model). Due to the large spatial domain, only the output of four GCMs was able to be corrected for bias; however, they were chosen using guidance from the Climate Change in Australia report (CCiA; CSIRO and Bureau of Meteorology, 2015), which suggested their inclusion for impact studies based on several criteria (see Box 9.2, p. 180 in CCiA). Furthermore, we explored the ensemble range of precipitation and surface temperature using the NHP subset of GCMs and found they spanned a similar range of wet and dry and cool and warm projections to the larger (40-member) ensemble of CMIP5 models used in CCiA. This resulted in the NHP having ensemble medians and a 10th–90th percentile ensemble spread of similar magnitude to the CCiA ensemble, although they showed slightly cooler and wetter median projections of surface temperature and precipitation. However, the application of the bias correction methods was found to apply corrections to the raw GCM data that were model and bias correction dependent and hence expanded the range of projections from consideration of the raw GCM data alone. Furthermore, the use of the CCAM RCM, which often had projections and temporal meanderings quite distinct from the host GCM, expanded the phase space of the original four host GCMs considerably.

To maintain some consistency with the CCiA report, we explored projections of precipitation, runoff, soil moisture, and potential evapotranspiration for the Natural Resource Management (NRM) regions. For this paper, we used the Southern and Southwestern Flatlands (SSWF) as an exemplar to present a range of analyses performed for the NHP. These included spatially averaged time series, spatially and temporally averaged bar plots of the change signal (future minus historical), and temporally averaged plots of the spatial change signal. To construct the change signal, we considered four future time slices of 30 years centred on 2030, 2050, 2070, and 2085 for RCP4.5 and RCP8.5 and analysed the magnitude and direction of change to provide detailed hydrological projections for each NRM. Readers are referred to the Australian Water Outlook (https://awo.bom.gov.au/about/overview/assessment-reports, last access: 9 April 2024) to explore the regional assessment reports in further detail. For further detailed regional analysis, guidance on the use of NHP data or further general information, please contact us via water@bom.gov.au.