## Abstract

As a consequence of the increase in atmospheric greenhouse gas concentrations, potential changes in both precipitation occurrence and intensity may lead to several consequences for Earth’s environment. It is therefore relevant to estimate these changes in order to anticipate their consequences. Many studies have been published on precipitation changes based on climate simulations. These studies are almost always based on time slices; precipitation changes are estimated by comparing two 30-yr windows. To this extent, it is commonly assumed that nonstationary processes are not significant for such a 30-yr slice. Thus, it frees the investigator to statistically model nonstationary processes. However, using transient runs instead of time slices surely leads to more accurate analysis since more data are taken into account. Therefore, the aim of the present study was to develop a transient probabilistic model for describing simulated daily precipitation from the Canadian Regional Climate Model (CRCM) in order to investigate precipitation evolution over North America. Changes to both the occurrence and intensity of precipitation are then assessed from a continuous time period. Extreme values are also investigated with the transient run; a new methodology using the models for precipitation occurrence and intensity was developed for achieving nonstationary frequency analysis. The results herein show an increase in both precipitation occurrence and intensity for most parts of Canada while a decrease is expected over Mexico. For the continental United States, a decrease in both occurrence and intensity is expected in summer but an increase is expected in winter.

## 1. Introduction

Precipitation plays a central role for ecosystems, human infrastructures, and freshwater supplies, among many other areas. Much infrastructure has been built to resist up to a particular water level. With climate change, precipitation-generating processes might evolve in time, leading to changes in the occurrence and intensity of precipitation. The potential impact on ecosystems and infrastructures is a major concern. To address this issue, temporal precipitation evolution as influenced by climate change must be described. Even though the general rainfall-generating mechanisms are fairly well known, the evolution of rainfall in a nonstationary climate is intractable. Indeed, the nonlinear and chaotic nature of the equations describing rainfall mechanisms prevents simple predictions for precipitation. Therefore, numerical simulations from climate models may be the best way to provide quantitative prediction of precipitations in a nonstationary climate (IPCC 2007). Regional climate models (RCMs) are possibly better tools for studying climate change impact on precipitation than global climate models (GCMs). The resolution of RCMs (from 50 to 10 km) enables a good representation of some important features in the formation of precipitation (mainly orography, surface processes, and hydrology). The spatial-scale characteristics of precipitation are reasonably represented in RCMs on a length scale of roughly 100 km (Rasmussen et al. 2012).

Most climate studies are based on the comparison of two standard time windows. In the IPCC Fifth Assessment Report (AR5; Christensen et al. 2013), GCM-projected climate changes are assessed for two 20-yr time slices: 1986–2005 for the present period and 2081–2100 for the future climate period. For RCM-projected changes, it is often two 30-yr time slices that are investigated, such as 1961–90 and 2071–2100. This is the case for the PRUDENCE experiment (Christensen and Christensen 2007; http://prudence.dmi.dk/main.html), the ENSEMBLES project (van der Linden and Mitchell 2009; http://www.ensembles-eu.org) over Europe, and NARCCAP (Mearns et al. 2009; http://www.narccap.ucar.edu) over North America. (Expansions of acronyms are available at http://www.ametsoc.org/PubsAcronymList.) This approach was used, for a certain time, because regional climate models were too expensive to run transient simulation. Even though transient runs are now available, time-sliced analyses are still largely performed, notably to avoid the modeling of nonstationary processes. Indeed, it is generally assumed that a 30-yr window is long enough for climate to dominate natural variability and short enough for nonstationary processes to be negligible. But from a statistical point of view, a large amount of information is lost in the omitted intermediate years, an unfortunate scenario given the high cost of such simulations. The intermediate years should be used in order to describe the nonstationary processes for the whole period 1961–2100. This may help to distinguish climate change from natural variability (Kay and Jones 2012). To our knowledge, only Goderniaux et al. (2011) have modeled daily RCM precipitation on the continuous time period 2010–85 on a probabilistic basis. Their probabilistic model was then used to supply a hydrological model with several proper precipitation scenarios.

In studies of extreme values, where only few observations per year are available, the need to consider transient simulations in order to increase the number of observations is crucial. For this reason, many recent papers used a transient run to adjust a time-varying generalized extreme value (GEV) distribution to describe a series of annual maxima (see, e.g., Kharin and Zwiers 2005; Hanel and Buishand 2011; Kay and Jones 2012) or a time-varying peaks-over-threshold (POT) model to describe a series of threshold exceedances (see, e.g., Sugahara et al. 2009; Mailhot et al. 2013). Such percentile estimation corresponding to a given return level is still challenging, even if the entire time period is used, because the nonstationary characteristics have to be assessed with only a few annual observations. For example, three nonstationary GEV parameters have to be estimated with only 140 annual maxima, those from the period 1961–2100.

On the other hand, regional frequency analysis stands as an alternative to assess extreme precipitation using 30-yr time slices. This approach compensates the lack of observational records by incorporating regional information given a set of homogeneous grid points. It has been applied for precipitation extremes (see, e.g., Fowler and Ekström 2009; Mladjic et al. 2011; Mailhot et al. 2012; Roth et al. 2014). However, finding a homogeneous set of points is the most challenging part of this approach that has not yet been completely solved for North America.

This paper focuses on the assessment of transient changes to precipitation occurrence, intensity, and 100-yr return level over North America for the continuous time period 1961–2100, using a climate simulation from the Canadian RCM (CRCM; Caya and Laprise 1999). Several papers have been published on CRCM precipitation changes between two subperiods (see, e.g., Mladjic et al. 2011; Mailhot et al. 2012; Monette et al. 2012; Paquin et al. 2014), but none on transient change. At this stage, only one simulation is investigated. The proposed methodology aims to extract more reliable information from a single simulation. Eventually, the methodology should be extended for an ensemble of simulations in order to give better estimates of climate change. The first step of the methodology is to fit several parametric models to CRCM simulated precipitation occurrence and intensity. Thereafter, only the best model, according to the Bayesian information criterion (BIC), will be kept. Finally, an approach to estimate nonstationary return levels will be described.

## 2. Data from the Canadian Regional Climate Model

### a. CRCM simulations description

The regional simulation under investigation was performed by the Canadian Regional Climate Model (CRCM 4.2.3, based on CRCM 4.2; Music and Caya 2007; de Elía and Côté 2010). It was run over the North American domain on a polar stereographic grid with a horizontal resolution of 45 km true at 60°N (totaling 200 × 192 grid points), using 29 vertical levels, a time step of 15 min, and an archival period of 6 h. The free domain, exclusive of sponge zone, is composed of 182 × 174 grid points. The regional model was driven by the outputs from the Third Generation of the Canadian Coupled Global Climate Model (CGCM3; Scinocca et al. 2008), which performed a transient simulation following the IPCC “observed 20th century” scenario (20C3M) for years 1961–2000 and SRES A2 (Nakicenovic et al. 2000) scenario for years 2001–2100. The simulation is driven by the fourth member from an ensemble of five twentieth-century and SRES A2 simulations.

CRCM convective precipitation is parameterized following Bechtold–Kain–Fritsch (Bechtold et al. 2001) while large-scale (stratiform) precipitation is parameterized as a simple supersaturation-based condensation scheme. CRCM’s sensitivity of precipitation simulations to the choices of parameterization scheme can be found in Paquin et al. (2002), Jiao and Caya (2006), and Christensen et al. (2008). CRCM, like many other RCMs, overestimates the occurrence of small precipitation amounts (Christensen et al. 2008). To avoid a large artificial occurrence probability, we assumed that under a threshold of 1 mm no precipitation occurs. The 1-mm threshold follows the recommendation of the Statistical and Regional Dynamical Downscaling of Extremes for European Regions (STARDEX; http://www.cru.uea.ac.uk/projects/stardex). We also reduced the corresponding amount of 1 mm when it occurs. Removing 1 mm frees us from the need to work with truncated statistical distribution for modeling daily intensity; otherwise, the distribution supports [1, ∞). If small-to-moderate precipitation amounts are important for an investigator, it would be possible to add 1 mm when postprocessing the results. Nevertheless, such a reduction does not influence the results for large precipitation events.

The CRCM simulation domain covers North America, as shown in Fig. 1. We restrict the study to the 12 570 grid points that represent land points. The sea grid points are not investigated.

### b. Notations

In this paper, the precipitation description is performed by seasons. We define the winter season as December, January, and February; the spring season as March, April, and May; the summer season as June, July, and August; and the autumn season as September, October, and November.

Let us introduce the notation used throughout the paper. For any given grid point and season, let *Y*_{kl} be the amount of daily precipitations in millimeters (after the bias reduction of 1 mm) on the *l*th day of year index *k*, where

The index *k* indicates the number of years elapsed since 1961. Thus, *k* = 0 stands for the year 1961 while *k* = 139 indicates the last year (2100). The index *l* stands for the day in a particular season of length *L*. For instance, assuming that we are interested in the summer season, the index *l* goes from 1 to 92; *l* = 1 stands for 1 June and *l* = 92 stands for 31 August.

To describe precipitation occurrence, we define a dichotomous random variable *X*_{kl} that indicates if precipitation occurs on the *l*th day of year *k*:

Thus, the number of days during which precipitation occurred in a given season and a given year *k* is

## 3. Methodology

Cooley and Sain (2010) claimed that “climate is what you expect and weather is what you actually get.” In this paper, we tried to describe climate with one of its possible realization: an RCM simulation. To infer climate from weather, parametric statistical distributions are fitted to precipitation data. Moreover, statistical models are needed for establishing the significance of the nonstationarity. Without such models, results can only be descriptive. Also, the model’s parameters are the best way to extract and summarize the information and to quantify the sampling uncertainty.

Precipitation is subject to seasonal fluctuations that are much larger than climate variations. Therefore, it is a better strategy to partition the calendar year into seasons and then study each part separately. Within a season, we assumed that both daily precipitation occurrence and intensity are independent and identically distributed. This assumption seems more suitable for convective than stratiform precipitation since the latter might be autocorrelated over days. Nevertheless, we made this assumption first to simplify the parametric models. Actually, this approach frees us from the needs to statistically model seasonal fluctuation, which could be very difficult and lead to large uncertainties.

In this paper, we developed several statistical models for describing daily precipitation occurrence, intensity and extremes. At every grid point, all these models have been estimated but only the best model, according to the BIC, was kept. The next sections are devoted to the description of the considered statistical distribution for modeling precipitation occurrence, intensity, and extremes.

### a. Parametric models for occurrence

For any given day *l* of year *k*, either precipitation occurs or it does not. The natural probability model for this dichotomous variable is the Bernoulli distribution. Let *p*_{k} be the probability that precipitation occurs a given day of year *k*, then

Assume that this probability does not evolve over time; that is, *p*_{k} = *p*, ∀*k* ∈ {0, 1, …, 139}. This leads to the first parametric model for precipitation occurrence:

Suppose now that this probability could evolve with time, this gives the second parametric model:

where

is the logit link function widely used in logistic regression. Note that we always have 0 < *p*_{k} < 1, with *p*_{k} close to 1 if *ζ*_{0} + *ζ*_{1}*k* is a large positive number and close to 0 if *ζ*_{0} + *ζ*_{1}*k* is a large negative number. Note that the first model is stationary while the second is nonstationary over years.

### b. Parametric models for intensity

For precipitation intensity, we only modeled strictly positive precipitation: the zero values are ignored. For a given year index *k* and conditional on days when precipitation occurred, eight different parametric densities were considered to model *Y*_{kl} | *X*_{kl} = 1:

where xp(*y*|*θ*) denotes the exponential density of parameter *θ* evaluated at *y*, amma(*y*|*λ*, *θ*) denotes the gamma density of parameters (*λ*, *θ*) evaluated at *y*,

for

Model _{3}, like models _{6}, _{7}, and _{8}, is a mixture model and more specifically a mixture of exponentials distributions. Such a model conditional on *X*_{kl} = 1 can be interpreted as follows: a particular realization *y*_{kl} of the random variable *Y*_{kl} has a probability (1 − *w*) of being a realization of the first mixture component Exp(*θ*_{0}) (where Exp denotes an exponential distribution) and a probability *w* of being a realization of the second mixture component Exp(*θ*_{1}). See Titterington et al. (1985) for more details.

The temporal functions in Eq. (8) have been chosen in order to insure that *θ*_{0}(*k*) and *θ*_{1}(*k*) are strictly positive for any value of *k*. Thus, the different models based on exponential and gamma distributions are always definite, no matter the value of *k*. Other functions like *θ*_{0}(*k*) = *α*_{0} + *α*_{1}*k* may be tricky to use in order to have *θ*_{0}(*k*) > 0 for *k* ∈ (0, ∞) because parameters constraints would be *α*_{0} > 0 and *α*_{1} ∈ ℝ. An exponential transformation could be appropriate though.

Tables 1 and 2 summarize the expectations and the variances as a function of the year index *k* respectively for the occurrence models and the intensity models. The expectations and variances of the occurrence models are *L* × *p* and *L* × *p* × (1 − *p*) for model _{1} and *L* × *p*_{k} and *L* × *p*_{k} × (1 − *p*_{k}) for model _{2}.

### c. Model estimation

Every model has been adjusted using maximum likelihood estimation. The maximum likelihood estimation of the parameters in models _{1}, _{1}, _{2}, _{4}, and _{5} is quite straightforward. For model _{2}, we used the well-known algorithms for estimating parameters of a logistic regression. For the mixture models _{3}, _{6}, _{7}, and _{8}, we had to use the expectation–maximization (EM) algorithm (Dempster et al. 1977) to maximize the likelihood function.

### d. Statistical model selection

For every grid points, two models were fitted for precipitation occurrence and eight for precipitation intensity. Then, we used the BIC to select the best model among those fit. The BIC is defined as

where is the likelihood of the parametric model, *κ* is the number of parameters of the model, and *n* is the number of observations. The best model, for a given grid point, is the one with the smallest BIC. The BIC has a penalty term for the number of parameters, which allows us to select a better model in the trade-off between adjustment and number of parameters.

### e. Return level

Extreme value theory (EVT) is generally considered to be the best tool for estimating return levels from a time series (Coles 2001). One of the central results of EVT tells us that, under appropriate conditions, the distribution of the maximum of a large number of random variables can be approximated by the GEV distribution. For our application, the seasonal maximum is extracted from a series of 90 observations. In general, it is unclear whether a sequence of 90 observations is sufficient to invoke the GEV distribution. However, in the context of precipitation maxima, a period of 90 days is usually assumed to be suitable (see, e.g., Fowler and Ekström 2009; Katz 2013).

Seasonal precipitation series are subject to nonstationarity due to climate changes. The standard approach for addressing nonstationarity in maxima series is to model the GEV parameters as a function of time (e.g., Coles 2001; Katz 2013). Examples of applications can be found for instance in Katz (2010), Gilleland and Katz (2011), and Cheng et al. (2014). A common pragmatic approach is to consider the shape parameter invariant in time because a large uncertainty is associated with its estimation even in the stationary case. Let be the seasonal maximum for the year *k* at a particular grid point. The usual approach consists in modeling as follows:

for different link functions *σ*_{k} and *μ*_{k}. For example, one could model those as

and estimate the parameter vector (*ξ*, *β*_{1}, *β*_{2}, *β*_{3}, *β*_{4}) with the vector of maxima . In most cases, these parameters have to be estimated with only a few maxima, leading to a large estimation variance.

To avoid estimating temporal trends with only a few maximum data, the information extracted from the daily series could be used instead. The idea is to model the entire precipitation distribution evolving over the whole period and then transform the data into a stationary series. This could be achieved with the intensity models described in section 3b. Let *ν*_{k} and be respectively the annual conditional mean (expectation) and conditional variance of the chosen model for a certain grid point:

We define *Z*_{kl}, the standardized amount of precipitation on day *l* of year *k*, as

Both the expectation and the standard deviation in the last equation can evolve with time depending on the model chosen for intensity. The standardized random variable *Z*_{kl} | *X*_{kl} = 1 has the following properties: 𝔼(*Z*_{kl} | *X*_{kl} = 1) = 0 and Var(*Z*_{kl} | *X*_{kl} = 1) = 1 for all *k* and *l*. The standardized series is thus second-order stationary while its two first moments are time invariant. However, let us assume that the transformed series is strictly stationary. We will see further in the present section that this assumption is equivalent to consider a constant shape parameter. The seasonal maxima series of the transformed series is defined as

where *L* denotes the set of days in a given season. Since we have assumed that after transformation the corresponding daily precipitation is stationary, the series of maxima may be considered as independent and identically GEV distributed:

A classical frequency analysis can thus be performed on the transformed series. Every theoretical result and definition, such as the return period, is directly applicable to the transformed series. Once the GEV is fitted, the quantile *z*_{(1−1/T)} corresponding to a return period of *T* years can be estimated. The quantile *z*_{(1−1/T)} is such that

The corresponding return level in the original precipitation space is obtained with the following inverse transformation:

Since **Y**_{k} could be nonstationary, a different return level is computed for each year in order to keep the exceedance probability constant in time. Katz et al. (2002) called these time-evolving quantiles “effective return levels.” In this paper, we define the return level *R*^{T} for a certain subset of years *K* as the maximum effective return level on this subset:

This definition of a *T*-yr nonstationary return level ensures that the estimate in the original precipitation space is at least a *T*-yr return level in all years of the period:

Despite the fact that there exist more appropriate alternatives to the definition of return period in a nonstationary context (see Rootzén and Katz 2013), we chose the maximum over the period because we believe that is easier to understand and to interpret. However, for practical purposes, we recommend using a more advanced approach.

Hereafter, the maximum distributions of the original series and the standardized series are compared. The usual approach consists in assuming that the seasonal maximum from the original series **Y**_{k} for a year index *k* at a certain grid point is approximately GEV distributed:

where *ξ* is assumed time invariant. Following this, the distribution of the standardized maximum *M*_{k} can be computed as

with

The standardized maximum is also GEV distributed and is characterized with the same shape parameter. Actually, the standardization is equivalent to the usual model [Eq. (11)] but using different link functions for both the scale and the location parameters. The difference lies in the fact that the parameters of these link functions are estimated with the daily precipitations and not only with seasonal maxima. Since more data are taken into account, the estimation variance should be reduced compared to the first approach.

## 4. Results

In the section 4a we show the results obtained for summer precipitation; section 4b shows the results for the winter, spring, and fall precipitation. We recall that we adjusted both occurrence models independently to every grid point and then we selected the best one according to the BIC. Every intensity model has been fit independently of the occurrence models and independently to every grid point. Again, the best model has been chosen according to the BIC. We show in this section the results for every studied grid point.

### a. Results for summer precipitation

Figure 2 shows the occurrence expectation of summer total precipitation for the 1961 climate, for the 2100 climate, and for the difference between both climates. We see that summer precipitation occurrence increases over time in the northern part of the domain (Canada and Alaska) while it globally decreases south of the Great Lakes.

Figure 3 shows the chosen model for the precipitation intensity in summer according to the BIC. Model _{1} was never chosen and model _{8} was only chosen for 0.7% of the grid points. The most popular models were _{2} (21%) and _{3} (46%). Since these two models are stationary, summer precipitation intensity may be stationary for 67% of the grid points. As shown in Fig. 4, intensity seems to increase in the northern part of Canada and Alaska and largely decreases in Mexico. For southern Canada and the continental United States, no changes are expected in summer. Increases in both precipitation occurrence and intensity in northern regions for the summer season come from an increase of convective precipitation, both in occurrence and intensity, as shown in Fig. 5.

To show that the chosen models for intensity fit the data very well, we present the quantile–quantile (QQ) plots for six grid points subject to different meteorological conditions (see Fig. 6). These grid points under consideration are the ones containing the cities of Montreal; Chicago; Churchill; Vancouver; Washington, D.C.; and Nashville (see Fig. 1). The chosen theoretical models (see Table 3) fit the data very well up to approximately the 99th quantile. For instance, the theoretical model fits the 6541 observations above the 99th quantile. Upon this point, the theoretical models underestimate intense precipitation values. This is why we did not use these models to estimate intense precipitation. Instead, the upper tail of the distribution is modeled by the GEV distribution. Nonetheless, it seems that the estimations of time-varying expectation and variance with the considered models are quite suitable.

On the other hand, we modeled intense precipitation with the GEV distribution. Figure 7 shows the QQ plots for the fitted GEV on the 140 transformed annual maxima (i.e., the vector **M**) for the six given grid points subject to different meteorological conditions. From this last figure, we see that the GEV fits the data very well, at least for the six grid points considered. We have no reason to believe that it could be different for the remaining grid points, except perhaps for a small number of them. In the potential scenario where the fit would be poor, the return level estimation variance would be larger but we have to live with it; when the sample size from which the maximum is extracted is large enough, the GEV distribution is the only distribution with reasonable mathematical foundation to model maxima.

For nonstationary models, the 100-yr return level estimation was the maximum of the 100-yr return level computed for the 140 years of the period. The estimated 100-yr return levels for the entire domain are shown in Fig. 8. We also plot in this last figure the difference between 100-yr return level for 2100 and 1961 climates. Note that the comparison does not concern two years, but the difference of model means evaluated at these specific years. It is rather a difference between the climates of 2100 and 1961. Besides, the transient statistical models that we developed provide a climatology for every years in the studied interval 1961–2100. It can be seen in Fig. 8 that 100-yr precipitation levels are stationary for continental United States, while they increase for Canada and decrease for Mexico.

### b. Results for the other seasons

Figure 9 shows the difference in precipitation between 2100 and 1961 climate for the winter, spring, and fall seasons. For the winter season, precipitation increases over most regions in Canada and in the United States for both occurrence and intensity. It is only in Mexico that winter precipitation occurrence largely decreases. The spring and fall precipitation are somehow a transition between summer and winter. Changes in winter are mostly due to an increase in both convective and stratiform precipitation in north of Mexico (see Fig. 10).

Figure 11 shows the 100-yr return level for winter, spring, and fall precipitation, as well as the differences between 2100 and 1961 climate. The increase of winter 100-yr precipitation level is mostly located in coastal areas. For spring and fall, the increases of 100-yr precipitation level are mostly on the Pacific coast and the eastern part of continental North America. In Mexico, no change in 100-yr precipitation level is expected.

## 5. Discussion and conclusions

### a. Results comparison

Our results using a CRCM transient simulation are consistent with those from de Elía and Côté (2010). Their study investigates changes in precipitation and temperature over North America for an ensemble of 18 RCM simulations. Changes were assessed using two 30-yr periods: 1961–90 and 2041–70. They showed that summer precipitation is expected to increase in northern Canada while it is expected to decrease in Mexico and the continental United States. In the same way, our results are consistent with those from Wehner (2013) for a time-sliced analysis on an ensemble of eight simulations. Moreover, the increase of winter convective precipitation over the Gulf of Mexico, the Great Lakes, and the Atlantic coast is also consistent with the study of Paquin et al. (2014).

For the 100-yr return level, our results over the United States are consistent with those from Bukovsky and Karoly (2011), who used two time slices, 1990–99 and 2090–99, from one simulation of the Weather Research and Forecasting Model. An increase of the 100-yr precipitation intensity is expected for almost all grid points in the United States. Increases in precipitation may come during any season but summer. For Canada, CRCM projections suggests a 100-yr precipitation increase, which is consistent with the studies of Mailhot et al. (2012) based on four simulations and Monette et al. (2012). Increases in precipitation are expected for all seasons but much larger for summer in northern Canada.

### b. Contributions and limitations

This study has provided a quantitative description of transient precipitation series generated by a RCM. We believe that using a transient simulation leads to a more accurate and complete description; change estimations are more accurate and nonstationary processes are described. Even if the description of daily precipitation is better, the major gain consists in achieving a nonstationary frequency analysis to assess the return level in a changing climate. Nonstationary processes are not estimated with only few data but rather with the entire dataset. This methodology lies in the assumption stating that nonstationary processes that drive extreme values are consistent with those that drive daily precipitation. Since we statistically modeled the distribution of precipitation, the last assumption was straightforward. Nonstationary statistical distributions were useful to transform the original precipitation series into a stationary series. Thenceforth, a classical frequency analysis was achieved on the transformed series. Thus, every theoretical result and definition, such as that for the return period, is directly applicable on the transformed series. We believe that we gain in interpretation since classical frequency analysis is well known and understood.

From a meteorological point of view, large-scale precipitation may be correlated on successive days. Nevertheless, the effects of such temporal dependence of large-scale precipitation are intractable for total precipitation. From a statistical point of view, the independence assumption that we have made is quite difficult to test for both occurrence and intensity of total precipitation because these quantities are nonnormal and nonstationary. A way to verify it, and to handle it if so, would be to introduce a Markovian dependence (see, e.g., Évin et al. 2011) in the occurrence and intensity models. If such models are chosen, dependence would be significant and would be taken into account. However, such Markovian modeling for all the considered statistical models is beyond the scope of the present paper. Considering the dependence, if so, could reduce the estimation variance. On the other hand, chosen independent models considered in the present manuscript fit the data very well as shown in Fig. 5. Therefore, we believe that the consequence of such assumption is quite limited.

However, we could further improve the methodology using regionalization. Such an approach takes advantage of spatial correlation to improve estimates, provided that spatial correlation can be modeled. For our approach, we are trying to develop regional models based on the occurrence and intensity models presented in this paper. This is an issue on which we are still working.

Despite the fact that the results of this study are consistent with those in the literature, they are based on only one CRCM simulation. The results could be substantially different for another simulation. It is generally recognized that considering a large number of models should minimize the structural uncertainty from which climate models suffer (Hagedorn et al. 2005). In that sense, our study does not cover a large range of possible outcomes for precipitation, which could be achieved by considering an ensemble of different climate models. We recommend that this be done. Nevertheless, we believe that the approach presented here is a significant contribution to the study of transient simulations.

## Acknowledgments

The CRCM data have been generated and supplied by Ouranos. We acknowledge Professor Alain Mailhot for his helpful comments and suggestions. We also acknowledge the two anonymous referees for their suggestions for improving the paper.

## REFERENCES

*Climate Change 2013: The Physical Science Basis.*T. F. Stocker et al., Eds., Cambridge University Press, 1217–1308.

*An Introduction to Statistical Modeling of Extreme Values.*Springer, 209 pp.

**39B,**

*Climatic Change 2007: The Physical Science Basis*. S. Solomon et al., Eds., Cambridge University Press, 996 pp.

*Extremes in a Changing Climate,*A. AghaKouchak et al., Eds., Springer, 15–37.

*J. Geophys. Res.,*

**117,**D13106, doi:.

*Special Report on Emissions Scenarios*. Cambridge University Press, 599 pp.

**52,**

*Statistical Analysis of Finite Mixture Distributions,*Wiley, 243 pp.