Articles | Volume 24, issue 11
Research article
17 Nov 2020
Research article |  | 17 Nov 2020

Estimation of rainfall erosivity based on WRF-derived raindrop size distributions

Qiang Dai, Jingxuan Zhu, Shuliang Zhang, Shaonan Zhu, Dawei Han, and Guonian Lv

Soil erosion can cause various ecological problems, such as land degradation, soil fertility loss, and river siltation. Rainfall is the primary water-driven force for soil erosion, and its potential effect on soil erosion is reflected by rainfall erosivity that relates to the raindrop kinetic energy. As it is difficult to observe large-scale dynamic characteristics of raindrops, all the current rainfall erosivity models use the function based on rainfall amount to represent the raindrops' kinetic energy. With the development of global atmospheric re-analysis data, numerical weather prediction techniques become a promising way to estimate rainfall kinetic energy directly at regional and global scales with high spatial and temporal resolutions. This study proposed a novel method for large-scale and long-term rainfall erosivity investigations based on the Weather Research and Forecasting (WRF) model, avoiding errors caused by inappropriate rainfall–energy relationships and large-scale interpolation. We adopted three microphysical parameterizations schemes (Morrison, WDM6, and Thompson aerosol-aware) to obtain raindrop size distributions, rainfall kinetic energy, and rainfall erosivity, with validation by two disdrometers and 304 rain gauges around the United Kingdom. Among the three WRF schemes, Thompson aerosol-aware had the best performance compared with the disdrometers at a monthly scale. The results revealed that high rainfall erosivity occurred in the west coast area at the whole country scale during 2013–2017. The proposed methodology makes a significant contribution to improving large-scale soil erosion estimation and for better understanding microphysical rainfall–soil interactions to support the rational formulation of soil and water conservation planning.

1 Introduction

Soil erosion plays a pivotal role in shaping the Earth's physical landscape; however, it can threaten both ecosystems and human societies (Alewell et al., 2015). Accurate quantification of soil loss impact at large spatial scales is therefore important for developing land-use planning and sustainable conservation practices (Bilotta et al., 2012). The soil erosion rate is driven by a combination of factors, including rainfall, topography, soil characteristics, land cover, and land management applications (Wischmeier and Smith, 1958; Panagos et al., 2015b). Among these, rainfall is a driving force that accounts for a large proportion of soil loss throughout most of the world (Panagos et al., 2015b). The erosive force of rainfall with consequent runoff is represented as erosivity of rainfall. This is a crucial factor for estimating soil loss in large-scale soil erosion models, for instance, the Universal Soil Loss Equation (USLE, Wischmeier and Smith, 1978; or RUSLE, Renard et al., 1997), Limburg Soil Erosion Model (LISEM) (De Roo et al., 1996), and USLE-M (Kinnell and Risse, 1998).

Rainfall erosivity estimation involves the microphysical properties of rainfall and rainfall–soil interactions on different time steps (Petan et al., 2010). The impact of rainfall, the main mechanism driving the splashing of soil particles from the soil mass, which leads to soil erosion through soil disintegration and mobilization, relies on the kinetic energy (KE) of raindrop motions (Wischmeier and Smith, 1958; Wang et al., 2014). Robust measurement of raindrop size and terminal velocity is vital for estimating and predicting rainfall erosivity. Many measurements can be used to obtain these two parameters, including the stained paper or flour pellet methods (Marshall and Palmer 1948; Wischmeier and Smith, 1958), high speed cameras (Jones, 1959; Kinnell, 1981; McIsaac, 1990), and disdrometers (Petan et al., 2010; Angulo-Martinez et al., 2012). Accurate measurements of raindrop size can be provided in all their methods, and terminal velocity of raindrops can be further measured by video cameras and disdrometers. Velocity can also be estimated as the function of raindrop diameter from the empirical relationship (Beard, 1976; Atlas and Ulbrich, 1977; Uplinger, 1981; Van Dijk et al., 2002). When using ground observations, rainfall KE can be estimated at a given site.

However, direct measurement of rainfall KE in a large area is difficult because it requires considerable effort, as well as a dense network of expensive instruments that provide accurate outputs (Fornis et al., 2005; Mikoš et al., 2006; Meshesha et al., 2016; Dai et al., 2017). Previous studies have, therefore, mainly employed more readily accessible records like rainfall intensity and attempted to estimate rainfall KE from the empirical relationship of unit KE (ke) with intensity (ke–I). Since Marshall and Palmer (1948) first observed a two-parameter exponential relationship between drop size and intensity, several forms of ke–I mathematical expressions for specific locations and climatic conditions have been proposed, including power-law (Park et al., 1982; Meshesha et al., 2016), linear (Sempere-Torres et al., 1998; Nyssen et al., 2005), polynomial (Carter et al., 1974), logarithmic (Wischmeier and Smith, 1978; Davison et al., 2005; Meshesha et al., 2014), and exponential (Kinnell, 1981; Brown and Foster, 1987) relationships. Among these, the exponential function is preferentially used currently (Van Dijk et al., 2002; Fornis et al., 2005; Petan et al., 2010; Sanchez-Moreno et al., 2012; Lim et al., 2015). Accurate raindrop size distribution (DSD) measured by disdrometers is widely used to derive ke–I relationships (Angulo-Martínez et al., 2016; Meshesha et al., 2016). However, such empirically derived formulas indicate that rainfall ke will increase infinitely with increasing intensity, whereas studies (Rosewell, 1986; Angulo-Martínez et al., 2016; Meshesha et al., 2019) have found that rainfall ke reaches a top value when intensity is around 70 mm h−1 (Hudson, 1963; Wischmeier and Smith, 1978). More importantly, such a ke–I relationship only represents local climate and precipitation microphysics and is valid for such regions. There is great uncertainty associated with rainfall erosivity estimation using this ke–I relationship in a large domain (Angulo-Martínez and Barros, 2015), especially due to the poor spatial and temporal predictability of the ke–I relationship. This has motivated researchers to directly calculate KE based on large-scale DSD measurements.

Ground- and space-based radar can be used to obtain DSD parameters (Atlas et al., 1973; Doelling et al., 1998). For example, the space-borne Dual-frequency Precipitation Radar (DPR) radar containing Ku- and Ka-bands in the Global Precipitation Measurement (GPM) satellite allows researchers to estimate the global three-dimensional spatial distribution of hydrometeors. Unfortunately, ground dual-polarization radars are available in limited areas (Prigent, 2010) with large uncertainties (Dai et al., 2019), and the GPM DPR instrument, which measures DSD with daily or longer temporal resolutions, fails to capture a full storm and meet the requirement for rainfall kinetic estimation. Mesoscale numerical weather prediction models, for instance, the WRF model, can simulate microphysical cloud processes and predict the evolution of particle size distribution through computationally feasible parametrization schemes (Dai and Han, 2014; Brown et al., 2016). DSD on the ground can be derived from the WRF model through consideration of various physical processes, types of hydrometeor, and free degrees of size distributions in hydrometeor. As such, a number of recent studies have investigated the retrieval and uncertainty of DSD parameters by WRF (Gilmore et al., 2004; Ćurić et al., 2009; Brown et al., 2016; Yang et al., 2019).

The WRF model runs with initial and boundary conditions using global reanalysis datasets, such as those of the European Centre for Medium-range Weather Forecasts (ECMWF) and National Centers for Environmental Prediction (NCEP). In other words, WRF-derived DSD can be obtained for any given area with fine spatial and temporal resolutions rather than traditional course linear interpolations. We therefore attempted to estimate rainfall erosivity for the entirety of the United Kingdom (UK) domain using WRF-derived DSD. For comparison, we also calculated interpolated traditional disdrometer-derived rainfall erosivity. To our knowledge, this work is the first attempt to take advantage of a numerical weather prediction model for estimating rainfall erosivity anywhere around the world. The current study contributes to the development of large-scale soil erosion estimation and provides a better comprehension of microphysical rainfall–soil interactions.

2 Methodology

2.1 Disdrometer-based rainfall KE estimation

KE dominates the ability of raindrops to separate soil particles. The KE (e, unit: J) of a raindrop with mass m (g) and terminal velocity v (m s−1) is defined by the following:

(1) e = 1 2 m v 2 .

Assuming a spherical volume for every raindrop shape, the mass of a drop can be calculated from the cube of the diameter D (mm). Because instruments (e.g., disdrometers) generally sample drop size, the mean radius and falling velocity of the corresponding sampling drop-size class are used to represent D and v, expressed as Di and vi, respectively. In such cases, the ei with any drop of a given class is given as follows:

(2) e i = 1 12 10 - 6 π ρ v i 2 D i 3 ,

where ρ is the water density (g cm−3). The sum of the KE of each individual raindrop within a given rain depth that hits a given area defines the total KE. The unit rainfall KE ket in the tth minute (MJ ha−1 mm−1) can be calculated as the sum of each drop KE in each size set, as follows:

(3) ke t = e sum A P t = 1 A P t i = 1 n i N i e i ,

where A represents the sample area of the sensor, Pt is rainfall depth at time t, and Ni is the drop number in class i. The instrument sums up the number of raindrops in each sampling class and produces the raindrop spectra for a time step. Here, we use the term ke to represent the disdrometer-based KE estimated by DSD directly measured every minute. The terminal velocity of a raindrop can be estimated from its power-law empirical relationship with raindrop diameter (Atlas and Ulbrich, 1977), with this considered more suitable for Chilbolton in the UK (Islam et al., 2012):

(4) v Atl = 3.78 D i 0.67

Thus, unit rainfall KE estimates per minute are obtained by replacing vi in Eq. (2) with vAtl.

The other form of rainfall KE is expressed at an event scale and represents the sum of the storm energy covering all time steps covering an event. The individual event energy (MJ ha−1) is calculated as follows:

(5) E = t = 1 nt ke t P t ,

where Pt is the rainfall amount (mm) in the tth minute and nt is the time step number. Historical rainfall data are divided into wet and dry periods. A string of erosive rainfall storms is first extracted through the predefined rules. A continuous 6 h dry-period interval was used to divide rainfall events (Hanel et al., 2016), following the “minimum dry-period duration” definition of a rainfall event (Bonta, 2004). Moreover, a rainfall amount of 12.7 mm was set as the threshold to filter effective rainfall events (Renard et al., 1997).

Rainfall KE is obtained for a given site based on size and velocity of raindrops. When disdrometer data are absent, energy can be estimated from empirical relationships using rainfall intensity I (mm). Five commonly used functions (including exponential, logarithmic, power law, and inverse proportion) have been mentioned in Section 1. Taking the exponential form as an example, the rainfall KE at any location can be estimated as follows:

(6) E max = e max 1 - a e - b I ,

where emax is the mean maximal value of energy measured under high rainfall intensity, and a and b are coefficients modelling the equation curve. Here, minimum KE can be determined by parameters a and emax together, while the overall shape of the curve is modelled by parameter b.

2.2 WRF-based rainfall KE estimation

Differing from disdrometer measurements, the complete DSD cannot be obtained from the WRF model. Instead, the DSD of the microphysical parameterization (MP) scheme is handled with a constrained-gamma distribution model, which is defined as follows:

(7) N ( D ) = N 0 D μ e - λ D ,

where N0, μ, and λ are the intercept, shape, and slope parameters of the DSD. In terms of double-moment bulk schemes, N0 and λ can be abstracted from the number concentration N and predicted mixing ratio q, as shown below:


c and d are the assumed power-law coefficients between diameter and mass (m=cDd), and Γ represents the function in gamma form (Morrison et al., 2009). The value of the shape parameter μ (μ=0) in double-moment schemes is fixed, except for the WRF double-moment 6-class (WDM6) schemes, following gamma distribution, which defined μ=1 (Jung et al., 2010; Johnson et al., 2016).

Because DSD retrieval is sensitive to MPs (Cintineo et al., 2014; Morrison et al., 2015), the WRF model in this study completely or partially incorporates three types of double-moment cloud MP schemes. The Morrison double-moment scheme involves the number concentrations and mixing ratios of multiple hydrometeors (Morrison et al., 2009). Moreover, the WDM6 scheme further considers a prognostic factor to estimate and predict the cloud condensation nuclei (CCN) number concentration (Hong et al., 2010; Lim and Hong, 2010). Finally, the Thompson aerosol-aware (TAA) scheme can predict both ice nuclei (IN) and CCN number concentrations (Thompson and Eidhammer, 2014).

The DSD parameters were thus obtained under the three WRF MPs. For theoretical DSD, ke estimates per minute were obtained by integration of the full raindrop size spectrum using the following:

(10) ke t = 1 A R t 0 N ( D ) 1 12 10 - 6 π ρ v i 2 D i 3 d D .

For the WRF-derived DSD covering the whole study area, there was no need to construct a ke–I relationship to interpolate KE in ungauged areas. The WRF-based rainfall KE under storm event scale is thus given as follows:

(11) E W = t = 1 nt ke t P t .

2.3 Rainfall erosivity estimation

Most storm events have relatively low intensities and KEs with occasional peaks, based on the disdrometer DSD data used to evaluate the rainfall ke–I function. Proper estimation of rainfall erosivity potential should consider total KE over a long period. The rainfall erosivity factor (or R factor) is calculated by a multi-annual average of the total storm erosivity index (Wischmeier and Smith, 1958; Van Dijk et al., 2002), while annual rainfall erosivity R can be obtained using the following:

(12) R = m = 1 M EI 30 m ,

where M is the total number of erosive events within a year. (EI30)m denotes total rainfall kinetic energy and maximum 30 min rainfall intensity recorded within 30 consecutive minutes (unit: mm h−1), respectively, for the mth event.

Wischmeier and Smith (1958) first proposed the use of EI30, as the rainfall erosivity for each event, based on research data from many sources. I30 was calculated to have higher relevance to soil erosion than maximum 5, 15, or 60 min rainfall intensities (Wischmeier and Smith, 1958). The calculation of EI30 initially uses recording-rain gauge data to divide continuous rainfall into time periods with equal rainfall intensity. Though rainfall measurements with high temporal resolutions are required, it is difficult to obtain them from general rainfall measurements. Therefore, short-time equal-interval rainfall data with higher accuracy over multiple years are preferred for estimating EI30. For example, Xie et al. (2016) used 1 min rainfall data instead of recording-rain gauge records. For coarse-resolution, equally spaced data, researchers have proposed a conversion factor to reduce bias error (Weiss, 1964; Williams and Sheridan, 1991).

The rainfall erosivity can be derived from rainfall KE. It plays a main dynamic role in USLE/RUSLE, representing the potential for soil erosion caused by rainfall. To distinguish the disdrometer- and WRF-derived rainfall erosivity in this study, we use the terms RD and RW, respectively.

2.4 Evaluation methods

Because there is no direct way to measure rainfall erosivity across a large area, it is difficult to validate outcomes using observations. However, RD is considered to be relatively accurate due to its specific measurement of raindrops. We therefore assumed that RW values were accurate if they closely matched RD of a given location. A long-term comparison of RW and RD at disdrometer stations was thus conducted to evaluate the validity of RW.

Three indicators were introduced for the evaluation: Pearson's correlation coefficient, mean absolute error (MAE), and coefficient of determination (R2) (Borrelli et al., 2017). The Pearson correlation coefficient is an index used to evaluate the linear correlation between two variables and is defined as follows:

(13) Pearson = n R D i R W i - R D i R W i n R D i 2 - R D i 2 n R W i 2 - R W i 2 ,

where n is the number of variable samples. Because this correlation cannot reveal the absolute bias of rainfall erosivity values, the MAE was also used; this is defined as follows:

(14) MAE = R W i - R D i n .

R2 is an indicator to assess the fit of the trend line, expressed as the ratio of the variance in the dependent variable predicted from the independent variable. It measures the extent to which the model replicates observations based on the proportion of the results interpreted by the model to the total change, written as follows:

(15) R 2 = 1 - SS res SS tot ,

where SSres is the sum of squares of residuals between two variables, and SStot is the total sum of squares.

3 Study area and data sources

The whole of the UK was set as the experimental area for investigating rainfall erosivity estimation. The UK consists of mostly lowland terrain, with a maximum elevation of 1345 m. Water and wind are the most significant forces of soil erosion in the UK. Together, they cause approximately 2.2 million tons of topsoil to be eroded annually, seriously affecting soil productivity, water quality, and aquatic ecosystems through siltation of watercourses (Environment Agency, 2004). According to the Environmental Agency, the total cost of soil erosion in the UK is approximately USD 88 million each year, including an agricultural production loss of USD 17.6 million (O'Neill, 2007). More importantly, the changing climate may exacerbate the degree of erosion. For example, hotter, drier climates make soils more susceptible to wind erosion, and intense storms increase rainfall erosivity (Defra, 2009). Studies of water erosion in England and Wales (Morgan, 1985; Evans, 1990) have found that loose soils (especially sand), such as the soils found in Shropshire and Herefordshire in Wales, are more susceptible to water erosion. In a study of rainfall erosion in Europe, Panagos et al. (2015a) found that the humid Atlantic climate results in highly variable rainfall erosivity, such as higher R factor values in western England and lower values in the eastern UK.

Figure 1Location of rain gauges, Joss–Waldvogel disdrometer (JWD) at Chilbolton Observatory, OTT Parsivel2 disdrometer (OPD) at Bristol Observatory and configurations of domain set-ups in the WRF model.

The gauge datasets used are from the land surface and marine surface measurements datasets (data availability: 1853–present) provided by the UK Met Office. A network of rain gauges covering 304 stations across the whole UK observes continuous rainfall data in hours (Fig. 1). The base data of most stations comprise the times of each tip (0.2 mm per tip), converted into 1 h rain accumulations. The rainfall observations are not always valid for each hour at each station. The hourly grid-based rainfall maps are then calculated based on ordinary kriging interpolation of rain gauge network data to obtain the spatial distribution of rainfall for each time step, as inputs for rainfall erosivity estimation. This wide-range-use geostatistical approach can account for both the distance and pairwise spatial relationship between points through variograms. The precipitation interpolation method uses sample gauge points taken at different locations and creates a continuous surface to achieve an accurate spatial variation estimation of rainfall patterns.

We used data from two disdrometers in southern England. The first was Chilbolton station (5108 N, 126 W), with an impact-type Joss–Waldvogel disdrometer (JWD) mainly used to compute rainfall erosivity. It can measure drop sizes from 0.3 to 5.0 mm in 127 bins. The sampling period and collector area were 10 s and 50 cm2, respectively. Data were available for April 2003 to July 2018. The second was the University of Bristol station (5127 N, 236 W), with an OTT Parsivel2 disdrometer (OPD). Data were available for November 2015 to December 2018. This disdrometer subdivides particles into appropriate classes and has a nominal cross-sectional area of 54 cm2. The 10 s period measurement data from the two disdrometers were averaged into a 1 min period to filter out time variations (Montopoli et al., 2008; Islam et al., 2012; Song et al., 2017).

Meteorological data come from the ERA-Interim dataset, a global atmosphere re-analysis product, generated by the ECMWF. For the scientific community, ERA-Interim is considered to be one of the most important atmospheric datasets, with its data-rich period available since 1979 and updated in current time (Dee et al., 2011). The Integrated Forecasting System released in 2006 contains a 12 h window-derived 4-D variational analysis, driving the data assimilation system to generate ERA-Interim. The dataset covers 60 vertical classes of approximately 80 km from the ground to 0.1 hPa. The gridded binary format is used to store data for 3 months in a separate file. A data processing scheme was established to collect and retrieve ERA-Interim data of each rainfall event.

The rain gauge and Chilbolton disdrometer datasets can be obtained from British Atmospheric Data Centre in the National Centre for Atmospheric Science research centre (Meteorological Office, 2012). ERA-Interim data can be obtained from the ECMWF Public Dataset website (, last access: 15 November 2020). Considering the availability of the above datasets and model requirements, we mainly used data covering the period 2004–2017.

4 Results

4.1 Empirically derived rainfall erosivity estimation

To evaluate the RW, the raindrop spectrum collected by the Chilbolton station disdrometer is used to estimate rainfall KE first. The key to estimating rainfall KE by disdrometer lies in building an empirical relationship between rainfall amount and KE. We used DSD measurements from 2004 to 2013 to establish five empirical relationships between unit rainfall kinetic energy (ke) and intensity (I) (Table 1) and used 2014–2017 data for the cross-validation. It can be seen from Table 1 that the inverse proportional relationship (Eq. III) had the worst performance, in that both the calibration and validation R2 values were < 0.3. The values of the other equations were > 0.48, among which the exponential formula (Eq. I) had the highest calibration R2 (0.50) and validation R2 (0.45), respectively. In addition, the power-law formula (Eq. V) showed a similar performance to the exponential formula at rainfall intensities < 5 mm h−1. However, the power-law formula also had a continuously increasing trend, which may not be suitable for high intensities. Figure 2 shows the ke–I relationship and five fitted curves at Chilbolton station. It can be seen that the two logarithmic curves (Eqs. II and IV) invariably overlap. The logarithmic form has been used for a long time in USLE (Wischemier and Smith, 1978). It describes ke well at both low and high I but does not have an upper limit. The power-law curve (Eq. V) can predict ke well at lower I but overestimates ke at high I. The exponent-based relationship (Eq. I) is widely used in the literature and in forecast models such as RUSLE (Renard et al., 1997), which fits the data particularly well in Fig. 2. Even though ke in an exponential curve has a minimum value at very low I, it also should be noted that higher rainfall intensities are much more important in determining overall storm energy than lower intensities. Therefore, we adopted it here as the empirical formula to estimate rainfall erosivity in the UK.

Table 1Relationship of ke–I at Chilbolton Station (2004–2013).

Download Print Version | Download XLSX

Figure 2Number of minutes per intensity class (x axis) and ke class (y axis) with five fitted ke–I curves at Chilbolton station (2004–2013), plotted on linear (a) and logarithmic (b) intensity scales.


Based on rainfall KE, the point RD can be obtained at a disdrometer location. In the current study, we established a method to estimate the R using 60 min rainfall data. EI30 obtained from 1 min DSD data was considered as the standard R at Chilbolton Station. Hourly rain gauge data at the same location were used to calculate (EI30)60, which refers to EI30 calculated from 60 min data. The regression relationship between EI30 and (EI30)60 was then established. The (EI30)60 of each month, obtained from the 60 min rainfall data of the Chilbolton Station rain gauge in 2004–2013, was calculated. The regression relationship between the monthly sum of (EI30)60 and the standard monthly EI30 from DSD was calculated to obtain a coefficient of 1.836. Rainfall erosivity can subsequently be calculated by multiplying (EI30)60 by the coefficient.

Figure 3Gauge-based interpolation maps of annual rainfall amount (Rain), rainfall kinetic energy (E), and rainfall erosivity (R) in 2013–2017.

Beyond assuming that the disdrometer-derived ke–I relationship can be applied to a whole study area, point rainfall measurements must be interpolated to obtain areal rainfall values in traditional rainfall erosivity estimation. We obtained 60 min rainfall data from 304 rain gauges around the UK from 2004 to 2017. Note that not all rain gauges were available for the whole period (available gauges each year are indicated in Fig. 3). We used the ordinary kriging interpolation method to obtain the spatial distribution of rainfall for each time step. This wide-range-use geostatistical approach can account for both the distance and pairwise spatial relationship between points through variograms. Figure 3 shows the results of annual rainfall (Rain), annual rainfall kinetic energy (E), and annual rainfall erosivity (R) for different years. The distribution trends of Rain, E, and R were similar and positively correlated except for certain locations or periods. For instance, in 2013, Rain in the northwestern UK decreased from west to east, while E and R decreased from south to north; furthermore, areas with large E and R values in the southeastern UK could not be directly observed from the rain map.

Figure 4Minutes number per intensity class (x axis) and ke class (y axis) with fitted ke–I curves at Bristol station (2015–2018), plotted on linear (a) and logarithmic (b) intensity scales.


The key concern in traditional rainfall erosivity estimation is the spatial predictability of the ke–I relationship. To verify the regional reliability of this relationship, we used data from a newer disdrometer located at the University of Bristol, approximately 87 km from Chilbolton Station. The validation data at Bristol Station discontinuously covered the period 2016–2019. Figure 4 shows the exponential relationship of ke–I at Bristol station, which differed substantially from that based on data from Chilbolton station. A comparison of the modelled and observed event rainfall erosivity is shown in Fig. 5. The modelled erosivity of rainfall event was not consistent with the observed event rainfall erosivity. The linear regression coefficient between these values was > 1.2, which was the result of the low ke for Bristol Station, and R2 was < 0.85, indicating considerable uncertainty associated with disdrometer-based rainfall erosivity estimation.

Figure 5Comparison of observed and modelled event rainfall erosivity at Bristol Station, covering the period of 2016–2019.


In summary, the point rainfall erosivity estimated by disdrometer is considered to be accurate compared to other methods. However, a large-scale rainfall erosivity through a simple interpolation of rainfall KE is subjected to a significant uncertainty. In the following analysis, the point RD is used to appraise the performance of the proposed WRF-based estimated method, and the RD in the whole UK is only used for a general comparison of spatial and temporal distribution of rainfall erosivity.

4.2 Rainfall and DSD estimation by WRF

We used the WRF model version 3.8, which has an Advanced Research WRF dynamical core, to downscale the ERA-Interim reanalysis data. The double-nested domain configuration used in the WRF model was centred at 5519 N, 221 W and applied at a downscaling ratio of 1 : 5, a finest grid of 5 km, and a temporal resolution of 1 h. Table 2 lists the detailed parameters used in this domain configuration. With the top pressure level set at 50 hPa in each, both domains include 28 vertical levels. To obtain favourable initial weather conditions, the model ran continuously to obtain 5 years of WRF simulation results.

Table 2The configurations of the WRF model for two nested domains.

Download Print Version | Download XLSX

Simulations were performed using three different bulk double-moment MPs: the Morrison (Morrison et al., 2009), WDM6 (Hong et al., 2010; Lim and Hong, 2010), and TAA (Thompson and Eidhammer, 2014) schemes. All three can predict the number concentration and hydrometeors mixing ratio for each time step. The WDM6 scheme also predicts the number concentration of CCN (Hong et al., 2010; Lim and Hong, 2010), while the TAA scheme is able to predict both IN and CCN number concentrations (Thompson and Eidhammer, 2014). Additionally, other physical parameterizations include the Dudhia shortwave radiation scheme (Dudhia, 1989), Mellor–Yamada–Janjic planetary boundary layer scheme (Janjić, 1994), RRTM longwave radiation scheme (Mlawer et al., 1997), Noah land-surface model (Ek et al., 2003), and Kain–Fritsch cumulus scheme (Kain, 2004).

The median volume diameter parameter (D0) and generalized intercept parameter (NW) are generally used in the DSD model of WRF (Islam et al., 2012).


where Dm is the mass-weighted mean diameter. The f(μ) is a function of the shape parameter μ. The parameter μ is assumed as zero or one (based on microphysical scheme configuration) in WRF. Figure 6 displays the spatial distribution of D0 and generalized intercept parameter NW for a given day with rainfall countrywide (January 10, 2013). D0 and NW had similar patterns and were mainly distributed across the southwestern and northeastern UK. The white strip in the middle of the maps in Fig. 6 represents an area that received no rain. However, the three MPs yielded large differences; D0 of MP-TAA was the highest among three MPs, whereas NW of MP-WDM6 was significant larger than the others. In addition, D0 and NW did not consistently show a positive correlation. The different MP estimation results underscore the complexity of the rainfall process, which is the reason we estimated rainfall KE using WRF schemes instead of traditional formulas.

Figure 6Map of average WRF DSD D0 and NW (10 January 2013).

4.3 Comparison of WRF- and disdrometer-derived rainfall erosivity at Chilbolton station

With the WRF-based rainfall intensity and DSD estimations, rainfall erosivity was derived using Eqs. (10)–(12). Hereafter, this is referred to as RW, which is further distinguished based on the three MP schemes used: RW-Morrison, RW-WDM6, and RW-TAA. Figure 7 compares disdrometer- and WRF-derived monthly rainfall erosivity estimations at Chilbolton station for the period 2014–2017. The general patterns of the four rainfall erosivity values were similar. RW-Morrison tended to be larger than RD in some months, whereas RW-TAA matched the RD value relatively well for smaller values. Because WRF data were taken from a 2 km × 2 km grid around Chilbolton station, there was a spatial error in addition to the systematic error of estimating rainfall erosivity. Based on the 4-year data, the study area is rainy throughout the year with little monthly R, or seasonal pattern changes (Fig. 8), influenced by the temperate oceanic climate. Figure 8 also indicated that through the perspective of monthly average results, RW-WDM6 values are low, RW-TAA has a good similarity with low RD, and RW-Morrison is the closest to RD in value.

Figure 7Comparison of disdrometer- and WRF-derived monthly rainfall erosivity estimations at Chilbolton station (2014–2017).


Figure 8Comparison of disdrometer- and WRF-derived average monthly rainfall erosivity estimations at Chilbolton station (2014–2017).


Table 3 shows the correlation indicator results between monthly RD and the three types of RW at Chilbolton station. The Pearson correlation coefficients generally exceeded 0.7, supporting the potential utility of WRF-based estimation. In terms of MAE, RW-TAA had the best performance (6.51), whereas RW-Morrison and RW-WDM6 showed slightly worse performance (approximately 8). Among the three schemes, RW-TAA had the best fit with RD. The indicators and comparison results suggest that the deviations in results need to be considered; therefore, a method of bias elimination is described in Sect. 4.4.

Table 3Indicator comparison between disdrometer-derived rainfall erosivity RD and three types of WRF-derived rainfall erosivity at Chilbolton station on monthly scale (2014–2017).

Download Print Version | Download XLSX

4.4RW estimation for the whole UK

The RW at Chilbolton station showed obvious systematic deviations compared with the disdrometer-derived results (see Sect. 4.2 and 4.3). Simple bias correction was therefore applied to adjust the individual storm KE estimations of RW. The biases from dividing average RW-Morrison, RW-WDM6, and RW-TAA by average RD during 2014–2017 were 0.55, 0.20, and 0.36, respectively.

Figure 9WRF-derived annual rainfall erosivity maps of the whole UK for different years.

The rainfall erosivity distribution for the whole UK was then obtained. Figure 9 shows the distribution of RW at the annual scale covering the period 2013–2017. The pattern of the rainfall erosivity maps showed a general regional-dominant characteristic. For example, it always decreased from west to east, predominantly shaped by orography. Affected by the prevailing westerly winds, there was abundant rainfall in the western and northern mountains, as indicated by high rainfall KE values in these regions. In addition, among the study years, 2014 and 2015 showed higher national rainfall erosivity, with a large range in the west coast area.

Figure 10The 5-year (2013–2017) average annual rainfall erosivity maps based on WRF grids and rain gauge interpolation.

Figure 10 shows the average R distribution for 2013–2017 estimated by rain gauges and WRF MPs. WRF grids could cover all regions in the UK evenly, offering more detailed erosivity results, especially in the mountainous northwestern region. Here, values of an average-R map calculated by rain gauges were much higher than three types of RW, although they all have R that decreased from west to east. Note that the ke–I empirical equation at Chilbolton station, used in the whole UK, will not always be accurate in regions with different rainfall characteristics. In terms of RW results, the three MPs obtained the same spatial pattern in rainfall erosivity, where RW-WDM6 yielded the greatest geographical difference. It is clear that the proposed WRF-based estimated method can capture more details of the spatial change of rainfall erosivity compared with the traditional disdrometer-based method.

The highest rainfall erosivity regions in the UK are concentrated in the mountainous areas along the western coast, related to their rainfall system. The moist air brought by the prevailing westerly wind from the Atlantic Ocean moves from west to east across the UK and rises when it encounters the mountains of western England. Therefore, the mountainous regions along the UK western coast have the highest rainfall amount and rainfall erosivity in the UK. In addition, western Scotland is under the subpolar oceanic climate, which enhances its humidity. On the contrary, eastern Scotland and northeastern England are more likely to be exposed to continental polar air mass, which brings dry and cold air and lower rainfall erosivity.

To evaluate the change in rainfall erosivity with time in the UK, the average value of all the WRF grids covering the whole UK was calculated over 2013–2017 (Fig. 11). The average RW trends of RW-Morrison and RW-TAA were similar, both increasing from a minimum in 2013 to a maximum in 2014 and then gradually decreasing from 2014 to 2017. The red line in Fig. 11 indicates a series of mean values of the three MPs results, which varied from 367.82 to 516.00 MJ mm ha−1 h−1 yr−1 (mean: 432.16 MJ mm ha−1 h−1 yr−1).

Figure 11WRF-derived average annual rainfall erosivity of all the WRF grids covering the whole UK (2013–2017).


The maximum values for RW-Morrison and RW-TAA occurred in 2014, whereas that of RW-WDM6 occurred in 2015. A sequence of extreme weather events occurred in the UK in 2014, including major winter storms in late January to mid-February, which caused widespread flooding and other economic losses and greatly increased rainfall erosivity that year. However, the gauge-based interpolation map shows the average annual rainfall amounts for the years 2013–2017 were 884.9, 1014.0, 1008.5, 894.9, and 937.3 mm, respectively. The large rainfall erosivity difference between 2014 and 2015 and the 2 years with similar rainfall amount indicate that much rainfall erosion occurs during the rainfall events of high intensity instead of simply high rainfall amount. A more notable variation pattern of rainfall erosivity may be found with longer simulation. The strength of the proposed method lies in its ability to estimate large covering and long-term rainfall erosivity.

5 Discussion

Compared to the previous large-scale rainfall erosivity studies based on empirical formulas and spatial interpolation, this study presents a WRF-driven approach directly using the simulated rainfall microphysical variables. As demonstrated in the literature, the relation between rainfall intensity and erosivity is not straightforward (Panagos et al., 2015a, 2017; Ballabio et al., 2017). However, although all works show that rainfall erosivity decreases from west to east in UK, previous studies (Panagos et al., 2015a; Naipal et al., 2015) using traditional methods lead to an overestimation of rainfall erosivity, which may be due to parameter a in the universal KE–I relationship being too high for the UK. Considering the 5 years (2013–2017) as a whole, the averaged RW-Morrison, RW-WDM6, and RW-TAA factor in each grid can be calculated. Nationally, the mean values of the three RW factors are 446.57, 640.92, and 416.35 MJ mm ha−1 h−1 yr−1, and their coefficients of variation (CVs) are 0.56, 0.81, and 0.59, respectively. Compared with the outcomes (mean R factor = 746.6 MJ mm ha−1 h−1 yr−1, CV = 0.81) of the Panagos et al. (2015a), the R factors of WDM6 scheme are quite similar, while other schemes have relatively low R factors and low CVs.

Although an acceptable rainfall erosivity estimation is obtained using the WRF model, some uncertainties associated with it cannot be ignored. For example, as the MPs of WRF were closely related to DSD, improper determination of MPs will introduce additional uncertainty. The marked discrepancy among the three schemes (especially between Morrison and the others) in this study underscored the possible uncertainty associated with RW. The reliability of the WRF model is heavily dependent on the model-driving initial data provided by mesoscale or global models and complicated scheme setting and parameter adjustment (Liu et al., 2012; Thompson and Eidhammer, 2014; Kumar et al., 2017). However, numerous uncertainties are observed in the parameterization of the WRF simulation, and the choice of microphysical schemes has a significant influence on the inverted DSD (Ćurić et al., 2009; Yang et al., 2019). Therefore, combining the DSDs obtained by an increasing number of disdrometers and the WRF model is valuable. For example, the Disdrometer Verification Network (DiVeN) in the UK (Pickering et al., 2019) started in February 2017 can be introduced to support and improve our estimation in future studies. Moreover, the measurement error by disdrometer may also contaminate the evaluation process. For example, when comparing the observed raindrop velocities based on the disdrometer at Bristol station with their empirical values, we observed dispersion of raindrops, with a number of drops showing significant deviations. This velocity distribution resulted in uncertainty in ke estimation.

Soil erosion in the UK is dominated by water erosion (10–30 t km−2 yr−1), especially in areas with abundant rainfall in Scotland, where the soil loss rate is approximately 5–10 times that of dry areas (McManus and Duck, 1996). Thus, it is significant to estimate rainfall erosivity to elucidate the microphysical characteristics of rainfall and rainfall–soil interactions. Benaud et al. (2020) collated empirical soil erosion observations from UK-based studies into a geodatabase. However, there is a limitation that this database does not cover the entirety of the UK, especially the limited records in northern Scotland. In our future work, we propose to compare the soil loss database with our estimated soil loss using WRF DSD-based rainfall erosivity and a soil erosion model (such as RUSLE). We believe that we can not only better analyse the impact of rainfall and rainfall erosivity on the UK soil loss, but also help to better understand microphysical rainfall–soil interactions to support the rational formulation of soil and water conservation planning.

6 Conclusions

This study presented a novel method for large-scale rainfall KE and erosivity estimation based on high-resolution, WRF-derived DSDs. Three microphysical parameterizations schemes (Morrison, WDM6, and TAA) were designed to obtain raindrop size distributions, rainfall KE and rainfall erosivity for the entire of the UK covering the period of 2013–2017. With validation from the long-term observations of a disdrometer, the WRF-based rainfall erosivity exhibited an acceptable performance at Chilbolton station. Among the three WRF schemes, TAA exhibited the most superior performance and was recommended for future investigation. The results revealed that high rainfall erosivity occurred in the west coast area of the UK. Compared with the traditional empirical method, the proposed method can explain rainfall erosivity from a microphysical perspective and reflect more spatial variation because of changes in rainfall KE at the whole-country scale. Therefore, the development of a numerical weather prediction model offers an opportunity to better understand rainfall erosivity directly from its true definition. More importantly, because the WRF model is able to be driven by the global reanalysis data to obtain large-scale rainfall kinetic information, the proposed scheme can be easily applied to other regions, especially in ungauged areas.

Some problems remain with the proposed scheme, as discussed in Sect. 5. Some of the problems, such as temporal downscaling of rainfall and point-to-area representative error by WRF, may introduce further uncertainty. This should be put in perspective for future work. It is expected that further exploration of research areas with different climatic and geographical characteristics would help us to establish a greater degree of accuracy on this matter.

Data availability

The rain gauge datasets and Chilbolton disdrometer data were sourced from the Met Office Integrated Data Archive System (MIDAS). Both datasets are available from the NCAS British Atmospheric Data Centre (, Met Office, 2006; and, Science and Technology Facilities Council et al., 2003). The ERA-Interim data driving the WRF model can be downloaded from the ECMWF Public Datasets web interface (, last access: 15 November 2020).

Author contributions

QD and JZ carried out the experiments, analysed the data, and prepared the paper with contributions from all the co-authors. ShuZ modified the text and provided financial support. ShaZ carried out quality checks of WRF. GL and DH principally conceived the idea and design of the study.

Competing interests

The authors declare that they have no conflict of interest.


This work was supported by the National Natural Science Foundation of China (nos. 41871299 and 41771424) and the National Key R & D Program of China (nos. 2018YFB0505500, 2018YFB0505502). The authors acknowledge the British Atmospheric Data Centre and the European Centre for Medium-range Weather Forecasts as the sources of data used in the study.

Financial support

This research has been supported by the National Natural Science Foundation of China (grant nos. 41871299 and 41771424) and the National Key R&D Program of China (grant nos. 2018YFB0505500 and 2018YFB0505502).

Review statement

This paper was edited by Rohini Kumar and reviewed by three anonymous referees.


Alewell, C., Egli, M., and Meusburger, K.: An attempt to estimate tolerable soil erosion rates by matching soil formation with denudation in Alpine grasslands, J. Soil Sediment., 15, 1383–1399,, 2015. 

Angulo-Martínez, M. and Barros, A.: Measurement uncertainty in rainfall kinetic energy and intensity relationships for soil erosion studies: An evaluation using PARSIVEL disdrometers in the Southern Appalachian Mountain, Geomorpholgy, 228, 28–40,, 2015. 

Angulo-Martinez, M., Beguería, S., Navas, A., and Machin, J.: Splash erosion under natural rainfall on three soil types in NE Spain, Geomorphology, 175, 38–44,, 2012. 

Angulo-Martínez, M., Beguería, S., and Kyselý, J.: Use of disdrometer data to evaluate the relationship of rainfall kinetic energy and intensity (KE-I), Sci. Total Environ., 568, 83–94,, 2016. 

Atlas, D., Srivastava, R., and Sekhon, R. S.: Doppler radar characteristics of precipitation at vertical incidence, Rev. Geophys., 11, 1–35,, 1973. 

Atlas, D. and Ulbrich, C. W.: Path-and area-integrated rainfall measurement by microwave attenuation in the 1–3 cm band, J. Appl. Meteorol., 16, 1322–1331,<1322:paairm>;2, 1977. 

Ballabio, C., Borrelli, P., Spinoni, J., Meusburger, K., Michaelides, S., Beguería, S., Klik, A., Petan, S., Janeček, M., Olsen, P., Aalto, J., Lakatos, M., Rymszewicz, A., Dumitrescu, A., Tadić, M. P., Diodato, N., Kostalova, J., Rousseva, S., Banasik, K., Alewell, C., and Panagos, P.: Mapping monthly rainfall erosivity in Europe, Sci. Total Environ., 579, 1298–1315,, 2017. 

Beard, K. V.: Terminal velocity and shape of cloud and precipitation drops aloft, J. Atmos. Sci., 33, 851–864,<0851:tvasoc>;2, 1976. 

Benaud, P., Anderson, K., Evans, M., Farrow, L., Glendell, M., James, M. R., Quine, T. A., Quinton, J. N., Rawlins, B., Rickson, R. J., and Brazier, R. E.: National-scale geodata describe widespread accelerated soil erosion, Geoderma, 371, 114378,, 2020. 

Bilotta, G., Grove, M., and Mudd, S.: Assessing the significance of soil erosion, T. I. Brit. Geogr., 37, 342–345,, 2012. 

Bonta, J.: Development and utility of Huff curves for disaggregating precipitation amounts, Appl. Eng. Agric., 20, 641,, 2004. 

Borrelli, P., Robinson, D. A., Fleischer, L. R., Lugato, E., Ballabio, C., Alewell, C., Meusburger, K., Modugno, S., Schütt, B., and Ferro, V., Bagarello, V., Oost, K. V., Montanarella, L., and Panagos, P.: An assessment of the global impact of 21st century land use change on soil erosion, Nat. Commun., 8, 1–13,, 2017. 

Brown, B. R., Bell, M. M., and Frambach, A. J.: Validation of simulated hurricane drop size distributions using polarimetric radar, Geophys. Res. Lett., 43, 910–917,, 2016. 

Brown, L. and Foster, G.: Storm erosivity using idealized intensity distributions, T. ASAE, 30, 379–386,, 1987. 

Carter, C. E., Greer, J., Braud, H., and Floyd, J.: Raindrop characteristics in south central United States, T. ASAE, 17, 1033–1037,, 1974. 

Cintineo, R., Otkin, J. A., Xue, M., and Kong, F.: Evaluating the performance of planetary boundary layer and cloud microphysical parameterization schemes in convection-permitting ensemble forecasts using synthetic GOES-13 satellite observations, Mon. Weather Rev., 142, 163–182,, 2014. 

Ćurić, M., Janc, D., Vučković, V., and Kovačević, N.: The impact of the choice of the entire drop size distribution function on Cumulonimbus characteristics, Meteorol. Z., 18, 207–222, 2009. 

Dai, Q. and Han, D.: Exploration of discrepancy between radar and gauge rainfall estimates driven by wind fields, Water. Resour. Res., 50, 8571–8588,, 2014. 

Dai, Q., Bray, M., Zhuo, L., Islam, T., and Han, D.: A scheme for raingauge network design based on remotely-sensed rainfall measurements, J. Hydrometeorol., 18, 363–379,, 2017. 

Dai, Q., Yang, Q., Han, D., Rico-Ramirez, M. A., and Zhang, S.: Adjustment of radar-gauge rainfall discrepancy due to raindrop drift and evaporation using the Weather Research and Forecasting model and dual-polarization radar, Water. Resour. Res., 55, 9211–9233,, 2019. 

Davison, P., Hutchins, M. G., Anthony, S. G., Betson, M., Johnson, C., and Lord, E. I.: The relationship between potentially erosive storm energy and daily rainfall quantity in England and Wales, Sci. Total Environ., 344, 15–25,, 2005. 

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., Monge-Sanz, B. M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.-N., and Vitart, F.: The ERA-Interim reanalysis: configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597,, 2011. 

De Roo, A. P. J., Wesseling, C. G., and Ritsema, C. J.: LISEM: a single-event physically based hydrological and soil erosion model for drainage basins. I: theory, input and output, Hydrol. Process., 10, 1107–1117,<1107::aid-hyp415>;2-4, 1996. 

Defra: Safeguarding our soils – a strategy for England, Department for Environment, Food & Rural Affairs, London, UK, available at: (last access: 16 October 2020), 2009. 

Doelling, I. G., Joss, J., and Riedl, J.: Systematic variations of Z–R-relationships from drop size distributions measured in northern Germany during seven years, Atmos. Res., 47, 635–649,, 1998. 

Dudhia, J.: Numerical study of convection observed during the winter monsoon experiment using a mesoscale two-dimensional model, J. Atmos. Sci., 46, 3077–3107,<3077:nsocod>;2, 1989. 

Ek, M., Mitchell, K., Lin, Y., Rogers, E., Grunmann, P., Koren, V., Gayno, G., and Tarpley, J.: Implementation of Noah land surface model advances in the National Centers for Environmental Prediction operational mesoscale Eta model, J. Geophys. Res.-Atmos., 108, 8851,, 2003. 

Environment Agency: The state of soils in England and Wales, Environment Agency, Bristol, UK, available at: (last access: 16 October 2020) 2004. 

Evans, R.: Soils at risk of accelerated erosion in England and Wales, Soil Use Manage., 6, 125–131,, 1990. 

Fornis, R. L., Vermeulen, H. R., and Nieuwenhuis, J. D.: Kinetic energy–rainfall intensity relationship for Central Cebu, Philippines for soil erosion studies, J. Hydrol., 300, 20–32,, 2005. 

Gilmore, M. S., Straka, J. M., and Rasmussen, E. N.: Precipitation uncertainty due to variations in precipitation particle parameters within a simple microphysics scheme, Mon. Weather Rev., 132, 2610–2627,, 2004. 

Hanel, M., Máca, P., Bašta, P., Vlnas, R., and Pech, P.: The rainfall erosivity factor in the Czech Republic and its uncertainty, Hydrol. Earth. Syst. Sc., 20, 4307–4322,, 2016. 

Hong, S.-Y., Lim, K.-S. S., Lee, Y.-H., Ha, J.-C., Kim, H.-W., Ham, S.-J., and Dudhia, J.: Evaluation of the WRF double-moment 6-class microphysics scheme for precipitating convection, Adv. Meteorol., 2010, 707253,, 2010. 

Hudson, N.: Raindrop size distribution in high intensity storms, Rhod. J. Agr. Res., 1, 6–11, 1963. 

Islam, T., Rico-Ramirez, M. A., Thurai, M., and Han, D.: Characteristics of raindrop spectra as normalized gamma distribution from a Joss–Waldvogel disdrometer, Atmos. Res., 108, 57–73,, 2012. 

Janjić, Z. I.: The step-mountain eta coordinate model: further developments of the convection, viscous sublayer, and turbulence closure schemes, Mon. Weather Rev., 122, 927–945,<0927:tsmecm>;2, 1994. 

Johnson, M., Jung, Y., Dawson, D. T., and Xue, M.: Comparison of simulated polarimetric signatures in idealized supercell storms using two-moment bulk microphysics schemes in WRF, Mon. Weather Rev., 144, 971–996,, 2016. 

Jones, D. M. A.: The shape of raindrops, J. Atmos. Sci., 16, 511–515,<0504:TSOR>2.0.CO;2, 1959. 

Jung, Y., Xue, M., and Zhang, G.: Simulations of polarimetric radar signatures of a supercell storm using a two-moment bulk microphysics scheme, J. Appl. Meteorol. Clim., 49, 146–163,, 2010. 

Kain, J. S.: The Kain–Fritsch convective parameterization: an update, J. Appl. Meteorol., 43, 170–181,<0170:tkcpau>;2, 2004. 

Kinnell, P. I. A.: Rainfall intensity-kinetic energy relationships for soil loss prediction, Soil. Sci. Soc. Am. J., 45, 153–155,, 1981. 

Kinnell, P. I. A. and Risse, L. M.: USLE-M: empirical modeling rainfall erosion through runoff and sediment concentration, Soil. Sci. Soc. Am. J., 62, 1667–1672,, 1998. 

Kumar, P., Kishtawal, C. M., and Pal, P. K.: Impact of ECMWF, NCEP, and NCMRWF global model analysis on the WRF model forecast over Indian Region, Theor. Appl. Climatol., 127, 143–151,, 2017. 

Lim, K.-S. S. and Hong, S.-Y.: Development of an effective double-moment cloud microphysics scheme with prognostic cloud condensation nuclei (CCN) for weather and climate models, Mon. Weather Rev., 138, 1587–1612,, 2010. 

Lim, Y. S., Kim, J. K., Kim, J. W., Park, B. I., and Kim, M. S.: Analysis of the relationship between the kinetic energy and intensity of rainfall in Daejeon, Korea, Quatern. Int., 384, 107–117,, 2015. 

Liu, J., Bray, M., and Han, D.: Exploring the effect of data assimilation by WRF-3DVar for numerical rainfall prediction with different types of storm events, Hydrol. Process., 27, 3627–3640,, 2012. 

Marshall, J. S. and Palmer, W. M. K.: The distribution of raindrops with size, J. Meteorol., 5, 165–166,<0165:TDORWS>2.0.CO;2, 1948. 

McIsaac, G.: Apparent geographic and atmospheric influences on raindrop sizes and rainfall kinetic energy, J. Soil Water Conserv., 45, 663–666, 1990. 

McManus, J. and Duck, R. W.: Regional variations of fluvial sediment yield in eastern Scotland, in: Erosion and Sediment Yield: Global and Regional Perspectives, edited by: Walling, D. and Webb, B., IAHS Press, Wallingford, UK, 157–161, 1996. 

Meshesha, D. T., Tsunekawa, A., and Haregeweyn, N.: Influence of raindrop size on rainfall intensity, kinetic energy, and erosivity in a sub-humid tropical area: a case study in the northern highlands of Ethiopia, Theor. Appl. Climatol., 136, 1221–1231,, 2019. 

Meshesha, D. T., Tsunekawa, A., Tsubo, M., Haregeweyn, N., and Adgo, E.: Drop size distribution and kinetic energy load of rainfall events in the highlands of the Central Rift Valley, Ethiopia, Hydrolog. Sci. J., 59, 2203–2215,, 2014. 

Meshesha, D. T., Tsunekawa, A., Tsubo, M., Haregeweyn, N., and Tegegne, F.: Evaluation of kinetic energy and erosivity potential of simulated rainfall using Laser Precipitation Monitor, Catena, 137, 237–243,, 2016. 

Met Office: MIDAS UK Hourly Rainfall Data, NCAS British Atmospheric Data Centre, available at: (last access: 15 November 2020), 2006. 

Meteorological Office: Met Office Integrated Data Archive System (MIDAS) land and marine surface stations data (1853-current), NCAS British Atmospheric Data Centre, available at: (last access: 16 October 2020) 2012. 

Mikoš, M., Jošt, D., and Petkovšek, G.: Rainfall and runoff erosivity in the alpine climate of north Slovenia: a comparison of different estimation methods, Hydrolog. Sci. J., 51, 115–126,, 2006. 

Mlawer, E. J., Taubman, S. J., Brown, P. D., Iacono, M. J., and Clough, S. A.: Radiative transfer for inhomogeneous atmospheres: RRTM, a validated correlated-k model for the longwave, J. Geophys. Res.-Atmos., 102, 16663–16682,, 1997. 

Montopoli, M., Marzano, F. S., and Vulpiani, G.: Analysis and synthesis of raindrop size distribution time series from disdrometer data, IEEE T. Geosci. Remote, 46, 466–478,, 2008. 

Morgan, R. P. C.: Assessment of soil erosion risk in England and Wales, Soil Use Manage., 1, 127–131,, 1985. 

Morrison, H., Milbrandt, J. A., Bryan, G. H., Ikeda, K., Tessendorf, S. A., and Thompson, G.: Parameterization of cloud microphysics based on the prediction of bulk ice particle properties. Part II: Case study comparisons with observations and other schemes, J. Atmos. Sci., 72, 312–339,, 2015. 

Morrison, H., Thompson, G., and Tatarskii, V.: Impact of cloud microphysics on the development of trailing stratiform precipitation in a simulated squall line: Comparison of one-and two-moment schemes, Mon. Weather Rev., 137, 991–1007,, 2009. 

Naipal, V., Reick, C., Pongratz, J., and Van Oost, K.: Improving the global applicability of the RUSLE model – adjustment of the topographical and rainfall erosivity factors, Geosci. Model Dev., 8, 2893–2913,, 2015. 

Nyssen, J., Vandenreyken, H., Poesen, J., Moeyersons, J., Deckers, J., Haile, M., Salles, C., and Govers, G.: Rainfall erosivity and variability in the Northern Ethiopian Highlands, J. Hydrol., 311, 172–187,, 2005. 

O'Neil, D.: The total external environmental costs and benefits of agriculture in the UK, Environment Agency, Bristol, UK, 2007. 

Panagos, P., Ballabio, C., Borrelli, P., Meusburger, K., Klik, A., Rousseva, S., Tadić, M. P., Michaelides, S., Hrabalíková, M., Olsen, P., Aalto, J., Lakatos, M., Rymszewicz, A., Dumitrescu, A., Beguería, S., and Alewell, C.: Rainfall erosivity in Europe, Sci. Total Environ., 511, 801–814,, 2015a. 

Panagos, P., Borrelli, P., Poesen, J., Ballabio, C., Lugato, E., Meusburger, K., Montanarella, L., and Alewell, C.: The new assessment of soil loss by water erosion in Europe, Environ. Sci. Policy, 54, 438–447,, 2015b. 

Panagos, P., Borrelli, P., Meusburger, K., Yu, B., Klik, A., Lim, K. J., Yang, J. E., Ni, J., Miao, C., Chattopadhyay, N., Sadeghi, S. H., Hazbavi, Z., Zabihi, M., Larionov, G. A., Krasnov, S. F., Gorobets, A. V., Levi, Y., Erpul, G., Birkel, C., Hoyos, N., Naipal, V., Oliveira, P. T. S., Bonilla, C. A., Meddi, M., Nel, W., Dashti, H. A., Boni, M., Diodato, N., Oost, K. V., Nearing, M., and Ballabio, C.: Global rainfall erosivity assessment based on high-temporal resolution rainfall records, Sci. Rep.-UK, 7, 1–12,, 2017. 

Park, S. W., Mitchell, J. K., and Bubenzer, G. D.: Splash erosion modeling: physical analysis, T. ASAE, 25, 357–361,, 1982. 

Petan, S., Rusjan, S., Vidmar, A., and Mikoš, M.: The rainfall kinetic energy–intensity relationship for rainfall erosivity estimation in the mediterranean part of Slovenia, J. Hydrol., 391, 314–321,, 2010. 

Pickering, B. S., Neely III, R. R., and Harrison, D.: The Disdrometer Verification Network (DiVeN): a UK network of laser precipitation instruments, Atmos. Meas. Tech., 12, 5845–5861,, 2019. 

Prigent, C.: Precipitation retrieval from space: An overview, C. R. Geosci., 342, 380-389,, 2010. 

Renard, K. G., Foster, G. R., Weesies, G., McCool, D., and Yoder, D.: Predicting soil erosion by water: a guide to conservation planning with the Revised Universal Soil Loss Equation (RUSLE), United States Department of Agriculture, Washington, DC, USA, 1997. 

Rosewell, C. J.: Rainfall kinetic energy in eastern Australia, J. Appl. Meteorol. Clim. , 25, 1695–1701,<1695:rkeiea>;2, 1986. 

Sanchez-Moreno, J. F., Mannaerts, C. M., Jetten, V., and Löffler-Mang, M.: Rainfall kinetic energy–intensity and rainfall momentum–intensity relationships for Cape Verde, J. Hydrol., 454, 131–140,, 2012. 

Science and Technology Facilities Council; Chilbolton Facility for Atmospheric and Radio Research; Natural Environment Research Council; and Wrench, C. L.: Chilbolton Facility for Atmospheric and Radio Research (CFARR) Disdrometer Data, Chilbolton Site, NCAS British Atmospheric Data Centre, available at: (last access: 15 November 2020), 2003. 

Sempere-Torres, D., Porrà, J. M., and Creutin, J. D.: Experimental evidence of a general description for raindrop size distribution properties, J. Geophys. Res.-Atmos., 103, 1785–1797,, 1998. 

Song, Y., Han, D., and Rico-Ramirez, M. A.: High temporal resolution rainfall rate estimation from rain gauge measurements, J. Hyydroinfom., 19, 930–941,, 2017. 

Thompson, G. and Eidhammer, T.: A study of aerosol impacts on clouds and precipitation development in a large winter cyclone, J. Atmos. Sci., 71, 3636–3658,, 2014. 

Uplinger, W.: A new formula for raindrop terminal velocity, in: Proceedings of the 20th Conference on Radar Meteorology, American Meteorological Society, 30 November–3 December 1981, Boston, MA, USA, 389–391, 1981. 

Van Dijk, A. I. J. M., Bruijnzeel, L. A., and Rosewell, C. J.: Rainfall intensity–kinetic energy relationships: a critical literature appraisal, J. Hydrol., 261, 1–23,, 2002. 

Wang, L., Shi, Z., Wang, J., Fang, N., Wu, G., and Zhang, H.: Rainfall kinetic energy controlling erosion processes and sediment sorting on steep hillslopes: a case study of clay loam soil from the Loess Plateau, China, J. Hydrol., 512, 168–176,, 2014. 

Weiss, L. L.: Ratio of true to fixed-interval maximum rainfall, J. Hydr. Eng. Div.-ASCE, 90, 77–82, 1964. 

Williams, R. G. and Sheridan, J. M.: Effect of rainfall measurement time and depth resolution on EI calculation, T. ASAE, 34, 402–406,, 1991. 

Wischmeier, W. H. and Smith, D. D.: Rainfall energy and its relationship to soil loss, Transactions American Geophysical Union, 39, 285–291,, 1958. 

Wischmeier, W. H. and Smith, D. D.: Predicting rainfall erosion losses-a guide to conservation planning, Department of Agriculture, Washington, DC, USA, 1978. 

Xie, Y., Yin, S., Liu, B., Nearing, M. A., and Zhao, Y.: Models for estimating daily rainfall erosivity in China, J. Hydrol., 535, 547–558,, 2016. 

Yang, Q., Dai, Q., Han, D., Chen, Y., and Zhang, S.: Sensitivity analysis of raindrop size distribution parameterizations in weather research and forecasting rainfall simulation, Atmos. Res., 228, 1–13,, 2019. 

Short summary
Rainfall is a driving force that accounts for a large proportion of soil loss around the world. Most previous studies used a fixed rainfall–energy relationship to estimate rainfall energy, ignoring the spatial and temporal changes of raindrop microphysical processes. This study proposes a novel method for large-scale and long-term rainfall energy and rainfall erosivity investigations based on rainfall microphysical parameterization schemes in the Weather Research and Forecasting (WRF) model.