Projecting Flood-Inducing Precipitation with a Bayesian Analogue Model

by   Gregory P. Bopp, et al.
Penn State University

The hazard of pluvial flooding is largely influenced by the spatial and temporal dependence characteristics of precipitation. When extreme precipitation possesses strong spatial dependence, the risk of flooding is amplified due to catchment factors that cause runoff accumulation such as topography. Temporal dependence can also increase flood risk as storm water drainage systems operating at capacity can be overwhelmed by heavy precipitation occurring over multiple days. While transformed Gaussian processes are common choices for modeling precipitation, their weak tail dependence may lead to underestimation of flood risk. Extreme value models such as the generalized Pareto processes for threshold exceedances and max-stable models are attractive alternatives, but are difficult to fit when the number of observation sites is large, and are of little use for modeling the bulk of the distribution, which may also be of interest to water management planners. While the atmospheric dynamics governing precipitation are complex and difficult to fully incorporate into a parsimonious statistical model, non-mechanistic analogue methods that approximate those dynamics have proven to be promising approaches to capturing the temporal dependence of precipitation. In this paper, we present a Bayesian analogue method that leverages large, synoptic-scale atmospheric patterns to make precipitation forecasts. Changing spatial dependence across varying intensities is modeled as a mixture of spatial Student-t processes that can accommodate both strong and weak tail dependence. The proposed model demonstrates improved performance at capturing the distribution of extreme precipitation over Community Atmosphere Model (CAM) 5.2 forecasts.



page 17

page 23


A semiparametric Bayesian model for spatiotemporal extremes

In this paper, we consider a Dirichlet process mixture of spatial skew-t...

Bayesian Modeling of Air Pollution Extremes Using Nested Multivariate Max-Stable Processes

Capturing the potentially strong dependence among the peak concentration...

A Hierarchical Max-infinitely Divisible Process for Extreme Areal Precipitation Over Watersheds

Understanding the spatial extent of extreme precipitation is necessary f...

Modeling Spatial Data with Cauchy Convolution Processes

We study the class of models for spatial data obtained from Cauchy convo...

Hierarchical Transformed Scale Mixtures for Flexible Modeling of Spatial Extremes on Datasets with Many Locations

Flexible spatial models that allow transitions between tail dependence c...

A flexible Bayesian hierarchical modeling framework for spatially dependent peaks-over-threshold data

In this work, we develop a constructive modeling framework for extreme t...

Spatial Modeling of Heavy Precipitation by Coupling Weather Station Recordings and Ensemble Forecasts with Max-Stable Processes

Due to complex physical phenomena, the distribution of heavy rainfall ev...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In this paper, we develop a mixture model for daily precipitation allowing for varying spatial dependence for different storm magnitudes. Instead of constructing a parametric formulation of the atmospheric dynamics regulating precipitation, we develop a Bayesian analogue method that can accommodate complex temporal dependence patterns, wherein precipitation analogues are established by identifying similar synoptic-scale atmospheric patterns from a historical library of observed climate states that are associated with precipitation outcomes.

Global climate models (GCMs) are our main tool for forecasting future climate conditions, providing a means to assess potential changes in the frequency and magnitudes of events such as heat waves, droughts and floods, all of which can have profound impacts on human health. Unfortunately, the coarse spatial resolution of GCMs cannot resolve fine scale hydro-meteorological processes associated with precipitation extremes (Boe06; Maraun10). However, GCMs also produce smoothly varying atmospheric variables such as atmospheric pressure and temperature whose spatio-temporal patterns are closely linked with precipitation outcomes (Xoplaki04; Raziei12).

Analogue methods try to address the problem of making forecasts in the presence of complex temporal dependence without a model parameterization of any geophysical dynamics. In their simplest form, analogue methods match the current climate state to observed climate states from a library of historical observationsc in order to forecast some future quantity (e.g. precipitation tomorrow) with the observed successor of the historical match (e.g. precipitation on the day following the historical match). While the atmosphere is known to be a chaotic dynamical system that is unstable under slight perturbations of initial conditions (Lorenz69), analogue approaches are justified by the tendency of that system to regularly revisit subsets of the phase space over time. Analogue methods were originally developed as empirical tools for short-term weather forecasting (Krick42) and climate modeling (Barnett78)

, but researchers have begun to recognize their utility in a variety of contexts, including machine learning

(Zhao16; Lguensat17), wind-speed modeling (Nagarajan15), and air quality monitoring (Delle14). Historically, analogue methods have been empirical, somewhat ad-hoc tools, but recently there have been attempts at putting these into a probabilistic framework (McDermott16; Mcdermott18).

Precipitation has been the subject of extensive study as it plays a central role in agriculture, flood-risk, and hydrology. Distributional forecasts of precipitation are critical for water management, infrastructure design, and developing disaster preparedness strategies. While modern numerical forecasting weather models have achieved considerable success at making short term forecasts by approximating solutions to the complex atmospheric processes governing the generation of precipitation, they typically do not account for uncertainty in their inputs. The function of stochastic weather generators, in contrast, is not to make short-term forecasts, but to accurately reproduce the distributional characteristics of precipitation on a fine spatial grid, while quantifying the uncertainty associated with those estimates

see Ailliot15 for a review.

As precipitation is an inherently spatial phenomenon, several approaches have been proposed to model the spatial dependence features of both occurrences (presence/absence of rain) and intensities (positive rain accumulations). A popular stochastic weather generator was proposed in a seminal paper by Wilks98 that models precipitation across multiple sites in two parts: (1) a two-state Markov process controlling the occurrence of precipitation at a given site, and (2) a precipitation intensity model that accounts for spatial dependence between sties but ignores temporal dependence. Berrocal08 and Kleiber12 develop similar two-stage spatio-temporal models for observations from rain-gauge stations, wherein a latent Gaussian process is thresholded to model precipitation occurrence, while another marginally transformed Gaussian process controls precipitation intensity. While most models assume the same spatial dependence structure across precipitation intensities, Bardossy09 raise the issue of varying spatial dependence types across different intensities and aim to address it with an empirical copula approach. In this paper, we address the issue of varying spatial dependence types by modeling precipitation as a mixture of Student-t processes with different spatial correlation structures. Our approach to this problem is similar to that of Gelfand05, who treat precipitation intensities as a Dirichlet process (DP) mixture of Gaussian processes. In the context of extreme value modeling, Fuentes13 have also explored DP mixtures of Gaussian processes to flexibly model spatial dependence, but in the context of modeling maxima, which support marginal transformations to generalized extreme value (GEV) margins. Gaussian processes and their mixtures, however, are characterized by weak tail dependence, and in the presence of strong spatial dependence among extremes, their application can lead to underestimation of flood risk. To allow for stronger tail dependence, similar approaches have been taken by Morris17 and Hazra18

who use skew-t processes and their mixtures, which possess strong tail dependence, to model ozone and fire risk extremes respectively.

Since extreme precipitation stands to do the most damage, accurately capturing the spatial dependence for high intensities necessitates special consideration. The last decade has produced many new methods for modeling spatial extremes. Two common approaches to modeling spatial extremes are based on liming results: max-stable processes see Davison12 for a review, which are appropriate for component-wise maxima over large temporal blocks, and generalized Pareto processes for high threshold exceedances (Ferreira14). However, full-likelihood inference for these models suffers from computational intractability when the number of spatial observation locations is large. While composite likelihood approaches have been proposed to remedy this issue (Padoan10), the scalability of full-likelihood methods remains an open problem. In addition to limiting models, several Bayesian hierarchical models that borrow ideas from classical geostatistics have been developed Cooley07; Sang09; Sang10; Reich12. Unlike in the classical setup to modeling extremes, which requires a somewhat arbitrary qualification of an extreme event (e.g. the threshold in a peaks-over-threshold model or block size in a model for block maxima), our proposed mixture model incorporates the entire distribution of data, while accommodating different dependence types for different storm intensities.

In addition to accurately representing the spatial dependence, capturing temporal dependence has also been central to the field of precipitation modeling. Several Hidden Markov models for unobserved weather states have been proposed that aim to capture temporal dependence at various scales, including those attributable to large, synoptic-scale weather patterns

(Bellone00; Ailliot09; Flecher10)

. Latent multivariate autoregressive models are also common approaches to modeling temporal dependence for both occurrence and intensity processes

(Bardossy92; Makhnin09; Rasmussen13). Methods based in physics have also appeared in the statistical literature; recently, Liu18 proposed a Lagrangian advection reference frames approach to modeling storm dynamics that couples radar reflectivity and wind field data to make short-term precipitation forecasts.

In the remainder of the paper, we take up a similar aim to that of Gao14 and Gao17 who use analogue methods to couple GCM forecasts of predictive atmospheric variables with the historical precipitation so as to understand the changes in the distributions of extreme precipitation under different climate forcing scenarios. Our method is different from these earlier efforts in that it is a hierarchical Bayesian model based approach that makes use of the full data likelihood and can easily account for multiple sources of uncertainty. The following sections are outlined as follows: in Section 2, we develop a Bayesian models for precipitation occurrences and intensities that capture temporal dependence with a probabilistic analogue formulation. Finally, in Section 3, we apply the model to precipitation data over the northeastern United States, and compare it with climate model and reanalysis distributional forecasts of extreme precipitation.

2 Model Definitions

One of the main challenges of modeling precipitation fields is that their distribution consists of a mixture of a preponderance of zeros and positive precipitation amounts. To account for this, the proposed spatio-temporal model for precipitation is made up of two parts: (1) an occurrence model for the presence versus absence of precipitation, and (2) an intensity model for positive precipitation amounts (Wilks90; Berrocal08). We begin by describing the occurrence model for precipitation.

2.1 Occurrence Model

A common approach to modeling spatially varying, binary data is via data augmentation, wherein a continuous latent process is thresholded into two categories (Albert93; Heagerty98; Collett02). To model the point referenced, binary occurrence of precipitation, we use a Gaussian process for the unobserved, latent component, such that it is positive at a location whenever there is precipitation at and negative otherwise. This commonly used model is referred to as a clipped Gaussian process or spatial probit model (DeOliveira00; DeOliveira18).

The spatio-temporal occurrence process , , consists of spatial random fields that encode the presence versus absence of precipitation (1: presence, 0: absence) at a location and time . It is defined in terms of a zero-thresholded, latent spatio-temporal Gaussian process :

The latent processes (by abuse of notation) are parameterized by a mean function and covariance function that induce spatial dependence in the occurrence process to reflect the fact that nearby locations are more likely to have common presence or absence of rain than distant ones. Conditional on the mean and covariance functions, the spatially dependent processes are assumed to be independent in time, each distributed as:

Temporal dependence is contained in the mean function. The covariance function can be expressed as a product of a variance parameter

and correlation function as . However, since the variance term of this model is unidentifiable (DeOliveira18), it is fixed at , making it sufficient focus on the spatial correlation function. Due to the flexibility provided by a parameter governing the smoothness of the process, we assume a correlation function from the general Matérn class (Stein99), which is both stationary and isotropic, making it expressible as a function of distance between locations : for range and smoothness parameters, where is a modified Bessel function of the second kind.

Using this construction, the marginal probability of occurrence at time

and site can be expressed in terms of the mean function of the latent process as , where

is a standard normal distribution function. To allow for spatially varying marginal occurrence probability, the mean function is further modeled as a linear combination of

spatial basis functions . For the spatial basis functions, we use Gaussian kernels although other choices are viable. For generic knot locations , the basis functions are defined as , where denotes a standard Gaussian density function and

is a bandwidth parameter. Denoting the vector of basis functions at location

as , we model the mean function at time as a sum of an offset term and a spatially varying term as . For the vector of spatial basis coefficients , we assume an independent normal prior .

Since the presence or absence of precipitation is determined by whether is positive or negative, the offset term can be thought of as governing the overall (non-spatially varying) probability of precipitation on day . We will use this term to capture the temporal dependence in occurrence of precipitation by leveraging synoptic-scale climate forcings (e.g. geopotential height and temperature) that are concomitant with precipitation by construction of a prior on . We defer discussion of the analogue prior on until Section 2.3 as it is also used in the model for precipitation intensities. The priors for the Matérn dependence parameters are taken to be and , where in subsequent sections and are taken to be the minimum and maximum distance between observation locations.

2.2 Intensity Model

In this section, we develop a Bayesian model for the positive spatial precipitation intensity process. The intensity process is treated as a continuous process defined on the entire spatial domain, , that is masked by the occurrence process. In other words, the intensity process is only observed at locations where the occurrence process is positive. We make use of the Student-t process because of its flexibility in capturing heavy tailed behavior, which is commonly observed in precipitation data (Vrac07; Naveau16).

Denote the precipitation intensity process by for times . While precipitation intensities are strictly positive, it is much more convenient to work with a spatial process defined on the whole real line. The softplus function , maps and is strictly increasing, preserving the natural ordering of the data. Both the softplus function and its inverse leave moderate to large values effectively unchanged. In the remainder of this section we define a model for the transformed precipitation intensities, .

A location-zero Student-t process can be expressed as a Gaussian process scale mixture see e.g., Shah14:


denotes a Inverse-Gamma distribution with shape

and scale . After marginalizing over ,

is a Student-t process with degrees of freedom

and scale . A Gaussian process is a limiting case of a Student-t process as .

To allow flexible spatial dependence types across different precipitation intensity levels, the transformed precipitation amounts are modeled as a finite mixture of Student-t processes. Finite mixtures can easily be accommodated via data augmentation. Let denote the latent mixture label for time , where is the total number of mixture components. A multinomial logistic model, also referred to as the Luce model in the probabilistic-choice econometrics literature when modeling latent classes (Luce59; Mcfadden73), is used for the mixture class membership. Denote covariates (e.g., containing synoptic scale atmospheric information) by for times , vectors of coefficients for each class by . Then, for linear predictor , the probability that the process at time belongs to mixture component class is modeled as:


For identifiability, the loadings for the class are fixed while the remaining loadings have independent normal priors:

Instead of using a common spatial dependence structure, each mixture class is free to have different spatial covariance, degrees of freedom , and scale parameters. Any spatial variation in the location surface at time is captured by an offset term and a linear expansion of spatial basis functions , just as was done in Section 2.1. Conditional on the mixture label at time , , the intensity process is modeled as follows:

For each mixture class, we assume an isotropic Matérn correlation function with potentially different smoothness and range parameters:

Just as in Section 2.1, we assume independent normal priors for the spatial basis coefficients , and we use an analogue prior for , which we discuss in Section 2.3. The following priors are used for the degrees of freedom and scale parameters: , and , where the Gamma distribution is parameterized with shape and scale respectively. Priors for the Matérn covariance parameters and are assumed to be indpendent across and the same as those in the occurrence model.

Because we use a mixture of Student-t processes, different mixture components are capable of capturing different spatial dependence types, allowing for varying dependence strengths across different intensities. The upper tail dependence of a generic stationary and isotropic spatial process can be characterized by the tail dependence function at level , defined as:


where denotes the marginal distribution function of . The limit determines the tail dependence class. The process is said to be asymptotically independent at distance is , and asymptotically dependent at distance otherwise. Gaussian processes are asymptotically independent processes with for all , making them suitable for modeling physical phenomena exhibiting weak tail dependence (Sibuya60). However, Student-t processes and their finite mixtures possess asymptotic dependence making them suitable in the case of strong tail dependence, but less so when the data exhibit weakening dependence at higher levels (Nikoloulopoulos09). Because a Gaussian process is a limiting case of the Student-t process, a Student-t process with large degrees of freedom is capable of capturing relatively weak dependence at sub-asymptotic levels.

2.3 Analogue Prior

Both models described in Sections 2.1 and 2.2 possess location offset terms and , which here we will denote generically as . The priors for the offset parameters have thus far been left unspecified. In this section, we describe how the similarities between climate forcings at different times can be used to model the temporal dependence in the offset term. Following McDermott16 and Mcdermott18, who take a similar approach to forecasting soil moisture and waterfowl settling behavior, we will refer to this as an analogue prior on .

The presence and intensity of precipitation is closely related to other atmospheric conditions such as atmospheric pressure, temperature, and water vapor. Rather than explicitly model the complex precipitation dynamics, we use closely related atmospheric variables to identify historical analogues (times ) of the atmospheric conditions at time . Because of concomitance of atmospheric conditions with precipitation, the precipitation conditions during identified historical analogue times can then be used to inform the precipitation for the reference time. In a classical analogue model (Barnett78; Sugihara90), identified historical analogue precipitation fields and their weighted averages would be used as the forecasts for the reference time. Instead, we make a weaker assumption, and only borrow information about the location offset terms from analogue precipitation fields.

To formalize this, for each time , define a vector of weights , quantifying the similarity between the atmospheric conditions at time and those at time . For identifiability, the elements of the weight vectors are restricted to , , and , for all . Combining the location parameters for other times into a single vector , we specify conditional normal priors on the location parameters , so that the prior mean is a weighted average of the precipitation location parameters from historical conditions with strong similarity to the reference conditions.

The weights are calculated using a kernel function applied to distances atmospheric conditions at different times, such that more similar historical atmospheric conditions receive larger weights. , we use an unnormalized, compact Gaussian kernel function: , where is a measure of distance, is a kernel bandwidth parameter, and is a threshold ensuring that the kernel is compact, which prevents irrelevant non-analogue days from receiving positive weight. In subsequent sections it is fixed such that on average (across times ) only the top nearest analogues receive weight. When the kernel bandwidth is large, many historical analogues will contribute a small weight, and when the bandwidth is small, fewer historical analogues will contribute a larger weight. Given distances between all times and , unnormalized weights are calculated as , and normalized to give weights .

Due to the high-dimension of geopotential height fields and temperature fields, some dimension reduction is needed before calculating distances between atmospheric covariates. To reduce the dimension while preserving the relevant variation in the data, each of the atmospheric variables are projected onto the first several components of the empirical orthogonal functions (EOFs) calculated from the data (Jolliffe02; Hannachi07; Demsar13). Rather than rely on snapshots of the atmospheric conditions on any given day, we employ the time lagging approach of Takens81. Combining lagged EOF loadings into trajectories of atmospheric conditions over time gives a fuller reconstruction of the state-space of the dynamical system (Sugihara90). Let denote the loadings for the first several components of relevant EOFs on day , we construct a new embedding matrix as , where is the number of lagged time steps. Doing this for all time points, the quality of as an analogue for is determined by the distance (e.g., Euclidean norm) between and .

For both models, posterior samples are made using Markov chain Monte Carlo (MCMC), the details of which can be found in Appendix


. Gibbs updates are available for some parameters, and the remaining are made with Metropolis-Hastings updates. Samples of the parameters are made for all observation times, and in-sample posterior predictive draws can be made directly, e.g.

for and non-observation location . However, to make out-of-sample posterior predictive draws, for example of for time , distances between atmospheric conditions on the future day and historical days must first be calculated, e.g. as well as their corresponding normalized weights from for , conditional on posterior samples of .

To assess the utility of the model two simulation studies are performed, the results of which are summarized in the Supplementary Material.

3 Precipitation Analysis

In this section, we apply our model to daily precipitation observed over the Susquehanna river basin in southern New York and Pennsylvania. By coupling daily precipitation accumulations observed at rain gauge stations with predictive atmospheric conditions, we can apply the proposed analogue model to make forecasts of precipitation over watersheds. The rain gauge data come from the National Oceanic and Atmospheric Administration (NOAA) ( and consist of daily precipitation accumulations (in mm) observed between 1986 and 2017 for gauge stations. To capture the temporal dependence, meteorological reanalysis estimates of geopotential height and surface temperature are used to construct the distance matrix for the analogue prior on location terms in the occurrence and intensity models. Reanalysis computer models infill meteorological fields on a spatial grid by assimilating historical, spatially varying atmospheric observations, which are treated as boundary conditions in a consistent model of the climate system. We consider 500 hPa geopotential height and surface temperature from the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) project (Gelaro17) for identifying analogues. The analogue prior distance matrix is constructed using the lagged EOF loading approach described in Section 2.3 with lagged time steps and the first 10 PCs of each of the two MERRA-2 atmospheric variables. An example of the 500 hPA and surface temperature fields for reference and nearest analogue days are shown in Figure 1.

Figure 1: The 500 hPa geopotential heights and surface temperatures for a reference day are show in the leftmost column. The middle and right columns show the corresponding fields for the top two most similar days based on the Euclidean distance between lagged vectors of the top 10 PC loadings of both variables.

We evaluate three models for daily precipitation: occurrence and intensity models with (M1) independence priors on and , (M2) analogue priors on and with distance matrix calculated using MERRA-2 atmospheric variables over the continental United States (CONUS), and (M3) analogue priors on and

with distance matrix calculated using MERRA-2 atmospheric variables over the a region surrounding Pennsylvania (local PA). The data is split into training (1986-2000) and holdout periods (2001-2017), during each year of which we consider only the months April, May, and June. During the training period, 20 stations are also held out for model evaluation of in-sample prediction (interpolation). Inference on each is performed with a Metropolis within Gibbs Markov chain Monte Carlo (MCMC) algorithm run for

iterations. The models are evaluated by two measures: by their ability to capture the distribution of precipitation when making (1) out-of-sample forecasts and (2) in-sample, kriging predictions. To evaluate (1), out-of-sample posterior predictive draws are made from each model during holdout period at each of the 174 gauge locations, and to evaluate (2), in-sample posterior predictive draws are made at holdout stations for each day of the training period. In addition to these three models, we also consider the distributional forecasts of precipitation from the LENS historical run during the holdout period and MERRA-2 estimates during the training period. LENS fields of 500 hPa geopotential height and surface temperature are projected onto the MERRA-2 EOF basis for comparability of fields when calculating distances. Since reanalysis estimates effectively condition on observed the atmospheric boundary conditions to infill precipitation, we also compare the in-sample, kriging predictions of the proposed model to MERRA-2 reanalysis estimates of precipitation. Both MERRA-2 and LENS precipitation forecasts are made on a grid. To resolve the spatial mismatch between grid cells and gauge locations, gauge station forecasts for the LENS and MERRA-2 models are made by assigning each gauge station the precipitation amount from its nearest LENS and MERRA-2 grid cell.

The models are evaluated based on how well they capture the distributional characteristics of extreme precipitation. For each day of the holdout period, the maximum daily precipitation accumulation across all stations is calculated. The distribution of maxima is then compared to the predictions from each model using tail weighted continuously ranked probability scores (TWCRPS) (Gneiting14). For a single sample , the TWCRPS is calculated as:


where refers to the model estimate of the target distribution from which is drawn, and is a weight function, which for this analysis we take to be where is the quantile of the simulated total precipitation data. In practice, as we do here, the TWCRPS is averaged over a sample as (where for in-sample prediction and is the number of holdout days for out-of-sample prediction). Lower TWCRPS scores correspond to better correspondence between the empirical and model-based distributions. We consider two TWCRPS weight functions: and where and

are the sample median and standard deviation of the observed precipitation. The TWCRPS results are summarized in Table

1. Both versions of the analogue prior model outperform the independence prior model in capturing the distribution of precipitation extremes during the holdout period. Interestingly, the model with analogue distance matrix constructed from locally defined EOFs around PA shows greater predictive skill than both the one constructed from CONUS EOFs and the LENS model forecasts. Moreover, the in-sample forecasts of precipitation using the Local PA analogue prior model are also competitive with the MERRA-2 model. As such, we focus the remainder of our analysis on the analogue prior model using local PA fields.

Predictions TW Fun. Indep. Prior CONUS Local PA LENS MERRA-2
In-sample 7.43 (7.35, 7.55) 7.41 (7.32, 7.54) 7.37 (7.28, 7.48) 7.30 (-,-)
In-sample 6.66 (6.59, 6.77) 6.64 (6.56, 6.77) 6.61 (6.53, 6.70) 6.54 (-,-)
Out-of-sample 20.9 (13.1, 25.2) 16.3 (15.2, 17.4) 11.8 (11.2, 12.5) 13.7 (-,-)
Out-of-sample 19.6 (12.0, 23.7) 15.0 (14.0, 16.1) 10.7 (10.2, 11.4) 12.7 (-,-)
Table 1:

Tail weighted continuously ranked probability scores (TWCRPS) for the Student-t process mixture, LENS model forecasts of the distribution of daily maximum precipitation among all rain gauge locations during the holdout period (2001-2017), and MERRA-2 model forecasts of the distribution of daily maximum precipitation among holdout gauge locations during the training period (1986-2000). The better value for each criteria is given in bold. Estimated standard errors are given in parentheses.

A map of the pointwise, marginal percentile estimates of the posterior predictive distribution of daily precipitation over the study region are given in Figure 2. Estimated precipitation levels tend to be highest over the northeastern part of the state as well as over parts of New Jersey and New York. To assess the correspondence between the empirical and model estimates of the spatial pattern of tail dependence, we also examine the F-madogram for monthly maxima Cooley06

. The F-madogram is analogous to the more traditional madogram from classical geostatistics, but is guaranteed to exist even when the first moment of the process under consideration is undefined. For a generic, stationary and isotropic spatial process

, , with marginal distribution function , and spatial lag , the F-madogram, defined as


and describes the dependence in the pairs that are a distance apart. When the pair are perfectly dependent, , and when the pair are independent, . The figure shows both good correspondence between the empirical and model estimates of as well as apparent non-zero dependence of monthly maxima even at spatial lags of km, which suggests that the model is capturing the spatial dependence properties of extreme precipitation well.

For illustration of the model based forecasts, the observed daily precipitation accumulations, a posterior predictive draw, and posterior predictive mean and standard deviation for a single day are shown in Figure 3. Observed and predicted zero accumulations are shown as gray. The figure shows good correspondence of the general spatial surface and smoothness characteristics in the observed and predicted precipitation amounts.

Figure 2: Left: Pointwise estimates of the

percentile of daily precipitation over the observation region based on the fitted analogue prior model overlaid with upper (0501), western (0502), and lower (0503) branch Susquehanna drainage basin boundaries. Right: Empirical F-madogram estimates (points) and 95% credible intervals (ribbons) for monthly maxima of daily precipitation accumulations. A dashed horizontal line corresponding to independence is plotted for reference.

Under the assumption that the library of observed historical climate states is sufficiently rich to match future states, our model can be used to make forecasts of future daily precipitation accumulations under various climate forcing scenarios. To do this, the 500 hPa geopotential height and surface temperature fields of the LENS RCP8.5 run are projected onto the first 10 principal components of their respective MERRA reanalysis fields. The distances between lagged trajectories of LENS EOF loadings and historical MERRA EOF loadings are calculated for all future days (e.g. for some future day , the distance vector consists of Euclidean distances between covariates on day and all historical days ). LENS RCP8.5 forecasts of geopotential height and surface temperature are coupled with the fitted analogue model to make out-of-sample posterior predictive draws of precipitation on a grid over three Susquehanna drainage basins from 2006-2100. Since the proposed model is defined on a continuous domain, predictions can be made on an arbitrarily fine grid. Summaries of the 3-month (April-May-June) maximum daily total precipitation volumes over basins are shown in Figure 4. The figure shows high uncertainty in the estimates, with 95% credible bands capturing slight increasing trends over time for all three basins. To quantify the flood-risk, the forecasts from this model could be input into a hydrological flow model that accounts for topography and land use among other factors.

Figure 3: Observed daily precipitation amounts at gauge locations (top left), posterior predictive draw (top right), posterior predictive mean (bottom left), and posterior predictive standard deviation (bottom right) for a single day during the observation period. Grey points and cells correspond to zero precipitation amounts. The general spatial mean surface as well as degree of smoothness in the observed precipitation data is well captured by that of the posterior predictive sample.
Figure 4: Global Historical Climate Network (GHCN) daily gauge station locations over Pennsylvania and surrounding states are overlaid with upper (0501), western (0502), and lower (0503) branch Susquehanna drainage basin boundaries (left). Forecasted yearly (during April-May-June) maximum total precipitation 95% credible intervals for each basin are plotted on the right.

4 Discussion

While analogue methods have a long history in meteorology, there have been relatively few attempts at using them in a fully probabilistic framework for precipitation modeling. The analogue prior is a very general approach to modeling temporal dependence that leverages climate model forecasts of atmospheric variables that are concomitant with precipitation. Since this model is developed in a Bayesian framework, it is possible to account for uncertainty in the predictive skill of concomitant atmospheric variables that are used to select analogues. Moreover, the Student-t process mixture is a flexible model that can accommodate a wide variety of spatial dependence types in both the bulk and tail of the distribution. This model could be extended further by considering non-linear dimension reduction techniques for identifying analogues such as Laplacian eigenmaps (Belkin02)

, self organizing maps

(Kohonen12), and diffusion maps (Coifman06). Similar to other analogue methods, the ability of the proposed model to produce accurate forecasts is predicated on the condition that there exist close historical analogues to future conditions see Vandendool94 for a discussion. As part of our ongoing and future work, we plan to apply our precipitation forecasts to a hydrological water-flow model to more directly interrogate changes in flood-risk over the Susquehanna drainage basins.

Appendix A MCMC Details

Metropolis-Hastings MCMC algorithms were implemented for making posterior draws of the parameters in both the occurrence and intensity models using the R programing language (

a.1 Occurrence Model Metropolis-Hastings Algorithm

The parameters , , and were updated using variable-at-a-time random walks. Truncated normal Gibbs updates are available for the latent Gaussian process at observation locations . Denoting occurrence indicators and mean function , the Gibbs updates for the latent Gaussian process are

where is the spatial covariance matrix for the observation locations, lower bounds have elements

and upper bounds have elements

Gibbs updates are also available for both the location parameters and . For the location parameter, the Gibbs update is

For the basis coefficients , the Gibbs update is

Inverse Gamma gibbs updates are used for and . Using inverse Gamma parameterized with shape and scale , having density . With prior , the Gibbs update for , is

and the Gibbs update for , assuming prior is

a.2 Intensity Model Metropolis-Hastings Algorithm

The parameters for , and were updated using variable-at-a-time random walks. The cluster labels were also updated variable-at-a-time, but with discrete uniform independence proposals, each on 1:. The location offset and basis functions Gibbs updates are similar to those in the occurrence model but are also dependent on the mixture labels.

where the covariance matrix is calculated using dependence parameters and corresponding to mixture class .

Similarly, the basis coefficients have Gibbs update


The Gibbs updates for the prior variances and are completely analogous to those in the occurrence model.