Significance Large accumulations of rainfall over a precipitation event can impact human infrastructure. Unlike precipitation intensity distributions, probability distributions for accumulations at first drop slowly with increasing size. At a certain size—the cutoff scale—the behavior regime changes, and the probabilities drop rapidly. In current climate, every region is protected from excessively large accumulations by this cutoff scale, and human activities are adapted to this. An analysis of how accumulations will change under global warming gives a natural physical interpretation for the atmospheric processes producing this cutoff, but, more importantly, yields a prediction that this cutoff scale will extend in a warmer climate, leading to vastly disproportionate increases in the probabilities of the very largest events.

Abstract Precipitation accumulations, integrated over rainfall events, can be affected by both intensity and duration of the storm event. Thus, although precipitation intensity is widely projected to increase under global warming, a clear framework for predicting accumulation changes has been lacking, despite the importance of accumulations for societal impacts. Theory for changes in the probability density function (pdf) of precipitation accumulations is presented with an evaluation of these changes in global climate model simulations. We show that a simple set of conditions implies roughly exponential increases in the frequency of the very largest accumulations above a physical cutoff scale, increasing with event size. The pdf exhibits an approximately power-law range where probability density drops slowly with each order of magnitude size increase, up to a cutoff at large accumulations that limits the largest events experienced in current climate. The theory predicts that the cutoff scale, controlled by the interplay of moisture convergence variance and precipitation loss, tends to increase under global warming. Thus, precisely the large accumulations above the cutoff that are currently rare will exhibit increases in the warmer climate as this cutoff is extended. This indeed occurs in the full climate model, with a 3 °C end-of-century global-average warming yielding regional increases of hundreds of percent to >1,000% in the probability density of the largest accumulations that have historical precedents. The probabilities of unprecedented accumulations are also consistent with the extension of the cutoff.

Occurrences of intense precipitation are projected to increase (1⇓⇓⇓⇓⇓⇓–8) associated with higher atmospheric moisture content (9, 10) under global warming. Measures of precipitation intensity, coarse-grained to 1- to 5-d intervals, exhibit end-of-century increases on the order of 20% for wettest annual 5-d rainfall (8) or 10% in average wet-day intensity (11) or 17% °C in the 99.9th to 99.999th percentiles of daily precipitation (12) in business-as-usual anthropogenic forcing scenarios. Associated with this, substantial increases in frequency of high-rain-rate events can occur (13⇓–15), and the return times of events exceeding a given threshold decrease (16).

Time-integrated accumulation, the amount of precipitation that falls during a single event, is of concern for many societal impacts (17). Because more intense precipitation could, in principle, yield shorter event durations (10), the expected change in accumulation probabilities is unclear. Here, we derive a stochastic prototype from a fundamental climate model equation. This leads to an explanation of key properties of the probability density function (pdf) of accumulations noted in station observations (18⇓–20)—why the pdf of accumulation size drops slowly with increasing size over many orders of magnitude before reaching a cutoff scale, after which the pdf drops rapidly for very large accumulations. From the theory, we show that physical balances creating this behavior regime imply extremely high sensitivity for the very largest events under climate change. This motivates an evaluation of the accumulation distribution and its changes under global warming in a global climate model, the Community Earth System Model [CESM1 (21)]. The fact that the integrated precipitation during the event plays a key physical role makes accumulation a natural variable for examining these impacts.

Moisture Equation and Dynamics of Accumulation The vertically integrated moisture equation for a climate model may be written ∂ t q = − P + C , [1]where ∂ t q is the rate of change of column water vapor q , and C is the vertically integrated moisture convergence into the atmospheric column, including evaporation at the surface and horizontal moisture convergence by atmospheric transport. This moisture convergence term can be split into a climatological component C ¯ and variations C ′ that have large fluctuations due to atmospheric internal variability. The precipitation P can depend on the vertical structure of both the water vapor and the temperature, but because moisture variations tend to be dominated by a deep vertical structure (22), the onset of precipitation (as measured by probability of precipitation occurring and by the conditional average precipitation) increases rapidly above a threshold value of q for a given tropospheric temperature in observations (23⇓–25) and climate models (26). Both the onset of conditional instability of atmospheric deep convection and large-scale saturation can contribute to this. The time-integrated accumulation s from the beginning of the precipitation event (when P first exceeds a small threshold value) at time t 0 to the end of the precipitation event at time t f is s = ∫ t 0 t f P ( t ′ ) d t ′ , [2]which gives the water lost from the atmosphere over the event. This moisture loss plays a role in terminating events, so it is useful to define a running accumulation s ∼ , i.e., the amount of water rained out up to time t within an event s ∼ ( t ) = ∫ t 0 t P ( t ′ ) d t ′ , i.e., d s ∼ = P d t . [3]The size s of the accumulation is simply the value of the running accumulation s ∼ when the event terminates. We now introduce two cases of a stochastic prototype model based on simplifications of Eq. 1. Case 1 is sufficiently simplified to permit analytic solutions that guide understanding of the physics, whereas Case 2 is integrated numerically in time and then diagnosed as for a climate model. For Case 1, we make use of the one-to-one monotonic relationship Eq. 3 between time and the running accumulation, using s ∼ instead of time during precipitation events. The increment d s ∼ gives the moisture loss P d t due to precipitation in a time increment d t , noting that s ∼ and s have units of vertically integrated moisture, commonly expressed in mm (equivalent to 1 kg m−2 of water). Precipitation during an event is an order of magnitude larger than climatological moisture convergence C ¯ , so for simplicity we neglect C ¯ in Case 1 (alternately, the definition of s ∼ could be altered to use P − C ¯ ). The fluctuations of moisture convergence C ′ are approximated as a noise D 1 / 2 d W s ∼ , yielding a simple stochastic equation for Case 1 d q = − d s ∼ + D 1 / 2 d W s ∼ [4]that serves as a prototype for moisture fluctuations within precipitation events. For Case 1, the stochastic process W s ∼ is assumed to have independent, Gaussian increments (a Wiener process, i.e., Brownian motion), and the coefficient D governs the amplitude of this moisture convergence noise term, here taken constant. For the sake of analytic solution, increments are taken to be independent in the s ∼ coordinate, an assumption that will be altered in the Case 2 numerical solution. Precipitation is taken to occur when q exceeds a threshold mimicking the observed onset above a critical value of column water vapor. For Case 2 of the stochastic approximation to Eq. 1, we approximate the precipitation dependence on q seen in observations and CESM (27), P = P ( q − q c ) , as a ramp function, zero below the threshold moisture q c and linear above, yielding d q = − P ( q − q c ) d t + C ¯ d t + D ∗ 1 / 2 d W t . [5]The climatological value of moisture convergence including evaporation, C ¯ , is taken constant. The variations C ′ are represented as a noise D ∗ 1 / 2 d W t , with W t again a Weiner process, but with independent increments in the conventional time coordinate. Although one can envision a hierarchy of sequentially more complex stochastic approximations to climate model equations, the cases here are sufficient to exhibit accumulation distributions that may be compared with the climate model. In both cases, the presence of a threshold for onset and termination of precipitation creates a first-passage problem (28, 29) in which trajectories in the precipitating regime above the threshold evolve by a competition between fluctuating moisture convergence and water loss by precipitation until the moisture first falls below the threshold that terminates the event. For Case 2, trajectories are computed numerically for a long time series that enters and exits events as for a climate model. For Case 1, one has a simple Fokker–Planck equation (30) corresponding to Eq. 4 for the evolution of the pdf p q ( q , s ∼ ) of the trajectories for q evolving as a function of s ∼ ∂ s ∼ p q ( q , s ∼ ) = ∂ q p q ( q , s ∼ ) + 1 2 D ∂ q 2 p q ( q , s ∼ ) . [6]This has a known solution (29) (reviewed in SI Text), and we return to the properties of the first-passage problem schematically in Physical Mechanisms and Implications for Robustness to discuss how the analytic solution relates to more complex cases. The key result for informing climate model analysis is the pdf p s ( s ) for events reaching an accumulation of size s at termination when q drops to the threshold value p s ( s ) ∝ s − τ exp [ − s / s L ] [7]with τ = 3/2 and s L = 2 D for Case 1, for s above a small-event cutoff (Eq. S3) that is too small to see in gauge observations (18, 19) or in the climate model. Features of this solution are seen in the top pair of curves of Fig. 1. There is a scale-free range over which the power law s − τ applies, within which the pdf p s ( s ) drops relatively slowly with each factor of 10 increase in s. A key transition occurs near s L , above which the exponential cutoff creates a sharp reduction in the pdf. The very largest events that would ever in practice be encountered are thus controlled by s L (i.e., by the coefficient D that measures the amplitude of moisture convergence variations due to atmospheric internal variability). Under global warming, the moisture content of the lower troposphere available for convergence tends to increase (9, 31). If the statistics of low-level wind convergence were to remain the same, many aspects of precipitation would tend to increase, known as the rich-get-richer or wet-get-wetter effect (10, 32). Here, the prediction of a change in s L , illustrated in the simplest prototype curve using s L = 145 and 210 mm, respectively, has substantial implications. The vertical arrow between these curves indicates where this modest change in the cutoff scale results in a factor of 10 change in the pdf for large events. This pdf increase grows rapidly with accumulation size above the historical cutoff value. Fig. 1. Probability density function (mm−1) of precipitation accumulation (in mm) for historical climate (1976–2005) and end-of-century (EoC) climate (2071–2100) from CESM and from two cases of a stochastic prototype. (Upper) Case 1 analytical solution, Eq. 7 with τ = 1.5 and cutoff scales s L = 145 and 210 mm, respectively (values as for E. China); Case 2 numerical solution to Eq. 5 (dashed) for 2 values of D ∗ , with comparison to the Case 1 analytical solution modified with τ = 1.25. (Lower) For 6 different regions (Materials and Methods) from CESM1 simulations (shifted vertically in 2-decade increments for readability). Bootstrap 25th to 75th and 5th to 95th percentiles are given by box and error bar. The pdfs for E. China and India have overlaid prototype solution curves with τ = −1.25 and cutoff scales as listed. Dashed curves show shifts of the historical case for comparison with EoC. The prototype makes clear that the changes are expected to occur disproportionately for the very largest events (i.e., approximately exponential changes in the rare portion of the size distribution above the regionally defined cutoff for the historical period). Dynamical feedbacks will affect this regionally, including the possibility of regional reductions in largest-event probability density if D decreases due to dynamical reductions in the variance of low-level convergence. Changes are small within the long power-law range because for τ > 1 , the pdf normalization is little affected by the changes in the large-event range. For Case 2 of the stochastic prototype, results (Fig. 1) exhibit close parallels to the analytic solution Eq. 7 of Case 1, but with the exponent τ adjusted to 1.25. A similar scale break occurs and is shifted to larger values when D ∗ is increased. Values in Fig. 1 are chosen to illustrate s L values relevant to the Eastern China (E. China) region of the CESM simulations in current climate and at the end of the century (see SI Text for details). In Case 2, deviations from the form of Eq. 7 as a function of scale may be noted in both the approximate power law and cutoff ranges, but the key point remains of an approximately scale-free range over which the pdf drops slowly, followed by a scale at which the pdf drops quickly. The Case 2 results can be understood as a modification of the Case 1 assumption of uncorrelated noise, as discussed in Physical Mechanisms and Implications for Robustness and SI Text.

SI Text Background on the First-Passage Process Simple Prototype. The simple stochastic prototype, Eq. 4, is used to approximate the moisture equation from the climate model, using the running accumulation s ∼ as a transformed temporal coordinate because it increases monotonically in time during a precipitation event but tracks an important physical variable, the integrated water loss due to precipitation. The steps implying the form of the accumulation distribution discussed in the text are expanded here. Once the problem is cast in the form of Eq. 4, and the noise term W is assumed to be from a Wiener process (Case 1 of the stochastic prototype), the corresponding Fokker–Planck equation has a classic form (29), except that s ∼ holds the place normally occupied by time. The evolution of the probability density p q = ( q , s ∼ ) as a function of column-integrated moisture q and s ∼ is given by this equation, Eq. 6, repeated here for convenience: ∂ s ∼ p q ( q , s ∼ ) = ∂ q p ( q , s ∼ ) + 1 2 D ∂ q 2 p ( q , s ∼ ) . [S1]For the first-passage problem, we consider that the precipitation event begins when the system first passes a threshold water vapor value q 0 and terminates when the water vapor drops below a slightly lower value q t = q 0 − ϵ , the threshold for precipitation termination. These conditions mimic the observed onset of precipitation at a critical value of column water vapor (for a given temperature profile), associated with conditional instability of entraining plumes or large-scale saturation in observations (22⇓–24) and in versions of the climate model considered here (26, 27). Within this simple prototype, the onset corresponds to an initial condition of a delta function of probability density at q 0 , and the termination corresponds to an absorbing boundary condition with probability density p q = 0 at q t . The evolution of the pdf for this case has a known solution (28, 29) p q ( q , s ∼ ) = A q σ ( s ∼ ) [ F ( q − q 0 + s ∼ σ ( s ∼ ) ) − B F ( q − q 0 + 2 ϵ + s ∼ σ ( s ∼ ) ) ] , [S2]where σ = ( D s ∼ ) 1 / 2 , A q is a normalization constant, and B = exp ( 2 ϵ / D ) is set by the boundary condition p q = 0 at the termination threshold. Note that the form F ( ⋅ ) is as in the scaling solution (Eq. 8), with the second occurrence shifted to maintain the boundary condition. For this case, F is simply a Gaussian. As the solution evolves forward in s ∼ , the probability of remaining above the termination threshold decreases by the mechanisms discussed in the main text. Evaluating the flux of probability across the termination threshold yields the probability that an event will terminate in a small increment surrounding a given s ∼ . This termination value of s ∼ defines the size s of the accumulation. This yields an inverse Gaussian distribution for the accumulation pdf p s ( s ) = A 0 exp [ − s S / s ] s − τ exp [ − s / s L ] , [S3]where A 0 is a normalization constant, τ = 3 / 2 , and s L = 2 D is the large-event cutoff, the key role of which is discussed in the main text. In the prototype full solution, there is also a small event cutoff s S = ϵ 2 / ( 2 D ) , which Doppler radar measurements suggest may be on the order of 10−3 to 10−2 mm (18). In gauge observations (19, 20) and in the climate model analyzed here, the small event portion of the range is not resolved. Thus, the distribution is simply considered over the range greater than the specified minimum accumulation s 1 ≫ s S , leading to the form p s ( s ) = A s − τ exp [ − s / s L ] , s > s 1 [S4]which corresponds to Eq. 7. The normalization constant A tends to be set by the smallest observable scale s 1 , because for τ > 1 , the low-s range dominates the integral when the power law range is long ( s 1 ≪ s L ) . As a result, increases in large-event probability density due to increases in s L create only a small adjustment in normalization constant, yielding small reductions in probability density over the power-law range. For instance, for s 1 = 0.2 mm, a change in s L from 145 to 210 mm, as in the top curve of Fig. 1, yields only a 0.04 % decrease in the normalization constant. Note that the time-mean rainfall is not closely related to the changes in the accumulation distribution. The time mean obeys long-term constraints from moisture and energy budgets that can affect the fraction of time spent precipitating. The dry-spell intervals between precipitation events have dynamics with strong parallels to that considered here for accumulation distributions, with the upward drift in moisture toward onset of precipitation driven by mean moisture convergence, including evaporation. For increasing moisture convergence variance, the probability density of the very longest dry spells in the prototype tends to increase consistently with the results shown here for the change in pdf of the largest precipitation accumulations. Considerations of Robustness and Relations with Other Systems. Despite the complexity of the climate model, the simple stochastic prototype was able to provide predictions of its behavior. While the analytic solution for Case 1 and the scaling solution aim at distilling physical insight, the numerical solutions for Case 2 of the stochastic prototype make clear that the slight adjustments to the exponent of the power law range seen in the climate model solutions are easily obtained for assumptions that are realistic in a climate-modeling context. Here, we elaborate briefly on how the form of the regime change at a cutoff scale is robust to these changes. The Case 2 model is solved in the time domain, as for a climate model (note that D * thus has different units than D). In Case 2, during nonprecipitating intervals, there is an upward drift in q due to C ¯ . Event onset is taken to begin for precipitation above a small threshold, such that P > C ¯ , and the drift is downward due to precipitation loss. The running accumulation s ∼ is then diagnosed within precipitation events in the Case 2 simulations. In this s ∼ coordinate, the Case 2 noise has temporal correlation and the spread of an ensemble of trajectories due to variations in moisture convergence increases less quickly than s ∼ 1 / 2 , consistent with the slower drop of probability density of accumulation size at termination of the event in the approximately scale-free range before the cutoff. This behavior corresponds to the subdiffusive case of anomalous diffusion, as commonly arises in tracer transport in complex flow (41). Because the time dependence of the spreading process is different from that of the drift process, as seen in the scaling solution, a regime change with a cutoff scale given by the competition of these processes must occur. The amplitude of the moisture convergence variations is the key parameter in this competition—an increase in this amplitude with increasing moisture necessarily extends the cutoff, and thus creates the selective increase in probability density of the largest events. The physics of this competition must occur even for more complex moisture convergence variations, as in the climate model. Prototypes for anomalous diffusion regimes and associated first-passage problems have been examined in a number of systems (42), so it is worth considering under what circumstances analogies drawn from these cases can be instructive. First-passage problem solutions exist for representation of these cases by fractional Fokker–Planck equations (43⇓–45). A key result is that the properties discussed in the Weiner process case are modified smoothly, with the exponent of the power law range adjusted. Properties are in some cases approximated by a generalized inverse Gaussian distribution as in Eq. 3 where τ can differ from 3/2 (46), similar to features noted numerically in both Case 2 and CESM1 here. Technical challenges involving divergent moments often arise in these representations, including failure of the method of images under non-Markovian conditions (42, 47). The system of interest here is Markovian with finite moments, and thus in a number of respects is simpler than these cases. The adjustment of the exponent is here associated with the observable of climate interest—the time-integrated loss term. Moment Ratio s M and Large-Event Cutoff Geographic Dependence in the Climate Model. In the main text, an event-size scale s M based on the ratio of second moment ⟨ s 2 ⟩ to the first moment ⟨ s ⟩ is used as an estimator proportional to the cutoff, as has been done to compare accumulation distributions for different regions estimated from the Department of Energy Atmospheric Measurement Program high-resolution observation sites (19) s M = ⟨ s 2 ⟩ / ⟨ s ⟩ . [S5]For the simple case (Eq. 3), the large-event cutoff s L is related to s M by s L = 2 ( ⟨ s 2 ⟩ − ⟨ s ⟩ 2 ) / ⟨ s ⟩ ≈ 2 s M . [S6]For Eq. 4, adjustments to the proportionality constant occur associated with s 1 and modifications to τ . A similar scaling also holds when the form of the cutoff exp [ − s / s L ] is replaced by amore general function G [ s / s L ] that approaches 1 for s ≪ s L and 0 for s ≫ s L (19, 20, 48). The form of a scale-free range followed by a cutoff is common to many systems. There is an underlying mathematical connection to the self-organized criticality literature that originally motivated observational work (18, 19)—certain first-passage processes can be put into one-to-one correspondence with simple models of self-organized criticality (29). Here, the derivation of the theoretical model from the equations of a climate model is key because we care foremost about the physical climate processes setting the cutoff in current climate and modifying it in future climate. Fig. S1 shows the estimate of the large-event cutoff s M in Eq. 6 as a spatial distribution for historical and end-of-century simulations, respectively. Fig. 3 shows the ratio of these. The full ensemble of 450 y is used, computing each of ⟨ s ⟩ and ⟨ s 2 ⟩ in Eq. 6. This first computation of such an estimate as a map shows spatial variations of s M are considerable within current climate (Fig. S1a), with values increasing substantially toward the tropics. In many regions, a tendency may be noted for s M to be larger over ocean regions than over land regions at comparable latitudes, likely associated with such factors as additional energy supply from the ocean surface and reduced impact of diurnal cycle. In the end-of-century case (Fig. S1b), large-scale spatial patterns remain similar, and the global warming increase in many areas can be seen primarily as a broad-scale increase on the order of 20 % , albeit slightly larger in some regions than others, as seen in the ratio in Fig. 3. This serves as a reminder that the changes in probability density of the largest events seen in Figs. 1 and 2 are increases in the largest events experienced for a given geographic region. Large-event changes should not be misconstrued as converting midlatitude accumulation cutoffs to values typical of the tropics. Nonetheless, because the s M changes translate into large changes in the pdf in the range of the very largest accumulations historically experienced in each region, these effects may be expected to play a key role for regional planners with regard to adaptation to potential flood magnitudes and other societal impacts.

Precipitation Event-Size Change Under Global Warming Motivated by the behavior predicted by the theory, we carried out an ensemble of 15 CESM1 simulations with the required high-resolution time output to assess these distributions under historical estimates of radiative forcing by greenhouse gases and aerosols and under Representative Concentration Pathway 8.5 (33) (Materials and Methods). Precipitation accumulation distributions for various regions in Fig. 1 consistently exhibit long power law ranges with exponent τ ≈ −1.25, a value in approximate agreement with prior observational estimates from shorter time series at various Department of Energy Atmospheric Radiation Measurement Program sites (19), and of slightly smaller magnitude than −1.5 (i.e., corresponding to a modestly subdiffusive case). Curves shown for China and India illustrate the simple form Eq. 7 with τ = −1.25. As expected, the power-law range is for each region followed by a large-event cutoff, above which the pdf of large accumulations drops steeply. Were the power law range to extend indefinitely, the mean and variance of the accumulation would diverge, an indicator of how important the physics of the cutoff is in limiting large events. Under global warming, the power-law portion of the distribution remains essentially unchanged, while, in many regions, the large-event cutoff changes, extending the power law range to slightly larger values on the s -axis. This extension increases the probability of these very largest events. To make this clearer, Fig. 2 shows the ratio of the pdf of each accumulation bin in the end-of-century climate to that under current climate. This corresponds to a risk ratio (34) for the s -intervals corresponding to each bin (i.e., for probability density rather than probability of exceedance). For events smaller than the cutoff, there is essentially no change in the pdf. Above the large-event cutoff, however, the pdf increases rapidly (roughly exponentially) for the very largest events. The extent of the increase for the various regions corresponds to differing degrees of change in s L . Changes for Midlatitude North America appear modest on this scale, but reach a 250% increase (150 to 450%, 25th to 75th) in its uppermost bin. Most of the regions have several bins in the large-event range, exhibiting increases of hundreds of percent relative to current climate. For E. China, the increase for the highest bin with historical precedent is >1,700%, with >1,000% for the 25th percentile bound. Fig. 2. Ratio of the pdf of a given accumulation under global warming to the pdf under current climate for the six regions (dots). Bootstrap 25th to 75th and 5th to 95th percentiles are given by box and error bar, respectively. The extension of the distribution tends to yield end-of-century occurrence of accumulation sizes unprecedented in the 450 y of the historical ensemble, indicated by ∞ symbols above the plot; the same effect yields long upward error bars in the highest bins (shifted slightly within bin for graphical clarity; bin boundaries are indicated by gray bars at bottom). Referencing future probabilities to current climate probabilities of a given accumulation, as in Fig. 2, is informative, but the end-of-century extension of the distribution yields accumulations that are unprecedented in the historical period. A direct measure of the change in distribution extent is thus useful for interpretation and for estimating statistical significance by geographic region. The change in distribution in Fig. 1 is very close to a simple rescaling, even for curves where the cutoff is more complex than exponential, as occurs for the lower four curves. This is indicated in Fig. 1 by overlaying curves interpolated for each historical distribution onto the corresponding end-of-century distribution, with s rescaled by a constant (a shift in log ⁡ s , with the pdf amplitude correspondingly rescaled). The match indicates that the end-of-century distribution is well captured by the historical distribution, rescaled by a single number for each region that typifies the change in the large-event cutoff s L . Probabilities of events that have no counterpart in the historical period (seen for E. China, India, and South American regions) are consistent with the rescaled distributions. Thus, the shift in the cutoff appears to be a useful paradigm for explaining, and potentially predicting, the occurrence of these unprecedented events. In addition to providing a simple hypothesis for future changes in corresponding observed distributions, the result that the distribution changes are governed by the shift in the cutoff can greatly aid in significance testing and displaying changes at a regional level. Fig. 3 uses an estimate s M from observational work (19) that is proportional to the large-event cutoff s L (Materials and Methods and Eq. 9) to indicate the simulated spatial distribution of these changes under future climate. The ratio of s M computed for 2071–2100 to the value from the historical period is displayed at each grid point, where this measure of the rescaling of the distribution differs significantly from 1 by a bootstrap test. As expected, there is regional variation, with some regions having little change, and even some regions of decrease. The predominant shift is toward increases (91% of grid points passing criteria for significance and number of events). Clausius–Clapeyron scaling of D (Materials and Methods), and thus of s L , in the prototype yields a factor slightly over 1.2 in s L , which is approximately the average over significant points in Fig. 3, although certain regions exhibit larger increases. Fig. 3. An estimate of the fractional change in the large-event cutoff scale s L . Accumulation moment ratio Eq. 9, s M , for the end-of-century (2071–2100) simulations divided by its value for historical climate simulations (1976–2005). Values > 1 imply increases in the cutoff scale, which yield large changes in probability of events larger than the cutoff. Increases are shown only where >95% of 1,000 bootstrap replications (Materials and Methods) exhibit an increase, and similarly for decreases. Points are masked (gray) where <4,000 events occurred in either 450-y ensemble. The regions over which accumulation distributions are shown in Fig. 1 are indicated as black outlines. The regions for which distributions are averaged in Fig. 1 are indicated in Fig. 3, chosen to sample regions of relatively modest change, such as the Midlatitude North American region, as well as regions of large change, such as a box in the Indian region and a box covering E. China and the neighboring ocean area including Taiwan. Factors rescaling the box-average distributions in Fig. 1 bracket the Clausius–Clapeyron scaling, from slightly below (1.14 for Australia and Midlatitude North America) to substantially above (1.46 for E. China). It is worth emphasizing that, even for regions where the rescaling in s is approximately consistent with Clausius–Clapeyron, the consequences for the frequency of the largest events can be substantial. A 20% increase in moisture convergence in Eq. 7 yields a fractional increase in the pdf of approximately exp ( 0.2 s / s L ) , which increases rapidly for the largest accumulation sizes. Thus, even the midlatitude North American region, with relatively modest increase in s M in Fig. 3, exhibits roughly a 250% increase in probability for the bin with the largest events that have historical precedents in Fig. 2.

Physical Mechanisms and Implications for Robustness This distinct behavior for the largest events is due to the time-dependent dynamics that yields the cutoff scale. Standard prototypes for changes in extreme events (17, 35) consider stationary distributions of a climate variable that in a warmer climate becomes shifted due to a change in mean or whose width changes due to a change in variability. For non-Gaussian distributions (Fig. 4A), such as precipitation rate (14, 15, 35, 36) or water vapor (27), the change in mean and variance are typically linked. In addition to affecting occurrences of extreme values (red), these mechanisms create changes throughout the probability distribution (decreases in blue; increases in light red). In the behavior for accumulations found here (Fig. 4B), changes are primarily in the probability of the most extreme accumulation range. How does this occur? Fig. 4. Mechanisms for changes in accumulation probabilities involve two evolution regimes. (A and B) Traditional schematics of a shift and change in width of a stationary distribution (A) (14,17, 35) differ from the effect for accumulations in which the pdf is susceptible to disproportional change for the largest events, above the cutoff (B). (C) The mechanism that creates this effect (schematized for Eq. 8) involves the evolution of the pdf of an ensemble of events (example trajectories shown in blue) before the first passage across an event-termination threshold at lower moisture. The power-law range of accumulations comes from spreading of the pdf due to internal variability, which dominates termination when accumulations are not too large. The large-accumulation regime occurs when drift toward termination by precipitation loss becomes important. A change in moisture-convergence variability changes the large-event cutoff that separates these regimes. A scaling argument from the theory here points to the essential features of the climate physics that create this sensitivity and provides a sense of the robustness. Probability distribution solutions of Eq. 6 during the event have the form p q ( q , s ∼ ) = 1 σ ( s ∼ ) F ( q − q 0 + s ∼ σ ( s ∼ ) ) , [8]where the width of the distribution σ spreads due to atmospheric variability as a function of the running accumulation s ∼ , which is acting as a transformed temporal coordinate. The prefactor σ ( s ∼ ) − 1 maintains normalization as the distribution spreads. For constant D and Gaussian noise, one obtains σ ∝ ( D s ∼ ) 1 / 2 , and F is Gaussian. The drift toward lower q is simply s ∼ in the numerator because this represents integrated loss by precipitation up to a given time. The event termination when moisture falls below a given threshold gives the first-passage process; the boundary condition at this threshold is met by a pair of such solutions (SI Text), so Eq. 8 is sufficient to show the key physical factors. Fig. 4C illustrates the collision of two temporal-dependences in Eq. 8 that yield the different behavior ranges. In early time, the spread of the distribution width like ( D s ∼ ) 1 / 2 dominates the flux of probability toward event termination at the threshold q t , yielding the power law range. At long time, the drift term s ∼ enters, causing a faster flux of probability across the event termination due to the integrated precipitation loss. The transition between these regimes occurs where these two terms are the same order, yielding a cutoff scale proportional to D . As illustrated by Case 2 of the stochastic prototype, the basic features resulting from this transition remain robust under circumstances for which analytic solutions are not available, but which capture the value of the exponent τ in the accumulation distribution of the climate model. The spread of the distribution per unit of water lost by precipitation increases at a rate slightly slower than s ∼ 1 / 2 , but the key transition to the cutoff regime occurs by the same physics as in the analytic solution of Case 1: namely, as accumulated water loss becomes important. This argument would apply to more complex spreading processes, explaining why these features carry over to the full climate model.

Discussion In the physical mechanism for changes in precipitation accumulation presented here, one of the essential ingredients, an increase in the variability of moisture convergence, is in common with the prevailing discussion of changes in precipitation intensity under warming. The other ingredients—an integrating variable and a threshold for event termination—dictate not a general broadening of the accumulation distribution, but a shift in the cutoff that limits very large events. These ingredients imply that, for integrated accumulations, a corollary of the rich-get-richer effect applies, especially to the largest events. This might be termed the biggest-get-more-frequent or the biggest-get-bigger, because successively larger increases in frequency occur for successively larger accumulation categories above the historical cutoff scale as the cutoff value increases. Geographic patterns of the changes should be viewed with caution from a single climate model, and caveats apply to climate model simulation of extreme precipitation events (6, 14, 27)—these results motivate evaluation over a wider ensemble of climate models, despite the need for high time-resolution output, and over observations more extensive than the set so far examined. The simplicity of the mechanism and of the resulting distributions for accumulation in Fig. 1 suggest that precipitation accumulation provides a natural coordinate for evaluating statistics of extreme precipitation change. The success seen in Fig. 1 in approximating changes in the pdf, even for unprecedented events, by a suitable rescaling of the historical distributions suggests that this can be exploited to evaluate scenarios for potential impacts.

Materials and Methods Precipitation accumulations in CESM are computed as the integrated precipitation from the first exceedance of a small threshold (0.4 mm/h) to the first drop below the threshold. These have not previously been computed in climate models—the high time-resolution data are not normally saved. Augmented output (for 30S to 50N) from an ensemble of runs with the fully coupled National Center for Atmospheric Research CESM1 (Version 1.0.5) (21) under historical estimates of radiative forcing by greenhouse gases and aerosols and under Representative Concentration Pathway (RCP) 8.5 (33) with precipitation saved at all time steps (30-min intervals) is used to compute the accumulation distributions in Fig. 1. The ensemble of 15 simulations of 30 y each are initiated from different atmospheric initial conditions for the historical and RCP8.5 simulations (37) to yield different instances of atmospheric and climate internal variability. CESM output data used in this study are available from the authors following guidelines in the CESM data plan (www.cesm.ucar.edu/management/docs/data.mgt.plan.2011.pdf). For the accumulation distributions and ratios in Figs. 1 and 2, the latitude–longitude ranges for the six averaging boxes are: Midlatitude North America (37N to 50N, 240E to 288E), Africa (20S to 18N, 10E to 44E), Australia (28S to 16S, 120E to 152E), South American region (20S to 0, 290E to 320E), India (15N to 25N, 70E to 90E), and E. China region (20N to 30N, 110E to 130E), which includes ocean regions surrounding Taiwan, the East China Sea, and the southwest islands of Japan. Accumulation distributions are computed for each spatial point and then averaged within these domains. As seen in Fig. 3, the latter two regions provide examples within the long band across Southeast Asia exhibiting substantial increases in s M ; the other boxes characterize broad continental areas while remaining within tropical/midlatitude ranges of historical s M (Fig. S1) . Shifts of the historical curves in Fig. 1 to approximately match end-of-century correspond to rescalings of s by a factor of 1.14, 1.22, 1.14, 1.35, 1.29, and 1.45, in the order above. Sensitivity tests on box choices (e.g., dividing the African box in two) give robust results consistent with s M changes in Fig. 3. The boxed regions are chosen primarily over land due to relevance for societal impacts; Fig. 3 provides a good indicator of where similar results apply over oceans. Fig. S1. The accumulation moment ratio, s M (mm), an estimator proportional to the large-event cutoff value s L . (A) For historical climate simulations (1976–2005). (B) For the end-of-century (2071–2100) simulations. The bin size in Figs. 1 and 2 has a value of 0.2 in log ⁡ s . The pdf is normalized by bin increment in s and by the total number of events (i.e., the pdf provides information about probability for different sizes of events, given that an event occurs). The frequency of all events over a given time interval provides complementary information that can be important in interpreting measures of precipitation change (38). The total number of events tends to increase in the RCP8.5 ensemble relative to historical (Fig. S2), but the changes are modest compared with the differential changes in pdf in the large-accumulation range. Fig. S2. The number of events in the 450-y ensembles evaluated. (A) For historical climate simulations (1976–2005). (B) For the end-of-century (2071–2100) simulations. Masked where <4,000 events occurred in either 450-y ensemble. For analysis of spatial distributions of the cutoff (Fig. 3), we use an estimator of the cutoff s M = ⟨ s 2 ⟩ / ⟨ s ⟩ (i.e., the ratio of second moment ⟨ s 2 ⟩ to the first moment ⟨ s ⟩ ). For the exact version of Eq. 7 (SI Text), this is related to the large-event cutoff s L by s L = 2 ( ⟨ s 2 ⟩ − ⟨ s ⟩ 2 ) / ⟨ s ⟩ ≈ 2 ⟨ s 2 ⟩ / ⟨ s ⟩ = 2 s M [9]The constant of proportionality cancels for ratios of s M between future and current climate. The Clausius–Clapeyron scaling of D to provide a sense of the simplest expected increase in s L and s M uses a 7.3% increase in vertical and global average moisture per degree of global average temperature increase (39), and a corresponding temperature increase of 3.0 K (2071–2100 minus 1976–2005) for this model. Small variations about this can occur through regional differences in base temperature or temperature change (31), but regional departures in Fig. 3 from large-scale Clausius–Clapeyron increases are likely due to changes in the dynamical portion of moisture convergence. This scaling includes all quantities with units of moisture [i.e., a rescaling of the pdf for rain intensities (5), which can reduce temporal duration, is implicitly included]. A bootstrap procedure (40) is used to establish confidence intervals in Figs. 1 and 2 and to mask small changes in Fig. 3. The fifteen 30-y runs for each of historical and end-of-century are broken into 5-y segments, for a total of 90 segments, which are considered independent for these precipitation statistics. A set of 1,000 bootstrap replications is constructed by random picks with replacement from the set of segments to create 1,000 artificial ensembles of 90 segments. Statistics for each replication are computed exactly as on the original ensemble of 90 segments, for historical and end-of-century respectively. For Fig. 3, ⟨ s ⟩ and ⟨ s 2 ⟩ are each computed for the 450 y averages of a given replication and the moment ratio for each of historical and end-of-century (Fig. S1) computed before taking the ratio of these. The bootstrap distribution provides an estimate of the variations in the computed statistic due to sampling error. The 1,000 instances of each statistic are ranked, and values of the 5th and 95th percentile bounds are computed. Increases in the moment ratio are shown only if the 5th percentile is >1 (i.e., only if at least 95% of the replications indicate an increase). Conversely, decreases are shown only if at least 95% of the replications show a decrease. In Fig. 4A, a gamma distribution pdf, commonly used for daily rainfall intensity distributions (36), is used to illustrate the typical occurrence of changes throughout the distribution (14) as mean and variance increase (linear axes; shape parameter = 2; scale parameter = 3 and 4, schematizing historical and EoC cases, respectively). In Fig. 4C, the pdf as a function of moisture and s ∼ corresponding to Eq. 8 is schematized as shading with black contours for an idealized case for which the effect of drift toward the event-termination threshold due to precipitation loss can be seen within the diagram. The no-drift case (dropping the drift term) is shown as a single contour for reference. The white arrow shows the difference to the corresponding contour in the full case, indicating the growing importance of the drift term. For the Case 2 stochastic prototype, solutions in Fig. 1 are for P = α ( q − q c ) , with α = 0.1 h−1 and q c = 60 mm. C ¯ = 0.1 mm h−1 with an increased value for q < 10 mm to limit negative excursions; precipitation events are computed over intervals for which P > C ¯ . Two values of D ∗ = 20 and 27 mm2 h−1 yield shorter/longer cutoff cases shown, chosen to illustrate cutoffs similar to E. China historical and end of century CESM1 results, respectively.

Acknowledgments We thank J. E. Meyerson for graphical assistance; O. Peters for discussions; and several colleagues, two reviewers, and the editor for comments on the manuscript. Computations were carried out at the National Center for Atmospheric Research with high-performance computing support from Yellowstone(ark:/85065/d7wd3xhc) provided by NCAR’s Computational and Information Systems Laboratory, sponsored by NSF. This work was supported in part by Department of Energy Grant DESC0006739; National Science Foundation Grants AGS-1102838 and AGS-1540518; National Oceanic and Atmospheric Administration Grant NA14OAR4310274; and Office of Naval Research Grant N00014-12-1-0744 (to S.N.S.).

Footnotes Author contributions: J.D.N. designed research; J.D.N., S.S., S.N.S., and D.N.B. performed research; S.S. carried out climate model analysis; S.N.S. collaborated on theory; D.N.B. ran the Community Earth System Model simulations; J.D.N., S.S., and D.N.B. analyzed data; and J.D.N. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1615333114/-/DCSupplemental.