Assessing the predictability of MLR models for long-term streamflow using lagged climate indices as predictors: a case study of NSW (Australia)

The geographic location and extensive topographic variations present high climatic variability in Australia resulting in even higher inter-annual streamflow variability across the country, which is almost twice that of the rivers in any other part of the world (McMahon et al. 1992). As a consequence, the irrigators, agricultural producers, water managers, and planners have to undergo many problems to allocate irrigation water and environmental flows, manage and operate reservoirs, supply municipal water, estimate future hydroelectricity supply, etc. One such severe impact was the overall reduction of gross domestic production (GDP) by 1.6% during the 2002 to 2003 Australian drought (Horridge et al. 2005). According to Dutta et al. (2006), streamflow forecast is more significant compared to rainfall forecast as it can be predicted with longer lead times, therefore it enables the water users to make decisions at an earlier stage of the year, which ultimately increase the potential of financial benefits.

Two main sources of streamflow forecasting are initial catchment condition (antecedent streamflow, antecedent rainfall, soil moisture, etc.) and climate variables, i.e., rainfall, climate indices, etc. (Robertson & Wang 2009). While comparing to initial catchment condition, remote climate drivers have better predictability of streamflow as the climate indices fluctuate at very low frequencies which can impact the streamflow easily. Moreover, developing streamflow forecasting models incorporating initial catchment conditions is more complicated.

Australia, being surrounded by the Pacific, Indian, and Southern Ocean, is greatly influenced by the climatic anomalies originating from the oceans. The impacts of climate indices including the sea level pressure (SLP) and sea surface temperature (SST) anomalies have spatial and seasonal variations. The climate of south-east Australia is influenced by four major climate drivers originating in the Pacific Ocean, the Indian Ocean, and the Southern Ocean, which are ENSO (El Niño–Southern Oscillation), IPO (Inter Decadal Pacific Oscillation)/PDO (Pacific Decadal Oscillation), SAM (Southern Annular Mode), and IOD (Indian Ocean Dipole) (Duc et al. 2017).

The ENSO phenomenon, which results from the large-scale interactions between ocean and atmospheric circulation processes in the equatorial Pacific Ocean, has direct influences on the climate variability over many parts of the world (Ropelewski & Halpert 1987; Kiladis & Diaz 1989; Nicholls et al. 1991). El Niño and La Niña events are responsible for the different climatic conditions around the Pacific including eastern Australia (Stone & Auliciems 1992; Latif et al. 1997; Nazemosadat & Cordery 1997; Hoerling et al. 2001; Chiew 2006; CPTEC 2006). Several studies have revealed the influences of ENSO on streamflow throughout Australia (Chiew et al. 1998, 2003; Piechota et al. 1998; Dettinger & Diaz 2000; Dutta et al. 2006; Sharma et al. 2015). Chiew et al. (1998) and Piechota et al. (1998) found that ENSO-based (SOI and SST) streamflow predictions in northeast Australia are better than the forecasts from climatology. ENSO anomalies are found to be the strongest predictors of seasonal (spring) streamflow and rainfall in some studies (McBride & Nicholls 1983; Robertson & Wang 2009). The study of Chiew & Leahy (2003) explained that spring rainfall and runoff had high correlation (0.3 to 0.5) against winter SOI throughout eastern Australia.

The El Niño–Southern Oscillation Modoki events have significant influences on the climate of many parts of the world including Japan, New Zealand, western coast of United States (Ashok et al. 2007), Australia (Taschetto & England 2009), and South China (Feng & Li 2011). According to Taschetto & England (2009), classical ENSO causes reduction in precipitation in north-east and south-east Australia, whereas EMI (El Niño–Southern Oscillation Modoki Index) is responsible for reducing precipitation in north-west and northern Australia regions.

Although the dominant source of inter-annual variability in Australian rainfall and streamflow is believed to be the ENSO phenomenon, some recent evidence shows that Eastern Australia is also influenced by IOD and interdecadal modulation of ENSO as a result of the low frequency variability in the Pacific Ocean, which is referred to as Pacific Decadal Oscillation or PDO (Westra & Sharma 2009). Cai et al. (2011) and Risbey et al. (2009) found that IOD has an impact on austral winter (June to October) in the southern part of Australia, whereas ENSO has a strong influence on austral spring rainfall as a result of the strong covariation of ENSO and IOD.

Many researchers (e.g., Power et al. 1999; Kiem et al. 2003; Kiem & Franks 2004) have demonstrated the influence of the IPO to be significant on rainfall and streamflow variation on a decadal to multidecadal timescale. King et al. (2013) suggested that IPO played a significant role in the frequency of major floods during the 1950s, 1970s, and 2010–2011. Verdon et al. (2004) explained the enhanced rainfall and streamflow in eastern Australia as being the consequences of the combined impact of ENSO (La Niña) and IPO negative phase. IPO and PDO indices are suggested to be highly correlated and useful in explaining various warming and cooling phases in both northern and southern hemispheres (Power et al. 1999; Franks 2002b).

Whiting et al. (2003) studied the rainfall in Sydney and demonstrated the existence of a greater correlation of annual rainfall in Sydney with the PDO index than with SOI. A combination of correlation and wavelet-based methods was applied to identify the principal sources of variation in reservoir inflows of Sydney (Westra & Sharma 2009) which found ENSO, PDO, and IOD to be influential.

To date, most research has been focused on identification of suitable predictor variable(s) for forecasting rainfall and streamflow on daily or monthly scales (Sharma 2000; Kiem & Franks 2001; Ruiz et al. 2007; Robertson & Wang 2008, 2009; Heller et al. 2009; Kirono et al. 2010), while very few established the seasonal relationship in different parts of Australia (Piechota et al. 1998; Owens et al. 2003; Wang et al. 2009). The majority of the previous studies investigated the concurrent relation of single climatic variable with daily, monthly, or seasonal streamflow. Even though some studies considered the lagged climate modes, none explored the combined impact of different climate anomalies on seasonal streamflow of eastern Australia. Although a few of the previous studies considered lagged and combined impact of different climate anomalies on seasonal streamflow of eastern Australia, most of those were probabilistic approaches.

Some recent attempts were made by Abbot & Marohasy (2015), Mekanik et al. (2013), and Rasel et al. (2016) to forecast seasonal rainfalls for Queensland, Victoria, and South Australia, respectively, exploiting linear (MLR) and non-linear artificial neural network (ANN) techniques and considering the combined influence of different climate modes. Kiem & Verdon-Kidd (2009) explored the relationships between large-scale climate drivers and Victorian streamflow and found that ENSO alone is able to explain a very small part of the variation in Victorian streamflow. Therefore, it is important to investigate the combined influence of multiple climate drivers on seasonal streamflow.

Various studies have intended to provide a probabilistic forecast (Piechota et al. 1998; Ruiz et al. 2007; Robertson & Wang 2009; Wang & Robertson 2011; Duc et al. 2017), whereas to solve the water management problems, a deterministic streamflow forecast is more useful as knowing the expected amount of future streamflow helps stakeholders to make more accurate decisions knowing the expected amount of future streamflow. It is to be noted that the Australian Bureau of Meteorology (http://www.bom.gov.au/water/ssf/index.shtml) provides seasonal streamflow forecasting using a Bayesian joint probability (BJP) method, which is again a probabilistic approach. Therefore, the present study aims to investigate the extent of interactions of large-scale multiple climate mode with seasonal streamflow of New South Wales (NSW) with a view to exploiting these relations to forecast seasonal streamflow with a deterministic output.

NSW, located in the south-eastern part of Australia, is affected by frequent droughts and floods, especially in the western and north-eastern regions of the state. Researchers (Erskine & Warner 1988; Franks 2002a; Franks & Kuczera 2002) have quantified the variability of flood risk and attributed these to ENSO and IPO. Climate variability has serious impacts on the yield of planted crops like wheat, rice, etc. in this agriculturally important region. Although efforts have been made to forecast streamflow and rainfall, none of the current practices provide reliable seasonal streamflow forecast, which can enable water stakeholders to make low-risk decisions at the early stage of the crop period (Khan et al. 2005).

This research is intended to provide a successful deterministic prediction of spring streamflow of NSW using the potential large-scale climate drivers with the help of the statistical MLR method. To accomplish this objective, NSW has been divided into four geographically distinct regions and for each region, at several streamflow measurement stations, investigations were conducted to identify the best predictor(s) for forecasting seasonal (spring) streamflow several months ahead in these regions. At this stage, current analysis is conducted for spring season only, considering the outcomes of the previous research works (McBride & Nicholls 1983; Robertson & Wang 2009). A preliminary study of this research on the concurrent analysis of seasonal streamflow and climate indices found spring streamflow provided the most significant correlations compared to other seasons (not shown in this paper). Thus, spring streamflow is expected to have better interactions with climate indices in the lagged correlation analysis.

NSW situated on the east coast of Australia covers a land area of 880,0642 km² and is the most populous state of Australia with a population of 7.5 million, two-thirds of which live in the Greater Sydney Area. The state is bordered on the north by Queensland, on the west by South Australia, on the south by Victoria, and on the east by the Tasman Sea. The two most important features of NSW are the Great Dividing Range (GDR) and Murray Darling Basin (MDB), which accounts for nearly 40% of the value of agricultural production in Australia and 65% of the irrigated land. NSW possesses almost 61% of the water resources plan area of MDB (Department of Industry, NSW Government) while all of Australia's irrigated rice is produced by Murrumbidgee and NSW Murray irrigation regions Murray-Darling Basin Authority (2010).

Economically, NSW is the most important state of Australia as it contributes most of Australia's agricultural production, which is spread throughout the eastern two-thirds of the state. According to ABARE, 80.92% of the state is agricultural land which contributed 23% of the total gross value of agricultural production in Australia in the year 2015–16. Considering the geographical location, regional climatic variation, and agricultural importance, for the current study NSW was divided into four regions: Northern New South Wales (NNSW), Southern New South Wales (SNSW), Central West New South Wales (CWNSW), and Western New South Wales (WNSW), as shown in Figure 1. To explore the spatial variation of influences of different climatic variables for each region several stations were selected based on their long data records and fewer missing values (Table 1). A total of 12 stations were chosen with data records considered to be of appropriate length for the statistical analysis carried out in this study and also considering the predominance factor of coastal rivers in eastern Australia (Verdon et al. 2004) as well as the consumptive water use of the locations, where streamflow predictions are important. It can be seen from Figure 1 that the locations of the streamflow stations provide a good spatial coverage of NSW.

For this study, historical streamflow data were collected from the Australian Bureau of Meteorology (http://www.bom.gov.au/waterdata/). Observed monthly streamflow (in cumec) was collected for 102 years, ranging from 1914 to 2015 for nine stations while 99, 101, and 88 years of data were collected for North Cuerindi, Wee Jasper, and Mittangang Crossing stations, respectively. These stations have less than 0.5% missing values, which are filled by the series mean of the streamflow data. Using these data, seasonal mean discharge data are derived for the spring (September–October–November) season.

Considering the aforementioned research works on rainfall and streamflow in this region and the outcomes of single concurrent and lagged correlation analyses, four climate drivers: ENSO-based SST anomalies NINO3.4, EMI, IPO/PDO, and DMI (IOD) were selected for the MLR analysis.

The ENSO phenomenon has two components, SST and atmospheric pressure, which are strongly correlated and can be represented by two types of indicators, the SLP indicator and the SST indicator (Duc et al. 2017). Various SST anomalies are available which can be derived from different areas of the equatorial Pacific Ocean (Kiem & Franks 2001). Generally, the SST anomalies are monitored in three geographic regions (Figure 2) of the equatorial Pacific and defined as NINO3 (5°S–5°N, 150°–90°W), NINO3.4 (5°S–5°N, 170°–120°W), and NINO4 (5°S–5°N, 160°–150°W) (Risbey et al. 2009) along with NINO1+2 (0–10°S, 90°–80°W). Hanley et al. (2003) compared the response of pressure-based and SST-based anomalies to ENSO extreme events and found that NINO3.4 and NINO4 indices are equally sensitive to El Niño events, whereas SOI is less sensitive to La Niña events than others.

The IPO is described as the Pacific ENSO-like pattern of SST which is found in the analysis of near-global inter-decadal SST (Folland et al. 1999). IPO has a cycle of 15–30 years and is characterized with two phases, namely, positive and negative (Salinger et al. 2001; Henley et al. 2015). While IPO is used for the whole Pacific Basin, PDO is used for the North Pacific, poleward of 20°N.

The five oceanic and atmospheric climate indices data were obtained from Climate Explorer website (http://climexp.knmi.nl), while the EMI data were collected from the website of JAMSTEC (http://www.jamstec.go.jp/frcgc/research/dl/iod/modoki) for a duration of 102 years (1914–2015). An overview of the used climatic variables for the current analysis is presented in Table 2.

In the current study, linear relationships between the selected climatic variables and spring streamflow of NSW were explored using a multiple linear regression technique. At first, single concurrent and lagged correlation analysis was performed to explore dominating indices on spring streamflow. Subsequently, MLR analysis was performed by developing MLR models incorporating combined influences of two indices.

There are several techniques for exploring relationships between two or more parameters. Regression analysis is one of the popular statistical approaches and is highly recommended for this kind of analysis (Ruiz et al. 2007; Mekanik et al. 2013; Pumo et al. 2016). The most commonly used form of linear regression is multiple linear regression analysis. MLR models establish a statistical relationship between two or more explanatory variables and a response variable and provide a linear equation as output which represents the significant correlation among the variables. In every equation, the value of every independent variable (NINO3.4, EMI, IOD, or IPO/PDO for the current study) is associated with the value of dependent (spring streamflow for current study) variable. In many studies, climate forecasting has been undertaken using the MLR model, due to the fact that this model comprises many regressors to deal with the time series data base.

In order to ensure independence of the residuals error of the model, the Durbin–Watson (DW) test was performed, which assesses the serial correlation between errors. DW parameter has a range of 0 to 4; a value of less than 1 or greater than 3 is certainly a matter of concern (Field 2009).

The performance of the developed MLR models has been assessed by several statistical performance measures which are widely used for the evaluation of regression models. Statistical measures, namely, root mean square error (RMSE), mean absolute error (MAE), Pearson correlation coefficient (r), and Willmott index of agreement (d) are exclusively chosen for this study.

A detailed study of past research works revealed that different climate anomalies have impacts on the seasonal streamflow of NSW which varies seasonally as well as spatially. Pearson correlation analysis was done for identifying the strength of relation between the seasonal climate anomalies. To analyze the extent of the influences of different climatic variables on seasonal streamflow of NSW, Pearson correlation analysis was applied in two different segments. The first segment took into account the concurrent relationship of climate indices and seasonal streamflow, while the second segment was conducted in order to investigate the lagged relationships of the climate indices and seasonal streamflow. Later on, the outcomes of the second phase served as the basis for selecting the suitable lagged climate indices to use as inputs while developing MLR models.

The linear relationships between spring, summer, autumn, and winter streamflow and climate anomalies of the same seasons were explored by doing the Pearson correlation analysis for each station of the four selected regions of NSW.

In Table 3, concurrent correlations between seasonal streamflow and climate indices are presented, where it is observed that spring (Table 3(d)) shows stronger correlations with most of the indices compared to other seasons (Tables 3(a)–3(c)). It is evident from Table 3(d) that all the climate indices (except IPO) have significant correlations with spring streamflow for all the selected regions. IPO shows significant correlations only for CWNSW and for two stations of SNSW regions which are geographically close to each other. Therefore, it can be anticipated that only in the central-west and southern parts of NSW IPO has strong influence on spring streamflow.

For each selected station of the four study regions of NSW, single lag correlation analysis was performed between spring streamflow at year ‘n’ and monthly (December_n-1 to August_n) values of each climate variable. The outcomes are presented in Tables 4 and 5.

It was observed from the single lagged analysis that different regions are influenced by different climatic variables. Lagged NINO3.4 and PDO have significant impacts on the spring streamflow of all four selected regions. Spring streamflows of CWNSW stations and two stations of SNSW (Gundagai and Wee Jasper stations, which are closely located to CWNSW stations) were also influenced by lagged IPO (Tables 4 and 5). Almost all the stations had significant correlations with lagged IOD indices, whereas lagged EMI had impacts on a very limited number of stations.

NINO3.4 shows statistically significant correlations up to a lagged period of four months (April to August) with the correlations ranging from −0.19 to −0.46 (Tables 4 and 5). These findings align with the study of Kirono et al. (2010), who found significant correlations for another two ENSO indices – NINO3 and SOI as 0.35 and 0.36, respectively, in the same study regions. SOI was found to have a correlation of 0.51 with spring streamflow in the study of Chiew & Leahy (2003). Wang et al. (2009), also explored the strong impact of NINO3.4 on spring rainfall up to a lag of two months. Duc et al. (2017) explained that ENSO indices have a strong impact on rainfall during austral spring. The lagged periods of IOD that had significant correlations were not consistent for all the stations, while most of the stations have significant correlations up to three months' lag. It is observed that EMI has significant correlation only up to two months' lag (Tables 4 and 5) for five stations (Wee Jasper, Gundagai, Wagga Wagga, Cowra, and Brewarrina) with a maximum significant correlation of −0.41 (Brewarrina station). The lagged relationship of PDO was quite different as for some stations it presents more significant correlations with longer lead times. The maximum lead times (up to eight months) with significant correlations for spring streamflow were also obtained for this climate variable of the Pacific Ocean (Tables 4 and 5). This is similar to the assessment of Whiting et al. (2003), who stated that PDO has greater correlation with annual rainfall of Sydney than that of SOI. Latif et al. (1997) also showed a strong relationship between PDO and Australian summer monsoon. Westra et al. (2008) evaluated correlation coefficients between seasonal inflows of a reservoir in Sydney and climate indices, where correlations of spring inflow with NINO3.4 and PDO were found to be −0.17 and −0.19, respectively.

Various models with different combinations of lagged months' indices were analyzed for all 12 stations in order to find out the best forecasting model for each of the four study regions. Eighty-five years (from 1914 to 1998) of streamflow data were selected for the calibration of the models, while the remaining 17 years' (from 1998 to 2015) data were selected for validation in order to assess the future streamflow predictability of the developed MLR models. For all the stations, the best models with lower errors while satisfying the statistical limits were selected. A similar approach for validating the results was applied by Hossain et al. (2018), Rasel et al. (2016, 2017), and Mekanik et al. (2013) while predicting rainfall using climate indices. The Pearson correlation (r) values for both calibration and validation periods were satisfactory, while the correlation values in the validation period were found to be higher than that of the calibration period. F-test was performed to evaluate the best model that fit the population of the sample data, while the t-test was conducted to identify the significance of the individual parameters. The best model for each station with their regression coefficients, VIF, and DW statistics are presented in Table 6.

It can be seen from Table 6 that VIF values for the selected models are close to 1, which means that there is no multi-collinearity problem between the predictors. According to Field (2009), values less than 1 or greater than 3 for DW test will indicate the presence of serial correlations between the model errors. Thus, it can be concluded from the results of Table 6 that the DW test of each selected model satisfies the statistical limits, which also establishes the goodness-of-fit of the models.

The results of the MLR analysis depict a clear view of regional variation in influence of combined multiple models throughout the study area. A good number of models combining PDO and NINO3.4 show statistically significant correlations with spring streamflow for the NNSW region which implies the strong impact of these two indices in the northern part of the state. Therefore, for this region the best model for forecasting spring streamflow was obtained with the combination of two months' lagged NINO3.4 and five months' lagged PDO (at Singleton station) showing a significant correlation of 0.41. PDO-NINO3.4 combined models were predominant in the southern part of the state also, where two stations, Kiosk and Mittagang Crossing in SNSW, showed the highest significant correlations (0.41 and 0.49, respectively) for the same combination. On the other hand, the other two stations (Wee Jasper and Mittagang Crossing) of SNSW had maximum significant correlations (0.57 and 0.51, respectively) with IOD_NINO3.4 and IPO_NINO3.4 combinations, respectively. This is similar to the findings obtained in the western part of the state because except for Cowra (in CNSW), all the other four stations (Corowa, Wagga Wagga, Barham, and Brewarrina) are evidently influenced by either IPO or IOD interaction with NINO3.4. Both the two stations, Barham and Brewarrina in WNSW, are influenced by the same combination which is IOD_NINO3.4 having maximum significant correlations of 0.44 and 0.56, respectively. IPO is found to be dominant in CWNSW for two stations, as Corowa and Wagga Wagga are influenced by IPO_IOD and IPO_NINO3.4 combined models, respectively. Compared to the other two stations of CWNSW (Corowa and Wagga Wagga), Cowra is located further north, where the PDO_NINO3.4 models were found to be dominant and satisfying this fact, Cowra also had maximum significant correlation (0.44) for the same combination. It is evident from the results that the IOD and IPO containing MLR models showed good performance with significant correlation for the southern and western parts (inland) of the state. However, in the coastal eastern part of the state IOD–NINO3.4 combinations were observed to be less effective, which align with the findings of Pepler et al. (2014), who stated that eastern seaboard rainfall is less influenced by tropical SST variability such as ENSO and IOD than inland because the effect of the IOD opposes ENSO during the cool season. In the western part of the country IOD–NINO3.4 and PDO–NINO3.4 combined models can be confidently used to forecast spring streamflow. Thus, interactions of NINO3.4, IOD, and IPO can be influential in the inland of NSW, whereas PDO and NINO3.4 interactions can be effective in the north-eastern and south-eastern coastal regions of the state.

In order to determine the accuracy of the developed MLR models, validation tests were performed. Table 6 shows the performance statistics of RMSE, MAE, and index of agreement (d) of the best developed models for the calibration and validation periods. It is evident from Figure 3 and Table 6 that there is significant increment (except for Kiosk station) of the correlation values from calibration to validation stage such as, for example, at Singleton station, where correlation value increased from 0.41 in the calibration stage to 0.65 in the validation stage for PDO_March–NINO3.4_June combination. Again, the highest correlation (r = 0.71) from the current analysis was also obtained for the Singleton station with PDO_Jan–NINO3.4_May combination in the validation period (not shown in Table 6). However, to get the best predictor model, a few unusual events which were outliers in a box-plot analysis were removed from the calibration and validation periods. Thereby, the ratio of the duration of calibration or validation period to the number of outliers may have affected the corresponding correlation (r) values.

The predictability of the best MLR models for each of the four regions has been explained through the time series plots of observed and simulated flows in Figure 4.

In Figure 4 some overestimations of the models can be observed during the validation stage which may result from the ‘millennium drought’ (Bond et al. 2008) periods that occurred from 1994 to 2010 over the continent. It was explored by Verdon-Kidd & Kiem (2009) that a combination of climate drivers in the Pacific Ocean (ENSO, PDO), IOD and SAM were responsible for the past three droughts in south-east Australia; the ‘Federation’ drought (1895–1902), the ‘World War II’ drought (1937–1945), and the ‘Big Dry’ (1994–2010). Again, some streamflow events (e.g., 1950–1960, 1970–1980) are underestimated by the developed models. This may happen due to the failure of capturing the unusual extreme flood events, for instance, the floods of 1954 and 1974. One analysis by Lismore City Council (https://www.lismore.nsw.gov.au/cp_themes/default/page.asp?p = DOC-SVI-55-40-21) revealed that the flood of 1974 is considered to be a once-in-70-year event which resulted from the simultaneous occurrence of La Niña condition and negative IOD phase. Apparently, a simple MLR model consisting of only two climate indices (e.g., NINO3.4 and PDO) is not likely to replicate an unusual phenomenon like ‘millennium drought’ or an extreme flood event like that of 1974. Another reason is that some other climate indices might have been more influential at that time rather than the selected indices for this study.

The capability of the developed models of forecasting spring streamflow with higher accuracy has been ensured as the values of RMSE, MAE, and d in the validation period showed good agreement with the calibration period. The index of agreement (d) for both the calibration and validation periods was close to 0.5, which ensured good forecasting ability of the models. Significant increase in the Pearson correlation values has proved that the combined models have greater skills for predicting streamflow than the single lagged indices. For instance, at Singleton station, during single correlation analysis, NINO3.4_June and PDO_March showed correlations of −0.43 and −0.29, respectively (Table 4), whereas during MLR analysis, NINO3.4_June–PDO_March combination (Table 6) showed a higher correlation of 0.65 (in validation period). Moreover, while comparing the outcomes of the present research with the previous research studies on forecasting streamflow, it is evident that MLR models of this study showed higher correlations than any of the analyses that used single lagged index. Even for single lagged correlation analysis, in general, the present study showed higher correlations than the previous research (Table 7).

The variations of influences of different climate indices on different study regions of NSW is comparable with the recent study outcomes of Duc et al. (2017). They have reported association of climate indices with NSW rainfall using Bayesian model averaging. Among their studied sites, outcomes of the sites which are within 160 km of our selected streamflow stations are similar to our findings. They have reported that single IPO cannot impact NSW rainfall significantly, however its association with ENSO is significantly influential on rainfall of almost the whole of NSW. The current study also found strong influence of PDO–NINO3.4 on spring streamflow across almost the whole state (it is to be noted that IPO and PDO are similar as IPO acts on the whole Pacific basin and PDO is active in the North Pacific, poleward of 20°N). This finding is strongly supported by the findings of many past studies that suggested IPO or PDO phases modulate the frequency and magnitude of ENSO events (Power et al. 1999; Folland et al. 2002; Franks 2004; Verdon et al. 2004) which is influential on the streamflow volumes of many parts of the world (Kahya & Dracup 1993; Moss et al. 1994; Piechota & Dracup 1996; Chiew et al. 1998; Piechota et al. 1998; Dettinger & Diaz 2000; Kiem & Franks 2001; Wooldridge et al. 2001).

In Wagga Wagga, evidence of strong IOD influence has been found in the study of Duc et al. (2017), which is consistent with our study, as IOD combined models performed significantly near this area. Similar outcomes have been obtained for Mittagang Crossing station, where Duc et al. (2017) found ENSO–IOD combined impact on spring rainfall to be very strong (posterior probability = 1) and the current study has found (not shown here) significant correlation (r = 0.42) with the same combination of indices, i.e., NINO3.4 and IOD.

In the current research, the multiple linear regression method was applied with a view to exploring the potential skills of combined multiple climate indicators to forecast the spring streamflow of NSW regions with a longer lead time than the usual practice. Before performing multiple regression analysis, first single correlation analyses were performed to identify potential climate predictors for the region. Through single correlation analysis, several indices (NINO3.4, IOD, EMI, PDO, and IPO) were found to have strong effects on spring streamflow of NSW with a lagged time of maximum three months. Some indices were found to give significant correlations with a lagged time of more than two months, however, in general, the correlation values decrease with the increase of lagged months. This study, through findings of five effective climate indices for the region, opened an opportunity to study with more than two indices which no one has ever done for this region. For the current study, to achieve better correlations (prediction capability) different combinations of two (out of the five significant) indices were tested in the multiple regression analysis. It was found that the same combination of indices did not turn out to be best for all the stations/regions, which is reasonable as we are dealing with a large region and the further the distance from a particular station the greater the likelihood of being influenced by other indexes. Also, the combined best models' lagged indices are not necessarily from the same month. In general, among the best combined dual indices, NINO3.4 is found to be significant for all the stations except one (Corowa), PDO is more significant towards the north-eastern and south-eastern coastal region, IPO is more significant towards the central-south, whereas IOD is more significant towards the west of NSW. The best correlation is obtained for Singleton station in NNSW for the PDO–NINO3.4 combined model with a correlation of 0.65 (in the validation period) for the prediction of spring streamflow with two months' lagged period. It is noteworthy that every time the combined model outperformed the models considering a single climate variable in terms of Pearson correlation (r), it was evidence of better predictive skills of the MLR models.

For this study, selections of the best models were based on the significant correlation values in both calibration and validation stages. However, while looking at the time series comparisons between the observed and simulated streamflow values, it is found that the developed models are unable to capture some unusual events like severe droughts or high floods. A simple multiple linear regression model consisting of only two climate indices is not expected to capture the complex relationships between streamflow and climate drivers very well, and thus is not anticipated to provide a very good match with observed values. Moreover, in fact, rainfall and streamflow are also influenced by some other local and/or regional factors (i.e., temperature, humidity, wind speed, soil moisture, etc.), which are not possible to consider in such regression models. Thus, the extension of this research work will include other non-linear techniques (i.e., ANN or Fuzzy logic) as some researchers (Mekanik et al. 2013; Abbot & Marohasy 2015; Rasel et al. 2016) successfully explained the non-linear relationship between rainfall and climate drivers using this technique, although they could not provide any output model which could be useful to stakeholders. Since the relationship between streamflow and remote climate drivers is likely to be non-linear (Piechota et al. 1998), a non-linear model is expected to give better results than a linear model. Also, an extension of the current research can be to explore the effectiveness of incorporating more than two indices in a multiple regression model. Nevertheless, the developed MLR models have the potential to provide indications on the possibility of getting increased or decreased amounts of streamflow and expected magnitude in the future season. Currently, water users in Australia get the seasonal predictions of streamflow just at the beginning of the season, which do not give them enough time for prudent decision making. Moreover, those predictions are stochastic, i.e., the users do not get any estimation of expected magnitude. The developed MLR models for the study area are expected to provide water users and planners with some insight which will enable them to take tactical cropping decisions three months in advance. This sort of study is mainly based on regional climate index/indices applicable for a region. However, a similar concept can be applied to other regions if any such index/indices are found to be effective for other regions.