This study uses observed rainfall (daily time series for 1980–2000) and erosivity (mean for 1971–2000) from the PRISM group on a 1/24∘grid (about 4 km). The rainfall data set has been extensively documented and used seewhich have more than a thousand combined citations, while the erosivity data set has been only partially described in the peer-reviewed literature and was obtained courtesy of C. Daly (Oregon State University). The station-based calculations of storm characteristics and erosivity on which the map is based are described by and as well as in a more detailed report to the funding agency . We summarize here the many steps involved.

First, 23 stations with continuous, breakpoint rainfall measurements (which record the time of significant change in rain rate, together with the value of the change) were used to estimate a correction to the maximum 30 min intensity I30 values obtained from discrete 15 min averages (I30)B=1.034⋅(I30)15. Second, for each storm the storm energy E was computed as the sum over all time intervals (covering the storm duration) of the product of rainfall depth and rainfall kinetic energy (an exponential function of intensity; see McGregor et al., 1995). Third, a transfer function was developed for 1842 quality-controlled 15 min stations, linking the storm erosivity values EI30 for each rainy day (calculated as above) to the daily rainfall accumulation P. The transfer function has the form EI30=aPb and the coefficients a and b were computed by least-squares fit using year-round data (r2 values for the regression were reported to be “typically above 0.5”). Note that the regression is based on all storms except those with a return period greater than 100 years – even though snow storms and storms with accumulation of less then half an inch of rain are traditionally not included in the calculation of annual erosivity R, but very intense storms are. Finally, the 1971–2000 mean R was calculated for daily-resolved gauge data for stations co-located or nearby the 15 min stations, according to the formula R=∑EI30≃aN∑n=1N∑d=1365Pdb, where N is the number of years (30 years) and Pd is daily precipitation.

The creation of the gridded erosivity map based on the 1842 point values of R is described in an unpublished EPA report . In general it follows the PRISM method outlined in , but using mean annual precipitation for the period 1961–1990 as guidance in the interpolation in regions with significant topography or near coastlines. In flat, homogeneous areas the interpolation acts instead to smooth out differences which were not deemed physical (but instead due to data issues, such as missing data, difference in measuring equipments or period of records, and random noise). The final R field obtained by this procedure has a spatial resolution of 1/24∘but, for comparison to model data, we regridded it to a resolution of 1/8∘(0.125∘). We refer to this observation-based estimate of R as our “target R”.

Climate change is portrayed in this study as the difference from the historical period and the decades 2079–2099 simulated under a business-as-usual scenario (the representative concentration pathway RCP8.5). This signal is estimated from a 21-model subset of the CMIP5 ensemble , which has been statistically bias corrected and downscaled to a resolution of 1/8∘ according to the method of . This data set, known as the Bias-Correction Constructed Analog version 2 data (BCCAv2), provides the daily precipitation data that is used to calculate present-day and future erosivity in the contiguous United States in the historical and RCP8.5 simulations. The steps that transform CMIP5 output into BCCAv2 are as follows: GCM output is first regridded to a 2∘ resolution and then a quantile matching procedure adjusts the rainfall values to match the observed climatology, then a constructed analog technique

An analog for a given daily weather pattern can be constructed by combining patterns from a library of previously observed patterns; using a companion library of high-resolution patterns creates the high-resolution downscaled estimate.

Figures  and  present, for observations and the model ensemble, respectively, a quick overview of the rainfall characteristics that are most likely to influence erosion: mean annual accumulation (in mm day-1), mean daily intensity (defined as the average rainfall accumulation on days when it rains), the modified Fournier index, and the burst factor. The observed characteristics are well known. The western coast ranges receive plenty of precipitation but otherwise the western states are quite dry; east of the Rockies precipitation is generally more abundant, especially in the southeast and in places along the Appalachian Range. The same pattern is repeated and intensified in the intensity field, with especially high intensities in the southeast. Unsurprisingly, both the burst factor B and the modified Fournier index F show a mixture of the accumulation and intensity patterns.

The model ensemble reproduces the annual accumulation and F very well, as expected from the fact that we are using bias-corrected data sets. Instead, the quantities that are based on daily, rather than monthly, data are not at all well captured. Maximum values in the intensity field are muted compared to observations by a factor of at least 4 (changing the threshold for a rainy day from 0 to 1 mm day-1 does not solve the problem, not shown), and the low-intensity also affects the burst factor, especially in the southeast. This can be explained in part by the lower resolution of the model data, but it seems that the bias correction does not fully address the notorious drizzle bias that plagues GCMs .

Because of concern with the unresolved bias in daily rainfall intensities, we proceed by estimating erosivity and erosivity changes using two predictors. One predictor is rainfall daily intensity and the method by which we estimate R (method 1) is more closely related to the original erosivity calculation (but is affected by the intensity bias). The second predictor is F, which is unaffected by the bias. For F, the formulation of the erosivity calculation is less well motivated and the relationship is more uncertain; thus, we try both a linear (method 2) and a non-linear (method 3) relationship between F and R.

Figure a shows the target (observed) 1970–2000 annual erosivity. Values are maximum along the gulf coast and Florida, and decline further to the north. Erosivity in the western US is extremely low, aside from the mountain regions of the west coast (note, though, that snow might be contaminating the observational estimate in these areas). Unsurprisingly, the R pattern is reminiscent of both annual mean P and intensity, but some details are different, such as the very high values of R in southern Florida. This pattern differs from that of erosion itself, as steep topography favors erosion in most of the western US and land-use practices enhance or mitigate the effect of erosive rainfall .

To compare climate models to observations, we would ideally calculate the model erosivity using the same formula that was used in the making of the observational data set. Unfortunately, this is not strictly possible. The main reason is that the observed coefficients a and b were not reported in the literature. A secondary reason is that we use gridded values of precipitation and, even though the grid of the downscaled data set is quite fine, coefficients estimated from gauge values may not necessarily be appropriate for a gridded data set. Still, we want to calculate R from daily rainfall using the same functional form and, in order to do that, we need an alternative way to estimate the coefficients a and b that do not require temporally resolved values of R. Our choice was to exploit the fine spatial resolution of the gridded R data set, and use that in place of the missing temporal resolution.

For every 1∘ × 1∘ square, we take the mean 1970–2000 erosivity and daily time series for the period 1980–2000 (PRISM daily rainfall is not available prior to 1980) and use the 625 vertex points of the 4 km resolution data sets that are included in this subdomain to determine the coefficients a and b in the formula R≃aN∑n=1N∑d=1365Pdb. This provides us with a 1∘ × 1∘ gridded map of a and b, which we then interpolate linearly to the same resolution as the PRISM data and smooth to the 1/8∘ resolution of the models. In the fitting routine we require that a be positive and b be greater or equal to 1 a choice consistent with the parameter values reported in. We want the formula linking rainfall and erosivity to be sufficiently robust to produce a good fit to the target R when applied to rainfall at somewhat different resolution, so we check its performance on the same observational rainfall data, but regridded at the coarser resolution of the downscaled data, 1/8∘. Although this check is not using fully independent data, using data at coarser resolution as input does introduce some degree of separation from the data used in the fit. In the eastern domain, grid points for which the fitting routine selected outlier values (very different from nearby points, but without a clear physical feature to justify the sharp gradient) yielded large errors when the regression was applied to the coarser data. To partially avoid this, we imposed as a search criterion in the least-square fitting procedure that the regression coefficients were bound to be close to the large-scale values over homogeneous areas, and outlier values of the regression coefficients were discarded. This is consistent with the choice made by the PRISM team to use a smooth interpolation in the eastern US.

Figure b compares to the target map the erosivity obtained by regression on the 1980–2000 1/8∘ daily PRISM rainfall. Figure  shows maps for the a and b coefficients. The exponent b is close to 1 over the eastern part of the US, but it rises above 1.5 in the west. The multiplicative constant a is most reminiscent of the R and P fields, with values close to 0 in the west (which offsets for the quadratic term in precipitation) and maximum values in the southeast. Figure  also shows the absolute and relative error between the two R estimates (between Fig. a and b). Overall, the regression produces a good match to the original field, but a close look reveals that there are some regional biases. In particular, the west coast has much lower erosivity in the target field than that calculated by the regression; in fact, when looking at the relative error, most of the west appears to be problematic: the errors are small scale and patchy, so the regional bias is not as strong as the local bias, but there is an overall tendency for the regression model to overpredict erosivity.

An increase in the intensity of rainfall, especially in the frequency of the most extreme events, is a robust expectation for future climate change, supported by theory e.g.,, modeling e.g.,, and observations e.g.,. The effect of such changes on erosivity can only be assessed by using daily rainfall data. For this reason, we have calculated erosivity from daily precipitation, as shown above. Nevertheless, the biases in daily intensity that persist even in this bias-corrected data set shed doubt on the ability of this method to accurately capture changes in erosivity and, therefore, we include estimates of erosivity obtained from monthly-averaged rainfall data – which is brought close to observations by the bias correction.

Specifically, we calculate R from the modified Fournier index F, which is one of the most widely used indices, with an established literature for the US. used high-frequency gauge data from 132 stations to derive a relationship between F and R for the continental US; the relationship is described by two equations: R=0.7397F1.847, and R=95.77-6.081F+0.4770F2, where Eq. () holds for F<55 mm and Eq. () holds for F>55 mm and F is defined as in Eq. ()

Note that these are the equations reported in the figures. The main text in Reinard and Freimund contains a graphical error in Eq. (), and subsequent citations in and contain a typo in Eq. ().

Other scholars have found different relationships for other locations see the many examples reported in, both linear and of the form R=aFb, and the literature repeats Arnoldus' (1980) recommendation that relations obtained using the modified Fournier index should be applied only to locations within homogeneous climatic regions. Therefore, we also apply the same method described above for daily precipitation and estimate local R to F relationships for points in each 1∘× 1∘ square, allowing for both a linear and a power law relationship. Thus, we have three ways of estimating R from F.

Applying Eqs. () and () to the gridded PRISM data leads to overestimating the erosivity everywhere in the US, with the only exceptions of eastern Texas and Florida (Fig. 3c). The fact that the error is of the same order of magnitude as R itself is not inconsistent with the errors shown by Reinard and Freimund (see their Fig. 3b), but the positive bias in the estimate is harder to explain. Nonetheless, we note that the relationship between F and R in the literature for many world regions is often linear and, as will be shown below, the local-regression methods suggest a much weaker dependence on F than the Reinard and Freimund formulas in all but a few mountainous regions of the US.

To determine the F to R relationship from local regressions, we have tried several functional forms: a linear relationship, a power law, the sum of the two, with and without subjective bounds on the coefficients. Here, we present results from the linear and power law functional forms. Nonetheless, we imposed an additional regional constraint to avoid small-scale noise in the coefficient fields: the noise was deemed unphysical and it gave rise to large local errors when the fit was applied to the 1/8∘ PRISM data. In particular, we have imposed that the relationship between F and R stays close to linear in the southeast. This trial-and-error fit without a posteriori verification would be problematic if we were to claim that it proves a physical connection between F and R. But our goal is much more modest: we assume that a relationship between F and R exists, as established in the literature, and simply obtain a scaling between the two quantities that performs better than the Reinard and Freimund formulas.

The end result of the local regressions are presented in Figs. 3d, e, , and . The case of a linear fit provides the best fit (Fig. ): the errors are small and very patchy, indicating that there are weak regional biases in the estimate of the erosivity field. The power law produces errors of larger magnitude, and their severity is worst in the western US. Here, the power law coefficient sets a stronger-than-linear relationship between F and R, and the end result is a regression that over-estimates erosivity. In both the linear and non-linear cases, the local regression allows us to calculate erosivity for the entire domain, including the west coast, but errors there are larger. Moreover, the issue of including snow in the erosivity calculation is not properly addressed by any of these methods, although changes from snow to rain might indeed be an important factor in determining future changes in erosivity.

R is better approximated using F than with the daily precipitation. This is somewhat surprising, given that it is the daily precipitation that goes into the original R calculation. On the other hand, it might be that the use of annual mean rainfall in the interpolation of R from point measurements to the gridded field introduces a bias in R, making the field more similar to rainfall accumulation than it would otherwise be.

Another possible choice for building the regression between rainfall (either at daily frequency or aggregated into F) and erosivity would be to base it on larger regions that encompass broader variations in the relevant fields than what is seen in a 1∘ × 1∘ square. Taking inspiration from the regional clusters of Palecki et al. (2005) – which select regions with fairly homogeneous mean storm characteristics – we have defined eight regions and regressed rainfall and rainfall characteristics across the grid points within each region. The regressions are more biased in this case than in the case of the local regressions, independently of the choice of rainfall variable. This is as expected, given that we use much fewer parameters, and again underscores the uncertainty in the estimates or erosivity that are based on rainfall accumulation. In the rest of the study, we proceed using the estimates obtained by local regressions, mindful that these regressions should be interpreted more as a simple scaling than as a robust physical relationship.

In this subsection, we describe the performance of the climate models in reproducing the observational estimates of rainfall erosivity R. All three methods of estimating R are evaluated.

When applied to the model ensemble daily rainfall data for the 1980–2000 period in the historical simulations, Eq. () with the coefficients of Fig.  yields a very reasonable pattern of erosivity (Fig. a), but with obvious biases. The absolute error (Fig. c) is large along the west coast, but comparable to what is obtained by using observed daily rainfall for the fit (see Fig. ), and also large in the southeast, where erosivity is overestimated by the regression by 10–40 %. This is somewhat surprising, given that in the same region the models underestimate rainfall intensity; yet, in this region the relationship between rainfall and R is very close to linear, so it is possible that the positive bias in the number of rainy days (given that the mean has been corrected, low intensity must go together with high frequency) translates to the overestimation of erosivity. While the relative error (Fig. d) in the eastern US is only prominent in Florida and coastal regions, it is large everywhere in the west. Nonetheless, here the error pattern is as patchy as was seen for observations in Fig. , so that the overall regional bias is not as severe. The western US is also the region for which there is more scatter across the model ensemble, as signified by the intra-ensemble coefficient of variation (the standard deviation of the inter-model spread, divided by the multi-model mean erosivity, Fig. b).

When we calculate erosivity from the models' F, we obtain somewhat better results (compare Fig. a to Figs. d and Fig. a to Fig. e). Similar to what is seen for the observations, the linear relationship captures R with the least bias (Fig. c, d), while the non-linear relationship (Fig. c, d) induces a worst overestimation of erosivity in coastal regions of the gulf and the southeast and in the western US. Erosivity in the Great Plains tends instead to be underestimated by the models. The biases in model-estimated erosivity are all qualitatively quite independent of the formula used to calculate the estimates, but they are smallest when R is linked linearly to F, as one would expect from the fact that monthly rainfall (and thus F) was the target of the bias correction.

Figure  shows the end-of-century projections for annual rainfall accumulation, mean daily intensity, the burst factor, and the modified Fournier index. All four fields show overall future increases over most of the North American continent. The annual precipitation signal in the northern part of the domain (Canada and US) is clearly positive, with anomalies replicated across the ensemble (stippling indicates at least two thirds of the models agreeing on the sign of the anomalies), but in the southern part of the domain the anomalies are weak and not very robust. Although barely significant (see the limited extent of the stippling), there is general drying in Mexico and Texas, extending into the southwest US, and more positive anomalies in the eastern gulf and Florida. The pattern of yearly anomalies shown here is consistent with the expectations reported in the literature for both CMIP3 and CMIP5 models see for example. It can be interpreted as the sum of a weak US-wide summertime drying and the wintertime pattern of increase in rainfall north of 40∘ N and decrease in precipitation in the southwest and Texas. Areas close to the wintertime zero-anomaly line are areas where the climate signal is projected to be small compared to natural variability, and model disagreement is not problematic. But in other areas, such as northern California and Texas, model disagreement is indicative of true uncertainty in the response . The positive weak anomalies in Arizona and central Mexico are not expected from the literature. They are a consequence of the bias-correction method: the projections for the same 21 runs in the original CMIP5 ensemble show that the negative rainfall anomalies do extend farther north than their downscaled counterparts. In general, the downscaling seems to dampen the negative anomalies. This is visible for both accumulation and the Fournier index when we compare the downscaled anomalies in Fig.  with the anomalies for the corresponding CMIP5 models at the native resolution (see Fig. ).

The intensity signal see also presents a similar pattern of stronger positive signal in the north and along the coasts, but in this case the region of consistency across models is more expansive across the US, with only the Colorado Plateau and the southern Great Plains not showing a robust increase in intensity (anomalies are nonetheless positive). The pattern can be interpreted as the superposition of the general (thermodynamic) tendency for more intense rainfall in a moister atmosphere , dampened in the middle of the continent by circulation-driven dry anomalies in summertime . Where these two effects are of comparable magnitude in the multi-model mean, the changes in intensity are less robust.

The modified Fournier index is more robustly positive than the annual accumulation, from the east coast to the Mississippi and from the west coast across the Rockies. The southern Great Plains show negative change in F, but the anomalies are inconsistent across the ensemble. Also, F increases at the gulf coast even though mean annual P declines here, a difference caused by an increase in monthly intensity that dominates over the mean drying. The strengthening of the robust pattern of anomalies in F compared to annual rainfall is an indication that the seasonal cycle of rainfall will be getting more peaked, i.e., more rain in the rainy season . As noted above, this behavior is more pronounced in the downscaled anomalies, while the same models at the native resolution would indicate negative anomalies extending over Texas, Louisiana, and Mexico (Fig. ).

Changes in the burst factor B are the most uniformly positive. B values reflect changes in both mean rainfall and rainfall extremes. This is especially true given that we have calculated the effective rainfall entering B using the observation-derived threshold of half an inch, which corresponds to higher rainfall percentiles in the models. Moreover, the maximum daily rainfall for the month enters the calculation directly. The domain-wide robust increase in B – extending farther than either mean rainfall or mean intensity – is consistent with the expectation that the intensity of extreme rainfall events will increase more than average rainfall and average intensity.

Given the pattern of robustness and uncertainty in the predictors of erosivity (that is, in mean annual rainfall, daily intensity, F and B), we can expect that erosivity changes will be generally positive but quite uncertain in the southern and central US and dependent on the method used to estimate erosivity from rainfall. In the following section we look at such changes in a more quantitative way.

We have noted how the downscaling method introduces its own uncertainty in the projected rainfall changes used in this study; another source of uncertainty comes from the scatter across the models itself and can be partially addressed by looking at the sensitivity of our results to the size of the climate model ensemble. Figure  compares the mean signal in annual mean accumulation and in the modified Fournier index when we use one realization of the 21 models used for downscaling from the native resolution CMIP5 ensemble to the projection from 76 runs from 38 models in the CMIP5 ensemble. The two patterns are very similar but they differ in details, such as on the sign of the anomalies over Louisiana. Finally, while here we focus on the climate change signal, it is important to remember that natural variability in rainfall is large over North America and capable of masking forced trends even over multi-decadal periods .

Keeping in mind all these sources of uncertainty in rainfall changes, we now proceed to estimate the projected changes in rainfall erosivity. These inherit all the uncertainty in the climate projections and compound those with the uncertainty in the relationship between accumulated rainfall and erosivity.

Figures , , and show the multi-model mean difference in erosivity between 1980–2000 and 2080–2100, both in erosivity units (MJ mm ha-1 h-1 year-1, top panels) and as percentages of the 20th century values (bottom panels), alongside the intra-ensemble scatter in the same quantities. In Fig. , R is estimated from the daily precipitation values (method 1), while in Figs.  and it is estimated from F, either through a linear (method 2) or through a non-linear regression (method 3), as described in Sect. .

As would be expected from the rainfall changes, the models suggest an overall increase in erosivity across the US, but the details of the projections change quite significantly depending on the method used to estimate R. These patterns match very closely the patterns of the erosivity predictors (rainfall changes and F), indicating that non-linearities have a weak role in setting the erosivity anomalies. For the same reason, the inter-model agreement pattern for erosivity matches that for the climatic predictors. When daily precipitation is used, the increase in erosivity covers the west, east and north-central US, with negative anomalies in Texas and north into the Colorado Plateau. Conversely, the estimates obtained from the modified Fournier index would indicate positive anomalies in Texas (as in most of the domain) and only a region of very weak decrease in erosivity further north around Colorado and Kansas. For R using daily data the changes are primarily driven by the changes in mean accumulation. For R calculated via F, changes in monthly intensity are more important and there are regions where, by this estimate, R increases even as accumulation declines, because an increase in monthly intensity dominates.

The fact that the pattern of change in rainfall accumulation differs from that of the modified Fournier index leads to increased uncertainty in erosivity changes, compared to what has been previously reported in the literature as a consequence of the climate uncertainty. For example, in eastern Texas, estimates based on F suggest a robust increase of erosivity across a majority of climate models (even if the inter-model spread is large, Fig. b, there is strong agreement on the sign of the change, Fig. a), but estimates based on daily precipitation accumulation paint the opposite picture. For this region, the uncertainty derived from the rainfall–erosivity relationship is greater than the uncertainty derived from scatter across climate projections.

In spite of the overall uncertainties and biases in the model estimates, it is nonetheless clear that the expectation of increased erosivity for most of the US is supported by this analysis and that strategies to prevent soil loss will need to be implemented in a variety of environmental conditions, from the mountains of the west to the northern Great Plains, to the Appalachian Range and the Eastern Seaboard.