Climate models robustly imply that some significant change in precipitation patterns will occur. Models consistently project that the intensity of individual precipitation events increases by approximately 6%–7% K−1, following the increase in atmospheric water content, but that total precipitation increases by a lesser amount (1%–2% K−1 in the global average in transient runs). Some other aspect of precipitation events must then change to compensate for this difference. The authors develop a new methodology for identifying individual rainstorms and studying their physical characteristics—including starting location, intensity, spatial extent, duration, and trajectory—that allows identifying that compensating mechanism. This technique is applied to precipitation over the contiguous United States from both radar-based data products and high-resolution model runs simulating 80 years of business-as-usual warming. In the model study the dominant compensating mechanism is a reduction of storm size. In summer, rainstorms become more intense but smaller; in winter, rainstorm shrinkage still dominates, but storms also become less numerous and shorter duration. These results imply that flood impacts from climate change will be less severe than would be expected from changes in precipitation intensity alone. However, these projected changes are smaller than model–observation biases, implying that the best means of incorporating them into impact assessments is via “data-driven simulations” that apply model-projected changes to observational data. The authors therefore develop a simulation algorithm that statistically describes model changes in precipitation characteristics and adjusts data accordingly, and they show that, especially for summertime precipitation, it outperforms simulation approaches that do not include spatial information.
Some of the most impactful effects of human-induced climate change may be changes in precipitation patterns (IPCC 2013). Changes in future flood and drought events and water supply may incur social costs and mandate changes in management practices (e.g., Rosenzweig and Parry 1994; Vörösmarty et al. 2000; Christensen et al. 2004; Barnett et al. 2008; Karl 2009; Nelson et al. 2009; Piao et al. 2010). Understanding and projecting changes in future precipitation patterns is therefore important for informed assessment of climate change impacts and design of adaptation strategies.
Future changes in spatiotemporal precipitation patterns are not well studied, but climate simulations under higher CO2 consistently show divergent increases in precipitation amount and intensity that, if true, imply that some changes must occur. Models robustly project uniform increases in precipitation intensity (amount per unit time and unit area when rain occurs) of approximately 6% K−1 temperature rise, following atmospheric water content governed by Clausius–Clapeyron (e.g., Held and Soden 2006; Willett et al. 2007; Stephens and Ellis 2008; Wang and Dickinson 2012). However, model total precipitation rates rise differently: in transient runs, in the global average, by only 1%–2% K−1 (e.g., Knutson and Manabe 1995; Allen and Ingram 2002; Held and Soden 2006; Stephens and Ellis 2008; Wang and Dickinson 2012). (Time-averaged precipitation rise can exceed Clausius–Clapeyron only in the deep tropics.) Nearly all latitudes then show a discrepancy between changes in rainfall intensity and total amount that must be “compensated” by some other change in precipitation characteristics (e.g., Hennessy et al. 1997). In the midlatitudes, the amount and intensity discrepancy generally resembles the global average (Fig. 1). Model midlatitude rain events must therefore experience changes in frequency (fewer storms), duration (shorter storms), or size (smaller storms).
Current approaches to studying future precipitation cannot determine changes in rainstorm characteristics because they do not consider individual events. Studies generally analyze precipitation at individual model grid cells (Knutson and Manabe 1995; Semenov and Bengtsson 2002; Stephens and Ellis 2008) or over broad regions (Tebaldi et al. 2004; Giorgi and Bi 2005). Effects from changes in rainstorm frequency, duration, and size are all confounded in local or spatially aggregated time series. To overcome this limitation, we need an approach that identifies, tracks, and analyzes individual rainstorms. Currently used rainstorm tracking algorithms are not appropriate for climate studies: they are designed for severe convective storms in the context of nowcasting and forecasting using weather forecasting models (Davis et al. 2006b), radar (e.g., Dixon and Wiener 1993; Johnson et al. 1998; Wilson et al. 1998; Fox and Wikle 2005; Xu et al. 2005; Davis et al. 2006a; Han et al. 2009; Lakshmanan et al. 2009), or satellite images (Morel and Senesi 2002a,b). Storm studies in climate projections have also focused only on very large and intense events (Hodges 1994; Assunção et al. 2012). These algorithms cannot efficiently handle lower-intensity precipitation events with more complicated morphological features and evolution patterns. Studying future precipitation patterns requires new storm identification and tracking strategies.
Such studies also require high-resolution model output. For a model to plausibly represent changes in real-world rain events it must at minimum explicitly represent those events. Computational constraints mean that typical general circulation model (GCM) runs used for climate projections have spatial resolutions on the order of 100 km, too coarse to represent morphological features of localized rainstorms. Fine-grained observational data can be generated from radar datasets; studies of future precipitation require dedicated model runs at similar resolution, at the price of more limited time series.
Finally, making practical use of insights into changes in rainstorm characteristics requires methods for combining model projections with observational data. Hydrological and agricultural impact assessments cannot use scenarios of future precipitation from even high-resolution models, since model precipitation can differ considerably from that in observations (e.g., Ines and Hansen 2006; Baigorria et al. 2007; Teutschbein and Seibert 2012; Muerth et al. 2013). The two primary current approaches to addressing these biases are bias-correcting model output based on observations (of means or marginal distributions) (e.g., Ines and Hansen 2006; Christensen et al. 2008; Piani et al. 2010a,b; Teutschbein and Seibert 2012) and “delta” methods that adjust observations by model-projected changes (in means or marginal distributions) (e.g., Hay et al. 2000; Räisänen and Räty 2013; Räty et al. 2014). Neither approach allows representing changes in spatiotemporal dependence. Some work has extended bias correction methods to also address biases in spatiotemporal dependence (Vrac and Friederichs 2015), but these methods do not allow for future changes in that dependence. An important objective of this work is therefore to extend delta methods to capture model-projected changes in rainstorm characteristics, while still ensuring the greatest fidelity to real-world precipitation statistics.
In this work we seek both to understand the changes in rainstorm characteristics in future model projections that compensate for the amount-intensity discrepancy and to develop an approach to transform observed rainstorms into future simulations to account for those changes. The remainder of this paper is organized as follows. In section 2 we describe the high-resolution regional climate model runs and the radar-based observational data products used in the study. In section 3 we develop an algorithm for identifying and tracking individual rainstorms, and in section 4 we develop metrics for analyzing spatiotemporal precipitation patterns. In section 5 we compare precipitation in observations and present-day and future climate projections, and in section 6 we develop a simulation method to generate future precipitation simulations that combines information from models and observations. Finally, in section 7 we discuss the implications of these results.
2. Climate model output and radar-based observational data
As discussed above, model output used to study rainstorm characteristics must be at high spatial resolution. In this study we use high-resolution dynamically downscaled model runs over the continental U.S. region, with a constant-spacing grid of 12 km × 12 km, described by Wang and Kotamarthi (2015). These runs use the Weather Forecasting and Research (WRF) Model (Skamarock and Klemp 2008) as the high-resolution regional climate model, nudged by a coarser simulation from the Community Climate System Model, version 4 (CCSM4; Gent et al. 2011), of the business-as-usual scenario [representative concentration pathway 8.5 (RCP8.5)]. [The CCSM4 run is an ensemble member from the CMIP5 archive; see Meinshausen et al. (2011) and WCRP (2010) for further information.] In the WRF simulation we apply spectral nudging, which is conducted in the interior as well as at the lateral boundaries, on horizontal winds, temperature, and geopotential height above 850 hPa. Because high-resolution runs are computationally demanding, WRF was run only for two 10-yr segments of the scenario, which we term the baseline (1995–2004) and future (2085–94) time periods. In this analysis we separately analyze summertime (June–August) and wintertime (December–February) precipitation because precipitation characteristics differ by season. U.S. precipitation is predominantly convective in summer and large-scale in winter (Fig. 2).
Both CCSM4 and WRF are widely used models for atmospheric and climate science. WRF has been extensively used for mesoscale convection studies, for dynamical downscaling of climate model projections (Lo et al. 2008; Bukovsky and Karoly 2009; Wang et al. 2015), and as a forecasting model for numerical weather prediction (e.g., Jankov et al. 2005; Davis et al. 2006b; Clark et al. 2009; Skamarock and Klemp 2008). Validation studies addressing these different contexts include that of Ma et al. (2015), who showed in a model–observation comparison that downscaling CCSM4 output using WRF improved the fidelity of modeled summer precipitation over China on both seasonal and subseasonal time scales. For forecasting, Davis et al. (2006b) studied intense and long-lived storms in 22-km resolution WRF forecasts over the United States and found that the model captured storm initialization reasonably well but showed some apparent biases in storm size, intensity, and duration.
We compare WRF and CCSM4 model output to an observational data product of similar spatial and temporal resolution: the NCEP Stage IV analysis (Stage IV data hereafter), which is based on combined radar and gauge data (Lin and Mitchell 2005; Prat and Nelson 2015). The Stage IV dataset provides hourly and 24-hourly precipitation at 4-km resolution (constant-spacing 4 km × 4 km grid cells) over the contiguous United States from 2002 to the present. (To match our model output, we use 10 years of data from 2002 to 2011 and aggregate the 1-hourly data to 3-hourly.) Stage IV data are produced by “mosaicing” data from different regions into a unified gridded dataset. This mosaicing leaves artifacts in the spatial pattern of average precipitation: the range edges of individual radar stations are visible by eye in Fig. 3. These artifacts do not compromise our analysis of rainstorm frequency, duration, and size, as mosaicing effects do not affect temporal evolution. Stage IV data have been shown to be temporally well correlated with high-quality measurements from individual gauges (see, e.g., Sapiano and Arkin 2009; Prat and Nelson 2015).
Comparing the present-day WRF run with Stage IV data shows that the model captures seasonal mean precipitation amounts well (Table 1, first row), but with a large bias in rainfall intensity. Rain rates in precipitating grid cells in WRF output are approximately 50% lower than in observations (Fig. 4). The bias is consistent across two orders of magnitude in intensity, nearly identical in both summer and winter, and unlikely to be a sampling artifact, since both model and observation gridcell sizes (144 and 16 km2) are well below typical areal sizes of precipitating events. Because the model matches observed total rainfall, the too-low intensity in individual rainstorms must be compensated by some other bias in precipitation characteristics: model precipitation events may be more frequent, more numerous, larger in size, or any combination of these factors. [Davis et al. (2006b) have noted a bias in WRF toward too-large rainstorm events, although in a study of only the most severe storms.] Understanding model–observation biases therefore becomes an analogous problem to understanding differences in present and future model projections. Both cases involve discrepancies in amount and intensity, which imply some change in spatiotemporal precipitation properties that can be understood only by identifying and tracking individual rainstorms.
Intensity differences do introduce one minor complication in analysis. It is typical in statistical analyses of precipitation to cut off all data with precipitation intensity below a particular threshold in order to remove spurious light precipitation and make analysis tractable. In our case the intensity distributions for model and observation are different, suggesting that different cutoffs might be warranted. For simplicity, we apply the same cutoff of 0.033 mm h−1 (black line in Fig. 4) to both model and observations. This choice should negligibly affect the overall results, since even in the worst case (model wintertime) the cutoff excludes less than 2% of seasonal precipitation.
3. Identifying and tracking individual rainstorms
As discussed previously, existing algorithms for finding precipitation events are appropriate only for severe storms. To decompose the entire precipitation field into a set of events, we develop here a more general approach. We define a rainstorm for this purpose as a set of precipitating grid cells that are close in location and that move together retaining morphological consistency. As in Davis et al. (2006a), we use only the precipitation field itself, allowing ready comparison to observations. (More complex definitions might also draw on information about the meteorological context.)
Defining rainstorm events in spatiotemporal data involves two tasks:
Rainstorm identification: Divide the precipitation field at each time step into individual storms.
Rainstorm tracking: Build rainstorm events evolving over time by tracking identified storms across consecutive time steps.
We describe each step below, and provide more detail on algorithms in section S1 in the supplemental material.
a. Identifying rainstorms at a single time step
We first identify rainstorms from the precipitation field at each time step t. The simplest approach would be to treat each contiguous region of grid cells with positive precipitation exceeding some threshold criterion as an individual rainstorm (Guinard et al. 2015). This algorithm, however, identifies excessive individual rainstorm events that are in reality meteorologically related: a single storm system can often have multiple separate precipitating regions (Fig. 5a, left). Rainstorm identification is properly a clustering problem that requires grouping-related but not necessarily contiguous phenomena. We group nearby regions into “rainstorm segments” using a technique termed “almost-connected component labeling” (Eddins 2010; Murthy et al. 2015). The basic idea is to identify precipitating regions close enough they would connect given some prescribed dilation of individual regions. This procedure provides a natural way of grouping regions based on their proximity and morphological features and has been previously used for rainstorms (Baldwin et al. 2005). Unlike other traditional clustering algorithms such as k-means clustering (Hartigan and Wong 1979), the approach does not require a prespecified number of clusters and hence enables quick and automatic clustering.
Almost-connected component labeling does often suffer from a “chaining effect” that produces overly large rainstorm segments: a small number of precipitating grid cells located between two large rain regions can cause awkward linkage between them (Fig. 5a, right). To avoid this, we use a four-stage procedure based on mathematical morphology that treats “large” and “small” regions separately. [A similar approach was used on radar reflectivity datasets by Han et al. (2009).] We first identify contiguous regions of precipitation [Fig. 6, top panel; our threshold is >0.1 mm (3 h)−1]. In the second stage, we form segments using only the large regions (Fig. 6, upper-middle panel). In the third stage we assign small regions to the closest existing segments from the previous stage if they are close to the existing segments (Fig. 6, lower-middle panel). In the fourth stage we form segments with the remaining small regions (Fig. 6, bottom panel).
b. Tracking rainstorms over different time steps
Once the rainstorm segments for all time points are identified, we link them in different time steps to form rainstorm events evolving over time and space. This task is complicated by the fact that rainstorms often split into multiple segments that drift in different directions (see Fig. 5b for illustration). Most existing algorithms are not designed to handle this possibility. Our new algorithm is partly inspired by the work of Hodges (1994) and Morel and Senesi (2002a), but it extends their work to more efficiently identify segments originating from the same rainstorm and track their movement over time.
The algorithm works sequentially in time. For each time step from t = 1, …, T, we assign each rainstorm segment to one of the existing rainstorm events based on the two criteria: 1) the shapes of linked rainstorms in two consecutive time steps are morphologically similar enough to be considered as the same rainstorm, and 2) the rainstorm location and the movement direction do not change too abruptly over time. If we cannot find any existing rainstorm events that meet these criteria for a rainstorm segment, we initialize a new rainstorm event starting with that segment. On the other hand, if more than two events satisfy the criteria, we assign the rainstorm segment to the largest event. Since the multiple rainstorm segments at the same time step can be assigned to the same rainstorm event, our algorithm can incorporate situations where a rainstorm splits into multiple segments. We denote the resulting rainstorm events S1, …, SN. Figure 7 shows example time steps for rainstorm tracking results from our algorithm, which appears to simultaneously track different rainstorms with various morphological features reasonably well.
4. Describing rainstorm characteristics
Once rainstorms are identified and tracked, the goal is to quantitatively describe them. There is little precedent in the literature for this task; most previous rainstorm studies focus only on analyzing storm trajectories or tracks (Hodges 1994; Morel and Senesi 2002b). We therefore develop a set of metrics for characterizing storms, and follow a four-step process for describing spatiotemporal precipitation patterns and comparing those patterns between different climate states (or between model output and observations):
Compute metrics for individual rainstorms: duration, size, mean intensity, and central location.
Extend those metrics to apply to aggregate precipitation. That is, we decompose the total precipitation amount into the product of average intensity, a size factor, a duration factor, and the number of rainstorms.
Estimate the geographic variation of those aggregate metrics. That is, we find spatial patterns of rainstorm properties.
Compare precipitation patterns across model runs or across models and observations.
We describe the algorithms for each of these steps below.
a. Metrics for individual rainstorms
We characterize each individual rainstorm event Si with four metrics: duration, size, mean intensity, and central location. For completeness, we describe below how each metric is computed. Note that size, location, and amount metrics are not scalars but vectors or matrices over the lifetime of each storm.
Duration: The storm identifying algorithm gives us the beginning and ending time steps of the lifetime of storm Si [denoted as b(Si) and e(Si), respectively]. Since all analysis involves 3-h time steps, the rainstorm duration is l(Si) × 3 h, where l(Si) = e(Si) − b(Si) + 1.
Size: Rainstorm size is a vector of length l(Si) time steps. The storm identification algorithm identifies a number s(Si, t) of grid cells as part of the storm at each time step t. (We do not need to estimate fractional coverage of individual grid cells when using a grid as fine as 12 km × 12 km or 4 km × 4 km.) The size at each time step is the s(Si, t) × 144 km2 for the model output and s(Si, t) × 16 km2 for the observational data.
Mean intensity: Also for a time series of length l(Si) at each time step t when the rainstorm Si exists over the contiguous United States, we compute its average precipitation amount a(Si, t) by taking the average of the precipitation intensity over all grid cells identified with the storm, that is, , where υi,t,k [where k = 1, …, s(Si, t)] is the precipitation intensity at each gridcell location of Si at time t.
Central location: The central location of each rainstorm event Si is an l(Si) × 2 matrix, where each row c(Si, t) is the center of gravity weighted by the precipitation amount in each grid cell. This location measure is invariant to rotation and translation of the rainstorm event on the surface of the globe. (See section S2 in the supplemental material.)
Application of this procedure does still result, in both model output and observations, in large numbers of small rainstorms with negligible precipitation, many with very short lifetimes. These tiny storms would not affect estimates of average storm intensity and size, which are weighted by precipitation amount and event size, but if included they could dominate estimates of number of storms and their mean duration, reducing the utility of these factors. Because our goal is to provide information about events that contribute significantly to overall rainfall, we therefore exclude in all subsequent analyses those rainstorms with the lowest individual precipitation amounts. That is, we remove those small events in the tail of the cumulative distribution whose values add up to 0.1% of total precipitation.
b. Factorizing total precipitation
Using the computed metrics for each rainstorm, we can summarize the spatiotemporal aspects of precipitation patterns over a specified region by decomposing the total precipitation amount into the following four factors:
Our definition of these factors (see below) is chosen to ensure that a given fractional change in any factor leads to the same fractional change in total precipitation amount. That is, the factors allow us to quantitatively compare the effects of changes in different rainstorm properties. In section 5, we use these factors to compare the model baseline period to Stage IV data to evaluate model–observation discrepancies in rainstorm characteristics. We also use the same approach to compare the model baseline and future periods to identify changes in spatiotemporal precipitation patterns and quantify how they contribute to the total precipitation amount change. Factorization results are presented in a series of tables (see Tables 1–5). Such an analysis necessarily involves aggregating events over some reasonably large region; to identify subregional differences we also conduct a spatial analysis described in section 4c below.
Each factor is defined as follows:
Duration factor (hours per rainstorm): Duration factor is the mean duration per each rainstorm event. The factor is computed as 3 h × , because the length of each time step is 3 h.
Number of rainstorms: The number of rainstorms is simply given by N.
c. Spatial analysis
We can also use the computed metrics to find and visualize the spatial distribution of rainstorm characteristics rather than aggregating regionally. Precipitation intensity can be simply computed at individual grid cells without making use of rainstorm identification, but size, duration, and number require rainstorm identification and some statistical method for assigning a value to an individual grid cell. We use for this purpose a nonparametric kernel estimator, a standard procedure in spatial point process analysis (see, e.g., Gelfand et al. 2010, chapter 21). The procedure is most simply illustrated for estimating the expected number of rainstorms n originating at a given spatial location s, as
Here k(s1, s2) is a kernel function that satisfies = 1 and smoothly decays as the distance between s1 and s2 increases, and is a set of all gridcell locations in the contiguous United States. (See section S3 in the supplemental material for details about the choice of the kernel function.) For rainstorm size and duration, we use a similar kernel estimator but now compute the spatially weighted average value of the rainstorm metric, with the spatial weighting given by the kernel function. For further details see section S4 in the supplemental material. (Note that the regional factors produced by this analysis will not sum exactly to the total changes in seasonal mean precipitation, especially when seasonal mean precipitation is determined as an unsmoothed quantity computed at each grid cell.)
d. Comparing rainstorm characteristics in two precipitation patterns
The approaches described in sections 4a–c above provide a framework to compare any two spatiotemporal precipitation datasets and to identify similarities and differences in rainstorm characteristics. In the analysis below, we decompose total precipitation amounts for our different model and observational datasets into their component factors using the methods of section 4b; to identify differences, we compute the ratios of the factors. This analysis helps us to understand, in an average sense, which characteristics of rainstorms drive differences in spatiotemporal precipitation. To identify geographic variations in the various rainstorm characteristics, we also use the spatial analysis of section 4c to map individual rainstorm metrics across the study area, and again compare the different datasets by taking ratios. To quantify the statistical uncertainty in estimating the ratios of factors, we construct 95% confidence intervals by bootstrap sampling (shown in the last column of Tables 2–5). See section S5 in the supplemental material for details of uncertainty analysis.
The methods described above allow us to compare rainstorm characteristics both between model output and observations and between model output for present-day and future climate states. In both cases we know a priori that there must be differences in spatiotemporal precipitation characteristics.
a. Model–observation comparison
We find that WRF output and Stage IV data show substantial differences in rainstorm characteristics, but with different discrepancies in summer (Fig. 8 and Table 1) and winter (Fig. 9 and Table 1). In both seasons, we find that intensity in modeled precipitation events is considerably too low (as was evident even without rainstorm identification) even though U.S. average total precipitation is similar in model and observations. In summer, the compensating factor that resolves this discrepancy is predominantly size. WRF summer rainstorms are roughly 30% too weak but also 40% too large (Table 1). The size bias is similar to that found by Davis et al. (2006b) for severe summer storms in a single-year WRF simulation. Whereas Davis et al. (2006b) report that the duration of severe storms is too long, we find no overall bias in summer rainstorm duration in WRF. In our model output, the total number of summer storms across the study area is similar to that in observations, but with geographic patterns of bias: too many storms in the north and too few in the south (Fig. 8).
In winter, when U.S. precipitation is predominantly large-scale rather than convective, model biases are both qualitatively and quantitatively different than in summer. The intensity and size effects seen in summer become still larger in winter, and other factors become discrepant as well. WRF winter rainstorms are substantially weaker, larger, fewer, and longer than those in Stage IV data (Fig. 9 and Table 1). Winter storms are half as intense in the model as in observations and more than twice as large, biases that roughly cancel. The number of winter storms is also half that of observations, and their duration twice as long. Changes in storm number are nearly uniform across the study area, but other biases do not cancel locally. They do approximately cancel across the study area, so that while local precipitation amounts can be discrepant from observations, model and observations again are in reasonable agreement for total U.S. average precipitation.
b. Model present–future comparison
While the model–observation comparison is sobering, WRF runs offer a self-consistent response to changed climate conditions that is useful to study. The simulations must involve some compensating change that reconciles the diverging increases in precipitation amount and intensity. Different models may in fact find differing solutions to these constraints, but the way to gain insight into potential future changes in precipitation characteristics is to carefully examine model responses.
We compare the baseline and future WRF runs using the same approaches described above. The precipitation responses to climate change again differ by season, and show similar distinctions between summer convective and winter large-scale precipitation as those seen in the model–observation comparison. In the summer, the main drivers of changing precipitation patterns are again rainstorm intensity and size (Table 2). Changes in the other aspects (duration and the number of storms) are relatively small and not statistically significant. For convective precipitation, at least in the WRF runs used here, the compensating mechanism that allows a 3% K−1 increase in total precipitation amount that differs from the 6% K−1 increase in precipitation intensity is a decrease in the size of individual rainstorms (23% K−1).
In the winter, changes in precipitation patterns in future climate conditions are more complicated (Table 3). Across the continental United States the winter intensity effect is about a third less than that in summer, but still significant (4.4% K−1) and larger than the change in total rainfall (2.5% K−1). In contrast to summer, reduction in the size, duration, and number of storms all play a role in compensating for that discrepancy. (Size is estimated as the most important factor, but these differences are not statistically significant.) Interestingly, while each of these factors shows different geographic patterns of change (Figure 10), the change in seasonal mean is much more spatially uniform. This regional variation in compensation mechanisms means that different regions show different spatiotemporal precipitation changes that can be significant for socioeconomic impacts.
Winter patterns of total precipitation changes are also less uniform than those in summer. In both seasons, a part of the U.S. Southwest shows decreased rather than increased precipitation under future climate conditions, but the area of drying becomes larger in the winter, extending to about one-third of the contiguous United States. This area provides a useful test area for examining the mechanisms driving changes in rainstorm characteristics, since it presents an even stronger discrepancy between changes in intensity and total amount of precipitation than does the U.S. average. In the drying region, the intensity response is muted (and even negative in places), but still on average positive (Tables 4 and 5, the first and second rows). That is, average storm intensity increases (from +2% to +3% K−1 in summer and winter) while total precipitation actually decreases (from −5% to −6% K−1). We therefore repeat the analysis of changes in precipitation characteristics on the region of drying. (The analyzed areas for each season roughly coincide with the red regions in the maps in the top panels of Fig. 10; see Fig. S4 in the supplemental material for the exact areas analyzed.) We find that reduction of rainstorm size remains the dominant compensating factor in both summer and winter, with storm shrinkage substantial enough (more than −6% K−1 in both seasons) to result in decreased total precipitation even though storms are stronger (Tables 4 and 5). The main drivers of changing precipitation patterns in drying areas again appear to be partially compensating changes in intensity and size.
Finally, we conduct a preliminary graphical exploration of potential changes in rainstorm trajectories in changing climate conditions. Once rainstorms are identified across time, the trajectory of each storm is easily found by linking its central locations at each time step. Figure S5 in the supplemental material shows rainstorm trajectories in the baseline and the future model runs. Although more detailed analysis would be useful, we see no clear sign of changes. That is, regional differences in total precipitation change do not appear to be driven by deviations on storm trajectories.
6. Simulating future precipitation patterns by changing size and intensity distributions
The WRF runs analyzed here suggest that rainstorm properties would change significantly under future climate conditions, but model biases are strong enough that model projections alone are unsuitable for impacts assessments. We therefore seek to produce simulations of future precipitation that reflect model-projected changes in rainstorm characteristics, but also incorporate information from observational data.
For simulating future rainstorms, attempting to correct underlying biases in model projections is not an appropriate strategy, because model biases in rainstorm characteristics are much larger than the expected changes from present to future climate states. A better approach is that of “data-driven simulations” (i.e., transforming present-day observations into a future simulation using information from climate model runs). The approach has been used in simulations of temperature (e.g., Ho et al. 2012; Hawkins et al. 2013; Leeds et al. 2015; Poppick et al. 2016) and in a few studies of precipitation (e.g., Hay et al. 2000; Räisänen and Räty 2013; Räty et al. 2014), but in no case allowing for changes in spatiotemporal dependence. (This weakness is shared with bias correction approaches applied to individual time series.)
In this section we propose an algorithm for precipitation simulations that modifies the two aspects of rainstorm characteristics that change most significantly in future conditions, size and intensity. This method should be highly appropriate for simulating future summer precipitation, when factors other than size and intensity are largely unchanged. For winter precipitation, it may be problematic because duration and number of storms also show large, generally compensating changes. (Accounting for those changes would be complicated, because the number of storms increases in future projections in certain regions, requiring the creation of rain events not present in observations.) Our simulation approach consists of two steps:
changing the size of individual rainstorms in the observational record, and
changing the distribution of precipitation intensity in the transformed data from step 1.
In both cases, we account for regional differences in projected rainstorm changes using location-dependent resizing factors and gridcell-specific intensity distribution transformations. (See sections 6a and 6b below for further details.) In the remainder of this section, we first describe the algorithm, and then validate it by testing whether it produces improved simulation results over a simple gridcell-level simulation that incorporates no spatial information. We describe only the basic steps here, and give details and equations in section S6 in the supplemental material.
a. Changing individual rainstorm sizes
To change the sizes of observed rainstorms, we first determine the location-dependent resizing factors by comparing the estimated size functions between the baseline and future model runs (computed using the approach in section 4c). We then resize each observed rainstorm at each time step using the resizing factor that corresponds to its central location. To avoid changing the shape of a storm, we follow a simple resizing procedure that involves changing the distances between the storm center and all individual grid cells in the storm by the same factor. We regrid the rainstorm mathematically to a grid whose spacing is the inverse of the resizing factor (Fig. 11a, center). We then resize these new grid cells to those of the original grid, while keeping the center of the resulting rainstorm as close as possible to its original location (Fig. 11a, right).
Note that in some cases a rainstorm event is split into two or more substorms that are far from each other. For resizing purposes, we consider that those two parts warrant individual treatment if the shortest distance between their edges is greater than 120 km. In this case we define separate central locations and separate resizing factors for each substorm in the scattered rainstorm event.
b. Changing intensity distribution at each gridcell location
In the second step of the algorithm, we change the time series at each individual grid cell to adjust the marginal distribution of intensities (see Fig. 11b for illustration). The basic idea is to change the observed time series using a transformation that turns the marginal distribution of the model baseline time series into that of the corresponding model future time series. Since the observational data are not on the same grid as the model output, we determine each transformation based on the closest model gridcell location. To ensure that the model baseline output is treated identically to the observational data when adjusting intensities, we begin by applying the same resizing procedure to the baseline model output as used on the data.
We then change the number of wet time steps (time steps with positive precipitation amount) at each location in the observational time series as suggested by the model at the same geographical location. If the model projection shows a decrease in wet time steps, we apply the same fractional change by turning the lowest observational intensities into zeros. If the number of wet time steps increases, we promote some dry time steps (time steps with no precipitation amount) to wet time steps to apply the same fractional change. To choose the dry time steps to promote, we use an idea similar to the protocol described in Vrac et al. (2007), which creates rainfalls close in space or time to existing rainfalls. We first promote as many time steps as possible based on the precipitation amounts in their spatially adjacent grid cells. If there are not enough grid cells with positive precipitation in spatially adjacent grid cells, we select the time steps to promote based on the precipitation amounts in their temporally adjacent time steps.
Once the number of wet time steps is changed, we transform the marginal distribution of positive intensities for each time series. For each value of an individual time series, we take as the rescaling factor the ratio between the corresponding quantiles of the baseline and future intensities. Multiplying each value in the time series by the appropriate rescaling factor produces a simulated time series for a single location. The process is then repeated across all grid cells to produce the simulation of future precipitation.
c. Comparing simulation results to gridcell-wise simulation
We evaluate the performance of our simulation approach by applying it to the baseline model run itself, whose future state is known. We then check if the resulting simulated precipitation patterns indeed reproduce statistical characteristics of the actual future run. This test can provide at least minimal validation that the algorithm accurately reproduces the changes in the WRF runs it was trained on, and that it does so better than simpler alternative approaches. If the algorithm were to be used as an emulator, it would be appropriate to conduct a more challenging test, applying the algorithm to different runs of the same model not included in the original training set. (Our model simulations preclude this, however, as they consist of a single baseline and future realization.)
As a metric for comparison, we evaluate the distribution of dry spell lengths (successive time steps without rainfall), a characteristic that is not itself explicitly adjusted in our simulation method. We compare the performance of our simulation, which incorporates spatial information about rainstorm events, to a simpler “gridcell-wise simulation” that does not use information from adjacent grid cells, but instead simply adjusts the intensity distribution at each location as in section 6b. Table 6 summarizes the results for five regions of the United States (Midwest, Great Plains, Northeast, Wet South, and Dry South; see Fig. S6 in the supplemental material for region definitions). In most regions, the rainstorm event-based approach offers advantages over the simpler gridcell-based approach in summer (when precipitation changes are dominated by intensity and size, the two factors our approach addresses). In winter, when rainstorm changes include features not captured by our approach (changes in number and duration), the event-based approach offers comparable performance but no clear advantage.
It may seem counterintuitive that incorporating spatial information is important for capturing dry spell lengths: at first glance, the event-based approach may seem relevant mostly for capturing the spatial characteristics of precipitation and not for improving the temporal distribution at an individual location. However, changing rainstorm events can lead to better simulation results even at individual grid cells by allowing more flexible changes. Figure 12 illustrates this point using an example time series of observational data at a location close to Chicago. Our rainstorm event-based approach, which alters sizes of events, often removes a whole precipitation event from the simulated time series if that event represents the edge of a larger system that is made smaller in the future simulation. The gridcell-level approach, by contrast, can adjust the magnitude of the precipitation rate and can reduce the duration of events, but cannot account for the changes in event occurrence that result when storms shrink.
Climate model projections robustly imply that the spatiotemporal characteristics of precipitation events must change in future climate states. To help understand those potential changes, we have developed a new framework for analyzing changes or differences in rainstorm characteristics, including metrics for average intensity, size, duration, and number. The analysis framework is applicable both for comparison of future to present-day model simulations and for characterizing and validating the performance of models against observations. The same metrics allow us also to construct a method for simulating future precipitation events that combines model-projected changes with observational data.
Using this framework, we compare rainstorm properties in present and future high-resolution (12 km) dynamically downscaled model runs (WRF driven by CCSM4), and between those runs and the Stage IV radar-based observational data product. In all cases, the largest factors driving differences in rainstorm properties are intensity and size. In the summer season, when U.S. precipitation is predominantly convective, intensity and size are virtually the only factors of importance. In the model–observation comparison, WRF summer storms are too large but also too weak (leaving total precipitation consistent with observations). In the present–future model run comparison, WRF summer storms become smaller and stronger under future climate conditions (allowing total precipitation to rise less steeply than storm intensity). The same size–intensity tradeoffs are apparent in winter, when U.S. precipitation tends to be large-scale, but differences in the duration and frequency of rainstorms also become important. In the model–observation comparison, WRF winter storms are not only too large and weak, but also too few and too long-lasting. In the present–future model comparison, WRF winter storms become smaller and stronger, as in summer, but also less numerous and of shorter duration.
These parallels may aid in understanding the underlying causes of model bias. In the WRF simulations analyzed here, model–observation biases are generally larger than the projected changes under nearly a century of business-as-usual CO2 emissions. (To compare model projections to model bias, see Tables S1 and S2 in the supplemental material, which reproduce Table 1 but for the Stage IV data area only.) The WRF future simulation projects that rainstorms become 27% and 20% more intense (in summer and winter, respectively), but those future model storms remain weaker than observed present-day storms (by 30% and 27%). Size biases are even more persistent: for example, future WRF winter rainstorms become 4% smaller, but remain 130% larger than in observations (Tables S1 and S2). The scale of these spatiotemporal biases, in one of the most commonly used models for regional simulations, suggests that identifying detailed rainstorm characteristics is essential for validating and improving the representation of precipitation in models.
Model bias in rainstorm characteristics, while a cause for concern, does not necessarily invalidate the utility of the present–future model comparison. The divergent changes in precipitation intensity and amount are well-understood consequences of physical constraints and are robust across models (Knutson and Manabe 1995; Allen and Ingram 2002; Held and Soden 2006; Willett et al. 2007; Stephens and Ellis 2008; Wang and Dickinson 2012). The rise in intensity is driven by the increased water content in warmer air; the rise in total precipitation is the consequence of a more infrared-opaque atmosphere, which mandates a greater export of energy from surface to atmosphere as latent heat. The robustness of these effects in models suggests that the same changes will manifest in the real world as climate evolves. Precipitation events in the real world would then also experience some compensating change in spatiotemporal characteristics. Model simulations provide a self-consistent system governed by these fundamental constraints and at least a potential analog for understanding the real-world response.
In the WRF runs analyzed here, compensation largely occurs through reduction in rainstorm size, with the size factor overwhelmingly dominant in summertime convective precipitation. These changes would have implications for the societal impacts of changing hydrology. The robust increases in precipitation intensity in model projections have prompted concern about increased severe flooding events. The WRF runs here suggest that risk may be lessened by reductions in storm size. At the scale of a drainage basin, summer precipitation per rain event in our simulations rises not as Clausius–Clapeyron but as the smaller rise in total rainfall amount. In winter, precipitation events show more complex and potentially more impactful changes than in summer. Intensity increases are somewhat less than in summer, but storm frequency (number of initializations) plays a larger role in compensating for the intensity–amount discrepancy. The combination means that at the scale of a drainage basin, precipitation per event would increase more in winter than in summer. In addition, wintertime precipitation changes in WRF are more spatially heterogeneous, implying that local impacts of future precipitation changes may be geographically diverse (which further limits the possibility of diagnosing aggregate effects through observational studies in limited regions; see, e.g., Berg et al. 2013). Given the societal importance of projected changes in spatiotemporal precipitation characteristics, a research priority should be to understand whether the response shown in these WRF simulations is robust across models. All model simulations must involve some change in precipitation characteristics that compensates for diverging changes in intensity and amount, but different models or model configurations (including explicit rather than parametrized convection) may lead to different solutions to those constraints. Analysis methods that identify individual precipitation events are essential for this purpose.
The societal importance of potential flood increases also motivates our development of precipitation simulation methods. Future projections must play a role in risk assessments, but the large observation–model biases we have demonstrated mean that even if model-projected changes are believed to be true, risk assessments should not use model output alone. Given the scale of the biases, simple bias correction methods are also not appropriate. We therefore advocate instead a data-driven approach that extracts only projected changes from model runs and use them to transform observations into future projections. We have shown that our event-based simulation of WRF outperforms simpler simulation methods for summer precipitation, where compensation occurs largely through reduction in rainstorm size. WRF simulation of winter precipitation (or precipitation in any model simulation with more complex changes) remains an outstanding research question. Since our WRF runs project a future increase in winter storm number in certain regions (e.g., parts of the East Coast and Midwest, and most of the Southeast), a data-driven simulation algorithm could not capture such an increase by simply modifying existing storms, and would have to create new rainstorms in the observational record. Creating physically plausible and spatially consistent new precipitation events would likely require consideration of many other variables, including temperature and relative humidity.
This analysis represents, to our knowledge, a first attempt to understand and simulate model-projected changes in precipitation characteristics such as size, duration, and frequency of rainstorm events. The methodology for identifying and tracking storms should open up many other potential areas of research. While we have analyzed only a single model, comparison across models and model configurations is a clear research priority. It would also be useful to subdivide precipitating events according to meteorological context, and to separately examine the behavior of individual convective cells. Another potentially important area of research is to understand future changes in storm tracks. While our preliminary analysis shows no obvious changes in rainstorm trajectories, further analysis would be useful. The existing studies that examine rainstorm trajectories do so rather informally (Hodges 1994; Morel and Senesi 2002b; Xu et al. 2005; Assunção et al. 2012) and formal statistical analysis of rainstorm trajectories is a largely unexplored area. Finally, studies at other latitudes may also provide additional insight. The changes in midlatitude precipitation characteristics shown here are not a simple function of temperature rise and would not be identical across geographic regions. In some parts of the tropics, total precipitation may rise greater than instead of less than the Clausius–Clapeyron relation, requiring a different compensating response. Understanding how precipitation characteristics change in response to geographically differing constraints on the hydrological cycle may provide new insight into convective organization and structure. All these studies are made possible only given a methodology for identifying and studying the physical characteristics of individual storms.
The authors thank Dongsoo Kim and Mihai Antinescu for helpful comments and suggestions. This work was conducted as part of the Research Network for Statistical Methods for Atmospheric and Oceanic Sciences (STATMOS), supported by NSF awards 1106862, 1106974, and 1107046, and the Center for Robust Decision-making on Climate and Energy Policy (RDCEP), supported by the NSF “Decision Making under Uncertainty” program award 0951576.
Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JCLI-D-15-0844.s1.
Current affiliation: Department of Mathematical Sciences, University of Cincinnati, Cincinnati, Ohio.