Rain erosivity map for Germany derived from contiguous radar rain data

Erosive rainfall varies pronouncedly in time and space. Severe events are often restricted to a few square kilometers. Radar rain data with high spatiotemporal resolution enable this pattern of erosivity to be portrayed with high detail. We used radar data with a spatial resolution of 1 km² over 452 503 km² to derive a new erosivity map for Germany and to analyze the seasonal distribution of erosivity. The expected long-term regional pattern was extracted from the scattered pattern of events by several steps of smoothing. This included averaging erosivity from 2001 to 2017 and smoothing in time and space. The pattern of the resulting map was predominantly shaped by orography. It generally agrees well with the erosivity map currently used in Germany (Sauerborn map), which is based on regressions using rain gauge data (mainly from the 1960s to 1980s). In some regions the patterns of both maps deviate because the regressions of the Sauerborn map were weak. Most importantly, the new map shows that erosivity is about 66 % larger than in the Sauerborn map. This increase in erosivity was confirmed by long-term data from rain gauge stations that were used for the Sauerborn map and which are still in operation. The change was thus not caused by using a different methodology but by climate change since the 1970s. Furthermore, the seasonal distribution of erosivity shows a slight shift towards the winter period when soil cover by plants is usually poor. This shift in addition to the increase in erosivity may have caused an increase in erosion for many crops. For example, predicted soil erosion for winter wheat is now about 4 times larger than in the 1970s. These highly resolved topical erosivity data will thus have definite consequences for agricultural advisory services, landscape planning and even political decisions.

1 Introduction

Soil erosion by heavy rain is regarded as the largest threat to the soil resource. Rain erosivity is rain's ability to detach soil particles and provide transport by runoff and thereby is one of the factors influencing soil erosion. The most commonly used measure of rain erosivity is the R factor of the Universal Soil Loss Equation (USLE; Wischmeier, 1959; Wischmeier and Smith, 1958, 1978) or the Revised Universal Soil Loss Equation (RUSLE; Renard et al., 1991), although other concepts also exist (Morgan et al., 1999; Schmidt, 1991; Williams and Berndt, 1977). The R factor is given as the product of a rain event's kinetic energy and its maximum 30 min intensity. Both components are usually derived from hyetographs recorded by rain gauges. Such rain gauge data are spatially scarce. For instance, in Germany only one rain gauge per 2571 km² was available for the currently used R map (Sauerborn, 1994; this map will be called “Sauerborn map” in the following). Hence, point information has to be spatially interpolated to derive an R map that enables us to estimate R for any location. Different interpolation techniques have been applied. Most often correlations (transfer functions) of R with other meteorological data available at higher spatial density were used (for an overview see Nearing et al., 2017). The Sauerborn map was based on correlations between R and normal-period summer rain depth or normal-period annual rain depth, which differ between federal states (Rogler and Schwertmann, 1981; Sauerborn, 1994, and citations therein).

Recent research has shown that the erosivity of single events exhibits strong spatial gradients (Fiener and Auerswald, 2009; Fischer et al., 2016, 2018b; Krajewski et al., 2003; Pedersen et al., 2010; Peleg et al., 2016). This is due to the small spatial extent of convective rain cells, which is typical for erosive rains. The resulting heterogeneity has two consequences. First, interpolation of erosivity between two neighboring rain stations will not be possible for individual rains because a rain cell in between may be completely missed. Second, even long records of rain gauge data may miss the largest events that occurred in close proximity to a rain gauge and thus underestimate rain erosivity. This is illustrated nicely by the data of Fischer et al. (2018b). They showed that the largest event erosivity, which was recorded by contiguous measurements over only 2 months, was more than twice as large as the largest erosivity recorded by 115 rain gauges over 16 years and the same area. Furthermore, this single event contributed about 20 times as much erosivity as the expected long-term average. Even in a 100-year record this single event would thus still change the long-term average erosivity. The large variability of erosivity in space and time then directly translates to soil loss. This may be illustrated by soil loss measurements in vineyards in Germany. Emde (1992) found a mean soil loss of 151 t ha⁻¹ yr⁻¹ averaged over 10 plot years while Richter (1991) only measured 0.2 t ha⁻¹ yr⁻¹, averaged over 144 plot years. The difference is due to the largest event during the study by Emde (1992), which obviously had too much influence on the mean compared to the size of his data set. Such an event was missing entirely in Richter's (1991) much larger data set. The inclusion of rare events when measured by chance by a rain gauge leads to statistical problems due to their extraordinary magnitude. They cause outliers in regression analysis and thus strongly affect transfer functions. To avoid an effect by single events to the transfer function, Rogler and Schwertmann (1981) excluded all events for which the estimated return period was more than 30 years (assuming that event erosivities followed a Gumbel distribution). In consequence, the largest event was replaced by zero erosivity and, in turn, soil erosion was underestimated.

The demand for contiguous rain data to create R factor maps was only recently met by satellite data (Vrieling et al., 2010, 2014) and by radar rain data with considerably larger spatial (presently up to 9-fold) and temporal resolution (presently up to 36-fold) (Fischer et al., 2016). Put simply, the measurements are based on the principle that radar beams are reflected by hydrometeors (Bringi and Chandrasekar, 2001; Meischner et al., 1997). The intensity of the reflection depends on rain intensity and the travel time of the reflected radar beam depends on the distance between the emitting and receiving radar tower and the hydrometeors within the measurement volume. Radars usually measure with a resolution of approx. 1^∘ azimuth and 125 to 250 m in the direction of beam propagation. The data are then typically transformed to grids of square pixels of 1 km² (Bartels et al., 2004; Fairman et al., 2015), 4 km² (Koistinen and Michelson, 2002; Michelson et al., 2010) or 16 km² (Hardegree et al., 2008) after many refinement steps.

An R factor (map) can serve two purposes with contrasting requirements. First, it can be used in combination with measured soil loss or reported damages (e.g., Mutchler and Carter, 1983; Vaezi et al., 2017; Fischer et al., 2018a). In this hindcast case, the highest possible spatial and temporal resolution is recommended. The second application of the R factor is for forecasting erosion, which is required, e.g., for field use planning (Wischmeier and Smith, 1978). In this case, the long-term expectation is of interest rather than the true R factor of the (near) past that was influenced by the stochastic location of individual rain cells. Thus, for future modeling, smoothing of the stochastic noise is necessary.

Existing R maps have also undergone a number of smoothing steps, although this is not explicitly stated in the corresponding reports. Most R maps use regressions between long-term averages of erosivity and long-term meteorological parameters, e.g., annual rain depth. For long-term averages, periods of more than 20 years are accepted (Chow, 1953; Wischmeier and Smith, 1978) to remove the stochasticity of individual events and leave the general pattern. In consequence, a temporal smoothing follows from using long-term averages, and a spatial smoothing follows from the transfer functions and their application to rainfall maps. These rainfall maps include a third step of smoothing because the meteorological recommendation is to use normal-period rainfall (30-year) data, and point data (meteorological stations) have to be extended to create a map. For example for the R map in Germany (Rogler and Schwertmann, 1981), the precipitation map of 1931 to 1960 was used, which was the last available normal period although rain erosivities were derived mainly from measurements in the 1960s and 1970s. This precipitation map was mainly based on educated guesses of the best meteorologists at that time as geostatistical tools were only developed later (Matheron, 1970). With the large increase in data availability by radar measurements and the development of (geostatistical) smoothing tools, this uncontrolled smoothing can be replaced by accepted statistical methods of smoothing. The general recommendation is to apply smoothing until the pattern that is intended to be shown can be seen (O'Haver, 2018; Simonoff, 1996; Quantitative Decisions, 2004) by using several smoothing steps in sequence (O'Haver, 2018; Wedin et al., 2008).

In this study, we used the new RW product from the radar climatology RADKLIM (RADar KLIMatologie) from the German Meteorological Service (Deutscher Wetterdienst, DWD). RW data provide gauge-adjusted and further refined precipitation for a pixel size of 1 km × 1 km (Winterrath et al., 2017, 2018). RW data of 17 years (2001–2017) were available as a contiguous source of rain information. Using these data to establish a new R factor map for Germany should be a major step forward compared to the Sauerborn map, which was derived from an inconsistent set of data compiled by different researchers (e.g., some had winter precipitation data available and used them while others did not; see Sauerborn, 1994) and with equations developed independently for 16 federal states. Our data set is much larger (by a factor of 2571 regarding locations) and, because of the contiguous data source, it does not require interpolation with transfer functions. Our first hypothesis was that there will be considerable changes in the pattern of erosivity due to the removal of transfer-function weaknesses. Our second hypothesis was that the R factor map will exhibit larger values than the Sauerborn map, for two reasons. Very large and rare events will no longer be missed, as occurred previously due to the large distances between meteorological stations, and there is no longer any need to remove these events to arrive at robust transfer functions. The second reason for larger R factors is due to global climate change, as Rogler and Schwertmann (1981) and Sauerborn (1994) mostly used data from the 1960s, 1970s and 1980s. Global climate change is expected to increase rain erosivity (Burt et al., 2016).

2 Material and methods

2.1 Radar-derived precipitation data

DWD runs a Germany-wide network of presently 17 C-band Doppler radar systems (Fig. 1). This network underwent several upgrades during the analysis period. At the start of the time period considered, five single-polarization systems (DWSR-88C, AeroBase Group Inc., Manassas, USA) were operated without a Doppler filter, the latter being added between 2001 and 2004. Between 2009 and 2017, DWD replaced the network of C-band single-polarization systems of the types METEOR 360 AC (Gematronik, Neuss, Germany) and DWSR-2501 (Enterprise Electronics Corporation, Enterprise, USA) with modern dual-polarization C-band systems of the type DWSR-5001C/SDP-CE (Enterprise Electronics Corporation), all equipped with Doppler filters. During this period, a portable interim radar system of the type DWSR-5001C was installed at some sites.

Figure 1 Coverage (blue circles) of the 17 German weather radars for a range (utilized radius) of 128 km and the tower locations in 2017 (locations of some radar towers have changed over time). Black lines denote federal states; the federal states of Bavaria (cross-hatched), Lower Saxony (hatched) and selected mountain ranges (light brown) are mentioned in the text. Axis ticks represent distances of 100 km. A detailed topographic map can be found in Fig. S1.

The radar systems permanently scan the atmosphere to detect precipitation signals. Every 5 min, the radars perform a precipitation scan, each with terrain-following elevation angle to measure precipitation near the ground. The resulting local reflectivity information over a range of currently 150 km in real time and a constant 128 km in the climate approach is combined to form a Germany-wide mosaic of about 1100 km in the north–south direction and 900 km in the west–east direction. The reflectivity information is converted to precipitation rates applying a reflectivity–rain rate (ZR) relationship (Bartels et al., 2004). An operational quality control system screens the radar data. To further improve the quantitative precipitation estimates, the radar-derived precipitation rates are summed to hourly totals and immediately adjusted to gauge data from more than 1000 meteorological stations resulting in RADOLAN (RADar OnLine ANeichung, i.e., online-adjusted, radar-derived precipitation), which provides precipitation data in real time, mainly for applications in flood forecasting and flood protection (Bartels et al., 2004; Winterrath et al., 2012).

Based on RADOLAN, the climate version RADKLIM is derived. Compared to the real-time approach, the data are additionally offline-adjusted to daily gauge data, combining a total of more than 4400 rain gauges measuring hourly and daily (equivalent to 1 rain gauge per 80 km²). The data are then reprocessed by new climatological correction methods, e.g., for spokes, clutter or short data gaps. Spokes result from permanent obstacles blocking the radar beam, while clutter is introduced by non-meteorological targets like windmills or birds. The final product (called RW data) has a temporal resolution of 1 h and a spatial resolution of 1 km × 1 km in polar stereographic projection. For more detailed information on RADKLIM the reader is referred to Winterrath et al. (2017). The RW data, restricted to the German territory, are freely available (Winterrath et al., 2018). For the first time, the RADKLIM data set provides contiguous precipitation data with high temporal and spatial resolution. It includes local heavy or violent precipitation events (for classification of heavy and violent see UK Met Office, 2007) that are partly missed by point measurements alone. Thus, it particularly improves the analysis of extreme precipitation events.

Two additional data sets were used to verify the validity of the approach and to examine effects of methodological details (see below). These data sets are erosivities derived from radar data at 5 min resolution taken from Fischer et al. (2016) and erosivities derived from rain gauge data of 115 stations in Germany from 2001 to 2016, which were taken from Fischer et al. (2018b).

2.2 Erosivity calculation procedure

According to Wischmeier (Wischmeier, 1959; Wischmeier and Smith, 1958, 1978), the erosivity of a single rain event (R_e in N h⁻¹) is the product of the maximum 30 min rain intensity (I_max30 in mm h⁻¹) and the total kinetic energy per unit area (E_kin in kJ m⁻²).

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{R}_{\mathrm{e}}=I_{\max \mathrm{30}}\cdot E_{\mathrm{kin}}\end{array}$$

An erosive rain event is defined to have at least a total precipitation amount (P in mm) of 12.7 mm or an I_max30 of more than 12.7 mm h⁻¹ that is separated from the next rain by at least 6 h. In order to scan and fulfil the 6 h criterion, we did not separate between days but used a continuous 17-year data stream. Specific kinetic energy e_kin,i per millimeter rain depth (in kJ m⁻² mm⁻¹) is given for intervals i of constant rain intensity I (in SI units according to Rogler and Schwertmann, 1981).

For 0.05 mm h${}^{-\mathrm{1}}\le I<\mathrm{76.2}$ mm h⁻¹,

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }{e}_{\mathrm{kin},i}=\left(\mathrm{11.89}+\mathrm{8.73}\cdot {\log }_{\mathrm{10}}I\right)\times {\mathrm{10}}^{-\mathrm{3}}.\end{array}$$

For I<0.05 mm h⁻¹,

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{e}_{\mathrm{kin},i}=\mathrm{0}.\end{array}$$

For I≥76.2 mm h⁻¹,

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{e}_{\mathrm{kin},i}=\mathrm{28.33}\times {\mathrm{10}}^{-\mathrm{3}}.\end{array}$$

We used the equation by Wischmeier and Smith (1978) to calculate specific kinetic energy although several others have also been proposed (van Dijk et al., 2002) with none being superior (Wilken et al., 2018). Our choice retained comparability with the Sauerborn map. Furthermore van Dijk et al. (2002) had shown that kinetic energy as obtained by the Wischmeier and Smith equation did not deviate from measured kinetic energy in Belgium neighboring Germany.

For all intervals i, e_kin,i is multiplied with the rain amount of this interval and then summed to yield E_kin for the entire event. The annual erosivity of a specific year is the sum of R_e of all erosive events within this year. The average annual erosivity (R) is then the average of all annual erosivities during the study period (17 years in this case). While in the USA and other countries the unit MJ mm ha⁻¹ h⁻¹ is often used for R_e, we use N h⁻¹ because it is the unit most often used in Europe and because of its simplicity. Both units can be easily converted by multiplying the values in N h⁻¹ with a factor of 10 to yield MJ mm ha⁻¹ h⁻¹. The unit for R is then N h⁻¹ yr⁻¹.

Rain erosivity strongly depends on intensity peaks. Fischer et al. (2018b) have shown that these peaks increasingly disappear the lower the spatial and temporal resolution becomes. This can be accounted for by scaling factors but these scaling factors can only adjust to an average behavior, while the factors may either be too large or too small for a specific event. A high spatiotemporal resolution should be used to determine R_e for individual events. This is not required to determine the long-term average pattern like an R factor map for planning and prediction purposes. In that case, data with lower resolution and the application of appropriate scaling factors are advantageous because this will reduce the noise introduced by large events of small spatial extent that would not be leveled out by averaging alone. We will use data in 1 h time increments as those are adjusted to rain gauge measurements and the number of data is reduced by a factor of 12 compared to 5 min increments. This is especially important when all calculations, including identification of rain breaks >6 h and periods of I_max30, have to be carried out for many years and many locations. In our case, roughly 7×10¹⁰ 1 h increments had to be processed.

According to Fischer et al. (2018b), the following modifications in the calculation of R_e had to be made to account for the temporal resolution of 1 h, the spatial resolution of 1 km² and the method of measuring rain by radar: (i) I_max30 was replaced by the maximum 1 h rain intensity and the threshold for I_max30 was lowered to 5.8 mm h⁻¹, while the total precipitation threshold remained at 12.7 mm. (ii) Five or more subsequent 1 h intervals without rain separated events, which assumed that rain events begin on average in the middle of the first nonzero rain interval and end again in the middle of the last nonzero rain interval, yielding a total rain break of at least 6 h. (iii) The temporal scaling factor was 1.9 and the spatial scaling factor was 1.13, to which 0.35 had to be added to account for the radar measurement instead of the rain gauge measurement. The total scaling factor [$(\mathrm{1.13}+\mathrm{0.35})\times \mathrm{1.9}$] was then 2.81.

Gaps in the time series were considered when calculating annual sums of erosivity by scaling the total sum of erosivity over the whole time series to 365.25 days. If the effective number of missing values exceeded 2 months per year, the respective year was excluded from the calculation for that pixel. If the effective number of excluded years was larger than one, the respective pixel was excluded. This was the case for 0.6 % of all pixels.

2.3 Steps to generate an R factor map

The reduction of noise by using 1 h increments and a 17-year mean was still not sufficient to level out the most extreme events. Two further smoothing steps were therefore applied. The first step was to winsorize the annual erosivities of the 17 years for each individual pixel by replacing the lowest value with the second-lowest value and the highest value with the second-highest value (Dixon and Yuen, 1974). Winsorizing is an appropriate measure for calculating a robust estimator of the mean in symmetrically distributed data, but it is biased for long-tailed variables like rain erosivity. Thus, the country-wide mean of all winsorized data (94 N h⁻¹ yr⁻¹) was smaller than the mean of the original data (96 N h⁻¹ yr⁻¹). In order to remove this bias, we binned all data in 26 bins of 20 N h⁻¹ yr⁻¹ width and calculated the mean R before and after winsorizing. For bins with R<180 N h⁻¹ yr⁻¹, comprising 95 % of all pixels, the bias increased linearly with R (r²=0.92; n=8) and amounted to 2.3 % of R. Above 180 N h⁻¹ yr⁻¹ there was no further increase in the bias (r²=0.01, n=18), which was, on average, 3.4 N h⁻¹ yr⁻¹. We removed the bias by adding 2.3 % to all values <180 N h⁻¹ yr⁻¹ and 3.4 N h⁻¹ yr⁻¹ to all values above.

The last smoothing step applied geostatistical methods. A semivariogram (over a range of 50 km) was calculated and ordinary kriging was applied. Geostatistical analysis was done using the program R (version 3.5.0; R Core Team, 2018) and gstat (Gräler et al., 2016). A block size of 10 km × 10 km was chosen to remove noise and to fill the pixels with data gaps, while the spatial resolution remained at 1 km. The missing information was obtained from neighboring pixels. The radar data outreached the German borders. In total, 452 503 pixels were used to ensure small kriging variances near borders or on islands, while the final map was restricted to the German land surface (357 779 pixels).

Using 1 h data instead of 5 min data reduced the effect of single extreme events at certain locations. Winsorizing reduced the effect of extreme years at a location, in addition to the effect of averaging 17 years. Finally, kriging used the information from neighboring pixels to reduce the effect of the extremes. This smooths among near neighbors (distance < 20 km) but does not affect the regional pattern (>20 km). To evaluate whether this was the case and to quantify the effect of all smoothing steps, we used the data from Fischer et al. (2016). They had calculated rain erosivity from 5 min resolution radar data for 2 years (2011 and 2012) and an area of 14 358 km² (yielding a total of 28 770 pixel years), which is called “test region” in the following. Using these data we calculated semivariograms from annual to biennial erosivities based on 5 min and 1 h resolution. These semivariograms were compared to those from 17-year average erosivities, 17-year winsorized average erosivities, and 17-year winsorized and kriged erosivities for the test region and for the entire area of Germany. Smoothing should reduce the influence of individual violent thunderstorm cells and reveal the regional pattern. In geostatistical analysis this decreases the sill of the semivariogram while the range increases as it changes from being dominated by thunderstorm cells to being dominated by the regional pattern. The regional trend was calculated as the difference between the square root of semivariances at distances of 40 and 20 km divided by the difference in distance of 20 km to examine whether it was influenced by the individual smoothing steps. The effect of violent rain cells was calculated as the square root of the semivariance at a distance of 20 km divided by the difference in distance of 20 km minus the regional trend.

2.4 Annual erosivity return periods

Rain erosivity usually follows long-tailed distributions. This leads to the question of how frequent years of extraordinarily large erosivity are. To answer this question, the development of cumulative distribution curves (for basic concepts see Stedinger et al., 1993) is required. A period of 17 years is not sufficient to reliably estimate a cumulative distribution curve for every pixel. Therefore, we combined all annual erosivities of the total data set (452 503 pixels and 17 years) after expressing each of them relative to the corresponding winsorized and bias-corrected mean of the pixel (in %). This enabled the cumulative distribution curves to be calculated from a large data set (n=7.7 million). The expected maximum relative annual erosivity for a given return period could then be estimated from the complementary cumulative distribution curve (exceedance). This was also done for the relative annual erosivities of the test region, calculated from 1 h rain data, to examine whether the general cumulative distribution curve also applies to smaller regions.

The erosivities, when calculated from 1 h rain data, are already smoothed and do not adequately reflect the extremes that result from data that are better resolved, such as the 5 min rain data. The cumulative distribution curve for the test region was also calculated using the 5 min rain data. Given that the cumulative distribution curves of the entire study area and the test region agree for the relative erosivities calculated from 1 h data, we expect that the relative erosivities calculated from 5 min rain data of the test region can serve as a first estimate for the entire study region. The cumulative distribution curve for the test region calculated from 5 min data will then be a fair estimate of the return periods anywhere in the entire research area.

2.5 Seasonal distribution of erosivity

The seasonal distribution of erosivity, calculated as the relative contribution of each day to total annual erosivity, is called the erosion index distribution or EI distribution (Wischmeier and Smith, 1978). It is required in erosion modeling to determine the influence of seasonally varying soil cover due to crop development. The convolution of the seasonal effect of soil cover with the seasonal EI distribution results in the so-called crop and cover factor (C factor) of the USLE. The EI distribution was calculated for each pixel and averaged over all 452 503 pixels. Seventeen years of data still did not suffice to show similar amounts of erosivity on subsequent days, despite the large number of pixels. There was still considerable scatter that required smoothing to illustrate the seasonal distribution. Smoothing between individual days during the year involved three steps (for details of the methods see Tukey, 1977): first a 13-day centered median was calculated for each day. A centered median smooths but preserves the common trend signal (Gallagher and Wise, 1981), which is also true for the two subsequent steps. A 3-day skip mean (leaving out the second day) was calculated from the results, followed by a 25-day centered Hanning mean (weighted mean with linearly decreasing weights). To account for the periodic nature of the EI distribution and to allow the smoothing methods to be applied at the start and the end of the year, the year was replicated and shifted by plus or minus 1 year.

Radar measurements tend to have larger errors during wintertime with snowfall. The reduced reflectivity of snow particles may lead to an underestimation of the precipitation rate, while the increased reflectivity of melting particles in the bright band may cause an overestimation. Moreover, the lower boundary layer promotes a potential overshooting of the radar beam with regard to the precipitating cloud (Holleman et al., 2008; Wagner et al., 2012). Such measurement problems, if relevant, should especially influence the EI distribution during winter months and cause a deviation from measurements at meteorological stations. Therefore, we also calculated the EI distribution using data from 115 rain gauges distributed throughout Germany and covering 2001 to 2016. These data were taken from Fischer et al. (2018b). This data set will also be used in the discussion for comparison of recent radar-derived erosivities with recent rain-gauge-derived erosivities and with historic rain-gauge-derived erosivities taken from the literature.

Figure 2 Annual average R factor (N h⁻¹ yr⁻¹) map of Germany from 17 years of radar rain data. The axes' ticks represent distances of 100 km. Color classes from yellow to dark blue comprise approximately 10 %, 20 %, 20 %, 25 %, 15 %, 4 %, 3 % and 3 % of the area, respectively. For a comparison with the Sauerborn (1994) map see Fig. S2 in the Supplement. For comparison with the map before winsorizing and before kriging see Figs. S3 and S4. Average R factors for the 401 local authority areas (average area 893 km²) are given in Table S1.

3 Results

3.1 The effects of smoothing

The effects of smoothing on the appearance of the maps were negligible (compare Fig. 2 with Figs. S3 and S4) because smoothing had only removed the extraordinarily large variability that exists on small temporal and spatial scales. However, the high data density revealed that even long-term averages were insufficient to remove all influence of erratic cells of violent rain, and further attenuating steps had to follow. Annual sums of rain erosivity from 5 min data for the test region varied most due to the dominance of individual cells of violent rain that did not overlap or fill the entire area (semivariogam I in Fig. 3a). This was indicated by the short range (20 km) and high semivariance for that range (2749 N² h⁻² yr⁻²) (Table 1). The standard deviation (SD) of two pixels separated by 20 km thus was 52 N h⁻¹ yr⁻¹ (square root of 2749 N² h⁻² yr⁻²), which is more than half of the average annual erosivity in Germany. After averaging both years (2011 and 2012), the semivariance for a distance of 20 km was only reduced to 1569 N² h⁻² yr⁻² and the range stayed the same at approximately 20 km (semivariogam III in Fig. 3a). Both findings indicated that, even after averaging 2 years, the individual cells of violent rain were still fully detectable and had not merged to form a larger pattern. In consequence, the regional trend, albeit detectable, appeared minor (Table 1).

Table 1 Influence of temporal resolution of rain data (5 min and 1 h), averaging (1, 2, and 17 years), winsorizing and kriging on the semivariance (γ) at three distances h. For complete semivariograms see Fig. 3a.

¹ The regional trend was calculated as the difference between the square roots of γ at distances of 40 and 20 km divided by the difference in distance of 20 km. ² The effect of violent rain cells was calculated as the square root of γ at a distance of 20 km divided by the difference in distance of 20 km minus the regional trend.

Figure 3 (a) Experimental semivariograms of annual erosivity of the test region for different temporal resolutions of rain data (5 min and 1 h), different averaging (1 year, 2011+2012; 2 years, mean of 2011 and 2012; 17 years, 2001 to 2017), winsorizing and kriging (for selected distance classes see Table 1). The line through semivariogram VI is a linear regression through the origin (r²=0.9889). (b) Comparison of semivariances for the 1 h, 17-year and winsorized data before kriging for the test region and for the whole of Germany.

The effect when using data with a resolution of 1 h was almost as strong as when 2 years were averaged. Semivariance at a distance of 20 km was only 1667 N² h⁻² yr⁻² for annual values (semivariogam II in Fig. 3a) and 953 N² h⁻² yr⁻² for biennial averages (semivariogam IV in Fig. 3a). Even more important, the regional trend became more visible due to smoothing of the extreme events by using 1 h instead of 5 min data. This regional trend is evident from the gradual increase in semivariance over the entire distance of 50 km shown in Fig. 3. Importantly, smoothing by using 1 h data did not change average erosivity because the difference was adequately compensated for by the temporal scaling factor. The biennial average for the test region was 115 N h⁻¹ yr⁻¹ when calculated from 5 min data and 114 N h⁻¹ yr⁻¹ when calculated from 1 h data.

Averaging annual erosivities of the test region over 17 years further reduced variability (semivariogram V in Fig. 3a). Semivariance strongly decreased to 197 N² h⁻² yr⁻² and the influence of individual cells of violent rain became small relative to the regional trend. This led to an almost linear increase in semivariance over distance. The influence of extreme years in individual pixels was further reduced by winsorizing, which slightly reduced semivariance at 20 km distance to 190 N² h⁻² yr⁻² (semivariogram VIII in Fig 3b). For all of Germany, winsorizing reduced the standard deviation of a pixel over time from, on average, 49 to 39 N h⁻¹ yr⁻¹, while bias correction left the mean of erosivity over all pixels unchanged at 96 N h⁻¹ yr⁻¹. The effect on the appearance of the map was small (compare Figs. S3 and S4) because only small erratic patches of extraordinarily high or low erosivity disappeared.

Finally, kriging reduced semivariance at 20 km distance to 121 N² h⁻² yr⁻², leaving mainly the regional trend (semivariogram VI in Fig. 3a). Thus, the step from 5 min to 1 h resolution reduced semivariance at 20 km distance by a factor of 1.6 while averaging 17 years reduced semivariance by a factor of 8.5. Winsorizing contributed a factor of 1.04 and kriging a factor of 1.6. In total, semivariance was reduced by a factor of 23, indicating a pronounced patchiness of erosive rains on the annual scale that could not be leveled out by averaging 17 years alone. The effect of each smoothing step decreased with increasing distance. For a distance of 10 km, the combined factor was 32 while it was only 13 for a distance of 30 km. This was due to the decreasing importance of thunderstorm cells relative to the regional trend. Independent of the degree of smoothing, the regional trend, extracted from the change in semivariance between distances of 20 and 40 km, remained practically unchanged at 0.2 N h⁻¹ yr⁻¹ km⁻¹ (Table 1). In contrast, the effect of violent rain cells decreased greatly using the smoothing steps from 2.4 to 0.3 N h⁻¹ yr⁻¹ km⁻¹. The effect on the appearance of the map was again small (compare Fig. S4 and Fig. 2) because only large contrasts between close neighbors disappeared, which are hardly visible due to the small pixel size. The main visible effect was the filling of the few gaps.

After winsorizing and kriging, the semivariances for the test region followed a linear regression through the origin almost perfectly (r²=0.9889, n=50; line through semivariogram VI in Fig. 3a). This indicated that the variation in erosivity over a distance of 50 km followed linear trends without any noise (nugget) or short-range structures that could be attributed to individual cells of violent rain. The semivariances, when calculated for the whole of Germany, were considerably larger (twice as large at a distance of 50 km; Fig. 3b, semivariogam VII) and close to a linear trend only for short distances (e.g., a linear regression through the origin for the first 15 km yielded r²=0.9905). For longer distances, the semivariogram followed an exponential model (nugget 4 N² h⁻² yr⁻², partial sill 970 N² h⁻² yr⁻², effective range 123 km). The larger semivariance and the exponential model were both caused by the inclusion of mountain areas with large erosivities and steep erosivity gradients that were missing in the test region.

3.2 R factor map

Erosivity was on average 96 N h⁻¹ yr⁻¹ but varied between 46 and 454 N h⁻¹ yr⁻¹. The regional pattern of erosivity (Fig. 2) was mainly determined by orography (for a detailed topographic map see Fig. S1 in the Supplement). The largest values (above 185 N h⁻¹ yr⁻¹) were found in the very south where the northern chain of the Alps reaches altitudes of almost 3000 m a.s.l. (above sea level). Lower mountain ranges are also characterized by larger mean annual erosivities than in their surrounding area (compare Fig. 1 or Fig. S1 with Fig. 2). For instance, the Bavarian Forest with elevations of up to 1450 m a.s.l. exhibited annual erosivities of above 155 N h⁻¹ yr⁻¹. The Ore Mountains with elevations of up to 1244 m a.s.l., had erosivities mostly between 125 and 155 N h⁻¹ yr⁻¹. Also mountain ranges like the Black Forest or the Harz mountains clearly shape the erosivity map. Additionally, upwind–downwind effects were detectable. For example, the areas west-northwest (upwind) of the Harz mountains had erosivities of between 70 and 80 N h⁻¹ yr⁻¹, while the areas east-southeast (downwind) received less than 65 N h⁻¹ yr⁻¹.

Figure 4 Cumulative distribution curve of the annual R factor relative to the long-term mean R factor of a pixel. Dashed black line applies for erosivities derived from 1 h data for the whole of Germany and 17 years (n=7.7 million). Solid green line applies for erosivities derived from 5 min data for the test region and 2 years (n=24 770). Straight vertical and horizontal lines indicate return periods between 2 and 100 years. The y axis is probability scaled; the x axis is log scaled.

3.3 Annual erosivity return periods

The cumulative distribution of the relative annual erosivities followed a straight line in a probability plot fairly well when the logarithm was used (Fig. 4). This indicated a log-normal distribution (log mean 1.96; log SD 0.19). A very similar cumulative distribution was found for annual erosivities derived from the 1 h data of the test region (log mean 1.97; log SD 0.18). The distribution based on the less-smoothed 5 min data was considerably wider (log mean 1.94; log SD 0.22). The annual expected erosivity was 88 %, 216 % and 273 % of the respective long-term mean for return periods of 2, 30 and 100 years when the 5 min data were used (Fig. 4). It is important to note that these values apply for averages of 1 km² pixels and include the smoothing that results from the radar measurement, the radar reprocessing and from using 5 min rain increments. Even more extreme years are expected to occur in reality.

3.4 Seasonal distribution of erosivity

There was a pronounced peak in the seasonal distribution of relative erosivity during summer months (Fig. 5). The daily erosion index increased rapidly from mid-April to mid-May and was 0.61 % day⁻¹ on average in June, July and August. From mid-August to September the daily erosion index declined rapidly. In winter months the daily erosion index was small (mean of December, January, February and March: 0.08 % day⁻¹). There was no detectable difference in the seasonal variation between different regions in Germany (see Fig. S5). The cumulative distribution functions of different regions correlated with at least r²=0.998 (n=365).

Figure 5 Measured (circles) and smoothed (solid blue line) daily erosion index derived from radar data. The daily erosion index calculated from measurements between 2001 and 2016 at 115 rain gauges distributed throughout Germany is given for comparison (dashed orange line). For C factor calculations the smoothed values can be taken from Table S2. Comparison of measured daily erosion indices separated for different regions in Germany and the respective cumulative distribution curves are depicted in Fig. S5.

Even more striking was the fact that this pattern required considerable smoothing to yield a continuous seasonal time course. The difference between subsequent days in the unsmoothed data was enormous (e.g., 1.5 % day⁻¹, 0.4 % day⁻¹ and 0.4 % day⁻¹ on 29, 30 and 31 July). This was despite the large number of measurements (17 years and 455 309 pixels) that were averaged for each day. It highlights the exceptional strength of some violent rains. Despite the rather small extent of individual erosivity cells, many of them occurred in the same day, making a large relative contribution to total erosivity for this day. While particular days of the year were influenced by heavy precipitation, during other days no erosive rainfall occurred anywhere within the research area. A period of 17 years was not sufficient to level out the contrast between subsequent days. The results of the smoothing procedure show that even 221 years (17 years multiplied by a moving-average window of 13 days) were not sufficient to level out these differences. Two additional smoothing steps had to be applied to arrive at a smooth time course. Despite the strong smoothing that was necessary for the probability density function, the smoothing did not change the cumulative distribution function (which is used for calculating C factors). The cumulative distribution functions of the original data and of the smoothed data correlated with r²=0.9998 (n=365; both functions are shown in Fig. S5).

The distribution of the daily erosion index calculated from rain gauge data (1840 station years) was very similar to the distribution calculated from the much larger radar data set (compare solid and dashed lines in Fig. 5). This was especially true during winter months, when values derived from both measurement methods were considerably larger than expected from previous analysis in the 1980s.

Figure 6 Comparison of past mean erosivities derived from rain gauge data of the 1960s to 1980s as reported by Sauerborn (1994) with recent mean erosivities of the 2000s to 2010s. Recent erosivities were either determined from rain gauge data at the same meteorological stations (mean of 2001 to 2016; taken from Fischer et al., 2019; n=33; filled circles) or from radar data (mean of 2001 to 2017 and all radar pixels at a distance of <1.5 km from the meteorological stations; n=101, open circles). Both axes are square root scaled to improve resolution at low erosivities. Dashed line denotes 1:1. Solid lines are regressions through the origin.

4 Discussion

4.1 Increase in erosivity

The most striking difference between the Sauerborn map based on data from the 1960s to 1980s and the radar-derived map is a pronounced increase in erosivity. A German average of 58 N h⁻¹ yr⁻¹ was derived from the Sauerborn map (Auerswald et al., 2009), while the radar-derived map suggests an average of 96 N h⁻¹ yr⁻¹. This increase will come along with an equal increase in predicted soil losses by 69 %. An almost identical increase resulted when the erosivity of meteorological stations, as reported by Sauerborn (1994), was compared with the erosivity derived from radar data at the same locations. This resulted in an increase of 63 % (open symbols in Fig. 6). Thus, the increase in erosivity is not an effect of the regression approach that was previously used or due to better capturing of extreme events by the contiguous radar data.

Fischer et al. (2018b) calculated erosivity for 33 of the Sauerborn stations from recent (2001 to 2016) rain gauge data. A comparison of these data with the Sauerborn data (1994) also showed a similar increase of 52 % (closed symbols in Fig. 6). The increase in erosivity between the Sauerborn map and the new radar-derived map is thus also not an artifact of using radar data but the result of a true change in erosivity over time. This is further corroborated by Fiener et al. (2013), who analyzed long-term records from 10 meteorological stations in western Germany. They found an increase in erosivity of 63 % between 1973 and 2007. Both independent findings leave little doubt that the pronouncedly higher values in the new erosivity map are a result of a change in weather properties and not a result of the difference in the applied methodologies, although we did expect the mean to increase due to the contiguous data set, which is better at recording rare extremes.

Figure 7 Average R factor relative to the 17-year mean radar-derived R factor depending on the mean year of data origin. Data below year 1990 are calculated from statewide averages determined from meteorological station records; year is the mean year of station records. Data above year 2000 are radar-derived R factors for all of Germany for individual years. The closed circle denotes the reference point (present map).

A time series of 17 years is regarded to be too short in meteorology for calculating temporal trends. The data in Sauerborn (1994) were derived from different periods for different states. If we calculate the statewide mean R factors from her transfer functions relative to the statewide mean R factors of the radar-derived map and plot this relative R factor against the mean year from which the state-specific data originated, a 23-year-long period can be covered by the means (Fig. 7; years < 1990; the total time period of individual years covers an even wider range, mostly about ±5 years around the mean year). During this period there was a slight but insignificant increase in erosivity with time. This increase smoothly leads over to the steeper increase in radar-derived Germany-wide annual R factors if we express them again relative to the 17-year mean (Fig. 7; years > 2000). Both data sets combined cover more than 60 years and yield a very highly significant regression (r²=0.7340, n=27) that indicates an accelerating increase in erosivity likely due to climate change. Furthermore, Fig. 7 indicates that at the end of the radar time series (2017) the R factor likely is already 20 % higher than the values depicted in Fig. 2.

4.2 Change in the regional pattern of erosivity

The regional patterns of the Sauerborn map and of the radar-derived map generally agree well but with two exceptions. First, the radar-derived map shows distinctly larger values southeast of the German Bight of the North Sea where the air masses coming from the North Sea are channeled by the Elbe river estuary and its Pleistocene meltwater valley and then hit the higher areas of the north German moraines. A large frequency of large rains is not unlikely in this situation. The reason that this was missed by Sauerborn (1994) using the data obtained by Hirche (1990) for Lower Saxony might be mainly due to the small data density and the regression with long-term rainfall. Only 18 stations were available for the whole of Lower Saxony and only five of them were in the area of large erosivity. Using the 18 stations in the state of Lower Saxony only, and ignoring the difference between landscapes, resulted in a rather poor regression with long-term annual rainfall (r² was only 0.32 for n=18), and therefore a large prediction error and considerable smoothing of the true erosivity pattern can be expected. For comparison, in Bavaria the regression with long-term rainfall yielded r² of 0.92 (for n=18; Rogler and Schwertmann, 1981).

The second difference in the pattern is that the radar-derived map reveals more detail than the regression-based map by Sauerborn (1994). This is especially evident in southern Germany where southwest–northeast-oriented structures seem to follow tracks of thunderstorm movement. In the northeast quarter of Germany, where the pattern is not shaped by mountain ranges, a rather patchy pattern resulted. Although Sauerborn (1994) had already found a patchy pattern in this area it appears to be patchier now. At present, it is difficult to decide whether this pattern is random due to large multicell clusters of rainstorms that will level out in the long term or whether landscape properties, e.g., the existence of large forests, cause a stable pattern in an area where other factors affecting the pattern are missing. More detailed variation may also be expected in mountainous areas but radar measurements cannot adequately show this variation. In the future, using data obtained by commercial microwave links as an additional source for retrieving precipitation (Chwala et al., 2012, 2016; Overeem et al., 2013) may improve high-resolution estimates, particularly in these areas.

4.3 Change in the seasonal distribution of erosivity

The third pronounced difference between past and recent erosivities was found for the erosion index distribution. This distribution is needed for C factor calculations (Wischmeier and Smith, 1978). A change in the seasonality of erosivity was already suggested by Fiener et al. (2013) analyzing an 80-year time series. However, Fiener et al. (2013) used data from April to October only, and their results therefore cannot be compared directly with our results that show the most pronounced changes for the period from December to March.

At present, the C factors for all of Germany (DIN, 2017) are based on the erosion index distribution developed for Bavaria by Rogler and Schwertmann (1981), although unpublished erosion indices are also available for other federal states (e.g., Hirche, 1990). The index distribution by Rogler and Schwertmann (1981) is characterized by very low values during winter months, which in turn causes a sharp increase during summer months. In contrast, the radar-based index, although still having a pronounced summer maximum, predicts a higher percentage of erosivity during winter. Rogler and Schwertmann (1981) found that only 1.5 % of the annual erosivity fell from January to March, while Fig. 5 indicates that these months contributed 6.9 % to annual erosivity. This deviation may be caused by a regional variation in the erosion index because the unpublished indices for other federal states also suggested a larger contribution by winter months (e.g., January to March contributed 7.5 % in Lower Saxony according to Hirche, 1990). However, restricting our data set to Bavaria led to a very similar index during winter months (e.g., 6.2 % for January to March) to the index for the whole of Germany, and the discrepancy with Rogler and Schwertmann (1981) remained. Furthermore we could not find significant differences when calculating the index distribution separately for different regions (Fig. S5).

A second explanation might be that the Rogler and Schwertmann (1981) data were too limited to capture enough erosive rains during periods of infrequent erosive events. This explanation is corroborated by the large scatter between individual days that still existed in our data set (Fig. 5), although our data set was more than 50 000 times larger than the data set used by Rogler and Schwertmann (1981).

A third explanation could again be climate change. In Germany the number of extreme wet months increased in winter by 463 % from the first to the second half of the last century, while summer and autumn remained unchanged (Schönwiese et al., 2003).

The change in erosion index distribution may be regarded as being rather unimportant at first glance because erosivity is still dominated by precipitation in summer. This small increase in erosivity during the winter months, however, could have important consequences for the C factor of crops that provide only small soil coverage during winter. As there is practically no growth during winter, these crops stay susceptible to erosion over a long period. Thus they experience a considerable amount of erosivity, even though erosivity per day is small. For example, the C factor for continuous winter wheat increases from 0.04 to 0.10 when using the soil loss ratios taken from Auerswald et al. (1986) that entered DIN (2017) and the new erosion indices instead of those from Rogler and Schwertmann (1981).

4.4 Stochasticity

Soil erosion is characterized by a large temporal variability at a small spatial scale due to the stochastic character of erosive rains. About 20 years are necessary, according to Wischmeier and Smith (1978), until this variability levels out and average soil loss approaches values predicted with the (R)USLE. Our data set covered 17 years but significant additional smoothing was still necessary. One of the smoothing steps was to use hourly data, although 5 min data would have been available. In one or two decades the data series may be long enough to remove some of the smoothing steps. In particular, it would be desirable to use data of 30 min or even 5 min resolution.

This pronounced stochasticity is due to the small size of convective rain cells. Just recently it has been shown by analyzing the radar-derived rain pattern of the largest rainfall events that on average the rain amount is halved within a distance of only 2 km around the central point of a rain cell (Lochbihler et al., 2017). Given that rain amount is squared in the calculation of rain erosivity, the R factor decreases to one fourth within this distance. Larger areas are only covered if there is movement of the rain cells. This small size of rain cells questions the use of sparsely distributed rain gauges to derive rain erosivity. The inconsistent transfer functions among German states to derive erosivity from rainfall maps likely originated in the high stochasticity of rain gauge measurements under such conditions. It was only the unintended but unavoidable smoothing that was inherent in previous approaches that allowed deriving such maps. Radar technology enables us to replace this unintended smoothing using clearly defined statistical protocols and to quantify the effect of smoothing.

Another implication of this large variability is that 20 years will still not be sufficient to level out extraordinary events. The largest event erosivity that Fischer et al. (2016) found in 2 years on ∼15 000 km² was 622 N h⁻¹. Even for a 20-year period, this event will add 31 N h⁻¹ yr⁻¹ to the average annual erosivity at the small location of only a few squared kilometers (km²) where it occurred. Most studies measuring soil erosion under natural rain use much shorter intervals that usually cover only a few years and rarely exceed 10 years (see Auerswald et al., 2009, for a meta-analysis of German studies and Cerdan et al., 2010, for European studies). The interpretation of such short-term studies and the applicability of the results are limited due to the pronounced variability of natural rains.

In addition, the erosion index distribution required considerable smoothing to improve representation of the seasonal variation. Without smoothing, the shift in a certain crop stage by only 1 day can cause large discrepancies in the resulting C factor, depending on whether a day of large erosivity in the past is included or excluded at the bounds of the crop stage period. Smoothing can prevent this. This is especially important for short crop stage periods, while the effect becomes small for longer periods. For instance, the monthly sums of the smoothed data correlated closely with the sums of the unsmoothed data (coefficient of determination: 0.995; Nash–Sutcliffe efficiency: 0.994).

5 Conclusions

Radar-derived rainfall data enable us to derive highly resolved and contiguous maps of erosivity with high spatial detail. This avoids errors in landscapes with insufficient rain gauge density. The analysis showed that present (2001 to 2017) rain erosivity is considerably higher than erosivity in the past (1960s to 1980s). Furthermore, the seasonal distribution of rain erosivity also deviates from that of the past period. Winter months contribute more to total erosivity than previously recorded. Considerably more erosion can be expected for crops that are at a highly susceptible stage of development in winter. In consequence, the predicted soil loss will change pronouncedly by using recent erosivity and the ranking of crops regarding their erosion potential will change. This will have definite consequences for agricultural extension and advisory services, landscape planning and even political decisions.