Hydrological effects of the temporal variability of the multiscaling of snowfall on the Canadian prairies

Daily historical snowfall data were analysed with the objective of determining the stability of their variability at short temporal scales. The data are weakly multifractal over periods shorter than one month, which controls their scaling properties and which can be used to statistically downscale monthly data to shorter-duration values. Although the daily snowfall values appear to be stationary, their multifractality displays much temporal variability, with most sites showing statistically-significant trends. Through use of a physically-based hydrological model, it is demonstrated that the variability of the multiscaling of snowfall can affect the timing and quantity of snow accumulation in catchments where the snowpacks are subject to wind redistribution. Therefore trends in scaling, based on multifractal characteristics, should be taken into account when downscaling climate model scenario outputs.


Introduction
The Canadian prairies are a cold, semi-arid region prone to prolonged droughts.Consequently, the potential effects of climate change on water resources are of great interest.Prairie streamflow is dominated by snowmelt runoff and the source snowfalls are subject to wind redistribution by blowing snow which transforms them into highly variable, redistributed snowcovers (Shook and Gray, 1997) and subjects them to substantial sublimation loss (Pomeroy et al., 1993;Pomeroy and Gray, 1995;Fang and Pomeroy, 2008).Predicting future streamflows therefore requires prediction of future snow accumulation and redistribution.Although physically-based hydrological models that incorporate blowing snow processes are required to understand the effects of Correspondence to: K. R. Shook (kevin.shook@usask.ca)climate change on prairie hydrology, these models generally need input data such as precipitation, wind speed, air temperature and humidity at hourly temporal resolutions which are coarser than can be simulated by climate models (Wilby and Wigley, 1997) and some form of temporal downscaling is needed.
The distribution of precipitation in space-time has long been known to display scaling, where the statistical properties of the distribution can be related over some (large) range of scales.For example, a simple scaling law can relate world-wide extreme precipitation to the time scale over ranges spanning from minutes to years (Galmarini et al., 2005).Data sets displaying scaling may be either monofractal or multifractal in their distribution, depending on whether or not their scaling behaviour can be characterised by a single fractal dimension (Lovejoy and Schertzer, 2006).
Precipitation can be modelled as being multifractal in both time and space (Lovejoy and Schertzer, 2006;Olsson, 1995).Molnar and Burlando (2008) found differences in the multiscaling of winter and summer precipitation in Switzerland, with the summer precipitation showing strong multifractality and the winter precipitation being close to being monofractal.By statistically reproducing the observed multifractality, precipitation has been statistically downscaled to shorter time scales, while conserving observed scaling behaviours (Gaume et al., 2007;Seuront et al., 1999).
Statistical downscaling methods require the existence of stationarity in the relationships between predictor and predictand variables (Fowler et al., 2007).In the case of multifractal downscaling, the requirement is for stationarity in the degree of multifractality over time.Therefore, the objectives of this research are to determine: -the existence of multiscaling in the temporal distributions of prairie snowfall, -the existence of temporal trends in the multiscaling which will affect downscaling of prairie snowfall, and -the effect(s) of any temporal trends in snowfall multiscaling on prairie hydrology caused by changes in the temporal frequency distribution of prairie snow accumulation.

Selection of data
Measurements of snowfall are highly influenced by wind induced undercatch Goodison (1978).Of equal importance for these analyses, the introduction of Nipher wind shields in Canadian use during the winter of 1960-1961 is known to have introduced a discontinuity in snowfall data (Yang et al., 1999).To remove the influence of wind speed and shielding on the snowfall data, rehabilitated snowfall data were obtained from the Historical Adjusted Climate Database for Canada.These values have been adjusted for the effects of wind as described by (Mekis and Hogg, 1999).
Daily snowfall data for Calgary, Medicine Hat, Saskatoon, Regina, Indian Head and Brandon were selected because they were widely dispersed throughout the prairies, had data preceding the 20th century, and the data collection sites were not moved during the period of record.The analysed records span the years 1895-2003, with the exception of the Regina data which begin in 1898.The locations of the six sites on the Canadian prairies are shown in Fig. 1.
All analyses were performed using the open source statistical program R, which is described by Ihaka and Gentleman (1996).The program may be downloaded at www.r-project.org.

Scaling analyses
The presence of scaling in a time series is shown by the power spectrum of the data, which will scale according to the relationship (Davis et al., 1994) 1), for wavelengths greater than 0.032. where: The value of β indicates the behaviour of the underlying time series.If β=0, then the time series is entirely random.If β <1, then the time series is stationary.If β >1 then the time series is non-stationary (Davis et al., 1994).
For periods shorter than one month, i.e. for wavelengths greater than 1/31=0.032,scaling behaviour (β >0) was observed for all of the data sets.As shown in Fig. 2, the power spectrum for the Brandon daily snowfall displays scaling for time scales shorter than one month; at longer time scales, the scaling behaviour disappears.All of the snowfall datasets had values of β <1 (mean=0.24,s.d.=0.06), indicating that the data are stationary.Although the mean of the snowfall frequency distribution may be stationary, it does not necessarily follow that the scaling behaviour of snowfall is constant over time.
The scatter about the best-fit line visible in Fig. 2 has more than one cause.Figure 3 shows that the magnitude of β for a given snowfall data set tends to be unstable over timescales on the order of decades, undergoing frequent statisticallysignificant shifts which indicate temporal variability in the scaling behaviour.In this analysis, and in subsequent analyses, the entire time series was divided into ten sections of equal length, and statistical properties were determined for each section.
At smaller scales (i.e.within each section) the scatter appears to be fairly consistent over time, as shown by the 5% and 95% confidence levels plotted in Fig. 3  temporal variation with respect to the value of β.The difference between the value of β and its 5% or 95% confidence levels varies by less than 10% over all ten sections of this dataset.
From a practical point of view, the variation in scaling within each section, which appears to be random, is of less concern than the potential existence of a trend.As is discussed below, datasets showing scaling can be downscaled stochastically, and random variability can be accomplished by adding random variation to the downscaling parameters.The presence of a trend in the scaling parameters would be more problematic as it would imply that scaling parameters used in the past may not be usable in the future.
Therefore, it is important to quantify the type and magnitude of the scaling present within snowfall data sets to determine the parameters required for downscaling snowfalls over time scales shorter than one month.Equally, it is important to determine the presence, if any, of trends in these parameters.

Multiscaling analysis
The type and magnitude of scaling present in a data set can be quantified in several ways; one of the simplest is through the moment scaling function K(q) which is related to the normalised values of the data set, ε(t), by (Olsson, 1995) (2) where: q = moment, and λ = scale factor (ratio of scale of interest to size of entire data set).
The symbol refers to ensemble averaging the results of the calculations over all possible values of ε(t).If the data are monofractal, then plots of K(q) vs. q will describe a straight line.If the data are multifractal, then plots of K(q) vs. q will exhibit curvature, which can be quantified by the universal multifractal parameters (constants) α and C1 as described by (Seuront et al., 1999) The parameters α and C1 are indices of the inhomogeneity of data.Small values of α, and large values of C1, imply a tendency to large fluctuations in value within a data set (Finn et al., 2001).The minimum possible value of C1 is zero, which corresponds to a dataset whose values are completely uniform.The maximum value of C1 is equal to the dimension of the observation space, which in the case of a time series is 1.The value of α is always between 0 and 2 (Seuront et al., 1999).

Scaling of extremes
The probability density functions of the snowfall data sets show extreme positive skewness, because most of the daily snowfall data values are zero.The upper tails of these distributions are very "fat", compared to Gaussian distributions, the largest values fitting the scaling relationship defined by Olsson (1995) as where: X = a given value, x = a threshold value, and q D = scaling exponent.
The scaling exponent q D has been identified as being the limiting moment beyond which moments will diverge from scaling relationships (Olsson, 1995).In practice, application of Eq. ( 4) is made more difficult by the fact that logarithmic plots of Pr(X > x) vs. snowfall typically describe a curve, rather than having the abrupt "sill" described by Olsson (1995), which requires the use of a threshold or restricted range for computing q D .
Using the 100 largest values from each data set, the mean value of q D for the snowfall datasets is 3.57 (s.d.=0.42).To ensure that Eq. ( 3) can be applied, the maximum value of q used for all calculations of K(q) in this research was set to 3.0.The values of K(q) determined from the snowfall datasets typically agree very well with the values derived from the universal multiscaling relationship for values of q between 1 and 3.
Values of C1 and α for each data set were determined using the Double Trace Moment (DTM) method, as described by Olsson (1995), for time periods up to 32 days, and for moments between 1.1 and 3 in increments of 0.1.The magnitudes of α (mean=0.62,s.d = 0.07) and C1 (mean=0.46,s.d.=0.03) computed for the snowfall datasets were similar to those found by Olsson (1995) of 0.626 and 0.44, respectively, for precipitation over time scales of 8 min to 11 days.

Multifractal stability
To test the temporal stability of the universal multifractal parameters, α and C1 were calculated for sub-sections, each consisting of one-tenth of each snowfall time series.As shown in Fig. 4, the plots of K(q) vs. q, and consequently values of α and C1, show frequent and significant changes throughout each of the datasets.Interestingly, the plots of K(q) vs. q are only slightly curved, indicating that the snowfalls on the Canadian prairies, like Molnar and Burlando's (2008) observed Swiss snowfalls, are weakly multifractal.
Although the weakness of the observed multifractality implies that monofractal methods might be used in practice for downscaling, this study uses multifractal methods which are more general and which allow the use of a simple downscaling method described below.
Maximum and minimum values for α and C1 for each data set are listed in Table 1, together with the results of Mann-Kendall analyses of the parameters for trends at the 5% level.Many of the datasets show the presence of significant trends in their α values.Although none of the C1 parameters showed significant trends, and their magnitudes showed less variability than did the values of α, the values of the parameters are inversely related, as shown in Fig. 5.

Snowfalls and snow accumulation
The hydrography of the Canadian prairies is very unusual because much of the land drains internally to small wetlands, rather than contributing to the flow of major rivers (Godwin and Martin, 1975).As described above, the hydrology of the region is also unusual as almost all surface runoff events are due to the spring melt of winter snowpacks, rather than to runoff from rainfall.Therefore, it is difficult to link the multifractality of rainfall to that of streamflow records in this region, as was done in France by Tessier et al. (1996).
As snowpacks accumulate over an entire season, they would seem to be insensitive to the distribution of individual snowfall events.However, blowing snow events transport snow horizontally among landscape units and sublimation of blowing snow removes snow from the landscape.The moisture fluxes are highly dependent on wind speed, sublimation varying as the 5th power, and horizontal transport varying as the 4th power (Pomeroy and Gray, 1995;Essery, 2001).The fluxes also depend on the state of the snowpack, old snow being more resistant to erosion than new snow (Li and Pomeroy, 1997).The difference between the probability Hydrol.Earth Syst.Sci., 14,[1195][1196][1197][1198][1199][1200][1201][1202][1203]2010 www.hydrol-earth-syst-sci.net/14/1195/2010/  .The horizontal distance between the two sites is approximately 3 km, and the difference in elevation is less than 200 m.
density function (PDF) of cumulative snowfall and that of the snow on the ground can be quite dramatic, as shown in Fig. 6a and 6b, which plots historical snow data surveyed by Alberta Environment.The site in Fig. 6a is near Banff, Alberta in a sheltered mountain valley and in this case wind redistribution is muted, but a shift is still evident possibly due to land use effects, though the PDF shape is conserved.The site in Fig. 6b is southeast of Edmonton, Alberta in the prairies where wind redistribution of snow can be significant and shows ablation of snowfall as well as a change in the shape of the PDF.It is unlikely that the snowfall ablation is strongly influenced by melt as significant melt events prior to 1 March are uncommon on the Canadian prairies; the change in the shape of the PDF is consistent with removal and redistribution of snow by wind.It is anticipated that the temporal and frequency distributions of snowfall events will significantly affect the accumulation of snow which in turn will affect the distribution of runoff events.The influence of snowfall disaggregation on snowcover generation was tested using the Cold Regions Hydro-  .The horizontal distance between the two sites is approximately 15 km, and the difference in elevation is less than 60 m.
logical Model (CRHM), a physically-based hydrological model developed at the Centre for Hydrology, University of Saskatchewan (Pomeroy et al., 2007).CRHM is able to accurately simulate the accumulation and sublimation of blowing snow of prairie landscapes and, being physically-based is well-suited to simulation of changed conditions (Fang and Pomeroy, 2008).The surface routines of CRHM use physical parameters, such a crop height, which can be measured rather than being calibrated.Although the parameters of CRHM's subsurface flow modules do require calibration, these modules were not used in this study.

Downscaling snowfall data
Simulated snowfall data were created using the random multiplicative cascade method (Gaume et al., 2007), where the monthly snowfall is repeatedly subdivided over several generations.In the first generation, the snowfall is divided into two portions, according to www.hydrol-earth-syst-sci.net/14/1195/2010/ Hydrol.Earth Syst.Sci., 14, 1195-1203, 2010 where: S 0 = monthly snowfall, S = snowfall assigned to first piece in the first generation, S 1,2 = snowfall assigned to second piece in the first generation, and η=random number in the range 0-1.The process is then repeated iteratively, doubling the number of pieces with each generation.As the number of pieces increases the data set becomes more heterogeneous, although the total value of all of the pieces does not change.The process is shown schematically for four generations in Fig. 7.
To produce daily values, 5 generations are used creating 32 pieces, which were reduced to the required daily values for each month by   S 5 = the 5th generation, and n = number of days in given month.By manipulating the distribution of η, data series having varying values of α and C1 can be generated.The values of η were drawn from a lognormal distribution whose parameters were manipulated by trial-and-error to produce the desired values of α and C1.The lognormal distribution was one of several evaluated by Gaume et al. (2007).Further research will be required to determine the optimal probability density function of η for snowfall disaggregation over the Canadian prairies.
As shown in Fig. 8, increasing the magnitude of α while decreasing the magnitude of C1 results in snowfall distributions which tend to produce more frequent, smaller events which are distributed more evenly in time.In both cases the total snowfall is 19.5 mm.Snowfall events having large values of α will result in snowpacks whose surfaces are relatively new and therefore more subject to blowing, at any given time.

Simulation
The basin simulated by the model is of a typical small prairie stream, consisting of upland fallow and stubble fields which drain runoff into a grassed intermittent stream channel.It is modelled on Creighton Tributary at Bad Lake, Saskatchewan, and all model parameters were taken from this location.The ability of CRHM to produce accurate simulations of snow accumulations at this location was demonstrated by Pomeroy et al. (2007).Because of the shortness of the recorded data set (less than 20 years), meteorological data from Saskatoon were used in this study.Although Saskatoon is approximately 150 km northeast of Creighton Tributary the intention is to simulate the effects of downscaling, rather than to accurately model the hydrology of a particular location.
CRHM uses Hydrological Response Units (HRU) to simulate the behaviour of each of the basins's sub-regions.Within each HRU all state variables, forcings, and parameters are assumed to be spatially uniform.An important feature of CRHM is that HRUs may be linked in a variety of ways.Although runoff drains to the stream channel, the action of blowing will transport snow from the fallow field to the stubble field, and from both the stubble and fallow fields to the grassed channel, based on their relative surface roughnesses.Traditional hydrological models, where connectivity is based on surface water drainage, cannot properly capture the redistribution of wind-blown snow.
Although CRHM uses daily snowfall values, the model is run on an hourly time step and requires hourly values of air temperature, relative humidity, wind speed, rainfall, and incoming solar radiation.Each simulation was run for the period 1960-2006, for which good values of the hourly forcing variables other than solar radiation were available.Daily incoming shortwave solar radiation was estimated from air temperature using the method of Annandale et al. (2001), and was downscaled to hourly values.
Because the snowfall downscaling process is stochastic, the set of downscaled data will depend on the particular set of random values of η generated from the specified distribution.To reduce the effects of the finite lengths of the simulations, fifteen realizations were run for each disaggregation using a single set of lognormal parameters to specify the frequency distribution of η.The downscaled datasets had fairly consistent values of α and C1 over each set of realizations, their mean coefficients of variation being 0.07 and 0.03, respectively.As shown in Fig. 9, the ranges of α and C1 values, and their relationship to each other, are similar for the synthetic data and for the observed Saskatoon daily snowfall values.

Results
The CRHM runs produced snow accumulations for each HRU for which mean daily snow water equivalent (SWE) values were computed over each set of realizations.As shown in Fig. 10, the variation in α and C1 caused variation in the timing and magnitudes of snow accumulation and melt for each HRU.The accumulation on the stubble HRU was least affected, probably because it both receives and transmits blowing snow.As a control, a CRHM run was executed using the monthly snowfall totals distributed evenly over all days in each month.
The effect of variation in the multifractal parameters on the mean peak accumulation of SWE is shown in Fig. 11.The plots show that increased values of α in the snowfall time series result in linearly increased values of mean peak SWE on the grass HRU, decreased values on the fallow field and largely unchanged values on the stubble HRU.These results are consistent with the assumption of increased blowing snow events resulting from increased values of α, which correspond to increased uniformity of snowfall.
When the control uniform snowfalls were used in the CRHM simulation, the magnitudes of the resulting mean peak SWE for the fallow, stubble and grass HRUs were 48 mm, 59 mm, and 92 mm, respectively.These values are consistent with the trends in modelled mean peak SWE with increasing uniformity in snowfall noted above.As would be expected, the control grass HRU peak SWE is greater than, the fallow HRU peak is smaller than, and the stubble HRU peak is approximately the same as, the values simulated using the maximum value of α.

Summary and conclusions
Analyses of daily snowfall on the Canadian prairies determined that the datasets exhibit multifractality at temporal scales shorter than one month, and that the multiscaling shows evidence of trends for decreasing intermittency of snowfall over the last century at the majority of sites tested.The transformation of snowfall into snow on the ground via redistribution and ablation processes means that frequency distributions of snowfall cannot be used directly to create distributions of snow on the ground.Blowing snow and snow ablation simulations driven by stochastic or measured snowfall estimates can estimate snow on the ground using physically-based calculations.Such simulations demonstrate that the timing and quantity of the spring snowmelt is affected by the observed multiscaling of daily snowfalls, probably because the effects of blowing snow on snowfall accumulation are strongly influenced by the timing of snowfall.Because the majority of annual runoff in this region is typically the result of spring snowmelt, and because climate models require downscaling to produce daily snowfalls, the variability in the multiscaling of prairie snowfalls represents an additional source of uncertainty in determining the hydrological effects of climate model outputs.This can be compensated for in historical analyses but at the present time there is no method to identify likely scaling trends for future scenarios.

Fig. 1 .
Fig. 1.Location of the study sites in Western Canada.The prairie ecoregion in Canada is outlined.

Fig. 4b .Fig. 5 .
Fig. 4b.Variation in of α and C1 by fraction of the Medicine Hat snowfall.Each section corresponds to one-tenth of the entire time series.The vertical bars represent the maximum and minimum α and C1 values computed over each interval.Temporal variability of multifractality of Medicine Hat daily snowfall data, 1895-2003.Each section corresponds to one-tenth of the entire time series.

Fig. 6b .
Fig. 6b.Density plot of 1 March cumulative seasonal snowfall at Edmonton, AB and snow on ground at Bellis snow course(1974- 2009).The horizontal distance between the two sites is approximately 15 km, and the difference in elevation is less than 60 m.

Fig. 7 .Fig. 8 .
Fig.7.Schematic illustration of four generations of the multiplicative cascade method.In each generation, the total value of all pieces is 1.

Fig. 9 .
Fig. 9. Mean values of α and C1 computed for all 10 section of the Saskatoon daily snowfall time series and from simulated daily data downscaled from measured monthly snowfalls.

Fig. 10 :Fig. 11 :
Fig. 10: Mean CRHM simulated snowcover accumulations for three HRUs, based on measured monthly snowfall downscaled to daily values.The curves represent accumulations from downscaled data having α and C1 values of 0.36 and 0.64, and 0.72 and 0.41, respectively, which are similar to the extremes of the values measured in sub-sets of historical data.

Fig. 10 .
Fig. 10.Mean CRHM simulated snowcover accumulations for three HRUs, based on measured monthly snowfall downscaled to daily values.The curves represent accumulations from downscaled data having α and C1 values of 0.36 and 0.64, and 0.72 and 0.41, respectively, which are similar to the extremes of the values measured in sub-sets of historical data.

Table 1 .
Temporal variability in values of α and C1, and the results of Mann-Kendall trend analyses.
Fig. 6a.Density plot of 1 April cumulative seasonal snowfall at Banff, AB and snow on ground at Bow River snow course