## Abstract

Precipitation has a high spatial variability, and thus some modeling applications require high-resolution data (<10 km). Unfortunately, in some cases, such as meteorological forecasts and future regional climate projections, only spatial averages over large areas are available. While some attention has been given to the disaggregation of mean areal precipitation estimates, the computation of a disaggregated field with a realistic spatial structure remains a difficult task. This paper describes the development of a statistical disaggregation model based on Gibbs sampling. The model disaggregates 45.6-km-resolution rainfall fields to grids with pixel sizes ranging from 3.8 to 22.8 km. The model is conceptually simple, as the algorithm is straightforward to compute with only a few parameters to estimate. The rainfall depth at each grid pixel is related to the depths of the neighboring pixels, while the spatial variability is related to the convective available potential energy (CAPE) field. The model is developed using daily rainfall data over a 40 000-km^{2} area located in the southeastern United States. Four-kilometer-resolution rainfall estimates obtained from NCEP’s stage IV analysis were used to estimate the model parameters (2002–04) and as a reference to validate the disaggregated fields (2005/06). Results show that the model accurately simulates rainfall depths and the spatial structure of the observed field. Because the model has low computational requirements, an ensemble of disaggregated data series can be generated.

## 1. Introduction

Many storm management issues (e.g., design of urban drainage systems) require high-spatial-resolution precipitation data to accurately estimate potential risks associated with specific events. For example, this is particularly the case for flash floods in small watersheds (less than 10^{2} km^{2}) where meteorological data are needed to simulate streamflows (e.g., Duan et al. 2003; Turcotte et al. 2004).

Unfortunately, a mismatch exists between the spatial resolution of forecasted or projected rainfall fields (e.g., meteorological forecasts and future climate projections being available on grid boxes larger than 10 km) and that needed for hydrological applications. For example, Environment Canada’s operational meteorological forecasts are provided on a 33-km grid (Velazquez et al. 2009) and regional climate models (RCMs) run typically on a 25–50-km grid (Maraun et al. 2010). With the increase in computer performance, the horizontal resolution of numerical models is improving. For example, the Canadian RCM (CRCM) now runs at a 15-km horizontal resolution. Nevertheless, this high resolution is still not used in current projections. While some RCMs or meteorological models may produce simulations at this resolution, there will always be a need for a higher spatial resolution with respect to simulated rainfall.

For future climate predictions, statistical methods have been developed, using for example multiple linear regression or neural networks to model local precipitation as a function of large-scale variables such as pressure, geopotential height, and temperature. These methods are referred to as statistical downscaling methods (for a review, see Fowler et al. 2007). Easily available data (e.g., precipitation data from meteorological stations and predictor variables from reanalyses) have led to the development of numerous statistical downscaling methods in recent decades. Statistical downscaling is generally used to downscale global climate model precipitation data, but it can also be used to downscale RCM data as well (Hellström and Chen 2003). A major drawback of these methods is that even if spatial variability is taken into account (e.g., see Hughes et al. 1999), they produce precipitation at specific points only. The points usually correspond to weather stations where observations are readily available to calibrate the model. Areal estimation based on direct weighted averages or surface-fitting methods can then be used to estimate regional precipitation from point values (e.g., Dingman 2002).

Another approach is to take the areal mean rainfall from a climate or meteorological model and spatially distribute it onto a higher-resolution grid (<10 km). This approach, which is referred to as statistical disaggregation, is used to create a high-resolution map depicting a realistic spatial structure that replicates the amounts of precipitation. To be “realistic,” a precipitation map should reproduce extreme values, spatial correlation (Pegram and Clothier 2001), anisotropy, and nonlinearity of the event (Lovejoy and Schertzer 2010b). These requirements exclude the use of classical linear interpolation methods, which are not appropriate to reproduce the nonlinear spatial structure of precipitation events.

Cascade*-*based models are the most widely used statistical disaggregation technique (e.g., Harris and Foufoula-Georgiou 2001; Harris et al. 2001; Bacchi and Ranzi 2003; Deidda et al. 2006). These models are usually used to downscale short-term meteorological forecasts (Ferraris et al. 2003), but they can also be used for the disaggregation of RCM rainfall fields (Sharma et al. 2007). They are based on the multifractal properties (spatial and/or temporal) of precipitation, assuming that the moment *q* of the precipitation field, calculated at the scale *λ*, is proportional to *λ ^{K}*

^{(q)}where

*K*(

*q*) is a scaling exponent. If

*K*is a linear function of

*q*, the field is monofractal; otherwise, it is multifractal (Deidda et al. 2006). While the theory exists to generate realistic fields using a continuous-in-scale multifractal model (Lovejoy and Schertzer 2010a,b), nearly all applications use a discrete-in-scale model. The common procedure is to divide the grid box into four identical squares and distribute the precipitation according to scale-dependent parameters. This step is repeated as many times as required to achieve the desired resolution. The parameters are related from one scale to another by a simple scaling function. This modeling procedure is simple in the sense that it is easier to conceptualize than a continuous-in-scale multifractal model, and it is not computationally intensive. The simplicity of discrete-in-scale, cascade-based models is probably the reason for their popularity. However, these models do not explicitly account for spatial correlation, and it is difficult to obtain a so-called realistic spatial variability as visible discontinuities between two adjacent squares often occurred (see, e.g., Ahrens 2003 and Sharma et al. 2007).

Statistical disaggregation methods based on autoregressive processes or precipitation cells (see Ferraris et al. 2003) have also been tested, but to a lesser extent. Alternatively, there are Markovian methods, where precipitation at each pixel (i.e., high-resolution subgrid box) is a random variable conditioned by neighboring pixels. Allcroft and Glasbey (2003) applied Markovian methods by conditioning precipitation at one pixel using the 24 nearest pixels. Wheater et al. (1999), Chandler et al. (2000), and Mackay et al. (2001) disaggregated rainfall fields in three steps: (i) generation of the number of wet pixels as a function of the mean rainfall on the grid box, (ii) spatial distribution of the wet pixels [the probability of a given pixel being wet is conditioned by the state (wet or dry) of the eight nearest pixels], and (iii) simulation of rainfall depth on the wet pixels. Since the spatial distribution of the precipitation at each pixel within a grid box is a priori unknown, the disaggregation must be initialized by selecting, arbitrarily or randomly, values for each pixel. Then, Gibbs sampling, a Markov chain Monte Carlo (MCMC) algorithm, is used to “update” the values at each pixel by considering precipitation values in the neighborhood and, in doing so, creating a more realistic spatial structure.

To develop a statistical disaggregation model, high-density rainfall data, from a high-density rain gauge network or meteorological radar, are required. However, it is difficult to find such datasets that cover a long period of time (Segond et al. 2006). This explains why many disaggregation models are developed based on only a few events and why statistical downscaling is more popular than statistical disaggregation of climate model outputs. In recent years, some long-term radar data products have become available at a high resolution over large territories. These datasets allow for the characterization of numerous events and the estimation of the relationships between parameters of disaggregation models and explanatory variables.

The objective of this study is to develop a disaggregation model that is capable of producing realistic rainfall spatial structures while being conceptually simple. The rationales behind developing a simple model are

easily adaptable to different regions and/or additional explanatory variables and

more appealing for end users who need fine-resolution precipitation fields.

The proposed model uses Gibbs sampling and is applied to hundreds of daily events. Sixteen 45.6-km grid boxes are disaggregated to various pixel sizes (3.8, 7.6, 11.4, 15.2, and 22.8 km). The rainfall depth at each pixel is defined by a lognormal distribution. To produce realistic spatial structures, the expected value is related to the rainfall depth of the eight nearest pixels. Thus, the spatial correlation, intra, and also inter grid box may be intrinsically considered. The standard deviation is related to the expected value and to the atmospheric variable convective available potential energy (CAPE), which is an indicator of atmospheric instability (Blanchard 1998). To our knowledge, no previous disaggregation model using Gibbs sampling has included atmospheric variables.

This paper is organized as follows: section 2 describes the proposed statistical framework; section 3, the Gibbs sampling algorithm; section 4, the case study, including a description of the study area and the datasets; section 5, the parameter estimation; and section 6, the simulation results and comparisons with observed data.

## 2. Statistical framework

Let *R _{i}*

_{,j}be the rainfall depth at pixel (

*i*,

*j*) for a given day. The rainfall depth

*R*

_{i}_{,j}is expressed as a random variable conditioned by the rainfall depths of neighboring pixels and potentially by the physiographic variables describing the area (e.g., variation of altitude), the atmospheric variables characterizing the event (e.g., CAPE), or the rainfall depths at the previous time step to account for temporal correlation (if the time step is small).

In the present study, the area is flat; therefore, no physiographic variables were retained for the analyses. Also, even if temporal correlation might be longer than one day (e.g., in the case of large-scale precipitation associated with extratropical cyclones), no temporal correlation term is retained to spatially disaggregate the daily rainfall depth. Since local events (<10 km) have a lifetime of about an hour (Bloschl and Sivapalan 1995), it is assumed that the previous day’s 3.8-km field does not provide any supplementary and additional relevant information about the current day’s 3.8-km field, given that the current day’s 45.6-km precipitation field is known. Preliminary analyses indicated that the CAPE variable has a significant impact on spatial variability. As such, it is the only atmospheric variable used in the present study. Note that Perica and Foufoula-Georgiou (1996) used CAPE in their multiscale disaggregation model.

The random distribution of *R _{i}*

_{,j}must be constrained to nonnegative values. Rainfall depth can be represented by a mixture distribution: either it does not rain and the distribution has a mass at zero with probability 1 −

*p*, or it rains and the distribution is continuous on positive values with probability

*p.*This approach has been applied in many studies (see Over and Gupta 1994; Guillot and Lebel 1999; Onibon et al. 2004; Chandler 2005). Alternatively, when

*p*is close to one, it may be easier to use a continuous distribution and set to 0 values close to zero. A mixture distribution may be more accurate to simulate close-to-zero values, but for daily events, these smaller values are not of interest. Therefore, a continuous distribution, without a mass at zero, was used in this study.

Gamma (Guillot and Lebel 1999; Wheater et al. 1999; Onibon et al. 2004) and lognormal distributions (Pathirana and Herath 2002; Kundu and Siddani 2007) are commonly used to characterize greater-than-zero daily rainfall values. In the present study, a lognormal distribution was used (section 5).

Parameters of the *R _{i}*

_{,j}distribution can be expressed in terms of the first two centered moments: the expected value and the standard deviation. To account for possible anisotropy, the expected value is a linear combination of the precipitation coming from four main directions (45°, 90°, 135°, and 0°), defined by the eight nearest pixels: , , , and (see Fig. 1c) with .

Preliminary analyses showed that the components of **A** are highly correlated. This implies that if the expected value is a linear function of the components of **A**, interpretation of the estimated parameters would be difficult. To find uncorrelated variables, principal component analysis (PCA) was applied to the pixels over the calibration period. Four physically meaningful and weakly correlated variables were thus defined. The expected value is then expressed as a function of these variables:

where is the mean of vector **A** (i.e., the mean of the eight neighboring values) and is a vector of parameters. The first parameter *β _{d}* accounts for the distance between neighboring pixels, all directions combined, and

*β*and

_{x}*β*account for anisotropy. The value of

_{+}*β*should be positive since

_{d}*R*

_{i}_{,j}should be closer to the average daily values of the four nearest neighbors than to the average daily values of the four next nearest neighbors. A positive

*β*value indicates a stronger correlation in the 45° direction than in the 135° direction, while a positive

_{x}*β*value indicates a stronger correlation in the 90° direction than in the 0° direction. It can be shown that Eq. (1) could be rewritten as a linear function of the four components of

_{+}**A**with the sum of the four parameters equal to 1. This shows that the mean of the neighboring pixels is stochastically preserved.

For the standard deviation of *R _{i}*

_{,j}, it is assumed that an increase of the expected value and of the CAPE value would lead to an increase of the standard deviation, but the exact relationship is unknown. Perica and Foufoula-Georgiou (1996) expressed the two parameters of their cascade-based model as linear functions of CAPE. Here, based on empirical analyses, the simplest standard deviation expression for

*R*

_{i}_{,j}consistent with the aforementioned assumption is the following:

where is a vector of estimated parameters, with *θ _{k}* > 0 for

*k*= (1, 2, 3);

*C*is the mean daily value of CAPE (J kg

^{−1}), which is a measure of atmospheric instability, over the whole area; [Eq. (1)]; and

*ɛ*is a random noise accounting for the variability unexplained by the model (see paragraph below). Thus, for a fixed value of CAPE, the standard deviation increases as the expected value of

*R*

_{i}_{,j}increases. Onibon et al. (2004) used the same relationship between the standard deviation and the expected value with

*θ*

_{3}= 0.60. Also, for a fixed value of the expected value, the standard deviation is assumed to increase linearly with the CAPE value.

The relationship between the expected value, CAPE, and the standard deviation may be good on average, but it is too simple to be the exact relationship between the standard deviation and the two predictor variables. As will be presented in section 5, there are some differences between the estimated standard deviation obtained from Eq. (2) (without *ɛ*) and that calculated from the reference pixels in the calibration period. This is why *ɛ*, a random term with null expected value, must be added. In the present application, *ɛ* is uniformly distributed and its boundaries are estimated in section 5.

## 3. Gibbs sampling

Gibbs sampling is an MCMC algorithm generating a multiple random variable too complex to be generated directly but for which the conditional distribution of each component is known. In the present case, for a given day, the rainfall depth for all pixels (the *R _{i}*

_{,j}s) is the multiple random variable, and the conditional distribution of a pixel is given in section 2.

The principle of Gibbs sampling is to arbitrarily allocate initial values to each pixel, and then to update these values one at a time according to the known conditional distribution of each pixel. When the number of updates of each pixel is large enough, Gibbs sampling converges to the target distribution. In the actual case, it means that the disaggregated rainfall fields would have a realistic spatial structure, assuming that the hypotheses introduced in section 2 adequately describe the spatial structure of the rain. For a more complete description of Gibbs sampling, the reader is referred to Geman and Geman (1984) and Roberts and Smith (1994).

In the present study, Gibbs sampling was used to disaggregate the known daily rainfall of *T*_{1} × *T*_{2} grid boxes (with *T*_{1} = *T*_{2} = 4) to obtain a higher-resolution map of *P*_{1} × *P*_{2} disaggregated pixels per grid box (for a 45.6-km grid box, *P*_{1} = *P*_{2}; from 2 for 22.8-km pixels, to 12 for 3.8-km pixels). Let *L*_{1} × *L*_{2} be the size, in number of pixels, of the whole area to be disaggregated (with *L*_{1} = *L*_{2} = *T*_{1} × *P*_{1} = 8 to 48). Figure 2 illustrates the Gibbs sampling procedure. Suppose that all parameters (** β**,

**, and**

*θ**ɛ*boundaries) are known. For each day, the

*L*

_{1}×

*L*

_{2}rainfall field is generated as follows:

The mean rainfall value of a grid box is used as the initial rainfall depth of each pixel (Fig. 2a). Then, a random value for

*ɛ*is generated from a zero-mean uniform distribution.Vector

**A**for pixel (1, 1) is calculated and a new rainfall value for the pixel is generated from the lognormal distribution according to Eqs. (1) and (2).Step (ii) is repeated for pixel (1, 2) and so on until the last pixel (

*L*_{1},*L*_{2}) of the area is reached. A multiplicative factor is applied to each grid box to ensure conservation of rainfall depths of each grid box. This is the first iteration (Fig. 2b).Steps (ii) and (iii) are repeated until a sufficient number of iterations have been generated to eliminate the effect of the initial condition. This is referred to as the burn-in period. This first simulated rainfall field is saved (Fig. 2d). A new random value for

*ɛ*is generated.Step (iv) is repeated until the correlation with the previous simulated rainfall field is minimized. This is referred to as the autocorrelated time. This simulated rainfall field is saved (Fig. 2e). A new random value for

*ɛ*is generated.Step (v) is repeated until the desired number of simulated rainfall fields for a given day has been generated.

At Steps (ii) and (iii), when estimating **A** for pixels along the boundaries, some neighboring pixels are outside the area. If a pixel needed to calculate a component of **A** is actually located outside the area, then the component is estimated using its only available pixel value. For example, for pixel (1, 2), the component *A*_{/} is equal to the rainfall value at pixel (2, 3) only (Fig. 1b). If the two pixels needed to calculate a component of **A** are outside the area, which happens at the four corners, then the component is defined by the mean of the two nearest 0°–90° neighbors. For example, for pixel (1, 1), the component *A*_{\} is the mean of the rainfall values at pixels (1, 2) and (2, 1) (Fig. 1a).

Note that the multiplicative factor applied to conserve the areal mean on a grid box, following each iteration, is necessary since Eq. (1) preserves the mean stochastically, rather than exactly. The disaggregation results showed that the multiplicative factor is generally close to 1, except for some grid boxes with low precipitation values. This may be caused by the fact that when only a small number of pixels are wet, there is not enough data to ensure the convergence to the mean gridbox value.

## 4. Case study

This section describes the study area as well as the data used for the parameter estimation and for the application of the disaggregation model.

### a. Study area

The study area (182 × 182 km^{2}) is located at the border of the states of Georgia and Florida, in the southeastern United States (Fig. 3). It is subdivided into 16 grid boxes (*T*_{1} = *T*_{2} = 4; resolution = 45.6 km). For disaggregation, each grid box is subdivided into *P*_{1} × *P*_{2} = 12 × 12 = 144 pixels (Fig. 3). The spatial resolution of each pixel is 3.8 km.

This area is chosen for the following reasons: 1) almost all precipitation is rain, 2) it receives a significant amount of rain with annual rainfall depth ranging from 845 to 1612 mm during the 2002–06 period, 3) there are large-scale to highly convective rainfall events in the area, and 4) the topography is relatively flat with altitudes lower than 150 m. While disaggregation might be more interesting for hydrological applications in highly variable topographic areas, the proposed model is developed on a flat area for simplicity.

### b. Precipitation data

Daily precipitation data are taken from the stage IV analysis, which is a mosaic produced at the National Centers for Environmental Prediction (NCEP). Precipitation fields are produced using radar and gauge stations. Data analysis is performed at the 12 River Forecast Centers (RFCs) across the continental United States. Automatic quality control (QC) is done soon after any rain depth is registered, followed by a manual QC. The QC may differ from one RFC to another. The study area is located within the Southeast RFC. For more details about the data, see Lin and Mitchell (2005).

Stage IV provides a 1121 × 881 polar stereographic grid. In this study, pixels 265–312 on the *y* axis and 897–944 on the *x* axis were used. Data were rounded to the nearest integer (in mm day^{−1}). The data from 2002 to present are available online (http://data.eol.ucar.edu/cgi-bin/codiac/fgr_form/id=21.093). The precipitation data from 2002 to 2006 are used to match the available atmospheric data. The 2002–04 period (1093 available days) is used to estimate model parameters (** β**,

**,**

*θ**ɛ*, burn-in duration, and autocorrelated time).

The 2005/06 period (729 available days) is used for model validation. The grid boxes are built by aggregating 12 × 12 stage IV pixels to create a mesoscale map of 4 × 4 grid boxes with a resolution of 45.6 km (larger squares in Fig. 3). A total of approximately 5000 wet gridbox values are available. To produce different resolutions of disaggregated fields (sections 5 and 6), disaggregation is performed on 3.8-km pixels (original stage IV resolution) and on 7.6-km (2 × 2 aggregated stage IV pixels), 11.4-km (3 × 3), 15.2-km (4 × 4), and 22.8-km (6 × 6) pixels.

The performance of the disaggregation model is evaluated by comparing the statistical properties of the disaggregated rainfall fields with the statistical properties of the stage IV observations over the 2005/06 period.

### c. Atmospheric data

CAPE values are available at 3-h intervals from the NCEP North American Regional Reanalysis (NARR) starting in 1979. This reanalysis is provided on a 249 × 277 northern Lambert conformal conic grid covering North America with a horizontal resolution of approximately 32 km in the southern portion of the map and with 29 pressure levels (for details, see Mesinger et al. 2004).

A unique daily CAPE value over the study area is obtained by averaging the data over 24 h for the 35 NARR grid points covering the area. For each precipitation grid box, a different daily CAPE value could have been used by taking the maximum of the 3-h CAPE values instead of the average over 24 h. Preliminary analyses showed that it did not improve the estimation of the standard deviation (not shown). Other atmospheric variables simulated by RCMs, meteorological models, or produced by the NARR could be of interest (e.g., wind speed and direction), but are not considered in the present study.

## 5. Parameter estimation

Only wet pixels (i.e., pixels with at least 1 mm on a given day) and pixels with at least one wet neighbor are used for model calibration during the 2002–04 period (referred to as the calibration pixels). Parameters are estimated on 3.8-, 7.6-, 11.4-, 15.2-, and 22.8-km pixels. For the 22.8-km pixels, about 22 000 pixel days are retained; for the 3.8-km pixels, more than 800 000 pixel days are retained. The very large number of data presents a challenge for parameter estimation.

The maximum likelihood (ML) method, assuming lognormal distributed data, was first used to estimate the parameters on the 3.8-km pixels. However, the ML method did not converge. A possible explanation is that the proposed statistical framework that describes the rainfall depth at a pixel (section 2) is too simplistic. Even if it satisfactorily explains the major part of the variability, it cannot precisely describe the rainfall depth distribution of a pixel. While the assumptions made in section 2 on the expected value, the standard deviation, or the statistical distribution may only be violated by a low percentage of pixels, it still represents a large number of pixels that could significantly reduce the ML value. Perhaps these pixels are weighted too heavily, which leads to unsatisfactory estimated parameters. It would be possible to reduce the weight of these problematic pixels or simply eliminate them from the estimation procedure. However, this is a daunting task considering the amount of data (i.e., even if only 0.1% of pixels are problematic, that represents about 1000 3.8-km pixels) and even more importantly, eliminating pixels would be quite arbitrary. Removing the pixels with rainfall depths far from what is expected would produce an artificial confidence in the proposed model. For all these reasons, ML estimation was not retained.

As an alternative, model parameters ** β** and

**are estimated in two steps based on the minimum sum of squared errors (SSE). In Eq. (1), estimated vector**

*θ***minimizes SSE between observed and expected rainfall depths for all wet calibration pixels. The estimated components of**

*β***are shown in Table 1. Anisotropy is observed, with a stronger 45° than 135° component (**

*β**β*> 0) for the five scales. This agrees with the NARR wind data (not shown), which are stronger on average in the NW–SE direction than in the SW–NE direction. The parameter

_{x}*β*slightly decreases as the pixel size increases. Figure 4 shows that the residuals (i.e., observed–expected 3.8-km rainfall depths) are quite well distributed around zero for each class of expected values. The lower classes show no extreme negative residuals since precipitation is bounded to zero. Note that the variability is larger for the larger classes. Since the estimation is made according to the SSE, this implies that for

_{d}**estimation, the weight of a pixel increases as its expected rainfall depth increases. This is a desired behavior since in daily event hydrological applications it is more important to accurately represent larger values than smaller values.**

*β*Depending on the daily values of the neighbors, Eq. (1) may return a negative expected value for rainfall. During the calibration period, this occurs approximately on 6% of the 3.8-km pixels. But these negative values are very close to 0. Indeed, less than 0.1% of the pixels have expected values smaller than −1. Still, a minimal threshold for expected value must be fixed. A threshold value of zero would imply that the simulated value would be 0 (a lognormal with mean zero implies a point mass at zero). Since some pixels in the calibration period with null or negative expected values from Eq. (1) have positive precipitation (1 mm or sometimes more), a positive threshold close to zero should be fixed. In the algorithm, this threshold is fixed at 0.2 mm since it is the smallest threshold that may produce positive (around 1 mm) simulated values. Note that disaggregation applications are more likely to focus on extreme local values than close-to-zero values; it is then assumed that the threshold value has a limited hydrological impact.

As mentioned in section 2, the 3.8-km pixels are first used to estimate the standard deviation expression and parameters [Eq. (2)]. The residuals are separated into 3925 classes, formed from 157 expected rainfall depth classes × 25 classes of CAPE values. These numbers are chosen to accurately describe the relationship between the standard deviation, the expected rainfall depth, and CAPE with the constraint that a sufficient number of pixels are needed in a class to produce a robust estimation of the standard deviation. Therefore, 3826 classes with at least 100 3.8-km pixels coming from at least five different days are retained. For each class, the standard deviation of the residuals is calculated. The estimated vector ** θ** minimizes the SSE between the calculated standard deviations for these classes and their expected values according to Eq. (2) (with

*ɛ*set to 0). Estimated

**values are presented in Table 1. The value of**

*θ**θ*

_{3}is very close to the value obtained by Onibon et al. (2004) for a larger spatial scale (100 × 100 km

^{2}grid).

To evaluate whether CAPE had a significant impact on the estimation of the standard deviation, parameters of Eq. (2) were estimated again, with *ɛ* and *θ*_{2} set at zero. For this investigation, the standard deviation was no longer related to CAPE—it was only related to the expected value. Figures 5a,b illustrate, for each of the 3826 classes, the observed standard deviation versus the estimated standard deviation when accounting for (*ɛ* = 0) and not accounting (*ɛ* and *θ*_{2} set at 0) for CAPE, respectively. The figures show that accounting for CAPE reduces the dispersion around the *x = y* line as the SSE is four times lower in Fig. 5a than that in Fig. 5b. It clearly demonstrates that accounting for CAPE significantly improved the estimation of the standard deviation.

To evaluate whether the linear relationship between the standard deviation and CAPE is adequate, a linear regression is performed for classes with the same expected values. Thereby, 156 regressions with at most 25 data points per regression are made. The last class of expected value is ignored since it did not contain enough different CAPE values. The 156 values of intercept and slope (defined by and ) are plotted and compared to their expected values based on the estimated ** θ** value. Figure 6 shows that the relationship is valid for the expected range of rainfall depths. As Perica and Foufoula-Georgiou (1996), who expressed their scaling exponent by a linear function of CAPE, the parameter

*θ*

_{3}was expressed as a linear function of CAPE to test whether an additional parameter in Eq. (2) would improve the standard deviation estimation. However, the results show no significant difference with

*θ*

_{3}constant. Equation (2) is then the simplest equation giving the best standard deviation estimation.

Figure 5a shows that the relationship is quite good between observed and expected standard deviation, but there is some dispersion around the *y = x* line. This dispersion can be seen as the variability unexplained by the model. The parameter *ɛ* is included to account for this behavior [Eq. (2)]. As most of the observed standard deviations are between 1 ±0.35 times the expected standard deviation, *ɛ* is defined as a uniform random variable within the interval (−0.35, 0.35).

The ** θ** parameters for the other pixel sizes are estimated in the same way, except that the number of estimation classes is reduced, since fewer pixels are available. The same criterion, that is to have at least 100 pixels in a class coming from at least five different days, is applied. Parameter

*θ*

_{2}, accounting for CAPE, is almost two times larger for 7.6–22.8-km pixels compared to the 3.8-km pixels (Table 1). Since less estimation classes are used for larger pixels, extreme pixel values may be underrepresented, possibly leading to an increase of CAPE influence. The parameter

*θ*

_{1}increases, and so does the standard deviation, as the pixel size increases, suggesting that the neighboring pixels are less correlated as the pixel size increases. The exponent parameter

*θ*

_{3}seems uncorrelated to the scale. The parameter

*ɛ*is not estimated for pixels larger than 3.8 km since too few estimation classes were available.

Now that the parameters of Eqs. (1) and (2) are estimated, it is possible to evaluate the cumulative distribution function (CDF) of the rainfall depths at each calibration pixel for the following common two-parameter distributions: lognormal, gamma, and normal. For each wet pixel, the expected value and standard deviation are calculated from the estimated parameters and then the CDF value is calculated for each distribution. Theoretically, if a distribution fits the data perfectly, CDF values should be uniformly distributed in the (0, 1) interval.

Figure 7a shows histograms (100 bins) of lognormal CDF values for all of the 3.8-km wet calibration pixels (about 700 000 pixels). The CDF is not uniformly distributed. Extreme low values (CDF < 0.01) occur about 2.5 times more often than they should have if rainfall depth truly followed a lognormal distribution, while extreme high values (CDF > 0.99) occur about two times more often than expected. This means the lognormal distribution is not skewed enough to reproduce the tails of the rainfall statistical distribution. Also, there are not enough moderately small CDF values (between 0.05 and 0.3), which indicates that the lognormal distribution would produce too many moderately small rainfall values. The same graph for gamma distribution shows a very similar pattern (not shown). A lognormal distribution has a positive skew (median < mean), therefore the normal distribution is also tested. CDF values are more symmetrically distributed, but the extremes values are more poorly represented (Fig. 7b). Since extreme values are important for hydrological applications, the lognormal distribution is retained.

The burn-in period is determined by disaggregating the calibration period, from 45.6-km grid boxes to 3.8-km pixels, using independent Gibbs sampling runs with different numbers of iterations. Then, the differences of a given statistic between two independent runs with the same number of iterations are compared with the differences between a run with *x* iterations and an independent run with a much larger number of iterations. Formally, let *n* be the total number of grid boxes that are disaggregated over the whole period and and be the values of a statistic (e.g., median, maximum rainfall, and standard deviation) on the *i*^{th} grid box (*i* = 1, … , *n*) for two independent runs with *x* iterations. The following ratio of differences is then calculated:

where *y* ≫ *x.* If the ratio is close to 1, then *x* is a sufficient burn-in period. This ratio is calculated for different statistics and different *x* values with *y* = 1000. Results suggest that 40 iterations are sufficient for the burn-in period. Using *y* = 500 led to the same conclusion.

To determine the autocorrelated time, all days from the calibration period are disaggregated with the parameter values in Table 1 (3.8-km pixel size), except that *ɛ* is fixed at 0. Simulated maps after the burn-in period are retained at every five (5) iterations until 50 simulated maps for each day (after 40, 45, … , 285 iterations) are produced. Then, autocorrelations of lag 5–100 are used to calculate different statistics in the grid boxes where at least half of their pixels are wet. For each of the statistics, the results suggest that mean autocorrelation functions become stable after 60 iterations.

As the autocorrelated time is longer than the burn-in period, if one wants more than one simulated map, it would be quicker to perform independent Gibbs sampling runs and stop each run after 40 iterations than to perform a single Gibbs sampling run and to retain a new map after every 60 iterations. However, the classical approach is applied; that is, the daily simulated maps shown in the next section come from a single Gibbs sampling run.

## 6. Disaggregation performance

Results from the disaggregation of 2005/06 daily rainfall depths over the study area are presented in this section. The sixteen 45.6-km grid boxes are disaggregated to 3.8-, 7.6-, 11.4-, 15.2-, and 22.8-km pixels. For each day and each pixel sizes, 10 disaggregated maps were retained (after 40, 100, 160, … , 580 iterations).

Disaggregation performance of various pixel sizes may depend at which spatial scale the results are evaluated. This is why the evaluation is made on various pixel sizes for all the disaggregations performed. For a given evaluation pixel size *λ _{e}* and a disaggregation pixel size

*λ*, the RMSE is calculated as follows:

_{d}where and are the mean precipitation volume of the 10 disaggregated maps and the observed precipitation volume (from stage IV) at the *k*^{th} day on the evaluation pixel (*i*, *j*), respectively; *N*(*λ _{e}*) is the number of evaluation pixels in a side of the study area (from 8 when

*λ*= 22.8 km to 48 when

_{e}*λ*= 3.8 km) and

_{e}*D*is the number of days (=729). Table 2 shows the RMSE values for each disaggregation and evaluation pixel sizes. The RMSE from disaggregation is always lower than the RMSE from nondisaggregated grid boxes. In general, the smaller the disaggregated pixel is, the lower is the RMSE. The fact that disaggregation produces pixel values closer on average to the reference value than the mean areal grid box shows that the model is able to account for the neighboring grid boxes that provide the mesoscale structure of the field.

The mean spatial performance (2005/06) of the model is shown in Fig. 8 for the 3.8-km and the 22.8-km disaggregated fields. These results indicate that disaggregation can preserve the global spatial structure, with the 3.8-km disaggregated field slightly better than the 22.8-km disaggregated field for both evaluation pixel sizes, as already introduced in Table 2. In the following figures, only the 3.8-km disaggregated fields are analyzed with respect to the 3.8-km grid (original resolution of stage IV data).

A spatial disaggregation model must provide reliable rainfall depth and realistic spatial structures at each time step. Thus, the performance of the model is evaluated using daily statistics. Observed and disaggregated median pixel values on a grid box are compared (Fig. 9). The large-scale (CAPE < 1000 J kg^{−1}) and more convective (CAPE > 1000 J kg^{−1}) daily events were well delineated. The simulated values should cover the observed values without too much dispersion. The median pixel value is well simulated, as shown in Fig. 9. Since the simulated mean rainfall depth over a grid box is forced to be equal to the observed mean, the median is relatively easy to simulate, especially for grid boxes with small rainfall depths. Nonetheless, the higher median values are accurately simulated as well.

Figure 10 shows that the daily maximum 3.8-km pixel values on a grid box are, in general, accurately simulated; that is, the highest and the lowest simulated values are always of the same order of magnitude. However, the highest values, which characterize large-scale events, are underestimated (Fig. 10a). A closer look at the six-highest maximum values, which occur on four different days, provides two potential explanations. Of these four days, one day had a strong 45° rainfall gradient and another day had a strong 135° gradient. The unique estimation of ** β** used for all days does not allow for this kind of spatial structure. For the other two maximum values, they are simply much larger than their neighbors. It suggests, as shown in Fig. 7, that the lognormal distribution can accurately describe rainfall depths for most pixels but that it is not skewed enough to accurately describe extreme values.

The main difficulty of any disaggregation model is not to provide correct rainfall depths within a grid box, but rather to accurately simulate daily spatial correlations at high resolution. To evaluate the spatial correlation, empirical semivariograms are calculated for each day for both the observed and the disaggregated maps. Since the maps are regular grids, a lot of pixels were available for all distance classes. Figure 11 presents the observed and simulated semivariance for nearest neighbor 3.8-km pixels (i.e., for pairs of pixels sharing a side). Once again, disaggregation generally simulated the semivariance well, except for higher values of the large-scale events that are mostly underestimated. Once again, the single estimation of ** β** that is used for all days and the use of a (lognormal) distribution that is not highly skewed for extreme values are the most likely causes of this underestimation. The same pattern is observed for the other distance classes (up to three pixels; not shown for conciseness).

## 7. Conclusions

The Gibbs sampling disaggregation model introduced in this paper accurately reproduced the high-resolution spatial structure of the majority of the daily observed rainfall depths. The model assumes that rainfall depth at a pixel is related to the rainfall depth of the eight nearest neighbors and to the atmospheric variable CAPE. One of the main advantages of the model is that iterations of the Gibbs sampling can generate a realistic (without discontinuity) spatial structure for an event, despite its conceptual simplicity (only six parameters). The algorithm used to disaggregate is easy to compute and is not computationally intensive. The model could be used to estimate local impacts of precipitation events (e.g., estimation of peak flows on a small watershed when used as input to a distributed hydrological model) when only mesoscale rainfall data are available.

To account for seasonality, it is common to estimate monthly parameter values (see, e.g., Chandler et al. 2000 and Wheater et al. 2005). Instead, the proposed model uses the atmospheric variable CAPE to characterize the type of each event. This variable is simulated by RCMs. By relating model parameters to atmospheric variables instead of months, the model is more directly related to the physics of the event. For example, if more convective events occur in the future, the proposed model will disaggregate the rainfall accordingly. Therefore, the model is more appropriate to face climate perturbations caused by climate change or by atmospheric phenomena such as the El Niño Southern Oscillation (ENSO) or the North American Oscillation (NAO). A drawback of the present study is that CAPE data do not come from the same source as the precipitation data. Precipitation data are used only on a 24-h time step. A better relationship between CAPE and spatial variability of precipitation would be expected if the temporal resolution of precipitation fields was finer.

For shorter time step applications, temporal correlation with the previous time step disaggregated map might be considered. For example, the resulting disaggregated map at time *t* could be a weighted mean of the disaggregated map at time *t* − 1 and a disaggregated map at time *t* from the proposed model without a temporal term. The weight of the previous time step map would be larger as the time step decreases.

Parameter estimation was performed using all 2002–04 days. This implies that estimated parameters are assumed to be representative of any type of events. The results shown in section 6 suggest that the model is accurate for most events, except for four extreme large-scale days. This may be caused by the fact that these events are underrepresented in the sample used for parameter estimation. Thus, parameter estimation must be performed according to the application being used. For example, if one wants to disaggregate only extreme events, parameter estimation should be done accordingly. In this case, a more skewed distribution than the lognormal may likely be more appropriate (e.g., the generalized extreme value distribution). However, since high-spatial-resolution precipitation data are not available for a long period of time, it is difficult to find enough extreme events to accurately estimate the parameters.

In the present application, CAPE was the only atmospheric variable used. Other atmospheric and physiographic variables may partially explain the spatial structure of the rainfall. For example, wind and/or storm direction may have a significant impact, especially in areas with complex topography. Even if these meteorological variables may have influenced the observations, they were not included in the disaggregation model since the study area is characterized by a relatively flat topography. Nonetheless, in future studies, the model will be applied to a region characterized by more complex topography and these explanatory variables could eventually be considered. Indeed, incorporating the topographic influence on the rainfall spatial structure is the next step toward improving the model. Preliminary analyses on a mountainous region showed that the mesoscale impact of mountains is reproduced after disaggregation with the same parameter values used in the present study (without accounting for topography).

Disaggregation of various pixel sizes showed that generally, the highest resolutions produce the smallest differences with observations, no matter the evaluation pixel size. Anisotropy parameter *β _{x}* slightly increases while influence of the four nearest neighbors, characterized by

*β*, slightly decreases as the pixel size increases. For the same expected value and CAPE value, the standard deviation is higher as the pixel size increases, meaning that the influence of neighboring pixels to predict the pixel value is weaker as the pixel size increases. It should be noted however that fewer pixels are available to estimate parameters when the pixel size increases, augmenting the uncertainty on parameters values. Variability not explained by the model was accounted for with the

_{d}*ɛ*parameter. In the present application,

*ɛ*was uniformly distributed in the (−0.35, 0.35) interval. Other distributions for

*ɛ*(e.g., Gaussian) or other ways to account for the extravariability not explained by the model could have been used. Application of the model to other areas would provide more information about the behavior of the parameters. The more datasets used to develop the model, the more accurate will be the parameter estimation. In particular, the

*ɛ*parameter could be removed if one can improve the relationship between the explanatory variables and the daily precipitation at the pixel scale.

Ultimately, this model will be applied to generate rainfall depth data that will be used as input to a hydrological model. The spatial resolution and the possibility to rapidly create an ensemble of simulated maps from a single mesoscale model run (from a climate or a meteorological model) represent the added values of disaggregated rainfall fields.

## Acknowledgments

This work was supported in part by the Fonds Québécois de la Recherche sur la Nature et les Technologies, the Natural Sciences and Engineering Research Council of Canada, and Ouranos—a consortium on regional climatology and adaptation to climate change, Montreal, Quebec, Canada. Stage IV and NARR data were provided by NCAR/EOL with sponsorship from the National Science Foundation (http://data.eol.ucar.edu) and NOAA (http://www.esrl.noaa.gov/psd/data/gridded/data.narr.monolevel.html). The authors thank Bronwyn Pavey for her attentive editing and Guillaume Talbot, research assistant at the Institut National de la Recherche Scientifique, Centre Eau, Terre et Environnement (INRS-ETE), for his computational expertise and for the generation of some figures. The authors would also like to thank René Roy of Hydro-Québec, Michel Slivitzky, emeritus professor at INRS-ETE, and the three anonymous reviewers for their helpful comments.