Climate models capture key features of extreme precipitation probabilities across regions

Global climate models (GCMs) are often the tool of choice to assess the precipitation response to climate change. However, despite their widespread use their ability to simulate daily precipitation probabilities in current climate, and consequently the confidence we should have in them regarding future projections, have often been called into question. For example, models are known to rain too often and too little [1–4], and have difficulty representing magnitudes of extremes [5–7], especially in the tropics [8, 9]. Indeed, precipitation is one of the climate variables that is most challenging to simulate, as it depends on parameterization of sub-grid scale processes, and changes can depend both on changes in moisture and changes in convergence due to feedbacks with heating, termed the thermodynamic and dynamic contributions, respectively [10–13]. The complexities and systemic problems of GCMs in simulating precipitation have led researchers to decry the 'Dreary state of precipitation in climate models' [14], and may lead to question how well founded are GCM projections of future changes in precipitation, including changes in the probability of the most societally impactful extreme events.

Previous assessments of changes in extremes have instinctively used fractional (i.e. percent) changes (instead of absolute), under the implicit assumption that despite biases in simulating the historical probability distribution of precipitation, fractional changes would still be useful. However, up to this point, there has not been a test that evaluates this assumption, and consequently how confident we should be in this approach to projecting future changes under global warming.

In this paper we focus on two key measures, which can be evaluated using historical simulations, that can decrease or boost confidence in projections of future fractional changes in precipitation extremes. These two measures are applied to the Coupled Model Intercomparison Project Phase 6 (CMIP6) Ensemble to assess

• (1)  
  CMIP6 models' ability to simulate the shape of the large daily precipitation probability tail in current climate (although not necessarily its magnitude), and
• (2)  
  CMIP6 models' ability to simulate fractional changes in the key scale that controls the extreme tail overall magnitude.

Good model performance in these two measures may be viewed as necessary conditions for reliable future projections but we underline that they are not sufficient conditions. Because our measures are designed to test performance against historical observations, they do not directly address all possible variants of model behavior under warming. Rather, they address whether the widespread existing practice of assessing changes in simulated probabilities relative to the simulated historical distribution [15, 16] can be given observational underpinnings for the large-event range.

2.1. Models and observational datasets

2.2. Probability distributions in usual and rescaled coordinates, and definition of the cutoff-scale

Probability density functions (PDFs) are estimated as normalized frequency histograms with bins approximately constant in log(P) space, with P denoting daily precipitation. In what follows we assume that daily precipitation PDFs are well characterized by a functional form:

$$\begin{equation} \textrm{PDF}(P)\propto P^{-\tau_{P}}G\left(-\frac{P}{P_{L}}\right), \end{equation} \tag{ 1 }$$

where G is a generic function and P_L is a precipitation scale. For observed distributions, this is a cutoff-scale because when $P\gt P_{L}$, the probability distribution drops much faster per order of magnitude precipitation increase compared to the slower power law range decay. The power law behavior ($\propto P^{-\tau_{P}}$) dominates the form of the observed precipitation PDF for low and moderate P values (for $P\lt P_{L}$). Theoretical distributions commonly employed to describe non-zero daily precipitation data and which follow this form include Exponential [24], Gamma [25, 26], Weibull [27, 28] and Generalized Gamma [29] distributions. To tease out common properties of the extreme tail (that is, for $P\gt P_{L}$), that are independent of scale, we assume that the cutoff-scale P_L is proportional to $\hat{P}_{L} = \sigma_{P}^{2}/\bar{P}$, where $\sigma_{P}^{2}$ and $\bar{P}$ are the variance and mean daily precipitation over wet days, which we define as days where P ⩾ 0.1 mm. This moment-based estimator is not the only way to estimate this scale, but works for the four aforementioned distributions and is reasonably robust for PDFs that depart from these. We also assume that the influence of the power law range is small in determining extreme event probability, which we find to be well justified in latter sections. In the case of a Gamma distribution $P_{L} = \hat{P}_{L}$ and $\tau_{P} = 1-\bar{P}/P_{L}$.

Based on these assumptions, we also calculate probability distributions in a scaled coordinate

$$\begin{equation} P^{*} = \frac{P}{\hat{P}_{L}}, \end{equation} \tag{ 2 }$$

which allows us to examine similarities in the shape (regardless of P_L value) between observed and modeled large daily precipitation probability tails (i.e. how close are modeled and observed G(P^∗) in equation (1)), as exemplified in section 3. Variants of this type of scaling have been used, for example, to study the large-event tail of precipitation event sizes (also known as event accumulations) [30–34]. To simplify notation, in what follows we drop the hat when making reference to $\hat{P}_{L}$.

The calculation of probability distributions can be noisy if short time series are used, especially for the extreme tail range. For this reason, to reliably calculate extreme tails we first divide the 50° S–50° N domain (the latitudinal extent of TRMM-3B42) into 120 regions of 10° latitude and 30° longitude and then calculate probability distributions in each one of them by pooling the time series from all grid points in a particular region. For consistency, we also use these regions for analysis in section 5.

Errors in the simulation of precipitation in climate models are difficult to quantify, due to substantial differences in observational datasets for daily precipitation [35]. This is illustrated in figure 1(a) which shows the daily precipitation PDF over the North Atlantic Storm Track region (30° N–50° N, 80 °W–40 °W), see figure 2(a)) calculated from four observational datasets. While PDFs for daily precipitation P calculated from reanalyses are well approximated by Gamma distributions of form PDF$(P)\mathop\propto\limits_{\sim} P^{-\tau_{P}}\textrm{exp}(-\frac{P}{P_{L}})$—an approximate power law range with exponent τ_P followed by a near exponential cutoff scale P_L limiting extremes [26, 36–40]—the probability distributions diverge sharply in absolute terms, as attested by their very different τ_P and P_L parameters. For example, probability distributions in GPCP and TRMM-3B42, which bracket the range of observational estimations in this region and most of the globe (figure S1), are characterized by τ_P  = 0.20, P_L  = 9.2 mm and τ_P  = 0.64, P_L  = 22.7 mm respectively. The effect of the cutoff scale P_L on extremes can be readily appreciated from figure 1(a), where a larger cutoff value in TRMM-3B42 corresponds to larger probability of extremes in absolute terms compared to GPCP.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As an example of model behavior, figure 1(a) also shows the probability distribution over this region calculated from output of GFDL-CM4, one of the models participating in CMIP6 [17]. Although models often simulate a more complex probability distribution shape than observed [41], this model and many others exhibit a relatively simple behavior at the extreme tail—an approximate exponential decay range, seen as an approximate straight line in the semi-log plot for values above the cutoff.

While it is clear that models do not capture observed daily rainfall probability distributions in absolute terms, the main quantity of interest for global warming projections of extremes is usually the relative change (i.e. a percent or fractional change) of extreme percentiles respect to a reference historical state [15, 16]. In this case, the benchmark for models to provide useful information is more modest. First, they should simulate the correct shape of the probability distribution extreme tail (for values larger than the range of validity of the power law range, i.e. above the cutoff scale) compared to observations. Second, they should correctly simulate the physical scale that governs the scaling of this tail across regions. Assuming a form of the extreme tail that is close to exponential, and that global warming induced changes in the extreme tail are completely described by a change in parameter δP_L , then the relative change of a extreme daily precipitation percentile P_q (over wet days) is given by:

$$\begin{equation} \frac{\delta P_{q}}{P_{q}}\approx\frac{\delta P_{L}}{P_{L}}, \end{equation} \tag{ 3 }$$

(see supplementary material), a relation we assess below. That is, to leading order a percent change in a extreme percentile is solely given by a percent change in the cutoff scale, and is approximately independent on how the power law range exponent τ_P adjusts.

As a baseline model, previous studies often assume that fractional increases in extreme precipitation percentiles follow a Clausius–Clapeyron (CC) scaling ≈7%/K [15, 42, 43]. However, from moisture budget considerations, the daily precipitation cutoff scale P_L scales with the amplitude of moisture convergence fluctuations within precipitating events [26, 44]. Thus, for a global warming scenario the fractional change of P_L ($\frac{\delta P_{L}}{P_{L}}$) incorporates a thermodynamical contribution given by increases in moisture (following approximately a CC scaling) [16], a dynamical contribution associated with changes in column convergence [10–13, 45] and changes in the covariance between these two terms. Regardless of whether contributions to fractional changes in P_L are dynamical or thermodynamical in origin, this fractional change is the same in usual P or scaled $P^{*} = \frac{P}{P_{L}}$ coordinates. That is $\frac{\delta P_{L}}{P_{L}} = \frac{\delta P^{*}_{L}}{P^{*}_{L}}$ (which can be easily derived from P^∗ definition. Note that $P^{*}_{L} = 1$), which implies that fractional changes in P_L are independent of the absolute value of P_L in current climate.

We illustrate this in the schematics in figure 1(c) showing a postulated change in the probability distribution under global warming given by a fractional increase in P_L , similar to what [26, 46, 47] have reported. Here the fractional change in the 99.9th percentile is the same whether this is assessed over daily precipitation in absolute terms (P, not shown) or over daily precipitation rescaled by the historical cutoff scale (P^∗, figure 1(c)).

Although we have shown that fractional changes in extreme daily precipitation percentiles do not depend on the absolute value of the current climate cutoff scale (for a given location), model performance still depends on an accurate simulation of the shape of the extreme tail. For this reason, we evaluate model representation of the extreme tail of daily precipitation after rescaling (that is, in P^∗ coordinates). To illustrate, figure 1(b) shows the effect of rescaling the observational datasets and GFDL-CM4 and several other CMIP6 models daily precipitation data in figure 1(a) by their cutoff scales on the respective probability distributions. This operation allows to evaluate similarities or differences in the shape of the large daily precipitation probability tail. Now that $P^{*}_{L} = 1$ in all cases, the rescaled probability distributions lie much more closely together, meaning that the shape of their respective extreme tails (above their respective cutoff scales) is similar. This implies that, although estimations of P_L differ in absolute terms between observational datasets and CMIP6 (e.g. figure S2 in supplementary material) (implying also a large difference in absolute probabilities), the same fractional change in P_L (as in figure 1(c)) will yield similar fractional changes in extreme daily precipitation percentiles in all cases.

As outlined in the Introduction, these considerations motivate a two-part test in which we evaluate (i) how well CMIP6 models capture the shape (after rescaling) of the daily precipitation extreme tail compared to observations, discussed in this section and (ii) the ability of CMIP6 models to simulate fractional changes in P_L , discussed in section 5. The latter point cannot be tested directly for global warming conditions, nor we can reliably test against past conditions, due to the insufficient length of quality global gridded datasets (sim20 years). For this reason, we test how well models simulate spatial fractional changes in current climate. While this is not a perfect substitute, it does test fractional changes over a larger range of meteorological conditions than the range projected for future changes. If both these aspects (extreme tail shape and P_L spatial distribution) are well simulated, then together they provide substantial underpinnings for the ability of models to simulate future fractional changes in extremes.

To assess how well models simulate the shape of daily precipitation probability distributions in their extreme tail we introduce a simple metric $e_{\textrm{tail}}$. This metric characterizes the distance between modeled and observed probability distributions above their respective cutoff scales. This metric is defined as

$$\begin{equation} e_{\textrm{tail}} = \int_{P_{L}^{*}}^{\infty}|\textrm{PDF}_{\textrm{model}}(P^{*})-PDF_{\textrm{TRMM-3B42}}(P^{*})|\textrm{d}P^{*}, \end{equation} \tag{ 4 }$$

where $\textrm{PDF}_{\textrm{model}}(P^{*})$ and $\textrm{PDF}_{\textrm{TRMM-3B42}}(P^{*})$ denote model and TRMM-3B42 observational estimates, respectively, of daily precipitation probability distributions rescaled by their respective cutoff scales (as in figure 1(b)). Under this rescaling, the cutoff scale occurs at $P^{*}_{L} = 1$ for both modeled and observed probability distributions. We have tested other definitions for $e_{\textrm{tail}}$ [49–51] with similar results (figure S3 supplementary material). In terms of model evaluation, a small value of $e_{\textrm{tail}}$ means that the shape of the extreme tail in the model closely follows TRMM-3B42 observational estimate, with a value of $e_{\textrm{tail}} = 0$ meaning that the modeled extreme tail perfectly tracks observations. For probability distributions normalized such that they integrate to one for values above $P^{*}_{L}$ the maximum value of $e_{\textrm{tail}}$ is equal to 2, however this is rarely found in practice. The measuring stick we use to assess the performance of the multi-model mean (figure 2(a)) or particular models (figure 2(b)) is to compare against a similarly calculated $e_{\textrm{tail}}$ distance between TRMM-3B42 and GPCP, which gives an estimation of the observational uncertainty. This uncertainty, in principle, could be tightened over North America or Western Europe where there is a dense network of stations (although issues in terms of the mismatch of scales between stations and models would apply), but here we are interested specifically in the near-global case.

Figure 2(a) shows a map of the CMIP6 multi-model mean $e_{\textrm{tail}}$, with stippling denoting regions where the multi-model mean is within observational uncertainty. This map can provide valuable information to decision makers to assess how much to trust model projections of changes in precipitation extremes at the regional level. The values of $e_{\textrm{tail}}$ are small in mid latitudes and modest in the tropics (typical values on the order of 0.03 and 0.2, respectively relative to a theoretical maximum of 2). The observational uncertainty for the rescaled distributions is also modest (figure S4 supplementary material), although larger in the tropics. In contrast with previous studies showing poor model performance in simulating tropical extremes in absolute terms [8, 16], in most of the tropics models perform within observational accuracy. Overall, the CMIP6 multi-model mean performs as well or better than the distance between observational datasets in 67% of the 50° S–50° N latitude band for which TRMM data are available, with exceptions occurring primarily in subtropical regions with little precipitation. At the individual model level, 28 out of 35 models perform better than this measure of observational uncertainty (figure 2(b)), when integrating $e_{\textrm{tail}}$ over 50° S–50° N. Note that the GPCP to TRMM-3B42 difference was chosen to span the typical range between the observational data sets. If we had some basis on which, e.g. GPCP could be eliminated from the range of observational estimates, reducing the comparison range to observational datasets more similar to TRMM-3B42 (figure 2(b), see also figure S1) then fewer (8 models) or no models perform better than the difference between TRMM-3B42 and PERSIANN v1 and CMORPH v1 respectively. However, we are not aware of clear evidence that would support such a reduction in the observational uncertainty. We also emphasize that the differences being evaluated here are the small differences in the rescaled distributions, typified by figure 1(b). We conclude that models generally capture the shape of the extreme upper tail of daily precipitation probability distribution as well, or better, than the extent to which two widely used observational data sets spanning current observational uncertainty agree with each other over most of the globe.

The second aspect we test is models' ability to simulate fractional changes in P_L as a function of spatial region and season. As illustrated in figure 1(a), models and observational datasets do not simulate in general the same magnitude for the cutoff-scale P_L but the shape of the extreme tail, after rescaling, is well simulated (section 3; see figure 1(b) for an illustration). However, regardless of the magnitude of the simulated cutoff-scale, modeled P_L fractional changes will compare favorably to observed fractional changes if modeled and observed P_L spatial patterns are well correlated (see supplementary text for details). So, we test spatial fractional changes by calculating the spatial correlation between modeled and TRMM-3B42 P_L spatial patterns. As before, the benchmark we use to judge models is to compare to P_L pattern correlation between TRMM-3B42 and GPCP, which typify the range among observational estimates.

The CMIP6 multi-model mean and TRMM-3B42 P_L spatial patterns are shown in figures 3(a) and (b) respectively. We note that the CMIP6 multi-model mean P_L pattern (figure 3(a)) tends to be smaller in magnitude compared to the TRMM-3B42 pattern (figure 3(b)), implying that extremes (in absolute terms) tend to be subdued in CMIP6 compared to this observational dataset (as occurs in figure 1(a)). However, despite differences in magnitude (which as long as each grid point has a similar bias are unimportant for fractional changes), the spatial pattern of the CMIP6 multi-model mean of P_L is reasonably correlated with the corresponding GPCP (r = 0.7) and TRMM-3B42 (r = 0.85) patterns at the regional level (figure 3(c)), with even larger correlation coefficients if one is less interested in the details of the regional pattern (r = 0.8 and r = 0.94 with GPCP and TRMM-3B42 zonal averages respectively). This relationship between P_L multi-model mean spatial pattern and TRMM-3B42 is tighter than the extent to which GPCP and TRMM-3B42 agree with each other (r = 0.67 at the regional level, r = 0.86 for zonal averages). This is also the case for most individual models (26 of 35 models correlate comparably or better than the two data sets), and also tends to hold for individual seasons (figure 3(f)).

Zoom In Zoom Out Reset image size

Despite absolute errors in the historical simulated PDFs (compare figures 3(a) and (b)), the good agreement between modeled and observed P_L spatial patterns boost confidence in models ability to simulate future P_L and extreme percentiles fractional changes. Although there are substantial regional variations in τ_P (not shown), the spatial patterns of extreme percentiles largely depend only on the P_L spatial pattern, in agreement with equation (3) (with δP_L and δP_q denoting spatial differences in this context). This is clearly the case in observations and also models, where there is a strong correlation between extreme percentiles P_q (over wet days) and P_L spatial pattern (r = 0.96 for 99.9th percentile P_99.9 in the CMIP6 multi-model mean, figure 3(d)). This is also true for individual models and seasons (figure S5 supplementary material), and also true for future global warming fractional changes in the CMIP6 ensemble (figure 3(e)).

Due caution is necessary in using behavior across spatial regions in current climate as a contribution to inferring confidence in behavior under warming. Both can involve changes in precipitation type, associated meteorological features and other complications [52]. However, the simple dependence of extreme daily precipitation percentiles on a single distribution parameter (P_L ) implies that within the specific context of this regime of the daily precipitation PDF, the effect of P_L changes due to spatial differences (figure 3(d)) map onto similar changes due to warming. Figure 3(f) illustrates this with the example of current and future climate PDFs in the US Midwest (40° N–48° N, 95 °W–80 °W) in the UKESM1-0-LL model. In current climate, the cutoff-scale P_L is 16.0 mm and 99.9th percentile (over wet days) P_99.9 is 78.2 mm (see dots in figure 3(d)) while in the global warming scenario at the end of the century, P_L has increased about 25% to 20 mm (horizontal red arrow in figure 3(d)). This future fractional increase is similar to the current climate difference between the Midwest and Florida (25° N–32° N, 84° W–80° W) in the model (figure 3(d)). Thus for the daily precipitation PDF, the current climate Florida PDF provides an example of future precipitation probabilities in the Midwest through the large event range (figure 3(f)). We underline that we are not advocating location to location comparisons, but rather illustrating that the simplicity of the large-event tail of the PDF suggests usefulness of the range covered in current climate for tests of properties relevant to future changes.

The insight gained here for fractional changes in extremes would be elusive if looking at the daily precipitation extreme tail in absolute terms. Currently, CMIP6 models exhibit a variety of different cutoff-scale magnitudes, with corresponding different extreme percentiles (e.g. figure S2, see also reference [53]), which tend to disagree with observational estimates. Assessments in absolute terms are also complicated because observational estimates themselves differ substantially (e.g. figure 1(a)). However, many studies look at relative measures of change, such as the change in probability of exceeding a threshold given by a particular percentile or return time in the models simulated historical climate [15, 16, 54–60]. Such measures make the implicit assumption that even though the absolute precipitation value corresponding to the model probability measure does not match that in observations, simulations of the fractional changes will nonetheless provide a reasonable estimate. This assumption has been difficult to test, but the results here isolate aspects of this assumption that can be evaluated from current climate. The majority of models simulate a satisfactory representation of the shape of the extreme tail after rescaling by the cutoff scale P_L , with most models providing results that are closer to TRMM-3B42 observational estimate than the current observational uncertainty. While future changes in P_L are not possible to evaluate observationally, we find that models simulate a current climate regional pattern of P_L that is well correlated with observations, and in the case of many models correlates even better that the extent to which GPCP and TRMM-3B42 agree. This suggests that physical processes (column water vapor and convergence spatial distribution within precipitating events) leading to the current spatial distribution of extreme precipitation are reasonably well simulated by the current generation of state-of-the-art models (albeit of different magnitude), and to the extent that model performance simulating current climate is a predictor of their ability to simulate future changes, this boosts confidence that this remains the case under global warming conditions. We note the caveats that this is a current-climate, not an emergent-constraint criterion, and that this test across regions does not guarantee changes will be well simulated for individual regions. Overall, despite issues in model simulations of daily precipitation probability distributions in absolute terms, the reasonable simulation of (i) the shape of the distribution, and (ii) the proportionality of the key scale of the distribution in models to that in observations across a wide range of conditions provide two current-climate measures that supports the usefulness of model projections of fractional changes in extremes.

This work was supported by National Science Foundation AGS-1540518, AGS-1936810 and National Oceanic and Atmospheric Administration NA18OAR4310280. We thank J Meyerson for graphical assistance. We acknowledge the World Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP6. We thank the climate modeling groups for producing and making available their model output, the Earth System Grid Federation (ESGF) for archiving the data and providing access, and the multiple funding agencies who support CMIP6 and ESGF. We also acknowledge high-performance computing support from Cheyenne [61] provided by NCAR's Computational and Information Systems Laboratory, sponsored by the National Science Foundation.

The data that support the findings of this study are openly available at the following URL/DOI: https://esgf-node.llnl.gov/search/cmip6/.