Assessment of Rainfall Kinetic-Energy–Intensity Relationships

Abstract

Raindrop-impact-induced erosion starts when detachment of soil particles from the surface results from an expenditure of raindrop energy. Hence, rain kinetic energy is a widely used indicator of the potential ability of rain to detach soil. Although it is widely recognized that knowledge of rain kinetic energy plays a fundamental role in soil erosion studies, its direct evaluation is not straightforward. Commonly, this issue is overcome through indirect estimation using another widely measured hydrological variable, namely, rainfall intensity. However, it has been challenging to establish the best expression to relate kinetic energy to rainfall intensity. In this study, first, kinetic energy values were determined from measurements of an optical disdrometer. Measured kinetic energy values were then used to assess the applicability of the rainfall intensity relationship proposed for central Italy and those used in the major equations employed to estimate the mean annual soil loss, that is, the Universal Soil Loss Equation (USLE) and its two revised versions (RUSLE and RUSLE2). Then, a new theoretical relationship was developed and its performance was compared with equations found in the literature.

1. Introduction

The phenomenon of soil erosion is due to the action of erosive agents that detach and transport individual particles from the soil mass [1]. Soil erosion is a major environmental problem worldwide that heavily affects the agricultural sector. Each year, around 75 billion metric tons of soil are detached from the original soil by wind and water [2], and this amount has an increasing trend. In Europe, the affected area is extremely large. Specifically, the Mediterranean region is more prone to erosion than northern Europe. The reason for this can be found in the long dry periods followed by heavy erosive rainfall events, falling on steep slopes with fragile soils. In contrast, in northern Europe, the rainfall is evenly distributed throughout the year, resulting in less soil erosion [3]. Soil erosion effects can be subdivided into on- and off-site effects [4]. On-site effects refer to the impoverishment of the soil directly affected by erosion. Croplands are the most susceptible terrains to erosion since their soil is repeatedly tilled and left without a protective cover of vegetation. The erosive agents remove the top soil and thus affect croplands, lowering their productivity and quality by reducing the infiltration rates, water-holding capacity, nutrients, organic matter, and soil depth [5]. In contrast, off-site impacts refer to the effects of soil transport where the eroded soil ends its journey. Off-site effects include the silting up of basins and rivers, which leads to a reduction of the capacity of cross sections, thus contributing to trigger flood events. Moreover, the impoverishment of the top soil layer of nutrients due to erosion causes the extensive use of fertilizers, which in turn are exposed to erosive agents. Fertilizers removed from the soil pollute rivers and any other receiving body of water [5]. In this framework, the Universal Soil Loss Equation (USLE), introduced by Wischmeier and Smith [6], was created to support soil conservation planning at the field scale to estimate the soil loss for the U.S. climate. The original formulation was later revised by Renard et al. [7] (i.e., (Revised) Universal Soil Loss Equation, RUSLE) to determine the mean annual soil loss A as the product of six parameters:

$$\mathrm{A}=\mathrm{R}\times \mathrm{K}\times \mathrm{L}\times \mathrm{S}\times \mathrm{C}\times \mathrm{P}$$

(1)

where R is the rainfall erosivity factor; K is the soil erodibility factor; L and S are the topographic factors slope length and slope steepness, respectively; C is the cropping-management factor; and P is the erosion-control practice factor.

The RK product, defined as potential erosion, is a reference factor in soil protection strategies, especially on a regional scale. This is rationalized by considering that the assessment of other factors, especially the topographic factors, is strongly dependent on the spatial scale of the process being studied. Under these considerations, RK information deduced on a regional scale can be a straightforward tool for identifying the area most subject to erosion processes, from the point of view of the rainfall force hazard (R-factor), and the intrinsic characteristics of soil vulnerability (K-factor)—that is, the assessment of potential erosion risk [8].

It is thus evident that the RUSLE is site specific and strongly depends on the characteristics of the cropland and the prevention activities adopted to limit the effect of erosion. Hence, to better understand the phenomenon of erosion and to quantify it, many authors (e.g., [9,10,11,12]) have focused on the triggering factor expressed by R, defined as

$$\mathrm{R}={\displaystyle \sum }\left({\mathrm{KE}}_{\mathrm{ev}}{\mathrm{I}}_{30}\right)/\mathrm{N}$$

(2)

where I₃₀ is the maximum 30 min rainfall intensity for an erosive event, N is the period of observation (i.e., year), and KE_ev is the total kinetic energy (KE) of precipitation for an erosive event [7]. According to Wischmeier and Smith [6], a rainfall event cannot be considered erosive if the corresponding rainfall event depth is lower than 12.7 mm and the minimum interevent time (MIT) is higher than 6 h, unless the rainfall event depth equals 6.35 mm in 15 min. This classification can be questionable, and some authors have provided different criteria (e.g., [1]); however, it represents a standard approach for soil erosion studies around the world [13].

Even though R depends on the total kinetic energy of precipitation, many authors have questioned what the best descriptor is of erosivity between kinetic energy and momentum (e.g., [14]). The former is defined as half the product of the mass and velocity squared of the falling drops, and the latter is the product of mass and velocity [10]. Many authors have investigated these two quantities to understand which one could better describe the soil detachment process. For instance, Rose [15] and Paringit and Nadaoka [16] concluded that momentum provides better results than KE. However, Hudson [17] showed that the momentum and kinetic energy of natural rainfall have a very similar relationship with intensity. Lal and Elliot [18] estimated kinetic energy as the major agent in starting the process of soil erosion. Several scholars proposed splash erosion formulas that include rainfall energy [19,20,21]. Other studies highlighted that the kinetic energy of rainfall alone has a tendency to overestimate the erosivity, as it does not take into account the role of the run-off process [8] For instance, Hudson [22] investigated the issue and introduced a rainfall threshold for rainfall intensity. Other authors, for instance, Wischmeier and Smith [23], addressed the issue by multiplying the storm kinetic energy with a measure of maximum rainfall intensity. Renard et al. [7] combined this method with a threshold value for the storm depth. As highlighted by van Dijk et al. [12], rainfall kinetic energy represents the total energy available to cause the detachment and then transport of soil.

In agreement with the RUSLE formulation, in this work, we analyzed kinetic energy as the major descriptor of soil erosion. Although several experiments have been conducted to assess rainfall erosivity, its direct measurements are spatially limited to a few locations [8]. The limitations of the instrumental series available are mostly attributable to the lack of direct KE measurements and long continuous records of rainfall intensity with finer resolution, which are typically available only for very few locations. Several studies carried out around the world have attempted to address the first issue by calibrating empirical relationships that link rainfall KE to the most commonly available rainfall quantities (i.e., rainfall intensity) (see van Dijk et al. [12] for an extensive review). However, there are sources of uncertainty; for instance, Angulo-Martínez et al. [24] found two sources of bias when using literature relationships. The first one is due to the use of theoretical raindrop terminal velocities instead of measured values, and the second one is due to time aggregation, as rainfall is usually recorded every given amount of time. Interestingly, Assouline [25] found that accounting for the variability of rainfall intensity during a storm has a strong impact on kinetic energy estimation.

Knowledge of the relationship between kinetic energy and rainfall intensity is of utmost importance for the prediction of soil erosion risk. There are many types of well-established relationships, and most of these equations, including the original formulation proposed by Wischmeier and Smith [23], are based on data from specific locations and thus may not perform equally satisfactorily in other locations [1]. The rationale can be found in the different origins and types of rainfall characterizing different locations [10]. The application of an empirical relationship calibrated on a different area may provide a kinetic energy value significantly different from the actual one. Therefore, R would be miscalculated, resulting in under- or overestimation of the soil loss.

On the other hand, a different approach has been introduced to evaluate hydrological and meteorological parameters (i.e., kinetic energy, rainfall rate, etc.) to model rainfall drop spectra (drop size distribution, DSD) retrieved by disdrometers by using theoretical probability distributions. Among others, Adirosi et al. [26], Cerro et al. [27], and Feingold and Levin [28], for the Italian, Spanish, and Israeli areas, respectively, concluded that the gamma distribution, which is considered to be the standard approach for DSD modeling [29], does not provide the best agreement with the observed drop spectra, with respect to other common distribution laws (i.e., Weibull and lognormal).

As soil erosion is a major issue in the Mediterranean area, a disdrometer located in Rome (central Italy) offers a unique opportunity to define a kinetic-energy–rainfall-intensity (KE–I) relationship for that area, which is well known for its cultivations. The use of a disdrometer allowed us to compare empirical values of kinetic energy, estimated from DSD, with those obtained from literature equations.

The present study aimed to (i) investigate the relationship between the properties of rainfall drops and rainfall intensity, (ii) assess the applicability of a literature kinetic-energy–rainfall-intensity relationship for the case study where an optical disdrometer is available, and (iii) establish a new formulation for a kinetic-energy–rainfall-intensity relationship which better represents the physical phenomenon being studied.

The paper is organized as follows. First, the case study is presented and the characteristics of the disdrometer used to collect the dataset are detailed. Then, the procedure to select erosive rainfall events is presented. Third, the contribution of each raindrop in terms of DSD is investigated; the most used formulations to estimate raindrop terminal velocity are analyzed and the most suitable is chosen. Rainfall kinetic energy is then estimated from the DSD and raindrop mass and terminal velocity. Second, the reliability of the formulations used in the USLE [6] and its successive updates (i.e., RUSLE [30] and RUSLE2 [31]) is evaluated. The equation proposed by Zanchi and Torri [32] is also assessed, as it is the only one available for central Italy. Finally, the results are discussed and a new kinetic-energy–rainfall-intensity relationship is proposed for the case study area.

2. Case Study

To estimate the kinetic energy of rainfall, we need to investigate the size distribution of raindrops. This can be obtained by direct measurements carried out with pressure transducers, acoustic devices, or disdrometers [33,34]. Performing disdrometric precipitation measurement means determining a pair of velocity–diameter values for each hydrometeor. The relief of dimensional distributions of raindrops was performed with an optical disdrometer, model OTT Parsivel² (PARticle SIze and VELocity), installed in Rome (central Italy) (Figure 1). The operating principle of the optical disdrometer is based on the measurement of the damping of the light produced by the passage of the drops through a control area of 54 cm² between two diodes.

The dataset covers a period of about five years, ranging from 20 June 2010 to 28 April 2015. Data gauging time resolution equals 1 min. The Parsivel² disdrometer classifies drops automatically into 32 separate size classes, ranging from 0.062 to 24.5 mm. However, in this study, D_max = 7.5 mm was considered as the upper bound for diameter classes, which is consistent with values observed in natural rain [35]. To filter out spurious drops, for each 1 min drop spectrum, the simultaneous presence of 10 drops and a minimum rain rate of 0.01 mmh⁻¹ were selected as thresholds for rain/no-rain events [26]. As a result, 105,433 effective minutes of rain were recorded, with the rainfall intensity varying between 0 and 318.742 mmh⁻¹.

These values were collected for the same time window (i.e., from 20 June 2010 to 28 April 2015). The disdrometer can also record rainfall intensity with a resolution of 0.5 mmh⁻¹. In this regard, as demonstrated in a deep investigation conducted by Adirosi et al. [36], rainfall intensity values collected by a disdrometer can be affected by bias. In this study, this issue was accounted for by filtering the rainfall intensity measured by the disdrometer with the intensity estimated from rainfall drops, as explained below.

2.1. Selection of Potential Erosive Events

Since this methodology aims at defining the relationship between rainfall intensity and kinetic energy which results in soil erosion, we selected erosive rainfall events out of those gauged by the disdrometer according to the classification proposed by Wischmeier and Smith [6]. This method allowed us to exclude events that were not representative for the modeling of soil erosion process. Moreover, raindrops larger than 7.5 mm were excluded from the analysis, since such large measurements can be due either to hail or two or more drops falling at the same moment.

These conditions, applied to our measured drop spectra, led us to identify 88 different events as potentially erosive. The season with the highest number of events (i.e., 28) was spring, while a lower number of events (i.e., 17) occurred in the summer (Figure 2). Based on the annual distribution, the highest number of events occurred in 2014 (Figure 2). In Appendix A, the main characteristics of the erosive events are reported by year of observation.

2.2. Preliminary Analysis of the Dataset

Most erosion control practices are not designed or expected to be proof against extreme events, because in most cases, failure does not cause devastating damages and the practices can be reinstalled without great costs. Thus, in the latest version of the soil loss equation (i.e., RUSLE2 [37]), erosive events with a return period higher than 50 years were removed from R-factor calculations. Therefore, a further investigation was carried out in order to classify the erosive events collected by the disdrometer on the basis of the return period. To this end, we used rainfall depth series recorded by the Frascati rain gauge (18 years of observation), which is the nearest recording station. As is well known, short record lengths can affect the estimation accuracy of high return periods [38]. However, this issue did not play a key role in this investigation. As described below, the relationship between return period and rainfall depth evaluated for the Frascati rain gauge was used only as a threshold, in agreement with Mineo et al. [39].

Rainfall depth precipitation data were collected from 1994 to 2017, which had a 5 min interval time. Then, by aggregation, cumulated rainfall values were obtained for short durations (i.e., 10, 15, 30, 45, and 60 min) and long durations (i.e., 3 and 6 h). For these durations, annual maximum rainfall depths were computed. For the selected station, it was assumed that the rainfall depth values of a specific duration d and return period Tr (i.e., h_d,Tr) follow the extreme value type 2 probability distribution law (EV2; [39,40,41]). The three-parameter, site-specific depth–duration–frequency (DDF) curve was adopted [42]:

$${\mathrm{h}}_{\mathrm{d},\mathrm{Tr}}=\frac{\mathrm{a}\left(\mathrm{Tr}\right)\times \mathrm{d}}{{\left(\mathrm{b}+\mathrm{d}\right)}^{\mathrm{m}}}$$

(3)

where a and b are positive site-specific parameters, m is an exponent ranging from 0 to 1, and d is the rainfall duration. The parameters were calculated by Mineo and Napolitano [43] by solving Equation (3) numerically. From the DDF curves estimated for the Frascati rain gauge (Figure 3), it was quite evident that potentially erosive events showed depth precipitation amounts lower than the h_d,Tr=50yr curve (i.e., no erosive events had a return period higher than 50 years). Therefore, no event was discarded in the subsequent analyses.

In addition, in order to carry out surveys on statistically significant samples [36,40], the following filtering criteria were applied to the measured drop spectra:

• the number of registered drops must be greater than 100;
• at least five raindrops are requested per diameter class;
• rainfall intensity must be greater than 0.5 mmh⁻¹.

Such filtering steps led to a dataset consisting of 28,287 1 min drop spectra.

3. Kinetic Energy Definition

Kinetic energy is the most widely used descriptor of soil erosion. It is estimated evaluating the contribution of all raindrops once their mass, terminal velocity, and DSD are known. The DSD refers to the number of drops present in a unit volume of air (i.e., N(D) (m⁻³ mm⁻¹)), with a diameter D (mm) ranging between D and D + dD that reaches a unit horizontal area during a unit time Δt (min) [29,44,45,46]. The empirical DSD, for each diameter class k, can be estimated using the following expression [26]:

$$\mathrm{N}{\left(\mathrm{D}\right)}_{\mathrm{E}\left(\mathrm{k}\right)}=\frac{1}{\mathrm{A}\mathsf{\Delta }\mathrm{t}}\frac{{\mathrm{n}}_{\mathrm{t}\left(\mathrm{k}\right)}}{{\mathsf{\Delta }\mathrm{D}}_{\mathrm{k}}{\mathrm{V}}_{\mathrm{t}}{\left(\mathrm{D}\right)}_{\mathrm{k}}}$$

(4)

where ΔD_k (mm) is the width of each diameter class, A (m²) is the gauging area of the disdrometer when the ith drop is falling, n_t(k) is the number of raindrops falling in the kth class, and V_t(D)_k (ms⁻¹) is the raindrop terminal velocity for each kth diameter class.

Kinetic energy is a function of both the mass of raindrops and their terminal velocity. The mass M_E (kg) is estimated from the diameters obtained from the DSD and is defined as

$${\mathrm{M}}_{\mathrm{E}}=\mathsf{\rho }\frac{\mathsf{\pi }}{6}{\displaystyle \sum }_{\mathrm{k}}\mathrm{N}{\left(\mathrm{D}\right)}_{\mathrm{E}\left(\mathrm{k}\right)}{\mathrm{D}}_{\mathrm{k}}^3{\mathrm{V}}_{\mathrm{t}}{\left(\mathrm{D}\right)}_{\mathrm{k}}$$

(5)

where $\mathsf{\rho }$ is the water density (1000 kg m⁻³) under standard conditions. Therefore, once the DSD is known, the mass is estimated from Equation (5). Then, the kinetic energy is half the product of the mass and the terminal velocity of drops at the power of two.

Two expressions of the kinetic energy have been discussed in the literature [9]: the rain kinetic energy rate, or time-specific kinetic energy (i.e., KE_time, expressed in J m⁻² h⁻¹), and the rain kinetic energy content, or energy per unit volume (i.e., KE_mm, expressed in J m⁻² mm⁻¹).

The empirical time-specific kinetic energy is defined as [9]

$${\mathrm{KE}}_{\mathrm{time}\left(\mathrm{E}\right)}=3.6\mathsf{\rho }\frac{\mathsf{\pi }}{12}{10}^{-6}{\displaystyle \sum }_{\mathrm{k}}\mathrm{N}{\left(\mathrm{D}\right)}_{\mathrm{E}\left(\mathrm{k}\right)}{\mathrm{D}}_{\mathrm{k}}^3{\mathrm{V}}_{\mathrm{t}}{\left(\mathrm{D}\right)}_{\mathrm{k}}.$$

(6)

The empirical rain kinetic energy content KE_mm(E) is directly proportional to KE_time(E) by the rain rate I (mmh⁻¹), as shown below:

$${\mathrm{KE}}_{\mathrm{time}\left(\mathrm{E}\right)}=\mathrm{c}\times \mathrm{I}\times {\mathrm{KE}}_{\mathrm{mm}\left(\mathrm{E}\right)}$$

(7)

where c is a constant depending on the unit used. For instance, c equals 1 if KE_time(E) is expressed in J m⁻² h⁻¹ and KE_mm(E) is in J m⁻² mm⁻¹. In Equation (7), empirically evaluating I (i.e., I_E, expressed in mm h⁻¹) leads to

$${\mathrm{I}}_{\mathrm{E}}=3.6\mathsf{\rho }\frac{\mathsf{\pi }}{6}{10}^{-3}{\displaystyle \sum }_{\mathrm{k}}\mathrm{N}{\left(\mathrm{D}\right)}_{\mathrm{E}\left(\mathrm{k}\right)}{\mathrm{D}}_{\mathrm{k}}^3{\mathrm{V}}_{\mathrm{t}}{\left(\mathrm{D}\right)}_{\mathrm{k}}.$$

(8)

It is thus evident that the most important descriptors of the kinetic energy rate are the diameter and the terminal velocity of the drops. The terminal velocity V_t (ms⁻¹) is a function of the drop diameter D (mm). Ulbrich [47] proposed a power law, here denoted as V_t(D)_U:

$${\mathrm{V}}_{\mathrm{t}}{\left(\mathrm{D}\right)}_{\mathrm{U}}=3.78{\mathrm{D}}^{0.67}.$$

(9)

V_t(D)_U is thus a monotonically increasing function of the drop diameter. Atlas et al. [48] proposed a formulation interpolating measurements performed by Gunn and Kinzer [49], who estimated the terminal velocity of thousands of artificial droplets with diameter D ranging from 0.2 to 5.8 mm. Atlas et al.’s expression V_t(D)_A (ms⁻¹) is as follows:

$${\mathrm{V}}_{\mathrm{t}}{\left(\mathrm{D}\right)}_{\mathrm{A}}=9.65-10.3\times \exp \left(-0.6\mathrm{D}\right).$$

(10)

Equation (10) provides a good fit for diameters larger than 0.4 mm but predicts negative terminal velocity values for very small drops. Finally, Ferro [50] obtained the following relationship considering that terminal velocity is an increasing function (or almost constant for D > 0.55–0.57 cm) of diameter D:

$${\mathrm{V}}_{\mathrm{t}}{\left(\mathrm{D}\right)}_{\mathrm{F}}={\mathrm{V}}_{\mathrm{n}}\times \left[1-\exp \left(-{\mathrm{a}}_{\mathrm{n}}\mathrm{D}\right)\right]$$

(11)

where V_n and a_n are constants, 9.5 (ms⁻¹) and 6 (cm⁻¹), respectively. Equation (11) was based on the measurements of several scholars [35,49,51,52,53,54]. V_t(D)_F equals V_t(D)_A for D < 0.55 cm and it is more reliable than V_t(D)_U for D > 0.55 cm. According to Ulbrich’s formulation [47], V_t(D)_U is monotonically increasing, while V_t(D)_A and V_t(D)_F tend asymptotically to a constant value (Figure 4).

It is worth noting that, according to the formulations proposed by both [48,50], V_t strongly depends on the diameter for D < 4 mm, while the dependency attenuates for diameters ranging from 4 to 6 mm. For D > 6 mm, the dependency of V_t on D is negligible. In this work, we estimated the terminal velocity according Ferro’s formulation, since Atlas et al.’s relationship inconveniently provides negative values of V_t for small droplet diameters.

KE–I Empirical Relationships

As the formulation of KE is estimated from the DSD, its direct measurement is not always a straightforward task. Therefore, the following approach was used to deal with this issue, that is, the indirect estimation of kinetic energy from another widely measured hydrological variable (i.e., rainfall intensity, I). The advantage of this approach is obvious, given that rain intensity has been measured by rain gauges all over the world [55]. The empirical expressions proposed by different authors reflect different pluviometric calibration regimes through different mathematical formulations that link kinetic energy to precipitation intensity. Henceforth, for the sake of brevity, the relationship between KE_mm and I is referred to as KE–I. Several KE–I mathematical relationships proposed in the literature were inspired by the following logarithmic model that implies no upper limit to kinetic energy [6,32]:

$${\mathrm{KE}}_{\mathrm{mm}}=\mathsf{\alpha }+\mathsf{\beta }\times {\log }_{10}\left(\mathrm{I}\right)$$

(12)

where $\mathsf{\alpha }$ and $\mathsf{\beta }$ are constants derived through regression. In a second version of USLE, Wischmeier and Smith [23] showed that kinetic energy does not increase with rainfall intensity indefinitely and proposed a conditional expression (henceforth called WS) on the basis of Equation (12):

$${\mathrm{KE}}_{\mathrm{mm}\left(\mathrm{WS}\right)}=\{\begin{array}{c}11.9+8.73{\log }_{10}\left(\mathrm{I}\right), \mathrm{I}\le 76 \mathrm{mm} {\mathrm{h}}^{-1}\\ 28.3,\hspace{1em} \mathrm{I}>76 \mathrm{mm} {\mathrm{h}}^{-1}\end{array}.$$

(13)

It implies that KE remains constant at 28.3 J m⁻² mm⁻¹ for rainfall intensity exceeding 76 mmh⁻¹. It was found for Washington (United States) and the DSDs were estimated according to [56].

On the other hand, several authors have proposed an exponential model [11,30,57,58,59,60,61]:

$${\mathrm{KE}}_{\mathrm{mm}}=\mathsf{\gamma }\left[1-\mathrm{f}\times \exp \left(-\mathrm{zI}\right)\right]$$

(14)

where $\mathsf{\gamma }$ (J m⁻² mm⁻¹), f, and z (h mm⁻¹) are regression parameters. According to Equation (13), Equation (14) leads to an asymptotic value (i.e., KE = γ) at high intensities. Kinnell [59] suggested that the exponential model performs better than the logarithmic one. Having one extra parameter, the exponential model is more flexible to fit the datasets. The developers of RUSLE [7] used Equation (14), with the parameters calibrated by Brown and Foster [30], as a standard method for rainfall energy calculation. Later, in RUSLE2 [37], the exponential model was confirmed, although the method adopts the equation proposed by McGregor et al. [31].

Figure 5 shows some of the expressions found in the literature calibrated for each corresponding pluviometric regime. It is possible to appreciate that the curves tend to overlap for low precipitation intensity values, while the difference between the predicted KE values becomes significant as the intensity increases.

In the 1980s, Zanchi and Torri [32] proposed a relationship (henceforth called ZT) calibrated with measurements performed in Florence (central Italy):

$${\mathrm{KE}}_{\mathrm{mm}\left(\mathrm{ZT}\right)}=9.81+11.25\times {\log }_{10}\left(\mathrm{I}\right).$$

(15)

McGregor et al. [31] proposed an equation (henceforth called MG) based on the model suggested by Kinnell [59], which results in

$${\mathrm{KE}}_{\mathrm{mm}\left(\mathrm{MG}\right)}=29\left[1-0.72\times \exp \left(-0.082\mathrm{I}\right)\right].$$

(16)

On the other hand, Brown and Foster [30] introduced an exponential relationship (henceforth called BF) for the United States as follows:

$${\mathrm{KE}}_{\mathrm{mm}\left(\mathrm{BF}\right)}=29\left[1-0.72\times \exp \left(-0.05\mathrm{I}\right)\right]\mathrm{for}\mathrm{I}=0\mathrm{up}\mathrm{to}250{\mathrm{mmh}}^{-1}$$

(17)

4. Results

As we used DSDs measured by the disdrometer to estimate the kinetic energy expenditure (i.e., KE_time(E)), the hydrological parameters of interest for the 1 min drop spectra were calculated by summation from the minimum to the maximum diameter collected [26,29]. The empirical DSDs were characterized by rainfall intensity, I_E (mmh⁻¹), mean diameter, m(D)_E (cm), standard deviation s(D)_E (cm), median value of volumetric size spectra D_0(E) (cm), and kinetic energy expenditure KE_time(E) (J m⁻² h⁻¹). The DSDs identified from the analysis of the 88 selected erosive events were characterized by a mean diameter m(D)_E ranging from 0.046 to 0.55 cm and a standard deviation s(D)_E from 0.0087 to 0.40 cm. According to what has been found by other authors (e.g., [44]), both the m(D)_E and s(D)_E of empirical DSDs weakly increase with intensity (Figure 6). In Figure 6a, it is possible to appreciate a range of values of m(D)_E that differs significantly from the others. This behavior is even more evident in Figure 6c, where an area that comprises unusually high values of KE_time(E) (even at low intensities) is clearly depicted.

The rationale for this behavior could potentially be due to the horizontal component of wind, as it is one of the main environmental factors that influences disdrometer measurements [1,62]. In fact, the structure of the disdrometer modifies the airflow around the instrument itself and thus affects the trajectories of the falling drops, which can miss the measurement area of the device and not be detected. Friedrich et al. [63] studied the effect of wind on Parsivel² measurements collected during a hurricane and found that artifacts caused by particles moving at an angle through the sampling area, most of the time, occurred for very strong wind speeds (i.e., more than 20 ms⁻¹), although they were also observed when the wind speed was lower than 10 ms⁻¹. The main effect of the wind is on the number of drops recorded. Another source of error that must be taken into account is splashing, that is, when one drop, hitting the instrument housing, breaks in two drops, which can then enter the measuring area with unrealistically low fall velocities. Wind and splashing effects increase in heavy rain conditions.

Environmental factors (i.e., wind effect and splashing) can therefore cause significant errors in the number of drops recorded by the device, which may not correspond to the intensity value measured by the device (i.e., I_M). In fact, the instrument calculates I_M regardless of the number of drops recorded, and it can therefore be assumed as a benchmark. Under this assumption, here, we dealt with the environmental factor effects by introducing a further filter on the 1 min drop spectra. Thus, the reliability of the number of drops recorded for each 1 min drop spectrum was assessed by comparing empirical rainfall intensities (I_E), calculated using Equation (8), and those registered by the disdrometer (I_M). A 1 min drop spectrum was omitted from the analysis if the following condition applied:

$$\left|{\mathrm{I}}_{\mathrm{E}}-{\mathrm{I}}_{\mathrm{M}}\right|\ge \mathsf{\Delta }\mathrm{I}$$

(18)

where ∆I is a set equal to 0.5 mmh⁻¹, consistent with device accuracy. This constraint led us to reduce the dataset from 28,284 to 24,330 1 min drop spectra, characterized by

• rainfall intensity, I_E, from 0.5 to 230.3 mmh⁻¹;
• mean diameter, m(D)_E, from 0.046 to 0.26 cm (Figure 7a);
• standard deviation, s(D)_E, from 0.0087 to 0.18 cm (Figure 7b);
• median value of volumetric size spectra, ${\mathrm{D}}_{0\left(\mathrm{E}\right)}$, from 0.042 to 0.49 cm;
• kinetic energy expenditure, KE_time(E), from 1.92 to 7463.8 J m⁻² h⁻¹ (Figure 7c).

The elimination of measurements distorted by potential environmental factor effects is more evident when observing the KE_time(E) plotted against the rainfall intensity (Figure 7c). As mentioned, disdrometrical measurements can be affected by different biases; however, it could indeed be that the drop velocity was drastically increased by the wind to cause the scattered KE_time(E) values.

The 24,330 rainfall minutes unfolded in a range of rainfall intensities between 0.5 and 231 mmh⁻¹. However, Figure 8 clearly shows that a high percentile (i.e., 99%) of the drop spectra occurred with a rainfall intensity lower than 15 mmh⁻¹.

Selecting the 1 min drop spectra by grouping per rainfall rate classes (width of 2 mmh⁻¹), we can observe that (i) the total number of drops falling each minute (i.e., n_t) increased with the rainfall intensity, (ii) the range of variability of n_t gradually increased with the intensity class, and (iii) the class corresponding to the mean intensity (i.e., 39.48 mmh⁻¹) contained 1 min drop spectra with the highest variability of n_t (Figure 9). This suggests that the lower intensity classes consisted of rainfall minutes with more homogeneous characteristics.

To investigate the link between the number of drops and rainfall intensity more deeply, in Figure 10, the box plots for j different intensity classes (with j = 8) are depicted. In each plot, the number of drops falling in eight out of the n_t(k) diametrical classes are shown for each 1 min drop spectrum. Figure 10 suggests the following considerations: (i) the median value of n_t(k)(j) assumes greater values as the jth intensity class increases; (ii) for intensity classes lower than 14 mmh⁻¹, the maximum median value n_t(k)(j)_max shifts towards greater diametrical classes as the jth class intensity increases; (iii) for intensity classes greater than 14 mmh⁻¹, the n_t(k)(j)_max shifts towards smaller diametrical classes as the jth intensity class increases.

The results suggest that, for low intensity classes, there was a prevalence of droplets (i.e., drops with a diameter of <1 mm). As the intensity increased, the number of large diameter drops increased, and the coalescence phenomenon was prevalent, while the percentage of droplets decreased progressively and then grew again at a specific threshold value. This trend can be explained by the fact that the larger drops, which become unstable, start to break into smaller drops (i.e., disaggregation), determining quasi-constant values of the empirical median volume diameter (Figure 11) and kinetic energy content (Figure 12) as well [11,60].

Regarding the investigated drop spectra, the maximum values of D_0(E) (Figure 11) and KE_mm(E) (Figure 12) were achieved at an intensity threshold of 32 mmh⁻¹. This result confirms the approach proposed by Wischmeier and Smith [23], which found a limit value of 76 mmh⁻¹, calibrated on the U.S. pluviometric regime. Different limit values were assessed for different case studies by other authors, for instance, Rosewell [60] found a value of 100 mmh⁻¹ in Australia, while Carollo et al. [44] found that this value is equal to 40 mmh⁻¹ for the Sicilian pluviometric regime.

5. Discussion: A New Analytical Relationship

Since this paper is focused on processes on the earth’s surface, it is of primary importance that proper estimates are obtained for KE values at higher rainfall intensities.

To determine the most suitable KE–I relationship for the case study under analysis, the relationship existing between kinetic energy and empirical drop size distributions (i.e., number of drops vs. rainfall intensity) was investigated. Using Equation (7) and the expression proposed by Ferro [50] to estimate terminal velocity (i.e., Equation (11)), we estimated KE_mm(E) from disdrometer measurements. KE_mm(E) was estimated and then plotted against rainfall intensity.

Figure 13 shows the comparison between pairs of KE_mm(E) mean values and I_M, and pairs of KE–I predicted by WS, ZT, MG, and BF relationships. We can clearly see the different trend of the ZT equation (green line), as it does not consider a threshold value for intensity.

Finally, recalibrating the rainfall energy equation proposed by Kinnel [59], a relationship relating KE_mm and I here was derived through regression analysis:

$${\mathrm{KE}}_{\mathrm{mm}\left(\mathrm{NF}\right)}=31.75\left[1-0.64·\exp \left(-0.05\mathrm{I}\right)\right]\mathrm{for}\mathrm{I}=0.5\mathrm{up}\mathrm{to}232{\mathrm{mmh}}^{-1}.$$

(19)

The equation better represents the KE–I trend over the area of the case study (Figure 14).

To estimate the goodness of fit of the analytical relations, the mean absolute percentage error (MAPE, %), the root-mean-square error (RMSE, J m⁻² mm⁻¹), and the normalized RMSE (NRMSE, dimensionless) were evaluated:

$$\mathrm{MAPE}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{N}}\left|\frac{{\mathrm{KE}}_{\left(\mathrm{E}\right),\mathrm{i}}-{\mathrm{KE}}_{\left(\mathrm{NF}\right),\mathrm{i}}}{{\mathrm{KE}}_{\left(\mathrm{E}\right),\mathrm{i}}}\right|\frac{100\%}{\mathrm{N}}$$

(20)

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{KE}}_{\left(\mathrm{E}\right),\mathrm{i}}-{\mathrm{KE}}_{\left(\mathrm{NF}\right),\mathrm{i}}\right)}^2}{\mathrm{N}}}$$

(21)

$$\mathrm{NRMSE}=\frac{\mathrm{RMSE}}{\overline{{\mathrm{KE}}_{\left(\mathrm{E}\right)}}}$$

(22)

where KE_(model),i is represented by the ith values of rainfall kinetic energy predicted by both the literature empirical formulations (WS, ZT, BF, and MG) and by Equation (19) (hereinafter KE_mm(NF) or new formulation). In Figure 15, the MAPE (Equation (20)) between the KE_mm(E) and the ones predicted by empirical formulations (WS, ZT, BF, and MG) are represented, depending on rainfall intensity.

In terms of lowest MAPE, among the expressions in the literature, the best performance was provided by MG, followed by ZT. The greatest percentage values of MAPE were obtained by comparing the values of KE_mm(E) with the expression from BF. Certainly, the best performance was provided by the relationship specifically recalibrated for the case study (Figure 16 and

Table 1. Percentile values of MAPE in terms of maximum and minimum values, root-mean-square error (RMSE), and normalized RMSE (NRMSE) between KE_mm(E) and_. KE_mm estimated according to the empirical relations analyzed here (WS, ZT, BF, and MG) and the formulation proposed (NF).

Metric	Unit	WS	ZT	BF	MG	NF
MAPE25%	%	6.00	4.63	9.90	9.90	2.13
MAPE50%	%	12.05	9.07	16.25	8.34	5.03
MAPE75%	%	18.44	14.49	23.28	14.02	11.08
MAPE95%	%	30.69	24.93	31.18	22.91	18.53
max KEmm	J m−2 mm−1	44.86	47.50	45.54	35.82	32.89
min KEmm	J m−2 mm−1	0.40	0.22	0.49	0.46	0.22
RMSE	J m−2 mm−1	21.33	20.14	23.82	19.58	17.48
NRMSE	-	0.44	0.41	0.49	0.40	0.36

). Therefore, even if a reasonable dispersion between KE_mm(E) and KE_mm(NF) is appreciable, Equation (19) offers an improvement compared with the other literature expressions. Therefore, for the case study, the results suggest that the new formulation is preferable to calculate the kinetic energy of precipitation.

However, it is important to highlight that, among the other formulations, the best performance was provided by MG´s expression. This is extremely interesting, since this relationship is used to estimate the RUSLE2, which is the most updated expression to estimate soil erosion and its use is widespread. The fact that MG´s expression could successfully describe the pluviometric regime under study highlights its robustness and its wide range of applicability. Once its suitability to a case study area is verified, it could be used in ungauged areas or in case of missing data. However, if data are available, the new relationship proposed here can be an effective replacement of MG´s expression.

6. Conclusions

In this paper, first, we estimated kinetic energy values from measurements of an optical disdrometer installed in Rome. Then, we assessed the performance of some literature relationships which have had a major role in the estimation of the R-factor in the RUSLE equation and its further updates: the relationships proposed by Wischmeier and Smith [6], used in USLE; Brown and Foster [30] (RUSLE); and McGregor et al. [31] (RUSLE2). The relationship proposed by Zanchi and Torri [32] was also assessed, as it was estimated from data collected in central Italy.

From analyzing the kinetic energy against the rainfall intensity obtained from data collected by the disdrometer, it was evident that a large portion of data was scattered and did not cluster with the rest of the data. While in usual practice these data are kept in the analysis, here, these measurements, likely distorted by the side wind effect, were eliminated to retrieve a proper KE–I relationship. From the measurements, for large values of rainfall intensity, the empirical kinetic energy content of rainfall had quasi-constant values. This was due to the fact that the larger drops, which became unstable, broke into smaller drops (i.e., disaggregation), determining quasi-constant values of the empirical median volume diameter. The maximum values of the energy of raindrops were quasi-constant at intensities exceeding 32 mmh⁻¹. This finding confirms what was found by many authors who established a limit value that varies with the pluviometric regime and, thus, with the case study.

From the analysis, it was thus clear that the KE–I relationship was strongly influenced by the pluviometric regime under study. Therefore, it was necessary to assess if KE–I literature expressions are suitable for the case study under analysis (i.e., Rome). Among literature expressions, the best performance was provided by McGregor et al.’s [31] expression, followed by Zanchi and Torri [32]. This has an important implication, since McGregor et al.’s [31] relationship is used in the framework of the RUSLE2, which is the most updated expression to estimate soil erosion and is used worldwide. Instead, the greatest percentile values of error were found by comparing the empirical values of KE_mm with the expression from Brown and Foster [30]. Finally, an exponential equation relating kinetic energy and rainfall intensity was recalibrated in order to better represent the KE–I trend over the case study. The application of a local kinetic energy expression resulted in, on average, the best performance in estimating the kinetic energy values, as compared with the equations proposed by the RUSLE manual, which were not developed specifically for the Lazio region. Even if a reasonable dispersion between empirical values of kinetic energy and those obtained with the new relationship is appreciable, the new relationship offers an improvement in the estimation of KE_mm. Since the results led us to verify that the McGregor et al. [31] expression is suitable to describe the pluviometric regime under study, in case of an ungauged area or missing data, it can also be used to estimate the kinetic energy for the case study area. However, if data are available, the new relationship proposed in this paper can be an effective replacement of McGregor et al.’s [31] expression for the case study area.