Projected changes in US rainfall erosivity

Downscaled rainfall projections from 21 climate models from the CMIP5 archive are used to estimate future changes in rainfall erosivity in the continental Unites States. To estimate erosivity from rainfall in the absence of sub-hourly data, we have used both daily rainfall values and the modified Fournier index–which is based on monthly rainfall accumulation–and derived the scaling 5 relationship between rainfall and erosivity from observational estimates of both. The expectation of overall increase in erosivity is confirmed by these calculations, but a quantitative assessment is marred by large uncertainties. Specifically, the uncertainty in the method of estimation of erosivity is more consequential than that deriving from the spread in climate simulations, and leads to changes of uncertain sign in parts of the Southwest and Texas. 10 We suggest that progress can be made by establishing a more reliable functional relationship between daily rainfall and erosivity.


Introduction
Soil erosion has significant consequences for the productivity of land and the health of the riverine systems that receive the erosion materials.In an eroded soil, nutrients are lost, the effective root depth decreases, and runoff increases.In rivers, suspended sediments increase turbidity and transport pollutants, reducing the health of algae and fish and degrading the quality of drinking water, while excess sediments affect spawning fish, lower reservoir storage, and interfere with navigation.Erosion by water already causes about 55 % of total global erosion (Bridges and Oldeman, 1999) threatening sustainability and the productive capacity of agriculture (Yang et al., 2003) A knowledgeable forecast of how ero-sion will change in the coming decades is thus necessary to plan for land stewardship and ecosystem preservation and it is made more urgent and complex by the general expectation that increases in rainfall intensity under global warming (Trenberth et al., 2003) will exacerbate erosion.There is widespread agreement (Soil and Water Conservation Society, 2003) on an increased risk of soil erosion under climate change, but the potential damage needs to be assessed against other drivers such as land cover (e.g., Nearing et al., 2005) and with better understanding of the uncertainty due to climate scenarios, spread in the climate projections, and methods of translating large-scale climate anomalies into forcings for erosion models (e.g., Mullan et al., 2012;Zhang, 2007;Zhang et al., 2011).
The forces that affect erosion are summarized in a widely used formula (the universal soil loss equation, USLE) and its derivates (RUSLE1 and RUSLE2) that were developed by the US Department of Agriculture.The USLE was developed for cropland in the early 1960s (Wischmeier and Smith 1965) and was later extended to other land uses (Wischmeier andSmith 1978, Dissmeyer andFoster 1980).Between the mid 1990s and the early 2000s, RUSLE1 and RUSLE2 expanded on the index-based, empirically derived USLE by use of hybrid models that add process-based equations, making these rules land-use independent and able to compute deposition.The empirical relationship at the basis of USLE describes soil loss as proportional to an erosivity factor (R, which describes the compound effects of the rainfall events that cause erosion), a soil erodibility factor (K, which takes into account the soil composition), a topographic factor (LS, basically the slope of the terrain), and two other quantities that describe the effect of vegetation, land use, and land management (C, the cover-management factor, and P, the support practices factor) (Wischmeier andSmith, 1965, 1978;Dissmeyer and Foster, 1980).
Published by Copernicus Publications on behalf of the European Geosciences Union.

M. Biasutti and R. Seager: Erosivity projections
While climate, vegetation, and soils can all be coupled in complex ways (for example, excessive heat and drought can kill the vegetation cover and induce more soil loss, even if rainfall erosivity declines; see Nearing et al., 2004), we can start to address the role of climate change on soil erosion by first isolating the changes in the erosivity of rainfall while setting aside the impacts of changes in vegetation, soil composition and topography, land use and management.That is, while we address changes in the direct climatic driver of erosion, namely erosive rainfall, future changes in erosion will depend on local susceptibility and thus on geography, land use and management, and the health of the vegetation cover.Addressing erosion changes is beyond the scope of this work.
The erosivity factor R is defined as the mean annual sum of erosivity from individual storms with effective rainfall 1 , i.e., rainfall greater than 0.5 in.(about 12 mm).In turn storm erosivity is computed (Angel et al., 2005) as the product of the rainstorm energy (E, an empirically derived function of intensity and depth of rainfall) and the maximum 30 min rainfall intensity for the storm (I 30 ): R = EI 30 . (1) The unit land plot in the USLE erosion calculation is assumed to be only about 20 m in length; accordingly, the formula for erosivity assumes that point-wise measurements of rainfall, i.e., rain gauge data, are used.The requirement that we know the maximum 30 min rainfall intensity adds the need for continuous (or at least very high frequency) measurements.These are stringent requirements that have made the calculation of erosivity a difficult task even in data-rich regions such as the continental US (Angel et al., 2005).In response, many approximations to the original formula have been developed in the literature.Typically, point-wise erosivity calculated from the original definition is used to derive an approximate formula that uses temporally averaged rainfall data (from 30 min, to annual averages), and the approximation is used for locations with sparser data.A field of erosivity values is then extrapolated from the available data points using other physical variables as added predictors (e.g., topography, annual mean rainfall, average rainfall intensity, or more derived quantities such as the modified Fournier index and the burst factor that will be introduced below).
The erosivity map for the continental United States that is in use at the Environmental Protection Agency was constructed in such a way (Hollinger et al., 2002;Daly and Taylor, 2002;Angel et al., 2005).Annual erosivity could be obtained from the original formula for 1505 stations with 15 min resolution.When aggregated via a cluster analysis into nine climatic regions, this erosivity data set was used as the basis for the storm characteristics and erosivity trends assessment published by Palecki et al. (2005) and Angel et al. (2005).Nevertheless, missing data and inhomogeneities in the length of record made the original R-values unfit as the 1 Snow is not considered an erosive event.
basis of a continental erosivity map.Instead, they were used to statistically determine an approximation for erosivity that used only daily data.The full daily precipitation time series in stations collocated or nearby the high-frequency gauges were then used to create consistent and complete time series of annual erosivity for the period 1971-2000.Finally, the PRISM spatial distribution method of Daly and Taylor (2002) was used to draw a complete erosivity map: the method combines the available point values for estimated R with the gridded mean annual precipitation at 2.5 arcmin (4 km) resolution, under the assumption that additional information on the spatial field structure of R can be inferred from the spatial structure of P .As explained below, the erosivity map thus produced is used in this study to build our model-based estimate of erosivity and erosivity changes.
Long and complete daily-resolved rainfall time series are still not too common, and many approximate estimates of erosivity that use rainfall records at much coarser temporal resolution have been pursued.Among the indices that use monthly-averaged rainfall data, we focus on the modified Fournier index, F , introduced by Fournier (1960) and Arnoldus (1980) as where P j is the monthly rainfall depth for the calendar month j .In this formulation, a month of intense rainfall is seen to create more erosion than the same rainfall spread over a longer time.Quantitatively, though, the link between F and erosivity is not well constrained by observations.Renard and Freimund (1994) proposed a power law for the continental US, but noticed that a quadratic formulation was a better fit for high F values; moreover, more recent estimates of the power-law exponent for the same region but using different data have found different values (USDA Agricultural Research Service, 2013).Outside the US, other formulations altogether have also been proposed (see for example Table 3 in Schönbrodt-Stitt et al., 2013).
Although less used for the US, we also present results using the burst factor (Smithen and Schulze, 1982), as it provides the example of an index that gives more weight to the most intense, and erosive, events.This index, alongside monthly rainfall accumulation, takes as input the maximum daily rainfall for the month.In our formulation, it is calculated as where M j is the maximum daily rainfall in month j and P e j approximates the effective rainfall as the sum of all daily rainfall with accumulation larger than 12 mm.The same requirements that make estimating erosivity from observations a complex task make it impossible when using climate models.Even the models with the finest resolution compute values for rainfall that are meant to represent averages over hundreds of square kilometers.Moreover, even though models use time steps that are less than 30 min, practical storage limitations restrict the temporal resolution archived to 3-hourly at the finest and, more commonly, 6hourly or daily.Indeed, the volume of data involved in looking at decades of climate simulated by dozens of climate models makes even daily data a challenge to work with.How can we then gain an estimate of the projected changes in erosivity and potential soil loss for the US?Some past studies have circumvented these restrictions by applying idealized changes in precipitation amounts and intensity (Dabney et al., 2012, for example), thus providing results more akin to a sensitivity assessment than a full projection (see also Pruski and Nearing, 2002).Others have looked at changes in F , as estimated from monthly mean precipitation changes in a handful of GCMs (global circulation models) (Nearing et al., 2004;Segura et al., 2014), and used the published relationship between R and F to deduce the change in erosivity.Others (e.g., Zhang et al., 2012) have temporally downscaled the GCM monthly output (by linking monthly rainfall totals to the transition probability between dry and wet days and using a weather generator to create daily time series) and used the downscaled precipitation in an erosion model to directly provide estimates of runoff and soil loss.
In contrast to the above studies, this study uses the daily rainfall output from 21 CMIP5 (Coupled Model Intercomparison Project Phase 5; Taylor et al., 2012) models, biascorrected and downscaled to 1/8 • resolution (the BCCDA data set Maurer et al., 2007Maurer et al., , 2014)), to estimate changes in erosivity in the continental US under the RCP8.5 (representative concentration pathway) emission scenario.Our approach is to duplicate for the models what the literature has done to estimate R where observations of high-frequency rainfall were inadequate.Namely, we estimate the parameters of a statistical relationship between erosivity and daily precipitation from the available data, i.e., from observations of the 1980-2000 period, and we then use the same relationship to relate the GCMs' future projections of daily and monthly rainfall to future projections of erosivity changes.Regrettably, the erosivity data available to us is limited to the climatological values; we do not have access to the full time series of erosivity.Ideally, we would estimate the relationship between daily rainfall and erosivity over a subset of the full time series, and use the remaining independent data to validate our results.We acknowledge that a lack of independent validation is a serious drawback.Thus, to overcome this limitation, we also estimate change in erosivity using the F and BF predictors as established in the literature, and we present the full range of uncertainty in our estimates.
The details of the methods used are summarized in Sect. 2. Changes in the characteristics of rainfall and in other predictors for erosivity are presented in Sect.3.1, followed by changes in erosivity itself, estimated by different methods, in Sect.3.2.Section 4 summarizes and discusses the remaining challenges.

Overview
This study uses observed rainfall (daily time series for  and erosivity (mean for 1971-2000) from the PRISM group (Daly and Taylor, 2002) on a 1/24 • grid (about 4 km).The rainfall data set has been extensively documented and used (see Daly et al., 2002Daly et al., , 2008, , which 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 Angel et al. (2005) and Palecki et al. (2005) as well as in a more detailed report to the funding agency (Hollinger et al., 2002).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 I 30 values obtained from discrete 15 min averages (I 30 ) B = 1.034 • (I 30 ) 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 EI 30 for each rainy day (calculated as above) to the daily rainfall accumulation P .The transfer function has the form EI 30 = aP b and the coefficients a and b were computed by least-squares fit using year-round data (r 2 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 where N is the number of years (30 years) and P d 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 (Daly and Taylor, 2002).In general it follows the PRISM method outlined in Daly et al. (2002), 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 (Taylor et al., 2012), which has been statistically bias corrected and downscaled to a resolution of 1/8 • according to the method of Maurer et al. (2007Maurer et al. ( , 2010Maurer et al. ( , 2014)).This data set, known as the Bias-Correction Constructed Analog version 2 data (BC-CAv2), 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 technique2 is used to downscale from coarse to fine resolution.Finally, another step was necessary to correct a dry bias that was reintroduced by the constructed analog downscaling method: a ratio scaling factor was uniformly applied to correct the annual mean precipitation to again match observation.Note that this rescaling is not expected to correct biases in higher moments of the rainfall distribution.Details can be found in Maurer et al. (2007Maurer et al. ( , 2010)).

Rainfall characteristics
Figures 1 and 2 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 (Dai, 2006).
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.

Erosivity calculations using observed daily rainfall (method 1)
Figure 3a 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 (Judson and Ritter, 1964) and land-use practices enhance or mitigate the effect of erosive rainfall (Gyssels et al., 2005;Montgomery, 2007).
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.d .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 Hollinger et al., 2002).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 3b compares to the target map the erosivity obtained by regression on the 1980-2000 1/8 • daily PRISM rainfall.Figure 4 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 4 also shows the absolute and relative error between the two R estimates (between Fig. 3a 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.

Erosivity calculations using the observed modified
Fournier index F (methods 2 and 3) 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., Trenberth, 1999;O'Gorman and Schneider, 2009), modeling (e.g., O'Gorman, 2012; Tebaldi et al., 2006;Sillmann et al., 2013), and observations (e.g., Alexander et al., 2006;Lenderink et al., 2011;Asadieh and Krakauer, 2015).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 biascorrected 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.Renard and Freimund (1994) 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: where Eq. ( 5) holds for F < 55 mm and Eq. ( 6) holds for F > 55 mm and F is defined as in Eq. (2)3 .Reinard and Freimund found that neither relationship could be applied to the stations of Washington, Oregon, and California because winter snow was coming into the calculation of F , but should not contribute to R. Their criterion to exclude such stations was that no month between October and April should contribute more than 15 % to the annual accumulation.We have followed the same, quite arbitrary, criterion.Other scholars have found different relationships for other locations (see the many examples reported in Schönbrodt-Stitt et al., 2013), both linear and of the form R = aF b , 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. ( 5) and ( 6) 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 localregression 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, 5, and 6.The case of a linear fit provides the best fit (Fig. 5): 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 overestimates 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.

Erosivity calculations using historical model data
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. ( 4) with the coefficients of Fig. 4 yields a very reasonable pattern of erosivity (Fig. 7a), but with obvious biases.The absolute error (Fig. 7c) is large along the west coast, but comparable to what is obtained by using observed daily rainfall for the fit (see Fig. 4), 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. 7d) 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. 4, 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. 7b).When we calculate erosivity from the models' F , we obtain somewhat better results (compare Fig. 8a to Figs. 3d and Fig. 9a to Fig. 3e).Similar to what is seen for the observations, the linear relationship captures R with the least bias (Fig. 8c, d), while the non-linear relationship (Fig. 9c, 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.

Rainfall
Figure 10 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 Seager et al., 2007Seager et al., , 2014;;Biasutti et al., 2011;Christensen et al., 2013;Maloney et al., 2014).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 zeroanomaly 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 (Maloney et al., 2014).
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. 10 with the anomalies for the corresponding CMIP5 models at the native resolution (see Fig. 11).
The intensity signal (see also Wuebbles et al., 2014) 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 (O'Gorman and Schneider, 2009), dampened in the middle of the continent by circulation-driven dry anomalies in summertime (Maloney et al., 2014).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 (Chou and Lan, 2012).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. 11).
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 Bextending 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 11 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 (Deser et al., 2010).

Erosivity
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.12, 13, and 14 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. 12, R is estimated from the daily precipitation values (method 1), while in Figs. 13 and 14 it is estimated from F , either through a linear (method 2) or through a non-linear regression (method 3), as described in Sect. 2.
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 precip- itation 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 intermodel spread is large, Fig. 13b, there is strong agreement on the sign of the change, Fig. 13a), but estimates based on daily precipitation accumulation paint the opposite picture.For this region, the uncertainty derived from the rainfallerosivity 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.

Conclusions
Changes in soil erosion have many drivers and a full evaluation of the effect of climate change on soil loss and sediment deposition can only be achieved by a truly integrated assessment that takes into consideration changes in vegetation and land use.At the moment, the coupling of climate, vegetation, and landscape dynamics at the scales relevant for erosion is not yet achievable, and an assessment of erosion risks under global warming must, by necessity, take a simplified approach.Here, we focus on the one driver of erosion most directly related to climate: rainfall erosivity.
The purpose of this study was to provide an assessment of the projected changes in erosivity for the continental United States, based on the most up-to-date projections of climate change.We have used downscaled rainfall projections from 21 climate models from the CMIP5 archive to estimate future erosivity changes.To estimate erosivity, we have used both daily rainfall values and the modified Fournier index -which is based on monthly rainfall accumulation -and derived the scaling relationship between rainfall and erosivity from observational estimates of both.While daily rainfall maps are available for the period 1980-2000, only the mean value for the 1971-2000 is available for erosivity.Such a dearth of data hinders a robust estimate of the relationship between rainfall and erosivity and introduces another large source of uncertainty in the estimate of future changes, comparable to that due to the spread in climate projections across the CMIP5 ensemble.
Overall, we confirm that the general expectation of worsening erosivity under climate change is correct, but we suggest that previous estimates have been too confident for large swaths of the United States, especially in the south and in the interior.The pattern of erosivity change estimated in our study from CMIP5 changes in F is consistent with what was found by Segura et al. (2014) using just three models (with three scenarios each) from the older CMIP3 archive.But the pattern changes when we estimate erosivity with a different method: when erosivity is estimated from daily rainfall values, the multi-model mean indicates a decrease in Texas and the southern Great Plains, while an increase in erosivity is simulated by a majority of models when the estimate is done using monthly precipitation.Within this region there are locations where the climate signal itself is not very robust (some climate models project wetting, others drying) but there are also locations where the climate models agree on the rainfall changes.There, uncertainty in the sign of the erosivity changes is solely a consequence of our poor knowledge of how best to link erosivity to rainfall accumulation.Thus, we conclude that, for some regions, uncertainty in the method of estimation of R can be more consequential than uncertainty derived from the spread in climate simulations.
A more quantitative assessment will remain a challenge, hindered both by model deficiencies and by lack of complete erosivity records from observations.Ameliorating the biases in the representation of rainfall intensity in climate models is bound to help, and the task is already on the modelers' agenda, as is the need for reducing uncertainty in the projections of regional climate change.But the uncertainty stemming from the relationship between rainfall and erosivity could be ameliorated without having to wait for progress in climate models.Progress could be made by recalculating erosivity time series from high-quality gauges, and highfrequency, high-resolution satellite estimates of rainfall, so that a more detailed and more reliable record of the connection between daily rainfall and erosivity can be established for the entire US.Outside the Unites States (and few scattered localities), the lack of high-frequency, high-resolution data and a worse uncertainty in the relationship between rainfall characteristics and erosion potential (van Dijk et al., 2002) combine to make the challenge even more formidable.

Figure 1 .
Figure 1.Observed rainfall characteristics for 1980-2000 in the PRISM data set (mm day −1 ).(a) Mean annual accumulation.(b) Mean daily intensity (i.e., accumulation on rainy days, where a rainy day is any day with rain above 0 mm day −1 ).(c) Modified Fournier index, the sum of the squared monthly rainfall, divided by the annual rainfall.(d) Burst factor, the sum of the product of monthly rainfall accumulation with the maximum daily rainfall for the month, divided by the annual rainfall.

Figure 2 .
Figure 2. As in Fig. 1, but for the average of 21 downscaled CMIP5 models (the BCCAv2 data set).

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

Figure 3 .
Figure 3. Overview of observational estimates of R. (a) Target erosivity field.(b) Erosivity calculated by regression from daily precipitation, using the formula R = a N N n 365 d P b d (details in the text).(c) Erosivity calculated from annual modified Fournier index F , using the published formula of Reinard and Freimund.(d) Erosivity calculated by regression from annual modified Fournier index F using the formula R = a • F + b (details in the text).(e) Same as in (d) but using the formula R = a • F b .

Figure 4 .
Figure 4. Evaluation of observational estimates of R. Method 1: erosivity calculated from daily precipitation, using the formula R = a N N n 365 d P b d (details in the text).(a) The multiplicative coefficient "a" in the above formula, calculated by regression on the grid points in each 1 • × 1 • box (b), as in (a), but for the exponent coefficient "b".(c) Absolute error between the estimate of R in Fig. 3b and the target field in Fig. 3a.(d) Relative error in percent.

Figure 5 .
Figure 5. Evaluation of observational estimates of R. Method 2: erosivity calculated from annual modified Fournier index F , using the formula R = a • F + b (details in the text).(a) The multiplicative coefficient "a" in the above formula, calculated by regression on the gridpoints in each 1 • × 1 • box (b), as in (a), but for the coefficient "b".(c) Absolute error between the estimate of R in Fig. 3d and the target field in Fig. 3a.(d) Relative error in percent.

Figure 6 .
Figure 6.Evaluation of observational estimates of R. Method 3: erosivity calculated from annual modified Fournier index F , using the formula R = a • F b (a) The multiplicative coefficient "a" in the above formula, calculated by regression on the grid-points in each 1 • × 1 • box (b), as in (a), but for the coefficient "b".(c) Absolute error between the estimate of R in Fig. 3e and the target field in Fig. 3a.(d) Relative error in percent.

Figure 7 .
Figure 7. Evaluation of climate model estimates of R. Method 1 (the method of Fig. 4).(a) Ensemble mean erosivity.(b) Intra-ensemble standard deviation, divided by the ensemble mean erosivity (%).(c) Absolute error of the ensemble mean, compared to the target R.(d) Relative error of the ensemble mean.

Figure 8 .
Figure 8. Evaluation of climate model estimates of R. Method 2 (the method of Fig. 5).(a) Ensemble mean erosivity.(b) Intra-ensemble standard deviation, divided by the ensemble mean erosivity (%).(c) Absolute error of the ensemble mean, compared to the target R.(d) Relative error of the ensemble mean.

Figure 9 .
Figure 9. Evaluation of climate model estimates of R. Method 3 (the method of Fig. 6).(a) Ensemble mean erosivity.(b) Intra-ensemble standard deviation, divided by the ensemble mean erosivity (%).(c) Absolute error of the ensemble mean, compared to the target R.(d) Relative error of the ensemble mean.

Figure 10 .
Figure 10.Projections of changes in rainfall characteristics: difference between 2079-2099 and 1980-2000, where the future period is simulated under the RCP8.5 scenario and the past is from the historical simulation.Shown are ensemble mean differences, with stippling indicating that more than two-thirds of the models agree on the sign of the change (stippling is plotted at coarser spatial resolution, for clarity).(a) Annual mean rainfall accumulation.(b) Rainfall mean daily intensity.(c) Modified Fournier index.(d) Burst factor.

Figure 11 .
Figure 11.Projections of changes in mean accumulated rainfall (top) and modified Fournier index (bottom): difference between 2079-2099and 1980-2000, where the future period is simulated under RCP8.5 scenario and the past is from the historical simulation.Shown are ensemble mean differences, with stippling indicating that more than two-thirds of the models agree on the sign of the change.In (a) and (c) the ensemble is limited to the runs included in the downscaled ensemble used in the rest of this study, but only regridded to a 2 • × 2 • resolution.In (b) and (d) the ensemble includes 38 coupled models from the CMIP archives; and for each model we have used the ensemble mean of all available realizations, for a total of 76 historical and RCP8.5 simulations.

Figure 12 .
Figure 12.Projections of changes in erosivity: difference between 2079-2099 and 1980-2000, where the future period is simulated under RCP8.5 scenario and the past is from the historical simulation.Erosivity is calculated from daily precipitation values following method 1.(a) Ensemble mean difference in erosivity units (MJ mm ha −1 h −1 year −1 ).(b) Intra-ensemble standard deviation in erosivity units.(c) Ensemble mean difference in percent of mean 20th century values.(d) Intra-ensemble standard deviation in percent.

Figure 13 .
Figure 13.As in Fig. 12, but for erosivity calculated by non-linear regression from the modified Fournier index F (method 2).

Figures
Figures 12, 13, and 14  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.12, R is estimated from the daily precipitation values (method 1), while in Figs. 13 and 14 it is estimated from F , either through a linear (method 2) or through a non-linear regression (method 3), as described in Sect. 2.

Figure 14 .
Figure14.As in Fig.12, but for erosivity calculated by non-linear regression from the modified Fournier index F (method 3).