Statistical Downscaling of General Circulation Model Outputs to Precipitation Accounting for Non-Stationarities in Predictor-Predictand Relationships

This paper presents a novel approach to incorporate the non-stationarities characterised in the GCM outputs, into the Predictor-Predictand Relationships (PPRs) in statistical downscaling models. In this approach, a series of 42 PPRs based on multi-linear regression (MLR) technique were determined for each calendar month using a 20-year moving window moved at a 1-year time step on the predictor data obtained from the NCEP/NCAR reanalysis data archive and observations of precipitation at 3 stations located in Victoria, Australia, for the period 1950–2010. Then the relationships between the constants and coefficients in the PPRs and the statistics of reanalysis data of predictors were determined for the period 1950–2010, for each calendar month. Thereafter, using these relationships with the statistics of the past data of HadCM3 GCM pertaining to the predictors, new PPRs were derived for the periods 1950–69, 1970–89 and 1990–99 for each station. This process yielded a non-stationary downscaling model consisting of a PPR per calendar month for each of the above three periods for each station. The non-stationarities in the climate are characterised by the long-term changes in the statistics of the climate variables and above process enabled relating the non-stationarities in the climate to the PPRs. These new PPRs were then used with the past data of HadCM3, to reproduce the observed precipitation. It was found that the non-stationary MLR based downscaling model was able to produce more accurate simulations of observed precipitation more often than conventional stationary downscaling models developed with MLR and Genetic Programming (GP).


Introduction
General Circulation Models (GCMs) are the main tools used for projection of global climate into the future [1], and they are based on the mathematical representations of physics of the atmosphere, ocean, ice caps and land surface processes [2]. For projecting the global climate into the future, GCMs are forced with various scenarios of future greenhouse gas (GHG) emissions. GCMs can reliably simulate the global climate at a coarse spatial resolution in the order of few hundred kilometres. However, owing to their coarse spatial resolution they cannot adequately simulate the catchment scale climate, which is largely influenced by the subgrid-scale features such as topography, land use, and convective processes [3]. The majority of the catchment scale applications need hydroclimatic data at a much finer spatial resolution than that of GCM outputs [4]. Therefore, to bridge the spatial scale gap between the coarse resolution GCM outputs and the fine resolution hydroclimatic information required for catchment scale applications, dynamic and statistical downscaling approaches are used [5].
In dynamic downscaling, GCM outputs are fed into a Regional Climate Model (RCM) as lateral boundary conditions [6]. The high computational cost is a major limitation in dynamic downscaling. Statistical downscaling is based on the concept of developing empirical statistical relationships between the GCM outputs and the catchment scale hydroclimatic variables [7]. Statistical downscaling depends on the assumption that the relationships between the GCM outputs (predictors-inputs to downscaling models) and the observations of the catchment scale hydroclimatic variables (predictands-outputs of downscaling models) in the past are valid for the future under changing climate [8]. This assumption is called the stationarity assumption of the predictor-predictand relationships (PPRs). Unlike dynamic downscaling, statistical downscaling is associated with much less computational costs [9] and hence it has gained wide popularity in the downscaling scientific community [10].
The projections of hydroclimatic variables produced in statistical downscaling studies are subject to uncertainties that arise from numerous sources (some of these uncertainties are also valid for dynamic downscaling). The sources of uncertainties in a downscaling study include: GHG emission scenarios used to force GCMs, GCMs used to provide inputs to downscaling models, various techniques used for developing downscaling models, non-homogeneity in inputs to downscaling models, specific methodologies employed in developing downscaling models (e.g. how predictors are selected and pre-processed, use of cross-validation instead of traditional calibration and validation, and how outputs are post-processed), stationarity assumption of PPRs and stationarity assumption of bias/bias-correction approaches (bias refers to errors in downscaling model outputs), and the quality of observations against which downscaling models are calibrated [11]. These uncertainties can largely influence the future hydroclimatic projections.
To the date some of the above sources of uncertainties have been investigated in detail and in certain instances potential solutions to reduce the uncertainties that arise from those sources have been developed. As examples: use of multiple GHG emission scenarios to quantify the impact of changes in the anthropogenic GHG emissions on catchment scale climate [12,13,14], scenario neutral approach for reducing the dependence on multiple GHG emission scenarios [15], identification of GCMs that better simulate the climate over the region of interest [16], use of ensemble techniques to combine the outputs of downscaling models produced with outputs of different GCMs [17], and derivation of homogeneous inputs to downscaling models [18]. In statistical downscaling the uncertainties arising from the stationarity assumption of PPRs and stationarity assumption of bias/bias-correction approaches, have not been investigated adequately. However, the validity of the stationarity assumption of the PPRs is questionable under non-stationary climate likely to occur in the future [19].
It has been found that the constants and coefficients in the PPRs developed using multi-linear regression (MLR) technique can be used to detect non-stationarities in the PPRs [20,21]. However, in these studies a methodology to account for the non-stationarities in the PPRs of statistical downscaling models under changing climate was not investigated. If non-stationarities in the PPRs are observed in the past, there is a high likelihood that such non-stationarities in the PPRs will be present in the future, under non-stationary climate [20]. Therefore, there is a need for developing downscaling approaches that can explicitly account for the non-stationarities in the PPRs under non-stationary climate [11]. To the date, there are only a few studies that have investigated potential approaches for handling non-stationarities in the PPRs in statistical downscaling [22]. In these studies, potential approaches for handling non-stationarities in the PPRs in statistical downscaling have been developed: considering the changes in the occurrence frequency of modes of natural variability of climate as an indicator of changes in the climate [23,19], using a moving window approach [22] and employing wavelet transformation [24]. This paper presents a novel approach based on the use of a moving window to incorporate the non-stationarities in the large scale climate in the PPRs in statistical downscaling models. The next sections of the paper provide the details of "Study Area and Data", "Methodology", "Results" and "Discussion" in relation to a case study.

Study Area and Data
The operational area (about 62,000 km 2 ) of Grampians Wimmera Mallee Water Corporation For the calibration and validation (calibration and validation together refer to model development) of the statistical downscaling models, observations of the predictand and reanalysis data pertaining to the predictors are required. For this investigation, daily observed precipitation for the observation stations at Halls Gap, Birchip and Swan Hill were obtained from the SILO database (http://www.longpaddock.qld.gov.au/silo/) of the Queensland Climate Change Centre of Excellence for the period 1950-2010. In the SILO database for infilling the missing daily precipitation values, ordinary kriging method has been used to spatially and temporally interpolate daily rainfall [25]. In order to analyze the temporal and spatial errors in the interpolated precipitation data, independent cross validation technique has been used. These daily observations were then added to produce monthly precipitation totals for the development of downscaling models. In general, water resources allocation modelling is conducted at a monthly time step and hence it needs monthly hydroclimatic inputs. Therefore, the scope of this study was limited to the investigation of downscaling models operated at a monthly time step.
The monthly reanalysis data pertaining to predictors were obtained for the period 1950-2010 from the National Oceanic and Atmospheric Administration / Earth System Research Laboratory (NOAA/ESRL) Physical Sciences Division. The outputs of HadCM3 (a GCM) pertaining to the 20 th century climate experiment (20C3M) were obtained from the Programme for Climate Model Diagnosis and Inter-comparison (http://www-pcmdi.llnl.gov/) for the period 1950-1999. The monthly 20C3M outputs of HadCM3 were used for reproducing the past observed precipitation as it allows the assessment of the ability of the downscaling models (developed with reanalysis outputs) in reproducing the observed precipitation with GCM outputs. HadCM3 is one of the GCMs that can properly simulate the precipitation over Australia and the El Niño Southern Oscillation [16]. Therefore, HadCM3 was used to provide the inputs to the downscaling models.

Overview of methodology
In the non-stationary downscaling approach, initially, a series of PPRs (i.e. downscaling models) were determined with the MLR technique by using a moving window on the past predictor data obtained from the NCEP/NCAR reanalysis data archive and observations of precipitation at each station. Then the relationships between the constants/coefficients in these PPRs and the statistics of past reanalysis data of predictors were determined. The non-stationarities in the climate are characterised by the variations in the statistics of the climate variables over time. Hence, the above process enabled linking the non-stationarities in the climate to the PPRs. These PPRs developed using reanalysis data were then modified according to the statistics of predictor data pertaining to the past climate simulated by the 20C3M outputs of HadCM3, yielding new PPRs corresponding to the past climate characterised by HadCM3. The ability of these PPRs to reproduce the past catchment scale precipitation with the outputs of a GCM was considered as an indication of the success of the non-stationary modelling approach. The series of PPRs developed using the above moving window approach resembles a non-stationary downscaling model and it is denoted as SDM 1(MLR) . The performance of this non-stationary modelling approach was compared with that of a stationary modelling approach. Under the stationary modelling approach, linear and non-linear downscaling models were developed using MLR and Genetic Programming (GP). These stationary downscaling models developed with MLR and GP are denoted as SDM 2(MLR) and SDM 2(GP) , respectively. The stationary downscaling models SDM 2(MLR) and SDM 2(GP) are denoted as SDM 2

Defining an atmospheric domain and selection of predictors
The atmospheric domain in a statistical downscaling study is the region of the atmosphere corresponding to which the large scale atmospheric information is obtained, for providing inputs to a downscaling model. In this investigation an atmospheric domain that covered the region bounded by the longitudes 135˚E-150˚E and latitudes 30˚S-42.5˚S was selected, and it is shown in Fig 3. It should be noted that for all downscaling models the same atmospheric domain was used.
In statistical downscaling, it is the common practice to select an initial pool of predictors called probable predictors that are the most likely to influence the predictand of interest. The set of probable predictors selected for this study consisted of; 500 hPa, 700 hPa, 850 hPa and 1000 hPa relative humidity; 500 hPa, 700 hPa, 850 hPa and 1000 hPa specific humidity; 500 hPa, 700 hPa, 850 hPa, 1000 hPa air temperature; surface air temperature; 200 hPa, 500 hPa, 700 hPa, 850 hPa and 1000 hPa geopotential heights; mean sea level pressure; surface pressure; surface precipitation rate; 850 hPa zonal wind speed; and 850 hPa meridional wind speed. These probable predictors were selected based on the previous downscaling studies over the same study area [17,26,27].
Potential predictors are the most influential predictors on the predictand of interest selected from the pool of probable predictors. The potential predictors vary spatiotemporally as the relationships between the large-scale atmospheric conditions and the catchment scale climate can vary over time and space [24]. In order to identify the potential predictors, initially, for each calendar month, the NCEP/NCAR reanalysis data of probable predictors pertaining to the 42 grid points (shown in Fig 3) and observed precipitation data at Halls Gap, Birchip and Swan Hill for the period 1950-2010 were split into three 20-year time slices in the chronological order. Thereafter, the Pearson correlations coefficients between the probable predictors and precipitation were computed for each time slice and the whole period (period which contained all 20-year time slices) for each grid point in the atmospheric domain. Then the probable predictors which displayed statistically significant correlations at the 95% confidence level with observed precipitation with a consistent sign (no changes in the sign of the correlation over time) in all three time slices and the whole period were selected as potential predictors for each calendar month for each station. For each calendar month the three potential predictors that showed the highest correlations with observed precipitation were selected as inputs to both stationary and non-stationary downscaling models. In the selection of the three best potential predictors, when a predictor showed statistically significant high correlations corresponding to multiple grid locations of the atmospheric domain, only the best correlated location was considered.
Above predictor selection procedure assisted in avoiding the influx of any redundant information to the downscaling models. These three potential predictors selected for each calendar month are called the best potential predictors. In Table 1 the three best potential predictors identified for each calendar month are shown in the descending order of the magnitude of the correlation coefficients which they showed with the observations of precipitation at Halls Gap, Birchip and Swan Hill over the period 1950-2010. As shown in Fig 2 defining an atmospheric domain and selection of predictors was a common step for both downscaling approaches.

Development of non-stationary downscaling models (SDM 1 )
For the development of SDM 1(MLR) , for each calendar month, for each station, two 20-year moving windows were defined in such a way that the 1 st and the 2 nd moving windows initially covered the 1 st and the 2 nd 20-years of the reanalysis data of the three best potential predictors (refer to Table 1) and the observations of precipitation over the period 1950-2010. Then for each calendar month, these two windows were moved at a 1-year time step. The 1 st window was moved until it swept over the whole set of data , and for that, 41 1-year movements were required. Once the 1 st window was moved to the position where it covered the data of the period 1971-1990, the 2 nd window would have covered the last 20-years of data which corresponded to the period 1991-2010. Thereafter, for each movement of the two windows, data were concatenated for the 2 nd window by considering data from 1950 onwards. (e.g. in determining the PPRs for period 1972-1991, the 2 nd window would have covered period 1992-2010 and year 1950).
For each 1-year movement of the two windows and also for their original positions, the reanalysis data pertaining to the three best potential predictors were standardised using their means and the standard deviations. Then using these standardised data, the constant and the coefficients of PPRs in SDM 1(MLR) were determined for each movement of the 1 st moving window by minimising the sum of squared errors between the model-reproduced values and observed precipitation. This phase of the model development is called model calibration. The PPRs of SDM 1(MLR) had the format given in Eq 1, where Y is the predictand and b i is the coefficient of predictor x i , D is the constant and ε is the Gaussian error. Since i = 3 there are three coefficients b 1 , b 2 and b 3 corresponding to the three best potential predictors x 1 , x 2 and x 3 respectively.
Then the PPRs developed for the 1 st moving window were used to reproduce the observed precipitation pertaining to the 2 nd moving window as a validation of the performance of SDM 1(MLR) . The calibration and validation process yielded a series of 42 PPRs for each calendar month for each station (total number of PPRs = 12 x 42) corresponding to the 1 st moving window. The use of moving windows at a 1-year step enabled the identification of fine changes in the PPRs even with a short record of data.
The   A test was then conducted to determine the statistical significance of the linear trends [28] of each coefficient and the constant in the PPRs for each calendar month for each station. The results of the statistical significance test are shown in Table 2. According to Table 2, it was realised that, the linear trends in the coefficients and constants in the PPRs were significant at the 95% confidence level in the majority of the calendar months at all stations. These significant trends in the coefficients and constants in the PPRs are indicative of the influence of non-stationarities in the large scale climate on the catchment scale precipitation. However, it should be noted that the changes in the constant and the coefficients of the series of PPRs could be due to other factors such as natural variability of the climate and data inconsistencies. The identification of such influences on the PPRs is beyond the scope of this study.

Relationships between PPRs in SDM 1 and statistics of reanalysis data
The changes in the statistics of global climate over time are regarded as indications of climate change. Hence, it was realised that relating the changes in the PPRs in SDM 1(MLR) and the statistics of the best potential predictors over time is a potential way to determine the influence of non-stationarities in the climate on the PPRs over time. For determining the relationships between the PPRs in SDM 1(MLR) and the statistics of the best potential predictors, initially, for  Downscaling Accounting for Predictor-Predictand Non-Stationarities each calendar month, for each 1-year movement of the 1 st window, the mean, the standard deviation and the 25 th , 50 th , 75 th , 95 th percentiles of the NCEP/NCAR reanalysis data of the three best potential predictors were computed, for each station. In order to achieve this, 42 PPRs for each calendar month were considered for each station. The 1-year movement of the 1 st window yielded the time series of statistics of the reanalysis data of the three best potential predictors. By computing the correlations between each of these statistics and each coefficient/ constant in the PPRs, the most influential statistic on each coefficient and the constant was determined. Table 3 shows the most influential statistic of the three best potential predictors on each coefficient and the constant in PPRs in SDM 1(MLR) for each calendar month, including their correlation coefficients (CC) for observation station at Halls Gap as an example. It was seen that the majority of the correlations were strong and statistically significant at the 95% confidence level for all stations.

Modification of PPRs in SDM 1 according to statistics of GCM outputs
The next step of the modelling process was to modify the coefficients/constant in the PPRs of SDM 1(MLR) developed with the NCEP/NCAR reanalysis outputs according to the statistics of 20C3M outputs of HadCM3 GCM. For this purpose, initially, the variations of each coefficient and the constant of the PPRs in SDM 1(MLR) and the most influential statistic of the three best potential predictors were visualized using scatter plots.   1990-1999, using linear interpolation as described in the next few paragraphs. In other words, the 42 x 12 PPRs derived using the NCEP/NCAR reanalysis data were used with the 20C3M outputs of HadCM3 for the derivation of 3 x 12 new PPRs (i.e. one PPR for each calendar month per time slice per station). It should be noted that the non-stationary downscaling models developed in this investigation are quasi-models, as the PPRs determined using the 20C3M outputs of HadCM3 for each period were stationary throughout that period. As a demonstration for the application of linear interpolation to derive new values for the coefficients/constant in the PPRs of SDM 1(MLR) , the following example is considered. In January, the 25 th percentile of surface precipitation rate at grid location {(4,4)} was the most influential statistic (correlation coefficient = -0.85) on the first coefficient (b 1 ) of the MLR based PPR for Halls Gap station (see Table 3). As shown in Fig 5a, let the 25 th percentile of the surface precipitation rate at grid location {(4,4)} computed from the 20C3M outputs of HadCM3 Table 3 for a certain period of interest (e.g. 1950-1969) be x i . To find the value of the first coefficient (b 1 ) corresponding to x i , initially, the two points in the scatter (see Fig 5a) which referred to the values x x and x z of the 25 th percentile of the surface precipitation rate at grid location {(4,4)} closest to x i on either sides of x i are found. Then, by following the solid arrow line shown in Fig 5a, the value of the first coefficient (b 1 ) pertaining to x i is found as y i . If the value of the most influential statistic of the potential predictor computed from the 20C3M outputs of HadCM3 were outside its range computed from NCEP/NACR reanalysis data, then the coefficient/constant determined using the NCEP/NACR reanalysis data was used without any modification. As an example, in Fig 5a, value of the 25 th percentile of the surface precipitation rate at grid location {(4,4)} x j is outside the range of the scatter and there is no known value of coefficient b 1 in the scatter pertaining to any x k (>x j ), hence interpolation is impossible. In such cases, extrapolation of the values of the coefficient/constant corresponding to the given value of the statistic pertaining to the past GCM data is seen as a solution. However, since extrapolation can introduce large errors to the estimated value of the coefficient/constant, it was not practised. Like any other approach such as fitting a linear or a non-linear curve to the scatter between a constant/coefficient and a statistic of a predictor derived from reanalysis data, the linear interpolation between points in the scatter can also introduce uncertainties to the estimation of the values of the constant/coefficients pertaining to the values of the statistic of the predictor derived from the GCM data. In certain instances the uneven scatter (e.g. scatter with higher data density in certain regions than others) between the constant/coefficients and the statistic of predictor derived from reanalysis data was seen (e.g. Fig 5a). In the regions of the scatter where the density of data points is relatively high, the uncertainties introduced to the estimation of the values of the constant/coefficient can arise from the relative spread of the points, and in the regions of the scatter where the density of data points is low, uncertainties can arise due the absence of data points.

Bias correction of GCM outputs and reproduction of observed precipitation
Once the coefficients and constant of the PPRs were determined corresponding to the statistics of the three best potential predictors pertaining to the 20C3M outputs of HadCM3 for the periods 1950-1969, 1970-1989 and 1990-1999, these new PPRs were used to reproduce the observations of precipitation for each calendar month, for each station using 20C3M outputs of HadCM3. Both reanalysis and GCM outputs contain bias, however, in general, GCM outputs tend to contain more bias than reanalysis outputs that are quality controlled and corrected against observations. The bias in the GCM outputs can propagate into the outputs of a downscaling model. Hence, it is an important task to correct the bias in the GCM outputs before any use. In this investigation, using the monthly bias-correction [29] the bias in the average and the standard deviation of 20C3M outputs of HadCM3 was corrected against the corresponding NCEP/NCAR reanalysis outputs for each calendar month.
In the application of the monthly bias-correction it is assumed that the bias in the variables over any period beyond the baseline period will remain the same as that of the baseline period. As the baseline period of the monthly bias-correction 1950-1969 was considered. Over the baseline period the monthly bias-correction was applied to the 20C3M outputs of HadCM3 by replacing their means and the standard deviations in each calendar month with the corresponding means and the standard deviations derived from the NCEP/NCAR reanalysis outputs pertaining to the same period. Beyond the baseline period (i.e. 1970-1989 and 1990-1999), the bias in the 20C3M outputs of HadCM3 were corrected by standardising these outputs with their means and standard deviations pertaining to the baseline period and rescaling with the corresponding means and standard deviations derived from the reanalysis outputs. Then the bias corrected 20C3M outputs of HadCM3 were used with those PPRs modified according to the statistics of 20C3M outputs of HadCM3 for the reproduction of observed precipitation for periods : 1950-1969, 1970-1989 and 1990-1999.

Development of stationary downscaling models (SDM 2 )
For each station, two different types of conventional stationary downscaling models (SDM 2 ) were developed: (1) a downscaling model based on MLR called SDM 2(MLR) and (2) a downscaling model based on GP called SDM 2(GP) . The precipitation observations at the Halls Gap, Birchip and Swan Hill stations and the NCEP/NCAR reanalysis data of the three best potential predictors for the periods 1950-1969 and 1970-1989 were considered for the calibration and validation of SDM 2(MLR) and SDM 2(GP) respectively. For the calibration and validation of SDM 2(MLR) and SDM 2(GP) , the above two periods were selected as it enabled the comparison of performance of SDM 2(MLR) and SDM 2(GP) with that of SDM 1(MLR) which was calibrated, validated and modified for the same periods. SDM 2(MLR) and SDM 2(GP) comprised of a set of 12 PPRs at a given station (each PPR pertaining to a specific calendar month).
In the development of SDM 2(MLR) , for each calendar month, constants and coefficients of MLR-based PPRs were determined for the calibration period by minimising the sum of squared errors between the outputs of the SDM 2(MLR) and observations of precipitation. Thereafter, these PPRs were run with the reanalysis data of the three best potential predictors pertaining to the validation period.
For the development of SDM 2(GP) for each calendar month, for each observation station, the GP [30] technique was employed with the attributes shown in Table 4. The GP algorithm started with the random generation of a pool of downscaling models called an initial population using the observations of precipitation and the reanalysis data pertaining to the three best potential predictors for the calibration period. Then the fitness of each downscaling model in the initial population was assessed. Thereafter, the downscaling models in the initial population were selected for the mating pool based on their fitness. In the mating pool, genetic operations (e.g. crossover) were performed on downscaling models to generate a new population of downscaling models. Then, again the fitness of each downscaling model in the new population was assessed. In the above manner new generations of downscaling models were evolved iteratively until a predefined number of generations was met. Finally, the fittest downscaling model was identified for each calendar month for each station and run with the reanalysis data of potential predictors pertaining to the validation period.
Once SDM 2(MLR) and SDM 2(GP) were calibrated and validated following the above procedure, they were used with the bias-corrected 20C3M outputs of HadCM3 as inputs for the reproduction of observed precipitation over the periods : 1950-1969, 1970-1989 and 1990-1999.

Results
A statistical comparison of the performances of SDM 1(MLR) , SDM 2(MLR) and SDM 2(GP) is presented in Tables 5, 6 and 7 for the three precipitation observation stations. The overall performance of each downscaling model was assessed with normalised mean square error (NMSE) in all three time slices : 1950-1969, 1970-1989 and 1990-1999. It should be noted that in the calculation of the normalised mean square error, the mean square error is normalised with the variance of the observed precipitation at a given station. This enabled the cross comparison of performance of downscaling models developed for stations pertaining to different climate regimes unlike the root mean square error or the mean square error which are sensitive to the order of magnitude of the predictand data. According to Tables 5, 6 and 7 it was seen that, when SDM 1(MLR) for the stations at Halls Gap, Birchip and Swan Hill were run with NCEP/NCAR reanalysis outputs of potential predictors, it displayed the lowest NMSE in all three time slices in comparison to NMSEs of SDM 2(MLR) and SDM 2(GP) . This was due to the fact that unlike SDM 2(MLR) and SDM 2(GP) which were only calibrated for period 1950-1969, SDM 1(MLR) was calibrated using a 20-year moving window moved at a 1-year time step that included the periods: 1950-1969, 1970-1989 and 1990-1999. As shown in Table 5 which refers to the Halls Gap station located in relatively wet climate, when SDM 1(MLR) was run with the 20C3M outputs of HadCM3, it displayed the smallest NMSE over periods 1970-1989 and 1990-1999 in comparison to that of SDM 2(MLR) and SDM 2(GP) . Also, SDM 1(MLR) was able to better reproduce the 50 th percentile of observed precipitation than that by SDM 2(MLR) and SDM 2(GP) in all three time slices. However, SDM 2(MLR) was able to capture the 75 th and 95 th percentiles of observed precipitation better than those by SDM 1(MLR) and SDM 2(GP) , in all three time slices when it was run with 20C3M outputs of HadCM3. Also, in all three time slices, the standard deviation of observed precipitation was better simulated by SDM 2(GP) than that by SDM 1(MLR) and SDM 2(MLR) with the 20C3M outputs of HadCM3.
According to Tables 6 and 7 which refer to the Birchip and Swan Hill stations located in intermediate and dry climate respectively, it was noted that, with the 20C3M outputs of HadCM3, SDM 1(MLR) was able to display the lowest NMSE over the periods 1950-1969 and 1990-1999. However, it was seen that with the 20C3M outputs of HadCM3 no downscaling model was able to outperform the others in better reproducing any of the statistics (e.g. percentiles) of observed precipitation consistently in all three time slices at Birchip or Swan Hill stations. Bold text refer to statistics of precipitation reproduced by downscaling models run with NCEP/NCAR reanalysis outputs which show the lowest bias in comparison to statistics of observed precipitation. Bold italicised text refer to statistics of precipitation reproduced by downscaling models run with 20C3M HadCM3 outputs which show the lowest bias in comparison to statistics of observed precipitation. Downscaling Accounting for Predictor-Predictand Non-Stationarities

Discussion
When run with the 20C3M outputs of HadCM3 (past GCM outputs) it was seen that no downscaling approach (stationary or non-stationary) was able to consistently outperform the other approaches in all three time slices at any of the three stations, in terms of normalised mean square error (NMSE). However, the downscaling models based on the non-stationary approach (SDM 1(MLR) ) were able to display the lowest NMSE at all stations in two of the three time slices with the 20C3M outputs of HadCM3. This indicated that the downscaling models based on the non-stationary approach were able to produce more accurate simulations of observed precipitation more often than linear and non-linear downscaling models based on the conventional stationary approach. For all stationary and non-stationary downscaling models when run with the 20C3M outputs of HadCM3 it was seen that the NMSE was relatively higher for the stations located in the dry and intermediate climate in comparison to that for the station located in the wet climate. This hinted that, irrespective of the downscaling approach, higher degree of error is associated with the simulations produced by the downscaling models pertaining to relatively dry climate.
In the application of the SDM 1(MLR) for future climate, the constants/coefficients of PPRs developed using the reanalysis data of predictors and the observations of precipitation should be updated according to the future climate simulated by the GCM for deducing new non-stationary PPRs pertaining to the future. For that purpose, scatter between the statistics of reanalysis data of predictors and constants/coefficients of PPRs for the past climate and statistics of data of predictors simulated by the GCM for the future are used. In such instances, there is a likelihood that the statistics of predictors derived from GCM data for the future lie outside the range of statistics of predictors derived from the past reanalysis data. This issue can make the model less non-stationary, as extrapolation of a value of a constant/coefficient outside the range of reanalysis data is discouraged. Such likelihood, will be higher in the distant future and smaller in the near future and may make the non-stationary downscaling model developed for the distant future to be more stationary rather than non-stationary. However, with time, the continuous updating of the scatter between the statistics of the predictors derived from reanalysis data which is expanding over time and the constants/coefficients in the PPRs, will minimize of likelihood of having the statistics of the predictors derived from GCM data for future to lie outside the range of the reanalysis data. Alternatively, a non-linear regression technique such as Genetic Programming can be used to develop a relationship between the values of a statistic of a predictor derived from the reanalysis data and the values of a coefficient/constant in the PPR of interest. Then this non-linear regression relationship can be used with the values of the statistic of the predictor derived from a GCM database which lie outside the range of past reanalysis data to determine the corresponding values of the coefficient/constant in the PPR. This regression based approach can even be used in conjunction with the already proposed continuous updating of the scatter between the statistics of the predictors derived from reanalysis data and the constants/coefficients in the PPRs.