Computation of Rainfall Erosivity and Surface Runoff

*Soil Erosion - Rainfall Erosivity and Risk Assessment*

Congress 2008. Honolulu, Hawaii:

Ahupua'A; 2008. pp. 1-10

[62] Kaffas K, Hrissanthou V. Computation of hourly sediment discharges and annual sediment yields by means of two soil erosion models in a mountainous basin. International Journal of River Basin Management. 2019;**17**(1):63-77. DOI: 10.1080/15715124.2017.1402777

[63] Mati BM, Morgan RPC,

esp.1347

Quinton JN. Soil erosion modelling with EUROSEM at Embori and Mukogodo catchments, Kenya. Earth Surface Processes and Landforms. 2006;**31**(5):579-588. DOI: 10.1002/

**12**

Chapter 2

Abstract

1. Introduction

activities [1].

15

and Chris Evangelides

Rainfall Erosivity and Its

Estimation: Conventional and

Rainfall erosivity concerns the ability of rainfall to cause erosion on the surface of the earth. The difficulty in modeling the distribution, the size, and the terminal velocity of raindrops in relation to the detachment of soil particles led to the use of more tractable rainfall indices. Thus, in the universal soil loss equation (USLE), the coefficient of rainfall erosivity, R, was introduced. This coefficient is based on the product of the rainfall kinetic energy of a storm and its maximum 30-minute intensity. An important problem in the application of USLE and its revisions in various parts of the world concerns the computation of R, which requires pluviograph records with a length of at least 20 years. For this reason, empirical equations have been developed that are based on coarser rainfall data, such as daily, monthly, or yearly, which are available on larger spatial and temporal extents. However, the lack of denser data is dealt more effectively by means of machine learning methods. Computational systems for this purpose were recently developed based on feed-forward neural networks, yielding significantly better results.

Soil erosion is the detachment, transport, and deposition of soil particles, and it takes place in the course of one or more processes. These processes are natural, such as rainfall with the energy that it carries, surface water runoff, wind, and gravity. The observed processes may be combined with human activities and works, such as intensive cultivations, deforestation, construction of public works, and mining

Earth terrain is strongly influenced by erosion, which progresses in geological time, in case it is a purely natural process. Otherwise, when human activities and works are involved, the phenomenon of erosion is accelerated. This accelerated erosion may cause uncontrollable soil loss with negative consequences for human

The universal soil loss equation (USLE) with its revisions (RUSLE, RUSLE2) can

Machine Learning Methods

Konstantinos Vantas, Epaminondas Sidiropoulos

Keywords: rainfall, erosivity, machine learning, erosivity density, universal soil loss equation, nonlinear regression, neural networks

health, the natural ecosystems, the climate, as well as the economy.

be used to predict the average rate of soil erosion by grouping the numerous parameters that affect soil loss into a set of factors and is the prediction equation

## Chapter 2

## Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods

Konstantinos Vantas, Epaminondas Sidiropoulos and Chris Evangelides

## Abstract

Rainfall erosivity concerns the ability of rainfall to cause erosion on the surface of the earth. The difficulty in modeling the distribution, the size, and the terminal velocity of raindrops in relation to the detachment of soil particles led to the use of more tractable rainfall indices. Thus, in the universal soil loss equation (USLE), the coefficient of rainfall erosivity, R, was introduced. This coefficient is based on the product of the rainfall kinetic energy of a storm and its maximum 30-minute intensity. An important problem in the application of USLE and its revisions in various parts of the world concerns the computation of R, which requires pluviograph records with a length of at least 20 years. For this reason, empirical equations have been developed that are based on coarser rainfall data, such as daily, monthly, or yearly, which are available on larger spatial and temporal extents. However, the lack of denser data is dealt more effectively by means of machine learning methods. Computational systems for this purpose were recently developed based on feed-forward neural networks, yielding significantly better results.

Keywords: rainfall, erosivity, machine learning, erosivity density, universal soil loss equation, nonlinear regression, neural networks

## 1. Introduction

Soil erosion is the detachment, transport, and deposition of soil particles, and it takes place in the course of one or more processes. These processes are natural, such as rainfall with the energy that it carries, surface water runoff, wind, and gravity. The observed processes may be combined with human activities and works, such as intensive cultivations, deforestation, construction of public works, and mining activities [1].

Earth terrain is strongly influenced by erosion, which progresses in geological time, in case it is a purely natural process. Otherwise, when human activities and works are involved, the phenomenon of erosion is accelerated. This accelerated erosion may cause uncontrollable soil loss with negative consequences for human health, the natural ecosystems, the climate, as well as the economy.

The universal soil loss equation (USLE) with its revisions (RUSLE, RUSLE2) can be used to predict the average rate of soil erosion by grouping the numerous parameters that affect soil loss into a set of factors and is the prediction equation

most widely used in the world [2]. According to Nearing et al. [3], "Soil loss refers to the amount of sediment that reaches the end of a specified area on a hillslope that is experiencing net loss of soil by water erosion." Development of soil loss equations began in 1940, and the universal soil loss equation was developed at the National Runoff and Soil Loss Data Center established in 1954 by the Science and Education Administration [4].

The USLE represents an erosion model developed for the prediction of soil losses in an average long-term sense. It is based on the knowledge of the physical characteristics of the field area under study, along with the prevailing cropping and management system. USLE has been widely tested in field conditions, and therefore its validity has been established.

USLE consists of the product of six factors, whose numerical values can be specified depending on a particular location. There is a considerable variation in the resulting erosion values, if observations are limited within short periods. However, long-term averages correspond more satisfactorily to predictions.

The soil loss equation is

$$A = R \cdot K \cdot L \cdot \mathbf{S} \cdot \mathbf{C} \cdot P \tag{1}$$

2. Computation of the R factor

DOI: http://dx.doi.org/10.5772/intechopen.85937

during year j, and EI<sup>30</sup> (MJ mm ha�<sup>1</sup> h�<sup>1</sup>

Wischmeier and Smith [4] was used:

The erosivity of a particular event is

the rainfall event:

a 30-minute duration.

intensity.

Figure 1.

et al. [61].

17

The R factor is computed for time periods greater than 20 years, so that wet or dry rainfall periods will be incorporated, eliminating any bias. It is set equal to the average of the sum of the erosivity values for every year's rainfalls. The R coefficient is defined as the product of the kinetic energy of an erosive rainfall event by the maximum intensity of a 30-minute duration rainfall, during

k¼1

where n is the number of years in the record, mj the number of erosive events

m r¼1

er � vr 

ð Þ EI<sup>30</sup> <sup>k</sup> (2)

� I<sup>30</sup> (3)

); vr is the

, where i is rainfall

) the rainfall erosivity for event k.

er ¼ 0:119 þ 0:0873log10ð Þi (4)

<sup>R</sup> <sup>¼</sup> <sup>1</sup> n ∑ n j¼1 ∑ mj

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods

EI<sup>30</sup> ¼ ∑

where er is the kinetic energy per unit of rainfall (MJ ha�<sup>1</sup> mm�<sup>1</sup>

with the upper limit of 0.283 MJ ha�<sup>1</sup> mm�<sup>1</sup> if i>76 mm h�<sup>1</sup>

rainfall depth (mm) for the time interval r of the hyetograph, which has been divided into r ¼ 1, 2, …, m subintervals; and I<sup>30</sup> is the maximum rainfall intensity for

For the computation of er, numerous empirical relations (Figure 1) involving rain intensity have been proposed [6]. Thus, in USLE, the empirical relation of

The three different kinetic energy equations used in USLE and its revisions. The points are actual data, showing the relation between intensity and the kinetic energy of the rainfall, converted to SI units and coming from Haan

where the meaning of the symbols is given in [5] exactly as follows:


Out of the six factors, only two have units, the erosivity R and the erodibility K. The remaining four, i.e., slope steepness, slope length, cropping management, and components of control practice, are dimensionless factors, because they represent ratios in relation to the unit plot [4].

Regarding R, it needs to be observed that its units stem from its definition as the energy multiplied by maximum 30-minute intensity term. In contrast, the erodibility factor values were determined empirically by calibration against measured erosion data. This is important due to the fact that if there is a change on the definition of the R factor, then the values of the K factor should be recalculated [3].

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods DOI: http://dx.doi.org/10.5772/intechopen.85937

#### 2. Computation of the R factor

most widely used in the world [2]. According to Nearing et al. [3], "Soil loss refers to the amount of sediment that reaches the end of a specified area on a hillslope that is experiencing net loss of soil by water erosion." Development of soil loss equations began in 1940, and the universal soil loss equation was developed at the National Runoff and Soil Loss Data Center established in 1954 by the Science and Education

The USLE represents an erosion model developed for the prediction of soil losses in an average long-term sense. It is based on the knowledge of the physical characteristics of the field area under study, along with the prevailing cropping and management system. USLE has been widely tested in field conditions, and therefore

USLE consists of the product of six factors, whose numerical values can be specified depending on a particular location. There is a considerable variation in the resulting erosion values, if observations are limited within short periods. However,

A ¼ R � K � L � S � C � P (1)

long-term averages correspond more satisfactorily to predictions.

where the meaning of the symbols is given in [5] exactly as follows:

• A is the computed soil loss per unit area, expressed in the units selected for K and for

• R, the rainfall and runoff factor, is the number of rainfall erosion index units, plus a factor for runoff from snowmelt or applied water where such runoff is significant.

• K, the soil erodibility factor, is the soil loss rate per erosion index unit for a specified soil as measured on a unit plot, which is defined as a 72.6-ft length of uniform 9%

• L, the slope length factor, is the ratio of soil loss from the field slope length to that

• S, the slope-steepness factor, is the ratio of soil loss from the field slope gradient to

• C, the cover and management factor, is the ratio of soil loss from an area with specified cover and management to that from an identical area in tilled continuous

• P, the support practice factor, is the ratio of soil loss with a support practice like contouring, strip cropping, or terracing to that with straight-row farming up and

Out of the six factors, only two have units, the erosivity R and the erodibility K. The remaining four, i.e., slope steepness, slope length, cropping management, and components of control practice, are dimensionless factors, because they represent

Regarding R, it needs to be observed that its units stem from its definition as the energy multiplied by maximum 30-minute intensity term. In contrast, the erodibility factor values were determined empirically by calibration against measured erosion data. This is important due to the fact that if there is a change on the definition

of the R factor, then the values of the K factor should be recalculated [3].

Administration [4].

its validity has been established.

Soil Erosion - Rainfall Erosivity and Risk Assessment

The soil loss equation is

the period selected for R.

fallow.

16

down the slope.

ratios in relation to the unit plot [4].

slope continuously in clean-tilled fallow.

from a 72.6-ft length under identical conditions.

that from a 9% slope under otherwise identical conditions.

The R factor is computed for time periods greater than 20 years, so that wet or dry rainfall periods will be incorporated, eliminating any bias. It is set equal to the average of the sum of the erosivity values for every year's rainfalls. The R coefficient is defined as the product of the kinetic energy of an erosive rainfall event by the maximum intensity of a 30-minute duration rainfall, during the rainfall event:

$$R = \frac{1}{n} \sum\_{j=1}^{n} \sum\_{k=1}^{m\_j} (EI\_{30})\_k \tag{2}$$

where n is the number of years in the record, mj the number of erosive events during year j, and EI<sup>30</sup> (MJ mm ha�<sup>1</sup> h�<sup>1</sup> ) the rainfall erosivity for event k.

The erosivity of a particular event is

$$EI\_{30} = \left(\sum\_{r=1}^{m} e\_r \cdot v\_r \right) \cdot I\_{30} \tag{3}$$

where er is the kinetic energy per unit of rainfall (MJ ha�<sup>1</sup> mm�<sup>1</sup> ); vr is the rainfall depth (mm) for the time interval r of the hyetograph, which has been divided into r ¼ 1, 2, …, m subintervals; and I<sup>30</sup> is the maximum rainfall intensity for a 30-minute duration.

For the computation of er, numerous empirical relations (Figure 1) involving rain intensity have been proposed [6]. Thus, in USLE, the empirical relation of Wischmeier and Smith [4] was used:

$$e\_r = 0.119 + 0.0873 \log\_{10}(i) \tag{4}$$

with the upper limit of 0.283 MJ ha�<sup>1</sup> mm�<sup>1</sup> if i>76 mm h�<sup>1</sup> , where i is rainfall intensity.

#### Figure 1.

The three different kinetic energy equations used in USLE and its revisions. The points are actual data, showing the relation between intensity and the kinetic energy of the rainfall, converted to SI units and coming from Haan et al. [61].

Later, in RUSLE [5], the exponential relation of Brown and Foster [7] was used:

$$e\_r = 0.29 \left( 1 - 0.72 \text{ } e^{-0.05i} \right) \tag{5}$$

3. Methods for the estimation of R

DOI: http://dx.doi.org/10.5772/intechopen.85937

geostatistical models.

The lack of dense time series from pluviographs gives rise to many difficulties in the application of USLE and of its revisions in many countries. In order to cope with these difficulties, many models have been developed, based on rainfall depth values for various time steps (daily, monthly, and yearly), under specific spatial parameters and climatological data. Thus, a number of researchers reported good correlation results between R and yearly rainfall in many countries of the world by using various schemes ranging from simple parametric equations to

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods

These methods were applied in the countries of West Africa [13]; Switzerland taking account of elevation, aspect, longitude, and altitude [14]; the USA [15]; India [16]; Spain [17]; Italy [18]; China taking account of maximal rainfalls of 1-hour and 1-day durations per year [19]; Korea [20]; Mexico [21]; Honduras in

combination with elevation [22]; Rhodesia [23]; and Hawaii [24] as well as for the

Finally, for the estimation of the yearly value of R, several researchers used parametric equations for the initial determination of the daily R value and the corresponding daily rainfall depths. Those equations were applied in the Eastern USA [37], Australia [38], Spain [17], Kenya [39], China [40], Nigeria [41], North

Indicatively, some of the empirical equations of the literature are given below in dimensionless form, representing various parametrizations and methodologies. In West African countries, Roose [13] developed a simple linear relation between

In the USA, Renard and Freimund [15] proposed the following nonlinear equations for continental regions of the country, in which no rainfall intensity data exist,

In Morocco, Arnoldus [32] used the modified Fournier index, which is equal to

<sup>i</sup>¼<sup>1</sup> pm,i

MFI <sup>¼</sup> <sup>∑</sup><sup>12</sup>

where P is the yearly rainfall and pm,i the monthly rainfall from January to

where α2, β2, α3, β3, and γ<sup>3</sup> are nonlinear regression parameters.

R ¼ α<sup>1</sup> þ β<sup>1</sup> � P (8)

<sup>R</sup> <sup>¼</sup> <sup>α</sup><sup>2</sup> � <sup>P</sup><sup>β</sup><sup>2</sup> (9)

<sup>R</sup> <sup>¼</sup> <sup>α</sup><sup>3</sup> <sup>þ</sup> <sup>β</sup><sup>3</sup> � <sup>P</sup><sup>γ</sup><sup>3</sup> (10)

<sup>P</sup> (11)

use of parametric equations, of the yearly values of R to the monthly values of rainfall depth. Those equations were applied in Venezuela [27], Germany [28], the USA [15], Italy [29], Iran in combination with the yearly and the maximum daily rainfall per year [30], North Africa [31], Morocco [32], Nigeria [33], South Italy and

A different model followed in several countries concerns the correlation, by the

development of a European [25] and Universal model [26].

Southeast Australia [34], Uruguay [35], and Sudan [36].

Africa [42], South Italy [43], Peru [44], and Slovenia [45].

where α<sup>1</sup> και β<sup>1</sup> are linear regression parameters.

yearly values of R and rainfall P:

or

19

also for yearly values of R and rainfall P:

December, and developed the formula

Upon comparison of the equations used in USLE and RUSLE, McGregor and Mutchler [8] found out that for rain intensities up to 35 mm/h, Eq. (5) yields results by 12% less than those of Eq. (4), and they proposed a modification in the value 0.05 that controls the rate of change of er with i. Thus, in RUSLE, the following relation has been adopted:

$$e\_r = 0.29 \left( 1 - 0.72 \text{ } e^{-0.082i} \right) \tag{6}$$

Eq. (5) was developed for an application concerning reordered rainfall intensity data and not natural rainfall data [6]. The systematic underestimation of R has been shown in many other studies [3, 7–9], so Eq. (5) should not be used for calculations. The rules that apply in order to single out the storms causing erosion and to divide rainfalls of large duration in RUSLE2 [10, 11] are:


The current revised version of USLE, RUSLE2 [10], introduced the erosivity density (ED), as a measure of rainfall erosivity per unit rainfall to develop erosivity values for the USA, because ED requires shorter record lengths, as 10 years leads to acceptable results, allows more missing data than R, and is independent of the elevation:

$$ED\_j = \frac{\sum\_{k=1}^{m\_j} (EI\_{30})\_k}{P\_j} \tag{7}$$

where mj is the number of storms during a time period j, ð Þ EI<sup>30</sup> <sup>k</sup> the erosivity of storm k, and Pj the total precipitation height for the period j (usually monthly or annual).

Vantas et al. [12] using a numerical scheme with data from Greece showed that ED is more robust against the presence of missing precipitation values, as reflected in the following results:


The use of ED permits the utilization of data consisting of daily rainfall values that are more abundant in comparison to data from pluviographs. The details of the method, which was employed in the USA in combination with kriging, may be found in the RUSLE2 Science Documentation [10].

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods DOI: http://dx.doi.org/10.5772/intechopen.85937

#### 3. Methods for the estimation of R

Later, in RUSLE [5], the exponential relation of Brown and Foster [7] was used:

Upon comparison of the equations used in USLE and RUSLE, McGregor and Mutchler [8] found out that for rain intensities up to 35 mm/h, Eq. (5) yields results by 12% less than those of Eq. (4), and they proposed a modification in the value 0.05 that controls the rate of change of er with i. Thus, in RUSLE, the following

�0:05<sup>i</sup> (5)

�0:082<sup>i</sup> (6)

(7)

er ¼ 0:29 1 � 0:72 e

er ¼ 0:29 1 � 0:72 e

Eq. (5) was developed for an application concerning reordered rainfall intensity data and not natural rainfall data [6]. The systematic underestimation of R has been shown in many other studies [3, 7–9], so Eq. (5) should not be used for calculations. The rules that apply in order to single out the storms causing erosion and to divide

1. A rainfall event is divided into two parts, if its cumulative depth for duration

3. All rainfalls with extreme EI<sup>30</sup> values and a return period greater than 50 years

The current revised version of USLE, RUSLE2 [10], introduced the erosivity density (ED), as a measure of rainfall erosivity per unit rainfall to develop erosivity values for the USA, because ED requires shorter record lengths, as 10 years leads to acceptable results, allows more missing data than R, and is independent of the

> <sup>k</sup>¼<sup>1</sup>ð Þ EI<sup>30</sup> <sup>k</sup> Pj

where mj is the number of storms during a time period j, ð Þ EI<sup>30</sup> <sup>k</sup> the erosivity of storm k, and Pj the total precipitation height for the period j (usually monthly or

Vantas et al. [12] using a numerical scheme with data from Greece showed that ED is more robust against the presence of missing precipitation values, as reflected

• Using only 5% of the data, annual R values are underestimated on average by

• In the presence of 50% of the data, R values are underestimated by 49%, while

The use of ED permits the utilization of data consisting of daily rainfall values that are more abundant in comparison to data from pluviographs. The details of the method, which was employed in the USA in combination with kriging, may be

85%, when the average estimation error of ED values is 50%.

at the same time, the estimation error of ED is 20%.

found in the RUSLE2 Science Documentation [10].

EDj <sup>¼</sup> <sup>∑</sup>mj

2. A rainfall is considered erosive if it has a cumulative value greater than

relation has been adopted:

12.7 mm.

are deleted.

elevation:

annual).

18

in the following results:

rainfalls of large duration in RUSLE2 [10, 11] are:

Soil Erosion - Rainfall Erosivity and Risk Assessment

of 6 hours at a certain location is less than 1.27 mm.

The lack of dense time series from pluviographs gives rise to many difficulties in the application of USLE and of its revisions in many countries. In order to cope with these difficulties, many models have been developed, based on rainfall depth values for various time steps (daily, monthly, and yearly), under specific spatial parameters and climatological data. Thus, a number of researchers reported good correlation results between R and yearly rainfall in many countries of the world by using various schemes ranging from simple parametric equations to geostatistical models.

These methods were applied in the countries of West Africa [13]; Switzerland taking account of elevation, aspect, longitude, and altitude [14]; the USA [15]; India [16]; Spain [17]; Italy [18]; China taking account of maximal rainfalls of 1-hour and 1-day durations per year [19]; Korea [20]; Mexico [21]; Honduras in combination with elevation [22]; Rhodesia [23]; and Hawaii [24] as well as for the development of a European [25] and Universal model [26].

A different model followed in several countries concerns the correlation, by the use of parametric equations, of the yearly values of R to the monthly values of rainfall depth. Those equations were applied in Venezuela [27], Germany [28], the USA [15], Italy [29], Iran in combination with the yearly and the maximum daily rainfall per year [30], North Africa [31], Morocco [32], Nigeria [33], South Italy and Southeast Australia [34], Uruguay [35], and Sudan [36].

Finally, for the estimation of the yearly value of R, several researchers used parametric equations for the initial determination of the daily R value and the corresponding daily rainfall depths. Those equations were applied in the Eastern USA [37], Australia [38], Spain [17], Kenya [39], China [40], Nigeria [41], North Africa [42], South Italy [43], Peru [44], and Slovenia [45].

Indicatively, some of the empirical equations of the literature are given below in dimensionless form, representing various parametrizations and methodologies. In West African countries, Roose [13] developed a simple linear relation between yearly values of R and rainfall P:

$$R = a\_1 + \beta\_1 \cdot P \tag{8}$$

where α<sup>1</sup> και β<sup>1</sup> are linear regression parameters.

In the USA, Renard and Freimund [15] proposed the following nonlinear equations for continental regions of the country, in which no rainfall intensity data exist, also for yearly values of R and rainfall P:

$$R = a\_2 \cdot P^{\beta\_2} \tag{9}$$

or

$$R = a\_3 + \beta\_3 \cdot P^{\gamma\_3} \tag{10}$$

where α2, β2, α3, β3, and γ<sup>3</sup> are nonlinear regression parameters. In Morocco, Arnoldus [32] used the modified Fournier index, which is equal to

$$\text{MFI} = \frac{\sum\_{i=1}^{12} p\_{m,i}}{P} \tag{11}$$

where P is the yearly rainfall and pm,i the monthly rainfall from January to December, and developed the formula

Soil Erosion - Rainfall Erosivity and Risk Assessment

$$\mathbf{R} = a\_{\ $} \cdot MFI^{\beta\_{\$ }} \tag{12}$$

data. Thus, through the use of weekly instead of daily values, divisions of storms due to time step were prevented, and a measure of rainfall duration in days has been

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods

Table 1 gives the concise statistics of the cumulative rainfall, the duration, and

1. The relation of the precipitation height of the erosive events to the EI<sup>30</sup> value in Figure 3, where a linear relation appears with a wide opening between the

2. The relation of EI<sup>30</sup> to the months in Figure 4, where it turns out that in July its maximum appears with an almost sine-like variation of the median values per

Variable Min Mean Median Max SD SW KR CV

Precipitation (mm) 12.7 29.4 22.2 322.4 21.9 3.4 19.8 0.7 Duration (h) 0.5 13.0 10.5 152.5 10.5 2.5 13.6 0.8 SD is an abbreviation for standard deviation, SW for skew, KR for kurtosis, and CV for coefficient of variation.

Average statistical properties of cumulative rainfall, the duration and EI<sup>30</sup> of the erosive events.

) 0.7 89.1 44.2 3845.0 157.4 7.6 101.0 1.7

the erosivity EI<sup>30</sup> of the erosive events that resulted following RUSLE2 rules.

In the sequel, a series of three figures is given that contain:

logarithms of rainfall values and the logarithms of EI30.

added to the observations.

month.

Figure 2.

Table 1.

21

EI30 (MJ mm ha<sup>1</sup> h<sup>1</sup>

Stations' location coming from Greece used in the analysis.

4.2 Exploratory data analysis

DOI: http://dx.doi.org/10.5772/intechopen.85937

where α<sup>5</sup> and β<sup>5</sup> are nonlinear regression parameters.

Richardson et al. [37] used the daily rainfall values for the estimation of the corresponding daily erosivity values Rd:

$$\mathbf{R}\_{\mathbf{d}} = a\_{m,s} \cdot P\_{d}^{\beta\_{m,s}} \tag{13}$$

where the parameters αm,s and βm,s correspond to station s and to month m and are evaluated by means of nonlinear regression.

A serious issue with the applicability of Eq. (13) is the difficulty to compute its coefficients in dry areas or months due to small number of erosive events. Another issue, even in wet periods, is the amount of missing values that also leads to a small set of observations and the same problems [46].

In order to reduce the computation burden of Eq. (13), Yu and Rosewell [38] proposed the following empirical equation:

$$R\_d = a\_t \cdot \left[\mathbf{1} + \beta\_s \cdot \cos\left(\frac{\pi}{6}m - \alpha\right)\right] P\_d^{\chi\_i} \tag{14}$$

where the parameters αs, βs, and γ<sup>s</sup> are different only for the different stations s and ω is the month with the highest average of daily R values.

It must be noted that the logarithmic transformation and the subsequent application of linear regression must be avoided in the above equations, because in that case the minimization is applied to the average of the differences of the logarithms of the aggregated daily EI30values (coming from pluviograph data) and the logarithms of the estimated Rd values (coming from daily rainfall data) and not of the corresponding absolute values themselves. Specifically, the application of the log transformation to Eq. (13) resulted in a �10% systematic error, when data from RUSLE2 was used, as noted on p. 59 of [10].

A different approach for the estimation of R, in the presence of missing values, was followed by Vantas and Sidiropoulos [47]. They compared empirical equations to machine learning methods, extracting for the latter methods better results, in terms of statistical significance, compared to the former methods.

#### 4. Utilization of neural networks for the estimation of R

#### 4.1 Data

The data utilized in the analysis were taken from the Greek National Bank of Hydrological and Meteorological Information [48] and came from 84 meteorological stations (Figure 2). By adding up the years of record registered on the various stations, the time series comprised a total of 2425 years of pluviograph records, with an average of 28.9 years per station. The time step was 30 min for the time period from 1953 to 1997, and the data coverage was equal to 56%.

For the above described data, it was deemed useful to change the time step and aggregate the data into weekly values, because 57% of the recordings were associated with storms occupying time periods covering parts of more than one calendar day, although only 17% of the storms had duration of more than 24 hours.

Under the time step of 1 week, it was found out that 80% of the values emanated from a single storm. The storms that were crossed temporally by 2 consecutive weeks were assigned to the first of the 2 weeks, and they comprised only 7% of the

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods DOI: http://dx.doi.org/10.5772/intechopen.85937

data. Thus, through the use of weekly instead of daily values, divisions of storms due to time step were prevented, and a measure of rainfall duration in days has been added to the observations.

## 4.2 Exploratory data analysis

<sup>R</sup> <sup>¼</sup> <sup>α</sup><sup>5</sup> � MFI<sup>β</sup><sup>5</sup> (12)

<sup>d</sup> (13)

where α<sup>5</sup> and β<sup>5</sup> are nonlinear regression parameters.

corresponding daily erosivity values Rd:

Soil Erosion - Rainfall Erosivity and Risk Assessment

are evaluated by means of nonlinear regression.

set of observations and the same problems [46].

proposed the following empirical equation:

RUSLE2 was used, as noted on p. 59 of [10].

4.1 Data

20

Richardson et al. [37] used the daily rainfall values for the estimation of the

Rd <sup>¼</sup> <sup>α</sup>m,s � <sup>P</sup><sup>β</sup>m,s

where the parameters αm,s and βm,s correspond to station s and to month m and

A serious issue with the applicability of Eq. (13) is the difficulty to compute its coefficients in dry areas or months due to small number of erosive events. Another issue, even in wet periods, is the amount of missing values that also leads to a small

In order to reduce the computation burden of Eq. (13), Yu and Rosewell [38]

h i � �

where the parameters αs, βs, and γ<sup>s</sup> are different only for the different stations s

It must be noted that the logarithmic transformation and the subsequent application of linear regression must be avoided in the above equations, because in that case the minimization is applied to the average of the differences of the logarithms of the aggregated daily EI30values (coming from pluviograph data) and the logarithms of the estimated Rd values (coming from daily rainfall data) and not of the corresponding absolute values themselves. Specifically, the application of the log transformation to Eq. (13) resulted in a �10% systematic error, when data from

A different approach for the estimation of R, in the presence of missing values, was followed by Vantas and Sidiropoulos [47]. They compared empirical equations to machine learning methods, extracting for the latter methods better results, in

The data utilized in the analysis were taken from the Greek National Bank of Hydrological and Meteorological Information [48] and came from 84 meteorological stations (Figure 2). By adding up the years of record registered on the various stations, the time series comprised a total of 2425 years of pluviograph records, with an average of 28.9 years per station. The time step was 30 min for the time period

For the above described data, it was deemed useful to change the time step and aggregate the data into weekly values, because 57% of the recordings were associated with storms occupying time periods covering parts of more than one calendar

Under the time step of 1 week, it was found out that 80% of the values emanated

day, although only 17% of the storms had duration of more than 24 hours.

from a single storm. The storms that were crossed temporally by 2 consecutive weeks were assigned to the first of the 2 weeks, and they comprised only 7% of the

<sup>6</sup> <sup>m</sup> � <sup>ω</sup>

Pγs

<sup>d</sup> (14)

Rd <sup>¼</sup> <sup>α</sup><sup>s</sup> � <sup>1</sup> <sup>þ</sup> <sup>β</sup><sup>s</sup> � cos <sup>π</sup>

and ω is the month with the highest average of daily R values.

terms of statistical significance, compared to the former methods.

4. Utilization of neural networks for the estimation of R

from 1953 to 1997, and the data coverage was equal to 56%.

Table 1 gives the concise statistics of the cumulative rainfall, the duration, and the erosivity EI<sup>30</sup> of the erosive events that resulted following RUSLE2 rules. In the sequel, a series of three figures is given that contain:


Figure 2. Stations' location coming from Greece used in the analysis.


#### Table 1.

Average statistical properties of cumulative rainfall, the duration and EI<sup>30</sup> of the erosive events.

Figure 3. Relation of cumulative precipitation of the erosive events to the EI30. The red line is the LOESS line.

30 mm, a very large range of EI<sup>30</sup> values, from about 10 to 700 MJ mm ha�<sup>1</sup> h�<sup>1</sup>

Relation of cumulative precipitation of the erosive events to the EI30. The red line is the LOESS line that is fitted

The outcomes, denoted as <sup>y</sup>ð Þ<sup>i</sup> , where ð Þ<sup>i</sup> is the index of a sample consisting of 19,688 values, represented weekly cumulative rainfall erosivity as calculated from pluviograph data. The features included the weekly cumulative rainfall of erosive events, the month to which the above individual m values are referred, the altitude of the meteorological stations, and the number of the days of the week for which rainfall is recorded. The matrix of the features is denoted as x(i), while h(x(i)) denotes the output vector corresponding to the hypothesis in question (i.e., as computed either from the neural network or the empirical equations). The subscripts test and train mark quantities that belong to the testing or the training set,

The erosivity estimation problem is set up as a scheme of machine learning. The data were split using 70% of them as the training set and 30% as a testing one. As

<sup>i</sup>¼<sup>1</sup> h xtestð Þ<sup>i</sup> � ytestð Þ<sup>i</sup> <sup>2</sup>

ð Þi

test <sup>2</sup> (15)

<sup>i</sup>¼<sup>1</sup> ^ytrain � <sup>y</sup>

• <sup>R</sup><sup>2</sup> <sup>¼</sup> 0 means that there is no improvement comparing to a simplistic model

• R<sup>2</sup> <0 means an algorithm that makes prediction worse than the base model.

to be assigned, corresponding to the erosive events connected to 30 mm.

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods

4.3 Methods

Figure 5.

respectively.

23

measures of the out-of-sample error were used:

:

• <sup>R</sup><sup>2</sup> <sup>¼</sup> 1 means a perfect algorithm.

<sup>R</sup><sup>2</sup> <sup>¼</sup> <sup>1</sup> � <sup>∑</sup><sup>m</sup>

that returns the average of the training set values ^ytrain.

∑<sup>m</sup>

(a) The coefficient of determination:

using the median values of EI<sup>30</sup> per duration width.

DOI: http://dx.doi.org/10.5772/intechopen.85937

with this form of R<sup>2</sup>

, is

#### Figure 4.

Relation of cumulative precipitation of the erosive events to the EI30. The red line is the LOESS line that is fitted using the median values of EI<sup>30</sup> per month.

3. The relation of the duration of erosive events to EI<sup>30</sup> in Figure 5, where it is shown that for the selected duration widths, the relation to the median values of EI<sup>30</sup> is parabolic.

In these three figures, local polynomial regression fitting (LOESS, [49]) was used as a nonparametric method to produce smooth curves through the plotted variables. In that method, the fitting is done locally for a given point x using all the points in the neighborhood of x and the distances of these points from x as weights.

From the above diagrams and from the statistics of the erosive events, it can be seen that the estimation of EI<sup>30</sup> values is a nonlinear problem with skewed data and with expected results of a relatively low accuracy. This point is particularly reinforced by Figure 3, in which it is shown that for a specific rainfall depth, such as Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods DOI: http://dx.doi.org/10.5772/intechopen.85937

#### Figure 5.

Relation of cumulative precipitation of the erosive events to the EI30. The red line is the LOESS line that is fitted using the median values of EI<sup>30</sup> per duration width.

30 mm, a very large range of EI<sup>30</sup> values, from about 10 to 700 MJ mm ha�<sup>1</sup> h�<sup>1</sup> , is to be assigned, corresponding to the erosive events connected to 30 mm.

#### 4.3 Methods

The outcomes, denoted as <sup>y</sup>ð Þ<sup>i</sup> , where ð Þ<sup>i</sup> is the index of a sample consisting of 19,688 values, represented weekly cumulative rainfall erosivity as calculated from pluviograph data. The features included the weekly cumulative rainfall of erosive events, the month to which the above individual m values are referred, the altitude of the meteorological stations, and the number of the days of the week for which rainfall is recorded. The matrix of the features is denoted as x(i), while h(x(i)) denotes the output vector corresponding to the hypothesis in question (i.e., as computed either from the neural network or the empirical equations). The subscripts test and train mark quantities that belong to the testing or the training set, respectively.

The erosivity estimation problem is set up as a scheme of machine learning. The data were split using 70% of them as the training set and 30% as a testing one. As measures of the out-of-sample error were used:

(a) The coefficient of determination:

$$R^2 = \mathbf{1} - \frac{\sum\_{i=1}^{m} \left( h\left(\mathbf{x}\_{\text{test}}^{(i)}\right) - \mathbf{y}\_{\text{test}}^{(i)} \right)^2}{\sum\_{i=1}^{m} \left( \hat{\boldsymbol{\mathcal{Y}}}\_{\text{train}} - \mathbf{y}\_{\text{test}}^{(i)} \right)^2} \tag{15}$$

with this form of R<sup>2</sup> :


3. The relation of the duration of erosive events to EI<sup>30</sup> in Figure 5, where it is shown that for the selected duration widths, the relation to the median values

Relation of cumulative precipitation of the erosive events to the EI30. The red line is the LOESS line that is fitted

Relation of cumulative precipitation of the erosive events to the EI30. The red line is the LOESS line.

Soil Erosion - Rainfall Erosivity and Risk Assessment

In these three figures, local polynomial regression fitting (LOESS, [49]) was used as a nonparametric method to produce smooth curves through the plotted variables. In that method, the fitting is done locally for a given point x using all the points in the neighborhood of x and the distances of these points from x as weights. From the above diagrams and from the statistics of the erosive events, it can be seen that the estimation of EI<sup>30</sup> values is a nonlinear problem with skewed data and

reinforced by Figure 3, in which it is shown that for a specific rainfall depth, such as

with expected results of a relatively low accuracy. This point is particularly

of EI<sup>30</sup> is parabolic.

using the median values of EI<sup>30</sup> per month.

Figure 3.

Figure 4.

22

(b) The root mean-squared error:

$$RMSE = \sqrt{\frac{1}{m} \sum\_{i=1}^{m} \left( h(\mathbf{x}\_{\text{test}}^{(i)}) - \mathbf{y}\_{\text{test}}^{(i)} \right)^{2}} \tag{16}$$

(c) The mean absolute error:

$$\text{MAE} = \frac{1}{m} \sum\_{i=1}^{m} \left| h\left(\mathbf{x}\_{\text{test}}^{(i)}\right) - \mathbf{y}\_{\text{test}}^{(i)} \right| \tag{17}$$

In Eqs. (15)–(17), m is the number of the samples in the test set, h xtestð Þ<sup>i</sup> � � is the estimations from a trained algorithm (i.e., either from a trained neural network or from a fitted empirical equation) used as inputs, the testing set xtest and ytest are the outcomes from calculations by means of pluviograph data, and ^ytrain is the average of the training set erosivity values coming also from pluviograph data.

In order to compare neural network performance to that of empirical equations, two alternative exponential models discussed in the previous section were tried, namely, those given by Eqs. (13) and (14) and leading to optimal adjustments of two respective hypotheses. The latter were determined by minimizing the following objective function by the trust-region-reflective method [50, 51]:

$$J(\theta) = \frac{1}{2m} \sum\_{i=1}^{m} \left( h\_{\theta} \left( \mathbf{x}^{(i)} \right) - \mathbf{y}^{(i)} \right)^{2} \tag{18}$$

The optimization of the Θð Þ<sup>j</sup> weights is executed in such a way as to minimize the

where m is the number of values h<sup>Θ</sup> xð Þ<sup>i</sup> obtained from the network using xð Þ<sup>i</sup> as input data and yð Þ<sup>i</sup> as the calculated values (i.e., the EI<sup>30</sup> values as computed from

The architecture of the neural network used in the analysis included two hidden layers that had 32 and 16 neurons, respectively. In order to keep h<sup>Θ</sup> xð Þ<sup>i</sup> nonnegative, the second hidden layer of the neural network used the rectifier activation

The training set values xtrain were used to compute the averages x and standard deviation sd xð Þ and normalize both xtrain and xtest using the normalizing transfor-

N xð Þ<sup>i</sup> <sup>¼</sup> <sup>x</sup>ð Þ<sup>i</sup> � <sup>x</sup>

The training of the networks was performed by the method of early stopping [55] by utilizing a random validation set consisting of 10% of the training data, so as

The data importing, calculation of EI<sup>30</sup> values, and analysis were done using the

R language for statistical computing and graphics [56] using the packages:

g xð Þ¼ max 0ð Þ ; x (22)

sd xð Þ (23)

This minimization is a difficult optimization problem, since there are no constraints related to these parameters. The parameters are usually initialized by random values, and in the sequel, specialized algorithms are used for their determination. A survey about neural networks and their forms as shallow and deep can be

<sup>h</sup><sup>Θ</sup> <sup>x</sup>ð Þ<sup>i</sup> � <sup>y</sup>ð Þ<sup>i</sup> <sup>2</sup>

(21)

sum of squares of the differences between computed and measured results:

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods

1 <sup>2</sup><sup>m</sup> <sup>∑</sup> m i¼1

Jð Þ¼ Θ

found in Schmidhuber [52] and Goodfellow et al. [53].

to avoid over fitting of the neural network.

the pluviographic data).

A neural network with a single hidden layer.

DOI: http://dx.doi.org/10.5772/intechopen.85937

function [54]:

Figure 6.

mation:

25

where θ denotes the vector of the respective parameters and y(i) are the outcomes, as defined in the beginning of this section.

Neural networks aim at mimicking the function of the brain. They constitute powerful nonlinear regression methods. The result or exit of the neural network is produced via a series of intermediate nodes arranged in successive layers, characterized as hidden. These layers mediate between an initial layer of input nodes and a final layer of exit nodes. The hidden nodes define linear combinations of the data contained in the initial input layer. These linear combinations are further transformed by a nonlinear function, which possesses a continuous derivative, such as the hyperbolic tangent function:

$$\mathbf{g}(\mathbf{x}) = \tanh(\mathbf{x}) \tag{19}$$

The example in Figure 6 displays three layers, of which the one to the left is the input layer, the middle one is the hidden layer, and the third one is the output layer. In the neural network, αð Þ<sup>j</sup> <sup>i</sup> denotes the activation of node i in layer j, and Θð Þ<sup>j</sup> denotes the matrix of weights that serve as coefficients of the linear transformation from layer j to layer j þ 1. The values at the output nodes in the example of Figure 6 are:

$$a\_1^{(2)} = \mathbf{g} \left( \Theta\_{10}^{(1)} \cdot \mathbf{x}\_0 + \Theta\_{11}^{(1)} \cdot \mathbf{x}\_1 + \Theta\_{12}^{(1)} \cdot \mathbf{x}\_2 + \Theta\_{13}^{(1)} \cdot \mathbf{x}\_3 \right)$$

$$a\_2^{(2)} = \mathbf{g} \left( \Theta\_{20}^{(1)} \cdot \mathbf{x}\_0 + \Theta\_{21}^{(1)} \cdot \mathbf{x}\_1 + \Theta\_{22}^{(1)} \cdot \mathbf{x}\_2 + \Theta\_{23}^{(1)} \cdot \mathbf{x}\_3 \right)$$

$$a\_3^{(2)} = \mathbf{g} \left( \Theta\_{30}^{(1)} \cdot \mathbf{x}\_0 + \Theta\_{31}^{(1)} \cdot \mathbf{x}\_1 + \Theta\_{32}^{(1)} \cdot \mathbf{x}\_2 + \Theta\_{33}^{(1)} \cdot \mathbf{x}\_3 \right)$$

$$a\_4^{(2)} = \mathbf{g} \left( \Theta\_{40}^{(1)} \cdot \mathbf{x}\_0 + \Theta\_{41}^{(1)} \cdot \mathbf{x}\_1 + \Theta\_{42}^{(1)} \cdot \mathbf{x}\_2 + \Theta\_{43}^{(1)} \cdot \mathbf{x}\_3 \right)$$

$$a\_1^{(3)} = \mathbf{g} \left( \Theta\_{10}^{(2)} \cdot a\_0^{(2)} + \Theta\_{11}^{(2)} \cdot a\_1^{(2)} + \Theta\_{12}^{(2)} \cdot a\_2^{(2)} + \Theta\_{13}^{(2)} \cdot a\_3^{(2)} + \Theta\_{14}^{(2)} \cdot a\_4^{(2)} \right) \tag{20}$$

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods DOI: http://dx.doi.org/10.5772/intechopen.85937

Figure 6. A neural network with a single hidden layer.

(b) The root mean-squared error:

Soil Erosion - Rainfall Erosivity and Risk Assessment

(c) The mean absolute error:

RMSE ¼

MAE <sup>¼</sup> <sup>1</sup>

1 <sup>m</sup> <sup>∑</sup> m i¼1

<sup>m</sup> <sup>∑</sup> m i¼1

of the training set erosivity values coming also from pluviograph data.

objective function by the trust-region-reflective method [50, 51]:

<sup>2</sup><sup>m</sup> <sup>∑</sup> m i¼1

contained in the initial input layer. These linear combinations are further

<sup>J</sup>ð Þ¼ <sup>θ</sup> <sup>1</sup>

comes, as defined in the beginning of this section.

as the hyperbolic tangent function:

αð Þ<sup>2</sup>

αð Þ<sup>2</sup>

αð Þ<sup>2</sup>

αð Þ<sup>2</sup>

<sup>10</sup> � <sup>α</sup>ð Þ<sup>2</sup>

<sup>1</sup> <sup>¼</sup> <sup>g</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>2</sup> <sup>¼</sup> <sup>g</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>3</sup> <sup>¼</sup> <sup>g</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>4</sup> <sup>¼</sup> <sup>g</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>0</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>2</sup>

In the neural network, αð Þ<sup>j</sup>

αð Þ<sup>3</sup>

24

<sup>1</sup> <sup>¼</sup> <sup>g</sup> <sup>Θ</sup>ð Þ<sup>2</sup>

� � �

s

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

h xtestð Þ<sup>i</sup> � �

In Eqs. (15)–(17), m is the number of the samples in the test set, h xtestð Þ<sup>i</sup> � � is the estimations from a trained algorithm (i.e., either from a trained neural network or from a fitted empirical equation) used as inputs, the testing set xtest and ytest are the outcomes from calculations by means of pluviograph data, and ^ytrain is the average

In order to compare neural network performance to that of empirical equations, two alternative exponential models discussed in the previous section were tried, namely, those given by Eqs. (13) and (14) and leading to optimal adjustments of two respective hypotheses. The latter were determined by minimizing the following

> <sup>h</sup><sup>θ</sup> <sup>x</sup>ð Þ<sup>i</sup> � � � <sup>y</sup>ð Þ<sup>i</sup> � �<sup>2</sup>

where θ denotes the vector of the respective parameters and y(i) are the out-

Neural networks aim at mimicking the function of the brain. They constitute powerful nonlinear regression methods. The result or exit of the neural network is produced via a series of intermediate nodes arranged in successive layers, characterized as hidden. These layers mediate between an initial layer of input nodes and a final layer of exit nodes. The hidden nodes define linear combinations of the data

transformed by a nonlinear function, which possesses a continuous derivative, such

The example in Figure 6 displays three layers, of which the one to the left is the input layer, the middle one is the hidden layer, and the third one is the output layer.

the matrix of weights that serve as coefficients of the linear transformation from layer

<sup>11</sup> � <sup>x</sup><sup>1</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>21</sup> � <sup>x</sup><sup>1</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>31</sup> � <sup>x</sup><sup>1</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>41</sup> � <sup>x</sup><sup>1</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>12</sup> � <sup>α</sup>ð Þ<sup>2</sup>

� �

� �

� �

� �

� �

<sup>2</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>2</sup>

j to layer j þ 1. The values at the output nodes in the example of Figure 6 are:

<sup>10</sup> � <sup>x</sup><sup>0</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>20</sup> � <sup>x</sup><sup>0</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>30</sup> � <sup>x</sup><sup>0</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>40</sup> � <sup>x</sup><sup>0</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>1</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>2</sup>

<sup>11</sup> � <sup>α</sup>ð Þ<sup>2</sup>

h xtestð Þ<sup>i</sup> ð Þ� ytestð Þ<sup>i</sup> � �<sup>2</sup>

� ytest ð Þi

� �

g xð Þ¼ tanhð Þ x (19)

<sup>i</sup> denotes the activation of node i in layer j, and Θð Þ<sup>j</sup> denotes

<sup>12</sup> � <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>22</sup> � <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>32</sup> � <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>42</sup> � <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>1</sup>

<sup>13</sup> � <sup>α</sup>ð Þ<sup>2</sup>

<sup>13</sup> � x<sup>3</sup>

<sup>23</sup> � x<sup>3</sup>

<sup>33</sup> � x<sup>3</sup>

<sup>43</sup> � x<sup>3</sup>

<sup>14</sup> � <sup>α</sup>ð Þ<sup>2</sup> 4

(20)

<sup>3</sup> <sup>þ</sup> <sup>Θ</sup>ð Þ<sup>2</sup>

� (17)

(16)

(18)

The optimization of the Θð Þ<sup>j</sup> weights is executed in such a way as to minimize the sum of squares of the differences between computed and measured results:

$$J(\boldsymbol{\Theta}) = \frac{1}{2m} \sum\_{i=1}^{m} \left( h\_{\boldsymbol{\Theta}} \left( \boldsymbol{\chi}^{(i)} \right) - \boldsymbol{\chi}^{(i)} \right)^{2} \tag{21}$$

where m is the number of values h<sup>Θ</sup> xð Þ<sup>i</sup> obtained from the network using xð Þ<sup>i</sup> as input data and yð Þ<sup>i</sup> as the calculated values (i.e., the EI<sup>30</sup> values as computed from the pluviographic data).

This minimization is a difficult optimization problem, since there are no constraints related to these parameters. The parameters are usually initialized by random values, and in the sequel, specialized algorithms are used for their determination. A survey about neural networks and their forms as shallow and deep can be found in Schmidhuber [52] and Goodfellow et al. [53].

The architecture of the neural network used in the analysis included two hidden layers that had 32 and 16 neurons, respectively. In order to keep h<sup>Θ</sup> xð Þ<sup>i</sup> nonnegative, the second hidden layer of the neural network used the rectifier activation function [54]:

$$\mathbf{g}(\mathbf{x}) = \max(\mathbf{0}, \mathbf{x}) \tag{22}$$

The training set values xtrain were used to compute the averages x and standard deviation sd xð Þ and normalize both xtrain and xtest using the normalizing transformation:

$$N\left(\mathbf{x}^{(i)}\right) = \frac{\mathbf{x}^{(i)} - \overline{\mathbf{x}}}{s d(\mathbf{x})} \tag{23}$$

The training of the networks was performed by the method of early stopping [55] by utilizing a random validation set consisting of 10% of the training data, so as to avoid over fitting of the neural network.

The data importing, calculation of EI<sup>30</sup> values, and analysis were done using the R language for statistical computing and graphics [56] using the packages:

### Soil Erosion - Rainfall Erosivity and Risk Assessment


R2 values are unitless and RMSE and MAE values are in MJ mm ha<sup>1</sup> h<sup>1</sup> . The metrics refer to the weekly erosivity time step.

Table 2.

Estimation of the out-of-sample error metrics for the two empirical equations and the neural network.


R2 values are unitless and RMSE and MAE values are in MJ mm ha<sup>1</sup> h<sup>1</sup> y 1 . The values are calculated using the estimation of the annual erosivity values.

#### Table 3.

Estimation of the out-of-sample error metrics for the two empirical equations and the neural network.

hydroscoper [57], hyetor [58], and ggplot2 [59]. The open source software library TensorFlow [60] was used in order to train the neural network.

### 4.4 Results

The comparison of the algorithms was based on their performance on the testing set using the three different error metrics and two time steps: weekly and annual.

The weekly and the annual calculated values were produced from the aggregated erosivity values coming from the pluviograph data. The annual predicted values are also the aggregated values coming from the output of the corresponding algorithm (the predicted weekly erosivity values in our case). The annual time step was used in order to examine if there is a systematic error on the estimations of the algo-

Calculated annual erosivity values coming from pluviograph data versus predicted values using the neural

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods

DOI: http://dx.doi.org/10.5772/intechopen.85937

network on the testing set. With red line, the identity function f(x) = x is symbolized.

These results are given in Tables 2 and 3, where it is shown that the neural network outperformed the two parametric methods, which both had similar results. Specifically, the machine learning method had better performance that ranged from

The universal soil loss equation with its revisions represents the established model and the means for estimating annual long-term average soil loss. An important factor of this equation is the rainfall erosivity, R, which is closely connected to the energy carried by rainfall. Rainfall data are essential for the estimation of R. In particular, data from pluviographs are more suitable for that task. However, such data are scarce, and various empirical methods have been presented for estimating

R from coarser data. These methods are reviewed in this chapter.

In Figures 7 and 8, the annual erosivity values are presented as calculated in the test set versus the values estimated from the 84 fitted Yu and Rosewell equations for each station and the single neural network. In these two diagrams, the neural network's estimations are closer to the calculated ones than those coming from the

rithms, which was not the case.

empirical equations.

Figure 8.

5. Conclusions

27

8 to 22%, depending on the error metric and time step.

#### Figure 7.

Calculated annual erosivity values coming from pluviograph data versus predicted values using the Yu and Rosewell model on the testing set. With red line, the identity function f(x) = x is symbolized.

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods DOI: http://dx.doi.org/10.5772/intechopen.85937

Figure 8.

hydroscoper [57], hyetor [58], and ggplot2 [59]. The open source software library

Estimation of the out-of-sample error metrics for the two empirical equations and the neural network.

y 1

Method RMSE MAE R<sup>2</sup> Richardson et al. 141.94 65.46 0.59 Yu and Rosewell 140.36 64.33 0.60 Neural network 116.83 54.33 0.73

Estimation of the out-of-sample error metrics for the two empirical equations and the neural network.

Method RMSE MAE R2 Richardson et al. 261.73 144.63 0.76 Yu and Rosewell 259.80 141.48 0.76 Neural network 223.01 124.24 0.82

. The metrics refer to the weekly erosivity time step.

. The values are calculated using the estimation of the

The comparison of the algorithms was based on their performance on the testing set using the three different error metrics and two time steps: weekly and annual.

Calculated annual erosivity values coming from pluviograph data versus predicted values using the Yu and

Rosewell model on the testing set. With red line, the identity function f(x) = x is symbolized.

TensorFlow [60] was used in order to train the neural network.

R2 values are unitless and RMSE and MAE values are in MJ mm ha<sup>1</sup> h<sup>1</sup>

Soil Erosion - Rainfall Erosivity and Risk Assessment

R2 values are unitless and RMSE and MAE values are in MJ mm ha<sup>1</sup> h<sup>1</sup>

4.4 Results

Figure 7.

26

annual erosivity values.

Table 3.

Table 2.

Calculated annual erosivity values coming from pluviograph data versus predicted values using the neural network on the testing set. With red line, the identity function f(x) = x is symbolized.

The weekly and the annual calculated values were produced from the aggregated erosivity values coming from the pluviograph data. The annual predicted values are also the aggregated values coming from the output of the corresponding algorithm (the predicted weekly erosivity values in our case). The annual time step was used in order to examine if there is a systematic error on the estimations of the algorithms, which was not the case.

These results are given in Tables 2 and 3, where it is shown that the neural network outperformed the two parametric methods, which both had similar results. Specifically, the machine learning method had better performance that ranged from 8 to 22%, depending on the error metric and time step.

In Figures 7 and 8, the annual erosivity values are presented as calculated in the test set versus the values estimated from the 84 fitted Yu and Rosewell equations for each station and the single neural network. In these two diagrams, the neural network's estimations are closer to the calculated ones than those coming from the empirical equations.

### 5. Conclusions

The universal soil loss equation with its revisions represents the established model and the means for estimating annual long-term average soil loss. An important factor of this equation is the rainfall erosivity, R, which is closely connected to the energy carried by rainfall. Rainfall data are essential for the estimation of R. In particular, data from pluviographs are more suitable for that task. However, such data are scarce, and various empirical methods have been presented for estimating R from coarser data. These methods are reviewed in this chapter.

Also, various difficulties are discussed that are associated with the issues of missing values and the separation of erosive events. In addition, indications are given of the adverse nature of the relation of erosivity to precipitation characteristics. All these facts raise objections to the suitability of empirical models and lead to the use of data-driven methods that will be able to handle more effectively the problematic nature of the data and also to discover more reliably the hidden relations between erosivity, precipitation, seasonality, and local characteristics. This is achieved, because in machine learning procedures, no human intervention or bias is involved in the selection or in the forms of these hidden nonlinearities.

More specifically, neural networks, as machine learning tools, are employed in this chapter for the estimation of erosivity, in comparison to two different empirical equations of the literature. The results demonstrate, first of all, the superior performance of the machine learning method and its ability for generalization, so as to perform equally well upon new data. Another advantage of the proposed learning algorithm is that it produces a single model accessing all available data, in contrast to empirical equations, where a researcher must fit and handle as many nonlinear equations as the number of available stations or the product of the number of the stations by the number of the used seasons, making the analysis cumbersome and prone to errors.

The latter fact is exemplified on the country-wide data from Greece, used in this chapter. The single neural network, which was produced, outperformed 84 different empirical nonlinear equations, each one of which was fitted using data from each station. The comparison of the performance was made using test data sets to which neither the empirical equations nor the neural network had any access during their training.

In conclusion, the present chapter gives an indicative and characteristic sample of the use of machine learning methods in the problems associated with erosivity. There are still many related problems, such as the spatial and temporal distribution of erosivity, with significant potential for the application of machine learning methods. The same is true of many areas of hydrology and water resources.

Author details

Macedonia, Greece

29

Konstantinos Vantas, Epaminondas Sidiropoulos\* and Chris Evangelides Faculty of Engineering, Aristotle University of Thessaloniki, Thessaloniki,

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods

DOI: http://dx.doi.org/10.5772/intechopen.85937

© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium,

\*Address all correspondence to: nontas@topo.auth.gr

provided the original work is properly cited.

## Conflict of interest

The authors declare no conflict of interest.

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods DOI: http://dx.doi.org/10.5772/intechopen.85937

## Author details

Also, various difficulties are discussed that are associated with the issues of missing values and the separation of erosive events. In addition, indications are given of the adverse nature of the relation of erosivity to precipitation characteristics. All these facts raise objections to the suitability of empirical models and lead to the use of data-driven methods that will be able to handle more effectively the problematic nature of the data and also to discover more reliably the hidden relations between erosivity, precipitation, seasonality, and local characteristics. This is achieved, because in machine learning procedures, no human intervention or bias is involved in the selection or in the forms of these hidden

Soil Erosion - Rainfall Erosivity and Risk Assessment

More specifically, neural networks, as machine learning tools, are employed in this chapter for the estimation of erosivity, in comparison to two different empirical equations of the literature. The results demonstrate, first of all, the superior performance of the machine learning method and its ability for generalization, so as to perform equally well upon new data. Another advantage of the proposed learning algorithm is that it produces a single model accessing all available data, in contrast to empirical equations, where a researcher must fit and handle as many nonlinear equations as the number of available stations or the product of the number of the stations by the number of the used seasons, making the analysis cumbersome and

The latter fact is exemplified on the country-wide data from Greece, used in this chapter. The single neural network, which was produced, outperformed 84 different empirical nonlinear equations, each one of which was fitted using data from each station. The comparison of the performance was made using test data sets to which neither the empirical equations nor the neural network had any access during

In conclusion, the present chapter gives an indicative and characteristic sample of the use of machine learning methods in the problems associated with erosivity. There are still many related problems, such as the spatial and temporal distribution of erosivity, with significant potential for the application of machine

learning methods. The same is true of many areas of hydrology and water

nonlinearities.

prone to errors.

their training.

resources.

28

Conflict of interest

The authors declare no conflict of interest.

Konstantinos Vantas, Epaminondas Sidiropoulos\* and Chris Evangelides Faculty of Engineering, Aristotle University of Thessaloniki, Thessaloniki, Macedonia, Greece

\*Address all correspondence to: nontas@topo.auth.gr

© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## References

[1] Boardman J, Poesen J. Soil erosion in Europe: Major processes, causes and consequences. In: Soil Erosion in Europe. Wiley-Blackwell; 2006. pp. 477-487

[2] Kinnell PIA. Event soil loss, runoff and the universal soil loss equation family of models: A review. Journal of Hydrology. 2010;385(1):384-397

[3] Nearing MA, Yin S, Borrelli P, Polyakov VO. Rainfall erosivity: An historical review. Catena. 2017;157: 357-362

[4] Wischmeier WH, Smith DD. Predicting Rainfall Erosion Losses— A Guide to Conservation Planning. USDA, Science and Education Administration; 1986

[5] Renard KG, Foster GR, Weesies GA, McCool DK, Yoder DC. Predicting Soil Erosion by Water: A Guide to Conservation Planning with the Revised Universal Soil Loss Equation (RUSLE). Vol. Agriculture Handbook No. 703. Washington, DC: United States Department of Agriculture; 1997

[6] Brown L, Foster G. Storm erosivity using idealized intensity distributions. Transactions of the ASAE. 1987; 30(2):379-386

[7] Ramon R, Minella JP, Merten GH, de Barros CA, Canale T. Kinetic energy estimation by rainfall intensity and its usefulness in predicting hydrosedimentological variables in a small rural catchment in southern Brazil. Catena. 2017;148:176-184

[8] McGregor KC, Mutchler CK. Status of the R factor in northern Mississippi. In: Soil Erosion: Prediction and Control. Soil Conservation Society of America; 1976. pp. 135-142

[9] Van Dijk A, Bruijnzeel LA, Rosewell CJ. Rainfall intensity–kinetic energy relationships: A critical literature appraisal. Journal of Hydrology. 2002; 261(1–4):1-23

Ebro Basin (NE Spain). Journal of Hydrology. 2009;379(1):111-121

DOI: http://dx.doi.org/10.5772/intechopen.85937

[27] Lujan DL, Gabriels D. Assessing the rain erosivity and rain distribution in different agro-climatological zones in Venezuela. Sociedade & Natureza.

[28] Sauerborn P, Klein A, Botschek J, Skowronek A. Future rainfall erosivity derived from large-scale climate models—Methods and scenarios for a humid region. Geoderma. 1999;93(3):

[29] Lastoria B, Miserocchi F, Lanciani A, Monacelli G. An estimated erosion map for the Aterno-Pescara River basin. European Water. 2008;21(22):29-39

[31] Smithen A, Schulze R. The spatial distribution in Southern Africa of rainfall erosivity for use in the universal soil loss equation. Water SA. 1982;8(2):

[32] Arnoldus H. Methodology used to determine the maximum potential average annual soil loss due to sheet and rill erosion in Morocco. FAO Soils Bulletins (FAO). 1977;34:8-9

[33] Igwe C, Mbagwu J, Akamigbo F. Application of SLEMSA and USLE models for potential erosion hazard mapping in south-eastern Nigeria. International Agrophysics. 1999;13:41-48

[34] Ferro V, Porto P, Yu B. A

comparative study of rainfall erosivity estimation for southern Italy and Southeastern Australia. Hydrological Sciences Journal. 1999;44(1):3-24

[35] Munka C, Cruz G, Caffera RM. Long term variation in rainfall erosivity in Uruguay: A preliminary Fournier approach. GeoJournal. 2007;70(4):

[30] Sepaskhah A, Sarkhosh P. Estimating storm erosion index in southern region of Iran. Iranian Journal of Science and Technology, Transaction B: Engineering. 2005;29(3):357-363

2005;1(1):16-29

269-276

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods

74-78

257-262

[18] Diodato N. Estimating RUSLE's rainfall factor in the part of Italy with a

Mediterranean rainfall regime. Hydrology and Earth System Sciences Discussions. 2004;8(1):103-107

[19] Zhang G, Nearing M, Liu B. Potential effects of climate change on rainfall erosivity in the Yellow River Basin of China. Transactions of the

[20] Lee J-H, Heo J-H. Evaluation of estimation methods for rainfall erosivity based on annual precipitation in Korea. Journal of Hydrology. 2011;409(1):30-48

[21] Millward AA, Mersey JE. Adapting the RUSLE to model soil erosion potential in a mountainous tropical watershed. Catena. 1999;38(2):109-129

[22] Mikhailova EA, Bryant R, Schwager S, Smith S. Predicting rainfall erosivity in Honduras. Soil Science Society of America Journal. 1997;61(1):273-279

[23] Stocking M, Elwell H. Rainfall erosivity over Rhodesia. Transactions of the Institute of British Geographers.

[24] Lo A, El-Swaify S, Dangler E, Shinshiro L. Effectiveness of EI30 as an

[25] Panagos P, Ballabio C, Borrelli P, Meusburger K, Klik A, Rousseva S, et al. Rainfall erosivity in Europe. Science of the Total Environment. 2015;511:

[26] Yang D, Kanae S, Oki T, Koike T, Musiake K. Global potential soil erosion with reference to land use and climate changes. Hydrological Processes. 2003;

erosivity index in Hawaii. Soil Conservation Society of America.

1976;2:231-245

1985;1:384-339

17(14):2913-2928

801-814

31

ASAE. 2005;48(2):511-517

[10] USDA-ARS. Science Documentation, Revised Universal Soil Loss Equation Version 2 (RUSLE2). USDA-Agricultural Research Service; 2013

[11] Yin S, Nearing MA, Borrelli P, Xue X. Rainfall Erosivity: An overview of methodologies and applications. Vadose Zone Journal. 2017;16(12):0

[12] Vantas K, Sidiropoulos E, Loukas A. Temporal and elevation trend detection of rainfall erosivity density in Greece. In: 3rd International Electronic Conference onWater Sciences (ECWS-3); 2018

[13] Roose E. Erosion et ruissellement en Afrique de l'Ouest, vingt années de mesures en parcelles expérimentales. O. R.S.T.R.O.M. 1977;78:30-33

[14] Meusburger K, Steel A, Panagos P, Montanarella L, Alewell C. Spatial and temporal variability of rainfall erosivity factor for Switzerland. Hydrology and Earth System Sciences. 2012; 16(1):167-177

[15] Renard KG, Freimund JR. Using monthly precipitation data to estimate the R-factor in the revised USLE. Journal of Hydrology. 1994;157:287-306

[16] Sudhishri S, Patnaik U. Erosion index analysis for Eastern Ghat high zone of Orissa. Indian Journal of Dryland Agricultural Research and Development. 2004;19:42-47

[17] Angulo-Martínez M, Beguería S. Estimating rainfall erosivity from daily precipitation records: A comparison among methods using data from the

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods DOI: http://dx.doi.org/10.5772/intechopen.85937

Ebro Basin (NE Spain). Journal of Hydrology. 2009;379(1):111-121

References

pp. 477-487

357-362

[1] Boardman J, Poesen J. Soil erosion in Europe: Major processes, causes and consequences. In: Soil Erosion in Europe. Wiley-Blackwell; 2006.

Soil Erosion - Rainfall Erosivity and Risk Assessment

[9] Van Dijk A, Bruijnzeel LA, Rosewell CJ. Rainfall intensity–kinetic energy relationships: A critical literature appraisal. Journal of Hydrology. 2002;

[10] USDA-ARS. Science Documentation, Revised Universal Soil Loss Equation Version 2 (RUSLE2). USDA-Agricultural

[11] Yin S, Nearing MA, Borrelli P, Xue X. Rainfall Erosivity: An overview of methodologies and applications. Vadose

[12] Vantas K, Sidiropoulos E, Loukas A. Temporal and elevation trend detection of rainfall erosivity density in Greece. In: 3rd International Electronic

Conference onWater Sciences (ECWS-3);

[13] Roose E. Erosion et ruissellement en Afrique de l'Ouest, vingt années de mesures en parcelles expérimentales. O.

[14] Meusburger K, Steel A, Panagos P, Montanarella L, Alewell C. Spatial and temporal variability of rainfall erosivity factor for Switzerland. Hydrology and Earth System Sciences. 2012;

[15] Renard KG, Freimund JR. Using monthly precipitation data to estimate the R-factor in the revised USLE. Journal of Hydrology. 1994;157:287-306

[16] Sudhishri S, Patnaik U. Erosion index analysis for Eastern Ghat high zone of Orissa. Indian Journal of Dryland Agricultural Research and Development. 2004;19:42-47

[17] Angulo-Martínez M, Beguería S. Estimating rainfall erosivity from daily precipitation records: A comparison among methods using data from the

261(1–4):1-23

2018

Research Service; 2013

Zone Journal. 2017;16(12):0

R.S.T.R.O.M. 1977;78:30-33

16(1):167-177

[2] Kinnell PIA. Event soil loss, runoff and the universal soil loss equation family of models: A review. Journal of Hydrology. 2010;385(1):384-397

[3] Nearing MA, Yin S, Borrelli P, Polyakov VO. Rainfall erosivity: An historical review. Catena. 2017;157:

[4] Wischmeier WH, Smith DD. Predicting Rainfall Erosion Losses— A Guide to Conservation Planning. USDA, Science and Education

Erosion by Water: A Guide to Conservation Planning with the Revised Universal Soil Loss Equation

(RUSLE). Vol. Agriculture

United States Department of

Agriculture; 1997

30(2):379-386

usefulness in predicting

1976. pp. 135-142

30

[5] Renard KG, Foster GR, Weesies GA, McCool DK, Yoder DC. Predicting Soil

Handbook No. 703. Washington, DC:

[6] Brown L, Foster G. Storm erosivity using idealized intensity distributions. Transactions of the ASAE. 1987;

[7] Ramon R, Minella JP, Merten GH, de Barros CA, Canale T. Kinetic energy estimation by rainfall intensity and its

hydrosedimentological variables in a small rural catchment in southern Brazil. Catena. 2017;148:176-184

[8] McGregor KC, Mutchler CK. Status of the R factor in northern Mississippi. In: Soil Erosion: Prediction and Control. Soil Conservation Society of America;

Administration; 1986

[18] Diodato N. Estimating RUSLE's rainfall factor in the part of Italy with a Mediterranean rainfall regime. Hydrology and Earth System Sciences Discussions. 2004;8(1):103-107

[19] Zhang G, Nearing M, Liu B. Potential effects of climate change on rainfall erosivity in the Yellow River Basin of China. Transactions of the ASAE. 2005;48(2):511-517

[20] Lee J-H, Heo J-H. Evaluation of estimation methods for rainfall erosivity based on annual precipitation in Korea. Journal of Hydrology. 2011;409(1):30-48

[21] Millward AA, Mersey JE. Adapting the RUSLE to model soil erosion potential in a mountainous tropical watershed. Catena. 1999;38(2):109-129

[22] Mikhailova EA, Bryant R, Schwager S, Smith S. Predicting rainfall erosivity in Honduras. Soil Science Society of America Journal. 1997;61(1):273-279

[23] Stocking M, Elwell H. Rainfall erosivity over Rhodesia. Transactions of the Institute of British Geographers. 1976;2:231-245

[24] Lo A, El-Swaify S, Dangler E, Shinshiro L. Effectiveness of EI30 as an erosivity index in Hawaii. Soil Conservation Society of America. 1985;1:384-339

[25] Panagos P, Ballabio C, Borrelli P, Meusburger K, Klik A, Rousseva S, et al. Rainfall erosivity in Europe. Science of the Total Environment. 2015;511: 801-814

[26] Yang D, Kanae S, Oki T, Koike T, Musiake K. Global potential soil erosion with reference to land use and climate changes. Hydrological Processes. 2003; 17(14):2913-2928

[27] Lujan DL, Gabriels D. Assessing the rain erosivity and rain distribution in different agro-climatological zones in Venezuela. Sociedade & Natureza. 2005;1(1):16-29

[28] Sauerborn P, Klein A, Botschek J, Skowronek A. Future rainfall erosivity derived from large-scale climate models—Methods and scenarios for a humid region. Geoderma. 1999;93(3): 269-276

[29] Lastoria B, Miserocchi F, Lanciani A, Monacelli G. An estimated erosion map for the Aterno-Pescara River basin. European Water. 2008;21(22):29-39

[30] Sepaskhah A, Sarkhosh P. Estimating storm erosion index in southern region of Iran. Iranian Journal of Science and Technology, Transaction B: Engineering. 2005;29(3):357-363

[31] Smithen A, Schulze R. The spatial distribution in Southern Africa of rainfall erosivity for use in the universal soil loss equation. Water SA. 1982;8(2): 74-78

[32] Arnoldus H. Methodology used to determine the maximum potential average annual soil loss due to sheet and rill erosion in Morocco. FAO Soils Bulletins (FAO). 1977;34:8-9

[33] Igwe C, Mbagwu J, Akamigbo F. Application of SLEMSA and USLE models for potential erosion hazard mapping in south-eastern Nigeria. International Agrophysics. 1999;13:41-48

[34] Ferro V, Porto P, Yu B. A comparative study of rainfall erosivity estimation for southern Italy and Southeastern Australia. Hydrological Sciences Journal. 1999;44(1):3-24

[35] Munka C, Cruz G, Caffera RM. Long term variation in rainfall erosivity in Uruguay: A preliminary Fournier approach. GeoJournal. 2007;70(4): 257-262

[36] Elagib N. Changing rainfall, seasonality and erosivity in the hyperarid zone of Sudan. Land Degradation & Development. 2011;22(6):505-512

[37] Richardson C, Foster G, Wright D. Estimation of erosion index from daily rainfall amount. Transactions of the ASAE. 1983;26(1):153-156

[38] Yu B, Rosewell C. An assessment of a daily rainfall erosivity model for New South Wales. Soil Research. 1996;34(1): 139-152

[39] Angima S, Stott D, O'Neill M, Ong C, Weesies G. Soil erosion prediction using RUSLE for central Kenyan highland conditions. Agriculture, Ecosystems & Environment. 2003; 97(1):295-308

[40] Xie Y, Yin S, Liu B, Nearing MA, Zhao Y. Models for estimating daily rainfall erosivity in China. Journal of Hydrology. 2016;535:547-558

[41] Salako F. Development of isoerodent maps for Nigeria from daily rainfall amount. Geoderma. 2010; 156(3):372-378

[42] Le Roux J, Morgenthal T, Malherbe J, Pretorius D, Sumner P. Water erosion prediction at a national scale for South Africa. Water SA. 2008;34(3):305-314

[43] Capolongo D, Diodato N, Mannaerts CM, Piccarreta M, Strobl R. Analyzing temporal changes in climate erosivity using a simplified rainfall erosivity model in Basilicata (Southern Italy). Journal of Hydrology. 2008; 356(1):119-130

[44] Elsenbeer H, Cassel D, Tinner W. A daily rainfall erosivity model for Western Amazonia. Journal of Soil and Water Conservation. 1993;48(5): 439-444

[45] Petkovšek G, Mikoš M. Estimating the R factor from daily rainfall data in

the sub-Mediterranean climate of Southwest Slovenia. Hydrological Sciences Journal. 2004;49(5)

Proceedings of the Fourteenth

Lauderdale; 2011. pp. 315-323

1994;6:303-310

International Conference on Artificial Intelligence and Statistics; FL, USA: Fort

DOI: http://dx.doi.org/10.5772/intechopen.85937

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods

[55] Wang C, Venkatesh SS, Judd JS. Optimal stopping and effective machine complexity in learning. Advances in Neural Information Processing Systems.

[56] R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2018

[57] Vantas K. Hydroscoper: R interface to the Greek National Data Bank for hydrological and meteorological information. Journal of Open Source

[58] Vantas K. hyetor: R package to analyze fixed interval precipitation time series. 2018. Available from: https://

[59] Wickham H. ggplot2: Elegant Graphics for Data Analysis. New York:

[60] Abadi M, Barham K, Chen J, Chen Z, Davis A, Dean J, et al. Tensorflow: A system for large-scale machine learning.

[61] Haan CT, Barfield BJ, Hayes JC. Design Hydrology and Sedimentology for Small Catchments. Academic Press;

Springer-Verlag; 2009. p. 213

In: OSDI. 2016. pp. 265-283

1990. pp. 238-246

33

Software. 2018;3(23):625

kvantas.github.io/hyetor/

[46] Hollinger SE, Angel JR, Palecki MA. Spatial Distribution, Variation, and Trends in Storm Precipitation Characteristics Associated with Soil Erosion in the United States. Prepared for the United State Department of Agriculture; 2002. p. 103

[47] Vantas K, Sidiropoulos E. Imputation of erosivity values under incomplete rainfall data by machine learning methods. European Water. 2017;57:193-197

[48] Vafiadis M, Tolikas D, Koutsoyiannis D. HYDROSCOPE: The new Greek national database system for meteorological, hydrological and hydrogeological information. In: WIT Transactions on Ecology and the Environment. Vol.3. WIT Press; 1994. pp. 1-8

[49] Cleveland WS, Grosse E, Shyu WM. Local regression models. Statistical models in S. 1992;2:309-376

[50] Coleman T, Li Y. On the Convergence of Reflective Newton Methods for Large-Scale Nonlinear Minimization Subject to Bounds. Vol. 67. Ithaca, NY, USA: Cornell University; 1994

[51] Coleman TF, Li Y. An interior trust region approach for nonlinear minimization subject to bounds. SIAM Journal on Optimization. 1996;6(2): 418-445

[52] Schmidhuber J. Deep learning in neural networks: An overview. Neural Networks. 2015;61:85-117

[53] Goodfellow I, Bengio Y, Courville A, Bengio Y. Deep Learning. Vol. 1. Cambridge: MIT Press; 2016

[54] Glorot X, Bordes A, Bengio Y. Deep sparse rectifier neural networks. In:

Rainfall Erosivity and Its Estimation: Conventional and Machine Learning Methods DOI: http://dx.doi.org/10.5772/intechopen.85937

Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics; FL, USA: Fort Lauderdale; 2011. pp. 315-323

[36] Elagib N. Changing rainfall, seasonality and erosivity in the hyperarid zone of Sudan. Land Degradation & Development. 2011;22(6):505-512

ASAE. 1983;26(1):153-156

139-152

97(1):295-308

156(3):372-378

356(1):119-130

439-444

32

[37] Richardson C, Foster G, Wright D. Estimation of erosion index from daily rainfall amount. Transactions of the

Soil Erosion - Rainfall Erosivity and Risk Assessment

the sub-Mediterranean climate of Southwest Slovenia. Hydrological Sciences Journal. 2004;49(5)

Agriculture; 2002. p. 103

2017;57:193-197

pp. 1-8

University; 1994

418-445

[47] Vantas K, Sidiropoulos E. Imputation of erosivity values under incomplete rainfall data by machine learning methods. European Water.

[48] Vafiadis M, Tolikas D,

Koutsoyiannis D. HYDROSCOPE: The new Greek national database system for meteorological, hydrological and hydrogeological information. In: WIT Transactions on Ecology and the Environment. Vol.3. WIT Press; 1994.

[49] Cleveland WS, Grosse E, Shyu WM. Local regression models. Statistical models in S. 1992;2:309-376

[50] Coleman T, Li Y. On the Convergence of Reflective Newton Methods for Large-Scale Nonlinear Minimization Subject to Bounds. Vol. 67. Ithaca, NY, USA: Cornell

[51] Coleman TF, Li Y. An interior trust

minimization subject to bounds. SIAM Journal on Optimization. 1996;6(2):

[52] Schmidhuber J. Deep learning in neural networks: An overview. Neural

[53] Goodfellow I, Bengio Y, Courville A,

[54] Glorot X, Bordes A, Bengio Y. Deep sparse rectifier neural networks. In:

Bengio Y. Deep Learning. Vol. 1. Cambridge: MIT Press; 2016

region approach for nonlinear

Networks. 2015;61:85-117

[46] Hollinger SE, Angel JR, Palecki MA. Spatial Distribution, Variation, and Trends in Storm Precipitation Characteristics Associated with Soil Erosion in the United States. Prepared for the United State Department of

[38] Yu B, Rosewell C. An assessment of a daily rainfall erosivity model for New South Wales. Soil Research. 1996;34(1):

[39] Angima S, Stott D, O'Neill M, Ong C, Weesies G. Soil erosion prediction using RUSLE for central Kenyan highland conditions. Agriculture, Ecosystems & Environment. 2003;

[40] Xie Y, Yin S, Liu B, Nearing MA, Zhao Y. Models for estimating daily rainfall erosivity in China. Journal of

isoerodent maps for Nigeria from daily rainfall amount. Geoderma. 2010;

[42] Le Roux J, Morgenthal T, Malherbe J, Pretorius D, Sumner P. Water erosion prediction at a national scale for South Africa. Water SA. 2008;34(3):305-314

Mannaerts CM, Piccarreta M, Strobl R. Analyzing temporal changes in climate erosivity using a simplified rainfall erosivity model in Basilicata (Southern Italy). Journal of Hydrology. 2008;

[44] Elsenbeer H, Cassel D, Tinner W. A

[45] Petkovšek G, Mikoš M. Estimating the R factor from daily rainfall data in

daily rainfall erosivity model for Western Amazonia. Journal of Soil and Water Conservation. 1993;48(5):

Hydrology. 2016;535:547-558

[41] Salako F. Development of

[43] Capolongo D, Diodato N,

[55] Wang C, Venkatesh SS, Judd JS. Optimal stopping and effective machine complexity in learning. Advances in Neural Information Processing Systems. 1994;6:303-310

[56] R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2018

[57] Vantas K. Hydroscoper: R interface to the Greek National Data Bank for hydrological and meteorological information. Journal of Open Source Software. 2018;3(23):625

[58] Vantas K. hyetor: R package to analyze fixed interval precipitation time series. 2018. Available from: https:// kvantas.github.io/hyetor/

[59] Wickham H. ggplot2: Elegant Graphics for Data Analysis. New York: Springer-Verlag; 2009. p. 213

[60] Abadi M, Barham K, Chen J, Chen Z, Davis A, Dean J, et al. Tensorflow: A system for large-scale machine learning. In: OSDI. 2016. pp. 265-283

[61] Haan CT, Barfield BJ, Hayes JC. Design Hydrology and Sedimentology for Small Catchments. Academic Press; 1990. pp. 238-246

Chapter 3

Abstract

Simulation of Surface Runoff and

Numerical simulation of surface runoff is used to understand and predict watershed sediment transport and water quality and improve management of agricultural watersheds. However, models currently available are either simplified or parameterized for efficiency. In this chapter, CCHE2D, a physically based hydrodynamic model for general free surface flow hydrodynamics, was applied to study watershed surface runoff and channel flows. Multiple analytical solutions and experimental data were used to verify and validate this finite element model systematically with good results. A numerical scheme for correcting the bilinear inter-

polation of the water surface elevation solutions from the cell centers to the computational nodes was developed to improve the model. The correction was found necessary and effective for the sheet runoff simulations over the irregular bed topography. The modified numerical model was then used to simulate storms in a low-relief agricultural watershed in the Mississippi River alluvial plain. This physically based model identified the channel networks, watershed boundary automatically, and helped to develop rating curves at the gage station of this complex watershed. The numerical simulations resolved detailed runoff and turbulent channel flows, which can be used for soil erosion and gully development analyses.

Keywords: overland flow, rainfall runoff, numerical modeling, physically based

Numerical simulation is increasingly used for studying overland flows. Since runoff drives soil erosion and landscape evolution, the runoff models provide a foundation for modeling soil erosion, rill erosion, and related processes at the watershed scale [1, 2]. Models involving different levels of abstraction have been proposed [3–5]. Two commonly used models are the diffusion wave (DW) and kinematic wave (KW) models [6–9]. The KW models set the friction slope to be equal to the bed slope and ignore the inertial terms [10]. The method has been successfully used to describe overland flows [11–14]. The governing equations are highly nonlinear and do not have general analytical solutions, so one has to solve them numerically for practical cases [15]. The models based on full Saint-Venant

(SV) equations have also been applied and produced better results.

Channel Flows Using a 2D

Yafei Jia,Tahmina Shirmeen, Martin A. Locke, Richard E. Lizotte Jr. and F. Douglas Shields Jr.

Numerical Model

model, model validation, soil erosion

1. Introduction

35

## Chapter 3

## Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model

Yafei Jia,Tahmina Shirmeen, Martin A. Locke, Richard E. Lizotte Jr. and F. Douglas Shields Jr.

## Abstract

Numerical simulation of surface runoff is used to understand and predict watershed sediment transport and water quality and improve management of agricultural watersheds. However, models currently available are either simplified or parameterized for efficiency. In this chapter, CCHE2D, a physically based hydrodynamic model for general free surface flow hydrodynamics, was applied to study watershed surface runoff and channel flows. Multiple analytical solutions and experimental data were used to verify and validate this finite element model systematically with good results. A numerical scheme for correcting the bilinear interpolation of the water surface elevation solutions from the cell centers to the computational nodes was developed to improve the model. The correction was found necessary and effective for the sheet runoff simulations over the irregular bed topography. The modified numerical model was then used to simulate storms in a low-relief agricultural watershed in the Mississippi River alluvial plain. This physically based model identified the channel networks, watershed boundary automatically, and helped to develop rating curves at the gage station of this complex watershed. The numerical simulations resolved detailed runoff and turbulent channel flows, which can be used for soil erosion and gully development analyses.

Keywords: overland flow, rainfall runoff, numerical modeling, physically based model, model validation, soil erosion

#### 1. Introduction

Numerical simulation is increasingly used for studying overland flows. Since runoff drives soil erosion and landscape evolution, the runoff models provide a foundation for modeling soil erosion, rill erosion, and related processes at the watershed scale [1, 2]. Models involving different levels of abstraction have been proposed [3–5]. Two commonly used models are the diffusion wave (DW) and kinematic wave (KW) models [6–9]. The KW models set the friction slope to be equal to the bed slope and ignore the inertial terms [10]. The method has been successfully used to describe overland flows [11–14]. The governing equations are highly nonlinear and do not have general analytical solutions, so one has to solve them numerically for practical cases [15]. The models based on full Saint-Venant (SV) equations have also been applied and produced better results.

Two-dimension models are generally used for cases with irregular domains. A distributed rainfall-runoff model using the KW approximation solved by an implicit finite difference scheme was developed [16], but channel flows are computed using a separate KW model. Fully two-dimensional shallow water equations are being utilized for modeling overland flows in late 1980s [17]. A twodimensional finite difference (FD) runoff model was developed by solving 2D SV equations [18]. Shallow water equation-based 2D models [19] were used for runoff over an irregular topography of experimental scale with infiltration processes considered and in rural semiarid watersheds for overland flows generated by storms [20].

and applied widely to simulations of channel flow, flooding, coastal flow, bed topographic change, and chemical contamination in aquatic environments [36–40]. The major objectives of the present paper are to assess the accuracy and the effectiveness of this FEM in predicting overland runoff processes, and its applicability to practical agricultural watersheds with ditches and natural stream channels. The approach of the study followed the recommendations of [41] for quality assurance that numerical models have to be verified and validated using analytical solutions, physical experimental data, and field data. The validated numerical model was used to simulate and characterize the hydrological processes of an agricultural watershed in the Mississippi River alluvial plain where farm fields are drained and separated by ditches and stream channels. A limitation was found in the interpolation method when it is applied to the water surface elevation of the sheet runoff. A numerical scheme was developed and implemented for improving the bilinear interpolation. The present study focused on watershed surface flow processes over bare soils; interception, evapotranspiration, and infiltration were not

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model

DOI: http://dx.doi.org/10.5772/intechopen.80214

Surface runoff due to precipitation is typically quite shallow and can be aptly represented by the 2D shallow water equations within the CCHE2D model [36, 38]. The water surface elevation of the runoff flow, η, is calculated by the continuity

in which h is the local water depth, t is time; R is rainfall intensity, which may vary in time and space, and u and v are depth-averaged velocity components in x and y directions, respectively. The depth-integrated 2D momentum equations for

> ∂x þ

∂y þ

in which g is the gravitational acceleration, ρ is water density, τxx, τxy, τyx, and τyy

The full Eqs. (1)–(3) are applicable for general flow conditions. In realistic cases where runoff and channel flow conditions coexist, a general flow model is necessary. Under the sheet flow condition, the advection and turbulence stress terms in

∂hτxx ∂x þ ∂hτxy ∂y

∂hτyx ∂x þ ∂hτyy ∂y

<sup>∂</sup><sup>y</sup> ¼ �gh <sup>∂</sup><sup>η</sup>

<sup>∂</sup><sup>y</sup> ¼ �gh <sup>∂</sup><sup>η</sup>

ology and techniques can be found in earlier publications [36, 38, 42].

are depth-averaged Reynolds stresses, and τbx, τby are bed shear stresses. In the overland runoff area, the Reynolds stress terms vanish, and Eqs. (2) and (3) become the shallow water equations. The Reynolds stress terms remain significant in the part of the domain with channel and concentrated flows. A special finite element method called the efficient element method is adopted in the model, in which a collocation approach is used to discretize the equations in a structured quadrilateral nonorthogonal mesh system. A partially staggered grid is used for solving these equations. A velocity correction method is used to couple the continuity equation and the momentum equations. More details about this model's numerical method-

<sup>∂</sup><sup>y</sup> <sup>¼</sup> <sup>R</sup> (1)

� <sup>τ</sup>bx ρ

� τby ρ

(2)

(3)

∂η ∂t þ ∂uh ∂x þ ∂vh

considered.

37

2. Mathematical model

turbulent flows are as follows:

∂uh ∂t þ ∂uuh ∂x þ ∂vuh

∂vh ∂t þ ∂uvh ∂x þ ∂uvh

equation in a Cartesian coordinate system

In addition to finite difference method (FDM), the two-dimensional finite element (FEM) and finite volume methods (FVM) have been used for overland flow simulations. A FEM KW model was developed by Liu et al. [21] for simulating runoff generation and concentration over an irregular bed and reproduced experimental results. Tests [15] indicated that the FVM-based 2D SV model performed better than that of FDM. Costabile et al. [22] solved the shallow water equations using the FVM and applied the resulting model to simulate a real event on a watershed of 40 km<sup>2</sup> . Nunoz-Carpena et al. [23] solved the KW equation using the Petrov-Galerkin method. Venkata et al. [24] developed a Galerkin DW FEM and applied it to a small watershed. Singh et al. [25] simulated runoff processes by solving the 2D shallow water equations with a shock-capturing scheme and the FVM. Shirmeen et al. [26] showed results of a validated, FEM 2D model in predicting runoff from a flat agricultural watershed.

In order to check numerical models' mathematical correctness and physical applicability, the developed computational models have been tested with analytical solutions, experimental, and field data. Iwagaki [27] studied runoff using analytical methods and experimental data; several specific solutions were developed based on the characteristic method. Govindaraju et al. [28] developed analytical solutions using KW and DW approximations. Comparisons of analytical solutions, numerical solutions, and experimental data were discussed. Singh [29] detailed the KW model's analytical and numerical solutions and their wide applications. Cea [30] tested FVM using an experimental watershed with a complex shape. These overland flow models use simplified equations and need to specify pre-existing channel networks, which make it difficult to simulate soil erosion cases with hill-slope evaluation and mixed sheet-channel flow conditions.

CCHE2D is a physically based model, which treats the entire watershed including the channels and ditches as one continuous domain. One does not need to differentiate overland sheet flow and channel flow calculation areas using grid cells and 1D channel networks as is done in GSSHA [31], WASH123D [32], NIKE-SHE [33], and SHETRAN [34]. It is also not necessary to employ arbitrarily shaped subwatersheds and 1D channel networks as is done in the CCHE1D model [35]. In these models, 2D DW equations or KW equations are solved for the overland flow using finite difference methods, and the 1D SV equation is solved in the prescribed channel networks. In contrast to these models, in CCHE2D, hydrodynamics over the entire watershed is simulated using only 2D equations discretized on an irregular quadrilateral finite element mesh, which is generated using digital elevation model (DEM) data. The simulated overland sheet flow and channel flow are seamlessly connected everywhere in the domain and the channel network is formed automatically. This method may be more applicable when sediment transport, rill erosion, or gully erosion processes in watersheds are considered.

In this study, the CCHE2D model is modified and applied to simulate watershed hydrological processes. CCHE2D is a general hydrodynamic model for unsteady, turbulent free flows, sediment transport, and pollutant transport. It has been validated

#### Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model DOI: http://dx.doi.org/10.5772/intechopen.80214

and applied widely to simulations of channel flow, flooding, coastal flow, bed topographic change, and chemical contamination in aquatic environments [36–40].

The major objectives of the present paper are to assess the accuracy and the effectiveness of this FEM in predicting overland runoff processes, and its applicability to practical agricultural watersheds with ditches and natural stream channels. The approach of the study followed the recommendations of [41] for quality assurance that numerical models have to be verified and validated using analytical solutions, physical experimental data, and field data. The validated numerical model was used to simulate and characterize the hydrological processes of an agricultural watershed in the Mississippi River alluvial plain where farm fields are drained and separated by ditches and stream channels. A limitation was found in the interpolation method when it is applied to the water surface elevation of the sheet runoff. A numerical scheme was developed and implemented for improving the bilinear interpolation. The present study focused on watershed surface flow processes over bare soils; interception, evapotranspiration, and infiltration were not considered.

#### 2. Mathematical model

Two-dimension models are generally used for cases with irregular domains. A distributed rainfall-runoff model using the KW approximation solved by an implicit finite difference scheme was developed [16], but channel flows are computed using a separate KW model. Fully two-dimensional shallow water equations are being utilized for modeling overland flows in late 1980s [17]. A twodimensional finite difference (FD) runoff model was developed by solving 2D SV equations [18]. Shallow water equation-based 2D models [19] were used for runoff over an irregular topography of experimental scale with infiltration processes considered and in rural semiarid watersheds for overland flows generated by

In addition to finite difference method (FDM), the two-dimensional finite element (FEM) and finite volume methods (FVM) have been used for overland flow simulations. A FEM KW model was developed by Liu et al. [21] for simulating runoff generation and concentration over an irregular bed and reproduced experimental results. Tests [15] indicated that the FVM-based 2D SV model performed better than that of FDM. Costabile et al. [22] solved the shallow water equations using the FVM and applied the resulting model to simulate a real event on a

Petrov-Galerkin method. Venkata et al. [24] developed a Galerkin DW FEM and applied it to a small watershed. Singh et al. [25] simulated runoff processes by solving the 2D shallow water equations with a shock-capturing scheme and the FVM. Shirmeen et al. [26] showed results of a validated, FEM 2D model in

In order to check numerical models' mathematical correctness and physical applicability, the developed computational models have been tested with analytical solutions, experimental, and field data. Iwagaki [27] studied runoff using analytical methods and experimental data; several specific solutions were developed based on the characteristic method. Govindaraju et al. [28] developed analytical solutions using KW and DW approximations. Comparisons of analytical solutions, numerical solutions, and experimental data were discussed. Singh [29] detailed the KW model's analytical and numerical solutions and their wide applications. Cea [30] tested FVM using an experimental watershed with a complex shape. These overland flow models use simplified equations and need to specify pre-existing channel networks, which make it difficult to simulate soil erosion cases with hill-slope

CCHE2D is a physically based model, which treats the entire watershed includ-

In this study, the CCHE2D model is modified and applied to simulate watershed hydrological processes. CCHE2D is a general hydrodynamic model for unsteady, turbulent free flows, sediment transport, and pollutant transport. It has been validated

ing the channels and ditches as one continuous domain. One does not need to differentiate overland sheet flow and channel flow calculation areas using grid cells and 1D channel networks as is done in GSSHA [31], WASH123D [32], NIKE-SHE [33], and SHETRAN [34]. It is also not necessary to employ arbitrarily shaped subwatersheds and 1D channel networks as is done in the CCHE1D model [35]. In these models, 2D DW equations or KW equations are solved for the overland flow using finite difference methods, and the 1D SV equation is solved in the prescribed channel networks. In contrast to these models, in CCHE2D, hydrodynamics over the entire watershed is simulated using only 2D equations discretized on an irregular quadrilateral finite element mesh, which is generated using digital elevation model (DEM) data. The simulated overland sheet flow and channel flow are seamlessly connected everywhere in the domain and the channel network is formed automatically. This method may be more applicable when sediment transport, rill

erosion, or gully erosion processes in watersheds are considered.

predicting runoff from a flat agricultural watershed.

Soil Erosion - Rainfall Erosivity and Risk Assessment

evaluation and mixed sheet-channel flow conditions.

. Nunoz-Carpena et al. [23] solved the KW equation using the

storms [20].

watershed of 40 km<sup>2</sup>

36

Surface runoff due to precipitation is typically quite shallow and can be aptly represented by the 2D shallow water equations within the CCHE2D model [36, 38]. The water surface elevation of the runoff flow, η, is calculated by the continuity equation in a Cartesian coordinate system

$$\frac{\partial \eta}{\partial t} + \frac{\partial uh}{\partial \mathbf{x}} + \frac{\partial vh}{\partial \mathbf{y}} = R \tag{1}$$

in which h is the local water depth, t is time; R is rainfall intensity, which may vary in time and space, and u and v are depth-averaged velocity components in x and y directions, respectively. The depth-integrated 2D momentum equations for turbulent flows are as follows:

$$\frac{\partial uh}{\partial t} + \frac{\partial uuh}{\partial \mathbf{x}} + \frac{\partial vuh}{\partial \mathbf{y}} = -gh\frac{\partial \eta}{\partial \mathbf{x}} + \left(\frac{\partial h\tau\_{\mathbf{x}\mathbf{x}}}{\partial \mathbf{x}} + \frac{\partial h\tau\_{\mathbf{xy}}}{\partial \mathbf{y}}\right) - \frac{\tau\_{h\mathbf{x}}}{\rho} \tag{2}$$

$$\frac{\partial vh}{\partial t} + \frac{\partial uvh}{\partial \mathbf{x}} + \frac{\partial uvh}{\partial \mathbf{y}} = -gh\frac{\partial \eta}{\partial \mathbf{y}} + \left(\frac{\partial h\tau\_{yx}}{\partial \mathbf{x}} + \frac{\partial h\tau\_{yy}}{\partial \mathbf{y}}\right) - \frac{\tau\_{by}}{\rho} \tag{3}$$

in which g is the gravitational acceleration, ρ is water density, τxx, τxy, τyx, and τyy are depth-averaged Reynolds stresses, and τbx, τby are bed shear stresses. In the overland runoff area, the Reynolds stress terms vanish, and Eqs. (2) and (3) become the shallow water equations. The Reynolds stress terms remain significant in the part of the domain with channel and concentrated flows. A special finite element method called the efficient element method is adopted in the model, in which a collocation approach is used to discretize the equations in a structured quadrilateral nonorthogonal mesh system. A partially staggered grid is used for solving these equations. A velocity correction method is used to couple the continuity equation and the momentum equations. More details about this model's numerical methodology and techniques can be found in earlier publications [36, 38, 42].

The full Eqs. (1)–(3) are applicable for general flow conditions. In realistic cases where runoff and channel flow conditions coexist, a general flow model is necessary. Under the sheet flow condition, the advection and turbulence stress terms in

the momentum equations vanish because the dominant forcing for the overland flow is the gravity and bed shear stress. The water depth is very small, and water surface slope and bed slope become almost the same:

$$
\frac{\partial \eta}{\partial \mathbf{x}} \approx \frac{\partial b}{\partial \mathbf{x}}, \qquad \frac{\partial \eta}{\partial y} \approx \frac{\partial b}{\partial y} \tag{4}
$$

in which b is the bed elevation. The general flow equations then become the KW equations. Under this condition, the flow is completely dominated by the bed slope. Shear stresses on the bed are evaluated in conjunction with the Manning equation as:

$$
\pi\_{b\mathbf{x}} = \frac{1}{h^{1/3}} \rho \mathbf{g} n^2 u U,\tag{5}
$$

A numerical scheme has been developed and implemented in CCHE2D to correct the interpolation error [43]. Figure 2 illustrates how the scheme is formulated in one dimension with an exaggerated vertical scale. Eq. (7) is the formulation to compute the correction value Δb for nonuniform meshes, and it is simplified to Eq. (8) if the mesh is uniform. It is straightforward to extend Eqs. (7) and (8) to two dimension. Water depth at the cell centers is positive, without this correction, the depth at the middle point would become negative because the interpolated water surface elevation is below the bed. This scheme is necessary and effective when

Definition sketch for the formulation of the correction (Eqs. (7) and (8)) to linear water surface elevation interpolation. b1, b2, and b3 are bed elevation. Δb is the interpolation correction and Δx1 and Δx2 are mesh

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model

<sup>2</sup> <sup>b</sup><sup>2</sup> � <sup>b</sup>1Δx<sup>2</sup>

where b1, b2, and b3 are bed elevation, Δb is the interpolation correction, b<sup>i</sup>

Two analytical solutions were obtained by solving a one-dimensional kinematic

linearly interpolated value, and Δx1 and Δx2 are mesh spacing (Figure 2). The

equation analytically for rain-generated runoff by [44, 45]. The solution of sustained rains for the runoff to reach a steady state [44] and the solution for rainfall that stops before the runoff becomes steady [45], including the tailing stage solution after rainfall stops, were provided. The governing one-dimensional kine-

> ∂h ∂t þ u ∂h

for a few simple cases: runoff due to steady rainfall intensity on a uniform planar area of 200 � 1 m with a slope of 1.0%. The rainfall intensity was <sup>R</sup> = 2.7 � <sup>10</sup>�<sup>5</sup> m/ s, and the Manning's coefficient was n = 0.02 m�1/3s. For comparison, the same case was simulated using CCHE2D and a 10 � 100 point 2D mesh with uniform spacing. The solutions were recorded at cross sections located at 50, 100, 150, and 200 m,

<sup>u</sup> <sup>¼</sup> <sup>α</sup>h<sup>k</sup>�<sup>1</sup>

in which q is the discharge of water per unit width (m<sup>2</sup>

Δx<sup>1</sup> þ Δx<sup>2</sup>

� <sup>b</sup>3Δx<sup>1</sup> Δx<sup>1</sup> þ Δx<sup>2</sup> (7)

<sup>4</sup> ð Þ <sup>b</sup><sup>1</sup> � <sup>2</sup>b<sup>2</sup> <sup>þ</sup> <sup>b</sup><sup>3</sup> (8)

<sup>∂</sup><sup>x</sup> <sup>¼</sup> <sup>R</sup> (9)

/s). These analytical solutions were realized

/s), k is an exponent

, q <sup>¼</sup> uh <sup>¼</sup> <sup>α</sup>hk (10)

<sup>2</sup> is the

cases with irregular topography are simulated

DOI: http://dx.doi.org/10.5772/intechopen.80214

<sup>Δ</sup><sup>b</sup> <sup>¼</sup> <sup>b</sup><sup>2</sup> � bi

matic equation for deriving these solutions is:

(=5/3), and α (=5) is a coefficient (m2�<sup>k</sup>

from the upstream end of the plane.

39

4. Analytical verification

Figure 2.

spacing.

<sup>2</sup> <sup>¼</sup> <sup>1</sup>

<sup>Δ</sup><sup>b</sup> <sup>¼</sup> <sup>1</sup>

interpolated water surface elevation needs to be corrected by Δb.

$$
\pi\_{\text{by}} = \frac{1}{h^{1/3}} \rho \text{gn}^2 vU \tag{6}
$$

in which <sup>n</sup> is the Manning roughness coefficient and <sup>U</sup> <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>u</sup><sup>2</sup> <sup>þ</sup> <sup>v</sup><sup>2</sup> <sup>p</sup> is the total velocity magnitude.

#### 3. Interpolation of water surface elevation

CCHE2D uses a partially staggered method: the velocities are solved at collocation points and the pressure (water surface) is solved at cell centers [36]. A bilinear interpolation method is used to interpolate the water surface elevation solution to the collocation nodes where the momentum equations are solved. The bilinear interpolation works well for general channel flow simulations because the water depth is large in comparison with the variation of bed surface and the mesh size. When overland sheet runoff is simulated; however, the water depth is very small; it is often less than the microelevation variation of bed topography represented in an element. In this case, the interpolated water surface elevation may be lower than the bed if the bed is concave down and vice versa. This is a limitation of the interpolation method. In the concave down case, dry nodes are created artificially; in the concave up case, artificial masses of water could be erroneously created. Figure 1 illustrates this problem in one dimension. The problem occurs whenever irregular bed topographies are encountered. A correction is therefore necessary to the interpolation over the surface runoff area.

#### Figure 1.

The error of underestimation and overestimation caused by linear interpolation of water surface elevation from cell centers to collocation nodes.

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model DOI: http://dx.doi.org/10.5772/intechopen.80214

Figure 2.

the momentum equations vanish because the dominant forcing for the overland flow is the gravity and bed shear stress. The water depth is very small, and water

> <sup>∂</sup><sup>x</sup> , <sup>∂</sup><sup>η</sup> ∂y ≈∂b ∂y

in which b is the bed elevation. The general flow equations then become the KW equations. Under this condition, the flow is completely dominated by the bed slope. Shear stresses on the bed are evaluated in conjunction with the Manning equation as:

<sup>h</sup>1=<sup>3</sup> <sup>ρ</sup>gn<sup>2</sup>

<sup>h</sup><sup>1</sup>=<sup>3</sup> <sup>ρ</sup>gn<sup>2</sup>

CCHE2D uses a partially staggered method: the velocities are solved at collocation points and the pressure (water surface) is solved at cell centers [36]. A bilinear interpolation method is used to interpolate the water surface elevation solution to the collocation nodes where the momentum equations are solved. The bilinear interpolation works well for general channel flow simulations because the water depth is large in comparison with the variation of bed surface and the mesh size. When overland sheet runoff is simulated; however, the water depth is very small; it is often less than the microelevation variation of bed topography represented in an element. In this case, the interpolated water surface elevation may be lower than the bed if the bed is concave down and vice versa. This is a limitation of the interpolation method. In the concave down case, dry nodes are created artificially; in the concave up case, artificial masses of water could be erroneously created. Figure 1 illustrates this problem in one dimension. The problem occurs whenever irregular bed topographies are encountered. A correction is therefore necessary to the inter-

The error of underestimation and overestimation caused by linear interpolation of water surface elevation from

(4)

uU, (5)

vU (6)

<sup>u</sup><sup>2</sup> <sup>þ</sup> <sup>v</sup><sup>2</sup> <sup>p</sup> is the total

∂η <sup>∂</sup>x<sup>≈</sup> <sup>∂</sup><sup>b</sup>

<sup>τ</sup>bx <sup>¼</sup> <sup>1</sup>

<sup>τ</sup>by <sup>¼</sup> <sup>1</sup>

in which <sup>n</sup> is the Manning roughness coefficient and <sup>U</sup> <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

surface slope and bed slope become almost the same:

Soil Erosion - Rainfall Erosivity and Risk Assessment

3. Interpolation of water surface elevation

polation over the surface runoff area.

Figure 1.

38

cell centers to collocation nodes.

velocity magnitude.

Definition sketch for the formulation of the correction (Eqs. (7) and (8)) to linear water surface elevation interpolation. b1, b2, and b3 are bed elevation. Δb is the interpolation correction and Δx1 and Δx2 are mesh spacing.

A numerical scheme has been developed and implemented in CCHE2D to correct the interpolation error [43]. Figure 2 illustrates how the scheme is formulated in one dimension with an exaggerated vertical scale. Eq. (7) is the formulation to compute the correction value Δb for nonuniform meshes, and it is simplified to Eq. (8) if the mesh is uniform. It is straightforward to extend Eqs. (7) and (8) to two dimension. Water depth at the cell centers is positive, without this correction, the depth at the middle point would become negative because the interpolated water surface elevation is below the bed. This scheme is necessary and effective when cases with irregular topography are simulated

$$
\Delta b = b\_2 - b\_2^i = \frac{1}{2} \left[ b\_2 - \frac{b\_1 \Delta \chi\_2}{\Delta \chi\_1 + \Delta \chi\_2} - \frac{b\_3 \Delta \chi\_1}{\Delta \chi\_1 + \Delta \chi\_2} \right] \tag{7}
$$

$$
\Delta b = \frac{1}{4}(b\_1 - 2b\_2 + b\_3) \tag{8}
$$

where b1, b2, and b3 are bed elevation, Δb is the interpolation correction, b<sup>i</sup> <sup>2</sup> is the linearly interpolated value, and Δx1 and Δx2 are mesh spacing (Figure 2). The interpolated water surface elevation needs to be corrected by Δb.

#### 4. Analytical verification

Two analytical solutions were obtained by solving a one-dimensional kinematic equation analytically for rain-generated runoff by [44, 45]. The solution of sustained rains for the runoff to reach a steady state [44] and the solution for rainfall that stops before the runoff becomes steady [45], including the tailing stage solution after rainfall stops, were provided. The governing one-dimensional kinematic equation for deriving these solutions is:

$$\frac{\partial h}{\partial t} + u \frac{\partial h}{\partial \mathbf{x}} = R \tag{9}$$

$$
\mu = ah^{k-1}, \ q = \imath h = ah^k \tag{10}
$$

in which q is the discharge of water per unit width (m<sup>2</sup> /s), k is an exponent (=5/3), and α (=5) is a coefficient (m2�<sup>k</sup> /s). These analytical solutions were realized for a few simple cases: runoff due to steady rainfall intensity on a uniform planar area of 200 � 1 m with a slope of 1.0%. The rainfall intensity was <sup>R</sup> = 2.7 � <sup>10</sup>�<sup>5</sup> m/ s, and the Manning's coefficient was n = 0.02 m�1/3s. For comparison, the same case was simulated using CCHE2D and a 10 � 100 point 2D mesh with uniform spacing. The solutions were recorded at cross sections located at 50, 100, 150, and 200 m, from the upstream end of the plane.

Figure 3 shows the comparisons of the simulated runoff and the analytical solutions for the sustained rain collected at the four cross sections. Hydrographs at each cross section indicate that equilibrium runoff (steady state) is reached before the rain stops at T = 1000 s. The runoff is always nonuniform, and the peak discharge increases in the downstream direction. At first, the flow is unsteady (rising limb), then becomes steady until T = 1000 s, and finally becomes unsteady in the falling limb. The runoff reaches equilibrium earlier at locations closer to upstream. The simulation is a little less than the analytical solution at the time approaching the peak discharge, particularly near the downstream. The solution can be improved by reducing the local mesh size effectively.

5. Experimental validation

DOI: http://dx.doi.org/10.5772/intechopen.80214

discharge.

5.1 Case 1

Figure 5.

41

CCHE2D model was validated using experimental data sets collected from the literature. All of these cases were carried out on impervious overland flow planes. The only quantity measured in these experiments was the downstream runoff

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model

Morgali and Linsley [46] obtained two sets of experimental runoff data. Their tests were carried out over a straight turf surface of 21.95 m long with a constant slope (0.04) and width. The Manning's coefficient, n, was found to be 0.5 m1/3s. The rains had two different intensities and were uniform along the slope for 1200 s (20 min). Figure 5 compares the experimental data and the numerical simulations. The analytical solution for this test condition [44] is also presented in Figure 5. It was found that these runoff experiments fit well with the analytical solution. A 110 10 points uniform mesh and 0.01 s time step were used for the numerical simulation. The CCHE2D numerical results showed good agreements with the analytical solution as well as the experimental results (Figure 5). The rising limb of the discharge hydrograph and the peak discharge were captured very well by the simulations. The processes of the two experiments, 1A (R = 92.96 mm/h) and 1B (R = 48.01 mm/h), look similar because the only difference in the experiments was rainfall intensity. The peak discharges of the experiments occurred at approximately 850 and 1100 s, respectively, for Case 1A and Case 1B. The numerical solutions of CCHE2D agreed well with the experimental data. The peak discharges

times to peak discharge for Case 1A resulted from the analytical solution, CCHE2D

Comparisons of measured data with analytical solution, results of CCHE2D, and other numerical models.

/s, respectively. The

for Case 1A and Case 1B are 5.67<sup>10</sup><sup>4</sup> and 2.93 <sup>10</sup><sup>4</sup> m3

Figure 4 shows a case in which the rainfall stops before runoff reaches steady state (T = 200 s); the hydrographs, thus, have a different pattern. The peak discharge is reached at the time the rain stopped and is the same for all cross sections. The peak discharge for the lower cross sections lasts longer because the flows at the lower locations are sustained by upstream contributions. The runoff recession is earlier for upstream locations. The shape of the two sets of simulated hydrographs at all cross-section locations corresponded well with the analytical solutions.

#### Figure 3.

Comparisons of the simulated runoff hydrographs and analytical solutions. The sustained rain stopped at T = 1000 s after the steady states have reached everywhere on the slope. Comparisons at four cross sections are shown. Δx is the mesh spacing in the runoff direction.

#### Figure 4.

Comparisons of the simulated runoff hydrographs and analytical solutions. The rain stopped at T = 200 s, before the flow at any of the four cross sections reached steady state. Δx is the mesh spacing in the runoff direction.

## 5. Experimental validation

CCHE2D model was validated using experimental data sets collected from the literature. All of these cases were carried out on impervious overland flow planes. The only quantity measured in these experiments was the downstream runoff discharge.

## 5.1 Case 1

Figure 3 shows the comparisons of the simulated runoff and the analytical solutions for the sustained rain collected at the four cross sections. Hydrographs at each cross section indicate that equilibrium runoff (steady state) is reached before the rain stops at T = 1000 s. The runoff is always nonuniform, and the peak discharge increases in the downstream direction. At first, the flow is unsteady (rising limb), then becomes steady until T = 1000 s, and finally becomes unsteady in the falling limb. The runoff reaches equilibrium earlier at locations closer to upstream. The simulation is a little less than the analytical solution at the time approaching the peak discharge, particularly near the downstream. The solution can

Figure 4 shows a case in which the rainfall stops before runoff reaches steady state (T = 200 s); the hydrographs, thus, have a different pattern. The peak discharge is reached at the time the rain stopped and is the same for all cross sections. The peak discharge for the lower cross sections lasts longer because the flows at the lower locations are sustained by upstream contributions. The runoff recession is earlier for upstream locations. The shape of the two sets of simulated hydrographs at all cross-section locations corresponded well with the analytical solutions.

Comparisons of the simulated runoff hydrographs and analytical solutions. The sustained rain stopped at T = 1000 s after the steady states have reached everywhere on the slope. Comparisons at four cross sections are

Comparisons of the simulated runoff hydrographs and analytical solutions. The rain stopped at T = 200 s, before the flow at any of the four cross sections reached steady state. Δx is the mesh spacing in the runoff

be improved by reducing the local mesh size effectively.

Soil Erosion - Rainfall Erosivity and Risk Assessment

Figure 3.

Figure 4.

direction.

40

shown. Δx is the mesh spacing in the runoff direction.

Morgali and Linsley [46] obtained two sets of experimental runoff data. Their tests were carried out over a straight turf surface of 21.95 m long with a constant slope (0.04) and width. The Manning's coefficient, n, was found to be 0.5 m1/3s. The rains had two different intensities and were uniform along the slope for 1200 s (20 min). Figure 5 compares the experimental data and the numerical simulations. The analytical solution for this test condition [44] is also presented in Figure 5. It was found that these runoff experiments fit well with the analytical solution. A 110 10 points uniform mesh and 0.01 s time step were used for the numerical simulation. The CCHE2D numerical results showed good agreements with the analytical solution as well as the experimental results (Figure 5). The rising limb of the discharge hydrograph and the peak discharge were captured very well by the simulations. The processes of the two experiments, 1A (R = 92.96 mm/h) and 1B (R = 48.01 mm/h), look similar because the only difference in the experiments was rainfall intensity. The peak discharges of the experiments occurred at approximately 850 and 1100 s, respectively, for Case 1A and Case 1B. The numerical solutions of CCHE2D agreed well with the experimental data. The peak discharges for Case 1A and Case 1B are 5.67<sup>10</sup><sup>4</sup> and 2.93 <sup>10</sup><sup>4</sup> m3 /s, respectively. The times to peak discharge for Case 1A resulted from the analytical solution, CCHE2D

and the experiment, are 760, 850 and 950 s, respectively. The differences among the three are less for the Case 1B (Figure 5).

Figure 5 also compares the simulation results of CCHE2D and the model results by Govindaraju et al. [28]; the two numerical solutions agree well for the case with the higher rain (1A), but the fit of their solutions based on the SV equations does not correspond well for the case with the smaller rainfall (1B). The results of CCHE2D also outperform the analytical solution of the DW approximation [28].

### 5.2 Case 2

Cea et al. [30] conducted three runoff experiments of complex topography and simulated these cases using a 2D unstructured FVM. The experimental watershed was a rectangle (2 2.5 m) made by three planes of stainless steel, each of them with a slope of 0.05 (Figure 6). Two dikes (1.86 and 1.01 m in length) were placed in the watershed to vary the topography. Rainfall intensity, duration, and runoff hydrographs were measured. As a result, the runoff direction, distribution, and pattern of the hydrograph were affected. The runoff was accumulated and became channel flows along intercepting lines of slopes and dikes. Since both overland flow and channel flow are involved, faithful simulation requires solving full governing Eqs. (1)–(3). The rainfall applied to each test case was different. In the first test (2A), the rainfall intensity was 317 mm/h for 45 s. In the second test (2B), the rainfall intensity was 320 mm/h for 25 s; then it was stopped for 4 s and restarted for an additional 25 s with the same intensity. In the third test (2C), rainfall intensity was 328 mm/h. The rainfall was applied for 25 s; then it was stopped for 7 s and then restarted for another 25 s.

In this study, CCHE2D was applied and the numerical results were compared with experimental data. The watershed was modeled using an irregular structured mesh with the cell size ranging from 0.034 to 0.009 m; the mesh was refined near the main channel and the outlet for improving results. The Manning's roughness coefficient was set equal to 0.009 m1/3s. The simulation time was 120 s for each case. The channel flow and runoff sheet flow coexisted: the runoff from the watershed surface was accumulated in the bottom of the watershed channel with a triangle-shaped cross section formed by the side slopes. Results of cases 2A and 2C are shown in Figures 7 and 8, respectively.

slightly overestimated, the shape of the hydrograph and the peak discharge were

Distributions of simulated (a) water depth contours (b) flow (unit discharge) distribution and (c) velocity

Comparison of measured and simulated hydrographs using rainfall with two peaks (Case 2C).

Comparison of measured and simulated hydrographs using rainfall with one peak (Case 2A).

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model

DOI: http://dx.doi.org/10.5772/intechopen.80214

Figure 8 shows the comparison between the numerical and experimental runoff hydrographs of Case 2C. The shape of the hydrograph was successfully predicted. The interval between the two rainfall peaks was 7 s. The first runoff peak discharge occurred at the time the rainfall stopped, at 25 s. The runoff discharge decreased for approximately 10 s and then increased. The second runoff peak discharge occurred at approximately 57 s. The simulated processes and the observed physical

aptly predicted.

Figure 8.

Figure 9.

43

vectors at t = 54 s for test Case 2C.

Figure 7.

Figure 7 shows the comparison between the numerical solution and experimentally observed runoff hydrograph of Case 2A. The solution of the CCHE2D model agrees very well with the experimental results. The flow discharge increased continuously once the rain started. The peak discharge occurred at the time the rainfall stopped (at 45 s). Although the rising and the falling limbs of the hydrograph were

Figure 6. Topography of the experimental watershed [30].

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model DOI: http://dx.doi.org/10.5772/intechopen.80214

Figure 7.

and the experiment, are 760, 850 and 950 s, respectively. The differences among the

Figure 5 also compares the simulation results of CCHE2D and the model results by Govindaraju et al. [28]; the two numerical solutions agree well for the case with the higher rain (1A), but the fit of their solutions based on the SV equations does not correspond well for the case with the smaller rainfall (1B). The results of CCHE2D also outperform the analytical solution of the DW approximation [28].

Cea et al. [30] conducted three runoff experiments of complex topography and simulated these cases using a 2D unstructured FVM. The experimental watershed was a rectangle (2 2.5 m) made by three planes of stainless steel, each of them with a slope of 0.05 (Figure 6). Two dikes (1.86 and 1.01 m in length) were placed in the watershed to vary the topography. Rainfall intensity, duration, and runoff hydrographs were measured. As a result, the runoff direction, distribution, and pattern of the hydrograph were affected. The runoff was accumulated and became channel flows along intercepting lines of slopes and dikes. Since both overland flow and channel flow are involved, faithful simulation requires solving full governing Eqs. (1)–(3). The rainfall applied to each test case was different. In the first test (2A), the rainfall intensity was 317 mm/h for 45 s. In the second test (2B), the rainfall intensity was 320 mm/h for 25 s; then it was stopped for 4 s and restarted for an additional 25 s with the same intensity. In the third test (2C), rainfall

intensity was 328 mm/h. The rainfall was applied for 25 s; then it was stopped for 7 s

In this study, CCHE2D was applied and the numerical results were compared with experimental data. The watershed was modeled using an irregular structured mesh with the cell size ranging from 0.034 to 0.009 m; the mesh was refined near the main channel and the outlet for improving results. The Manning's roughness coefficient was set equal to 0.009 m1/3s. The simulation time was 120 s for each case. The channel flow and runoff sheet flow coexisted: the runoff from the watershed surface was accumulated in the bottom of the watershed channel with a triangle-shaped cross section formed by the side slopes. Results of cases 2A and 2C

Figure 7 shows the comparison between the numerical solution and experimentally observed runoff hydrograph of Case 2A. The solution of the CCHE2D model agrees very well with the experimental results. The flow discharge increased continuously once the rain started. The peak discharge occurred at the time the rainfall stopped (at 45 s). Although the rising and the falling limbs of the hydrograph were

three are less for the Case 1B (Figure 5).

Soil Erosion - Rainfall Erosivity and Risk Assessment

and then restarted for another 25 s.

are shown in Figures 7 and 8, respectively.

Topography of the experimental watershed [30].

5.2 Case 2

Figure 6.

42

Comparison of measured and simulated hydrographs using rainfall with one peak (Case 2A).

#### Figure 8.

Comparison of measured and simulated hydrographs using rainfall with two peaks (Case 2C).

slightly overestimated, the shape of the hydrograph and the peak discharge were aptly predicted.

Figure 8 shows the comparison between the numerical and experimental runoff hydrographs of Case 2C. The shape of the hydrograph was successfully predicted. The interval between the two rainfall peaks was 7 s. The first runoff peak discharge occurred at the time the rainfall stopped, at 25 s. The runoff discharge decreased for approximately 10 s and then increased. The second runoff peak discharge occurred at approximately 57 s. The simulated processes and the observed physical

#### Figure 9.

Distributions of simulated (a) water depth contours (b) flow (unit discharge) distribution and (c) velocity vectors at t = 54 s for test Case 2C.

processes showed a good general agreement; it also matched well with the model results of [30].

connected to the streams and ditches with drainage culverts, which often convey water from one sub-watershed to another. The locations of culverts in the study watershed were identified in a field survey and incorporated in the numerical mesh. Soil data were obtained from the Soil Survey Geographic (SSURGO) database [47]. The watershed is covered mostly by soils with high clay content, which is typical of the region [48]. Infiltration is, therefore, negligible and was not considered in the simulation. Precipitation and flow stage data were measured by field instrumentation. Because the stream instrumented with the gage station has complex conditions, it was difficult to collect reliable velocity data during a rain event. Only stage data were available. As a result, the gage station does not have a discharge-stage rating curve. Development of a rating curve using simulations and measured data for this site would be helpful for understanding the hydrologic

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model

Because the Howden Lake watershed is of low-relief, it was often difficult to determine the boundaries between sub-watersheds in field surveys or on topographic contour maps. For example, the runoff from a piece of field may flow in two directions into two sub-watersheds, and the location of the divide line might be identified only from the runoff flow distribution during a simulation. Normally, the outline of a watershed is a given condition for a hydrology study. In this study, the exact boundary outline was not firmly established even after field surveys. A larger area containing the studied watershed was simulated, and the watershed boundary and area were finally defined by the simulated runoff and channel flow patterns. The boundary outline of the studied watershed (Figure 11) contributing to the gage was identified by visually checking the simulated overland flow directions of

In the simulations, the streams and ditches between farming plots were represented using DEM elevations like flat surface areas. No channel networks were prescribed, but the simulated surface runoff flowed logically to the ditches and to the stream channel. No other watershed analysis tools were needed. Although the study results presented later are for this identified watershed, the spatial domain of numerical simulations was several times larger (Figure 10). The northern side of the stream channel had been blocked by farmers, so the overland flow from the watershed entered the stream in the middle and flowed in a southwesterly direction (Figure 11). The water from this identified watershed pasted the gage, while runoff from the region outside this watershed was discharged from the simulation domain via other ditches and streams. The area of this watershed, including cultivated land,

drainage ditches and a stream segment, was found to be 973,700 m<sup>2</sup>

the topographic elevation ranges from approximately 46.77–47.49 m in one plot and from 47.27 to 48.09 m in another. The mean slope of the fields is 0.0097 and

Several observed storm events were selected for the model application. To reduce minor losses of water due to evaporation, soil wetting and infiltration, etc., only large rain events were considered. The rainfall event in April 2011 (Table 1) was first used for simulation. Figure 12a shows the detailed ground elevation contour of a small simulation area (dashed rectangle area in Figure 11). The elevation of this area ranges from about 46.8 to about 47.4 m. Figure 12b shows the direction vectors of the runoff near the end of the simulation. Because the water is very shallow, the flow direction is highly affected by the ground topography. Figure 12c and d shows the direction vectors and water depth distribution at the

Although the variation of the bed surface topography is very small, the simulation shows how the runoff is controlled by microtopography (Figure 12a). At the peak time of the rainfall, the overall water depth in this area (Figure 12d) is much

. In this area,

processes in these watersheds.

DOI: http://dx.doi.org/10.5772/intechopen.80214

CCHE2D.

0.0098, respectively.

peak time of the rainfall.

45

Figure 9 shows the simulation results at t = 54 s (the peak of the second rainfall) for Case 2C: (a) simulated water depth contour distribution, (b) simulated flow unit discharge pattern and (c) velocity vector distribution in the watershed. The distributions indicate how the overland sheet flow, under the influence of dikes and topography, concentrates into channels and flows out of the watershed. The flows over the slopes are sheet runoff, but complex recirculations are developed in the main channel. The water surface is no longer parallel to the bed surface. These flows cannot be represented by KW, DW, and SV models.

## 6. Application to a real watershed mathematical model

This section presents the application of CCHE2D to a sub-watershed of the Howden Lake watershed, an 18 km<sup>2</sup> agricultural watershed in the Mississippi River alluvial plain (Figure 10). In this region of low relief, watersheds are configured by farm fields drained by culverts, ditches, and intermittently flowing streams called bayous. During periods between runoff events, the channels contain standing water. The studied sub-watershed was upstream of a gaging station on an intermittently flowing bayou. The average annual precipitation in this region is about 1440 mm. Precipitation occurs as intense thunderstorms or low-intensity rains associated with major frontal movements. The latter type of events may stretch over several days of drizzle and sporadic showers. During growing seasons, channels experience some flow and stage fluctuation due to irrigation withdrawals and return flows.

Watershed topography was surveyed by airborne LiDAR with a 1.5 m horizontal resolution. The vertical accuracy was 15.0 cm RMSE or better. The watershed elevation ranges from approximately 43.89–48.99 m. A nearly uniform fine mesh (mesh spacing = 3.76–4.98 m) was generated for the simulation with the ditches and small streams between the plots further refined locally. Cultivated fields are

#### Figure 10.

Location and topography of the Howden Lake watershed. Dashed curve encloses the runoff simulation area, and the dark closed curve is the gaged watershed.

#### Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model DOI: http://dx.doi.org/10.5772/intechopen.80214

connected to the streams and ditches with drainage culverts, which often convey water from one sub-watershed to another. The locations of culverts in the study watershed were identified in a field survey and incorporated in the numerical mesh.

Soil data were obtained from the Soil Survey Geographic (SSURGO) database [47]. The watershed is covered mostly by soils with high clay content, which is typical of the region [48]. Infiltration is, therefore, negligible and was not considered in the simulation. Precipitation and flow stage data were measured by field instrumentation. Because the stream instrumented with the gage station has complex conditions, it was difficult to collect reliable velocity data during a rain event. Only stage data were available. As a result, the gage station does not have a discharge-stage rating curve. Development of a rating curve using simulations and measured data for this site would be helpful for understanding the hydrologic processes in these watersheds.

Because the Howden Lake watershed is of low-relief, it was often difficult to determine the boundaries between sub-watersheds in field surveys or on topographic contour maps. For example, the runoff from a piece of field may flow in two directions into two sub-watersheds, and the location of the divide line might be identified only from the runoff flow distribution during a simulation. Normally, the outline of a watershed is a given condition for a hydrology study. In this study, the exact boundary outline was not firmly established even after field surveys. A larger area containing the studied watershed was simulated, and the watershed boundary and area were finally defined by the simulated runoff and channel flow patterns. The boundary outline of the studied watershed (Figure 11) contributing to the gage was identified by visually checking the simulated overland flow directions of CCHE2D.

In the simulations, the streams and ditches between farming plots were represented using DEM elevations like flat surface areas. No channel networks were prescribed, but the simulated surface runoff flowed logically to the ditches and to the stream channel. No other watershed analysis tools were needed. Although the study results presented later are for this identified watershed, the spatial domain of numerical simulations was several times larger (Figure 10). The northern side of the stream channel had been blocked by farmers, so the overland flow from the watershed entered the stream in the middle and flowed in a southwesterly direction (Figure 11). The water from this identified watershed pasted the gage, while runoff from the region outside this watershed was discharged from the simulation domain via other ditches and streams. The area of this watershed, including cultivated land, drainage ditches and a stream segment, was found to be 973,700 m<sup>2</sup> . In this area, the topographic elevation ranges from approximately 46.77–47.49 m in one plot and from 47.27 to 48.09 m in another. The mean slope of the fields is 0.0097 and 0.0098, respectively.

Several observed storm events were selected for the model application. To reduce minor losses of water due to evaporation, soil wetting and infiltration, etc., only large rain events were considered. The rainfall event in April 2011 (Table 1) was first used for simulation. Figure 12a shows the detailed ground elevation contour of a small simulation area (dashed rectangle area in Figure 11). The elevation of this area ranges from about 46.8 to about 47.4 m. Figure 12b shows the direction vectors of the runoff near the end of the simulation. Because the water is very shallow, the flow direction is highly affected by the ground topography. Figure 12c and d shows the direction vectors and water depth distribution at the peak time of the rainfall.

Although the variation of the bed surface topography is very small, the simulation shows how the runoff is controlled by microtopography (Figure 12a). At the peak time of the rainfall, the overall water depth in this area (Figure 12d) is much

processes showed a good general agreement; it also matched well with the model

This section presents the application of CCHE2D to a sub-watershed of the Howden Lake watershed, an 18 km<sup>2</sup> agricultural watershed in the Mississippi River alluvial plain (Figure 10). In this region of low relief, watersheds are configured by farm fields drained by culverts, ditches, and intermittently flowing streams called bayous. During periods between runoff events, the channels contain standing water. The studied sub-watershed was upstream of a gaging station on an intermittently flowing bayou. The average annual precipitation in this region is about 1440 mm. Precipitation occurs as intense thunderstorms or low-intensity rains associated with major frontal movements. The latter type of events may stretch over several days of drizzle and sporadic showers. During growing seasons, channels experience some flow and stage fluctuation due to irrigation withdrawals and

Watershed topography was surveyed by airborne LiDAR with a 1.5 m horizontal

resolution. The vertical accuracy was 15.0 cm RMSE or better. The watershed elevation ranges from approximately 43.89–48.99 m. A nearly uniform fine mesh (mesh spacing = 3.76–4.98 m) was generated for the simulation with the ditches and

small streams between the plots further refined locally. Cultivated fields are

Location and topography of the Howden Lake watershed. Dashed curve encloses the runoff simulation area, and

cannot be represented by KW, DW, and SV models.

Soil Erosion - Rainfall Erosivity and Risk Assessment

6. Application to a real watershed mathematical model

Figure 9 shows the simulation results at t = 54 s (the peak of the second rainfall) for Case 2C: (a) simulated water depth contour distribution, (b) simulated flow unit discharge pattern and (c) velocity vector distribution in the watershed. The distributions indicate how the overland sheet flow, under the influence of dikes and topography, concentrates into channels and flows out of the watershed. The flows over the slopes are sheet runoff, but complex recirculations are developed in the main channel. The water surface is no longer parallel to the bed surface. These flows

results of [30].

return flows.

Figure 10.

44

the dark closed curve is the gaged watershed.

#### Figure 11.

Numerical simulation identified watershed for the gage station. Simulation results in the dashed rectangle area are shown in Figure 12.

measurement, and therefore, water discharge was not measured. In order to better

Information and simulation results in an area indicated in Figure 11 in a dashed rectangle: (a) bed elevation contours, (b) velocity direction distribution near the end of the simulation, (c) velocity direction, and (d) the

was developed using simulated discharge, in which L is the measured water surface elevation, r and z are parameters, and L0 is the initial water surface elevation prior to a rainfall event. Eq. (11) has two unknown parameters, but there is only one relationship available for determining their values. The total volume of runoff, obtained by numerically integrating Eq. (11) in time, is equal to the rain

r Lð Þ <sup>i</sup> � <sup>L</sup><sup>0</sup> <sup>z</sup>

in which Li is the measured water surface elevation at the gage station. With

hydrograph computed using Eq. (11), and that of the numerical simulation, were

and a mean base stage L0, but the result showed unacceptable discrepancies. L0 varied due to antecedent precipitation, downstream discharge control, sedimentation, and water usage between events. The range of L0 for the studied events is 0.33 m (Table 1). Given the complexities of the hydraulic regime in the water body, varying from standing to moving state and with varying downstream controls, variable rating curve parameters are sensible. Event-specific rating curve parame-

Manning's roughness coefficient (n) is a major factor in the determination of watershed runoff characteristics and generally reflects ground cover and management. The event on April 27–28, 2011 was used for initial calibration of Manning's coefficient. The studied watershed is cultivated with soybeans (Glycine max L. Merr.), corn (Zea mays L.), cotton (Gossypium hirsutum L.), and rice (Oryza sativa L.).

Attempts were made to fit all simulated curves using a single set of values for r, z

Eq. (12) satisfied, values of r and z that best fit the shape of the discharge

<sup>Q</sup> <sup>¼</sup> r Lð Þ � <sup>L</sup><sup>0</sup> <sup>z</sup> (11)

Δt ¼ VR (12)

understand the watershed hydrology, a rating curve of the form:

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model

DOI: http://dx.doi.org/10.5772/intechopen.80214

water depth distribution at the peak of the April 2011 rainfall.

∑ i

ters are not ideal but are useful in a research context.

determined for each event by trial and error.

QiΔt ¼ ∑ i

volume, VR:

47

Figure 12.


#### Table 1.

Parameters of selected runoff events for numerical simulations.

deeper, and the flow directions (Figure 12c) are less affected by the local microtopographic features. The flow on the right side of the domain is still sheet runoff under the KW condition; while on the left side, the water depth is more than 0.2 m, and the flow is no longer governed by the KW condition. This model provides the outflow hydrograph as well as the temporal and spatial distribution of the water depth and flow velocity, which can be used for studying soil erosion, agropollutant transport, and water quality.

The gage station (Figure 11) recorded the channel water surface elevation at regular time intervals, but velocities were generally too low for accurate

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model DOI: http://dx.doi.org/10.5772/intechopen.80214

Figure 12.

Information and simulation results in an area indicated in Figure 11 in a dashed rectangle: (a) bed elevation contours, (b) velocity direction distribution near the end of the simulation, (c) velocity direction, and (d) the water depth distribution at the peak of the April 2011 rainfall.

measurement, and therefore, water discharge was not measured. In order to better understand the watershed hydrology, a rating curve of the form:

$$Q = r(L - L\_0)^x \tag{11}$$

was developed using simulated discharge, in which L is the measured water surface elevation, r and z are parameters, and L0 is the initial water surface elevation prior to a rainfall event. Eq. (11) has two unknown parameters, but there is only one relationship available for determining their values. The total volume of runoff, obtained by numerically integrating Eq. (11) in time, is equal to the rain volume, VR:

$$\sum\_{i} Q\_{i} \Delta t = \sum\_{i} r (L\_{i} - L\_{0})^{x} \Delta t = V\_{R} \tag{12}$$

in which Li is the measured water surface elevation at the gage station. With Eq. (12) satisfied, values of r and z that best fit the shape of the discharge hydrograph computed using Eq. (11), and that of the numerical simulation, were determined for each event by trial and error.

Attempts were made to fit all simulated curves using a single set of values for r, z and a mean base stage L0, but the result showed unacceptable discrepancies. L0 varied due to antecedent precipitation, downstream discharge control, sedimentation, and water usage between events. The range of L0 for the studied events is 0.33 m (Table 1). Given the complexities of the hydraulic regime in the water body, varying from standing to moving state and with varying downstream controls, variable rating curve parameters are sensible. Event-specific rating curve parameters are not ideal but are useful in a research context.

Manning's roughness coefficient (n) is a major factor in the determination of watershed runoff characteristics and generally reflects ground cover and management. The event on April 27–28, 2011 was used for initial calibration of Manning's coefficient. The studied watershed is cultivated with soybeans (Glycine max L. Merr.), corn (Zea mays L.), cotton (Gossypium hirsutum L.), and rice (Oryza sativa L.).

deeper, and the flow directions (Figure 12c) are less affected by the local microtopographic features. The flow on the right side of the domain is still sheet runoff under the KW condition; while on the left side, the water depth is more than 0.2 m, and the flow is no longer governed by the KW condition. This model provides the outflow hydrograph as well as the temporal and spatial distribution of the water depth and flow velocity, which can be used for studying soil erosion, agro-

regular time intervals, but velocities were generally too low for accurate

The gage station (Figure 11) recorded the channel water surface elevation at

Numerical simulation identified watershed for the gage station. Simulation results in the dashed rectangle area

4/27–4/28/2011 88.39 85,817 2.4 1.223 0.45 10/30–11/4/2013 53.59 52,182 1.9 4.24 0.78 11/21–25/2011 62.99 61,333 1.4 1.613 0.45 5/20–24/2013 48.77 47,483 1.0 1.436 0.59 9/25–27/2011 52.32 50,946 1.8 5.211 0.48

) zr L0 (m)

Event Measured rainfall (mm) Runoff volume\* (m<sup>3</sup>

pollutant transport, and water quality.

\*Computed from the main bulk of the rain event.

Parameters of selected runoff events for numerical simulations.

Soil Erosion - Rainfall Erosivity and Risk Assessment

Figure 11.

Table 1.

46

are shown in Figure 12.

The sensitivity of the CCHE2D model in Howden Lake watershed to Manning's n was examined using a wide range of values from 0.030 to 0.30 m1/3s. Smaller Manning's n results in a higher runoff peak discharge and an earlier peak flow arrival time. A visual comparison of discharge hydrographs based on stage measurements and numerical simulation (Figure 13) indicates that n = 0.3 m1/3s is the most appropriate choice for the overland runoff area because the peak times of these runoff events are consistent. Considering that the depths of the sheet runoff are much smaller than the microtopographic irregularities over the fields, the calibrated n represents not only the bed resistance but also form drags due to the microbed forms, crop residue, and vegetation. This n value agrees with the recent runoff studies [25, 31, 49] in cases of overland flows, including those in the Goodwin Creek Experimental Watershed in Northern Mississippi. There are numerous trees, bushes, and weeds growing along and within the channel, thus, n = 0.16 m1/3s was used for the channel and kept unchanged for other rain event cases.

The total observed rainfall volume for the April 27–28, 2011 event (Figure 13) was approximately 86,000 m<sup>3</sup> (88.32 mm). The total simulated runoff volume is about 80,600 m<sup>3</sup> (83.78 mm), which is reasonable because the hydrograph recession limb extended past the simulation termination at 47 h. There were several small rain events that occurred before the event shown in Figure 13, so the runoff volume based on the observed water surface elevation may include recession of the earlier events.

Figure 14 compares the discharge hydrographs of several additional runoff events computed using Eq. (12) and that of the numerical simulations. The identified parameters for these events, r and z, are listed in Table 1. Events 9/2011 and 10/2013 have one major peak, while those of 11/2011 and 5/2013 each have two major peaks. The simulated hydrographs fit well with those computed using Eq. (11). The two rain peaks of the 5/2013 event were separated by about 2 h, but those of the 11/2011 event were separated by 15 h. The runoff of the 5/2013 event showed only one peak because the two rain peaks were very close, and the runoff peaks were superimposed. However, the temporal separation of the two peaks of the 11/2011 event was much longer. Therefore, the superimposed hydrologic response also displayed two peaks. These watershed responses were reproduced by the numerical simulations.

As noted above, the watershed has multiple field ditches that convey runoff into the channel (Figures 10 and 11). Ditch and channel flow were simulated together

with the overland sheet runoff. Figure 15 shows the simulated flows in the channel network of the watershed. The contours represent the distribution of the unit flow discharge. The vectors in the ditches and in the stream formed a channel network indicated by the large velocity vectors; those in the runoff area are too small to be seen. The flows in the stream are turbulent when the rainfalls are large. Because no velocity data were acquired, the simulated velocity results in the channel were not

Simulated flow in the network of drainage ditches and the stream in the watershed.

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model

DOI: http://dx.doi.org/10.5772/intechopen.80214

The numerical model CCHE2D was used to model sheet runoff from watersheds, large and complex enough to include both overland and channel flow processes. The model was systematically verified and validated using analytical solutions and

validated.

49

Figure 15.

Figure 14.

Comparisons of simulated runoff and Eq. (11).

7. Conclusions

Figure 13. Sensitivity of simulated hydrograph to Manning's coefficient.

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model DOI: http://dx.doi.org/10.5772/intechopen.80214

Figure 14. Comparisons of simulated runoff and Eq. (11).

The sensitivity of the CCHE2D model in Howden Lake watershed to Manning's n was examined using a wide range of values from 0.030 to 0.30 m1/3s. Smaller Manning's n results in a higher runoff peak discharge and an earlier peak flow arrival time. A visual comparison of discharge hydrographs based on stage measurements and numerical simulation (Figure 13) indicates that n = 0.3 m1/3s is the most appropriate choice for the overland runoff area because the peak times of these runoff events are consistent. Considering that the depths of the sheet runoff are much smaller than the microtopographic irregularities over the fields, the calibrated n represents not only the bed resistance but also form drags due to the microbed forms, crop residue, and vegetation. This n value agrees with the recent runoff studies [25, 31, 49] in cases of overland flows, including those in the Goodwin Creek Experimental Watershed in Northern Mississippi. There are numerous trees, bushes, and weeds growing along and within the channel, thus, n = 0.16 m1/3s was used for the channel and kept

The total observed rainfall volume for the April 27–28, 2011 event (Figure 13) was approximately 86,000 m<sup>3</sup> (88.32 mm). The total simulated runoff volume is about 80,600 m<sup>3</sup> (83.78 mm), which is reasonable because the hydrograph recession limb extended past the simulation termination at 47 h. There were several small rain events that occurred before the event shown in Figure 13, so the runoff volume based on the observed water surface elevation may include recession of the earlier

Figure 14 compares the discharge hydrographs of several additional runoff events computed using Eq. (12) and that of the numerical simulations. The identified parameters for these events, r and z, are listed in Table 1. Events 9/2011 and 10/2013 have one major peak, while those of 11/2011 and 5/2013 each have two major peaks. The simulated hydrographs fit well with those computed using Eq. (11). The two rain peaks of the 5/2013 event were separated by about 2 h, but those of the 11/2011 event were separated by 15 h. The runoff of the 5/2013 event showed only one peak because the two rain peaks were very close, and the runoff peaks were superimposed. However, the temporal separation of the two peaks of the 11/2011 event was much longer. Therefore, the superimposed hydrologic response also displayed two peaks. These watershed responses were reproduced by

As noted above, the watershed has multiple field ditches that convey runoff into the channel (Figures 10 and 11). Ditch and channel flow were simulated together

unchanged for other rain event cases.

Soil Erosion - Rainfall Erosivity and Risk Assessment

the numerical simulations.

events.

Figure 13.

48

Sensitivity of simulated hydrograph to Manning's coefficient.

Figure 15. Simulated flow in the network of drainage ditches and the stream in the watershed.

with the overland sheet runoff. Figure 15 shows the simulated flows in the channel network of the watershed. The contours represent the distribution of the unit flow discharge. The vectors in the ditches and in the stream formed a channel network indicated by the large velocity vectors; those in the runoff area are too small to be seen. The flows in the stream are turbulent when the rainfalls are large. Because no velocity data were acquired, the simulated velocity results in the channel were not validated.

## 7. Conclusions

The numerical model CCHE2D was used to model sheet runoff from watersheds, large and complex enough to include both overland and channel flow processes. The model was systematically verified and validated using analytical solutions and

experimental data due to steady and unsteady rainfall intensity, and applied to a real world watershed. Good agreement between the analytical solutions, experimental data, and numerical simulations were obtained. For the experimental cases involving complex watershed shapes, the numerical model has the ability to simulate runoff over the slope surfaces and the channel flows.

A numerical scheme was developed to correct the bilinear interpolation of the water surface elevation from its solutions at the staggered cell centers to the collocation nodes. The scheme was necessary and effective for obtaining good sheet runoff simulation results in watersheds with irregular topography. One would have to smooth the ground topography if a model requires the interpolation of water surface solution under this condition.

The model was applied to an agricultural watershed in the Mississippi River alluvial plain. It was useful to identify the boundary of the monitored watershed and develop the rating curve at the gage station of the watershed. Several significant runoff events were selected for simulation. Each of the simulated runoff hydrographs and the rating curves agreed well with those observed in the field. The sensitivity of the model to overland sheet flow friction was studied. An increase in the bed surface friction coefficient significantly diminishes the peak of runoff discharge, delaying its time of arrival. Values of n = 0.2–0.3 m1/3s for overland flow were found to be adequate to best fit the numerical simulations and the observed data in the studied watershed. With a high-resolution mesh, the model can predict the complex surface runoff pattern over the agricultural land. Ditch and stream channels in the domain are a connected channel network. The model is able to simulate sheet runoff, turbulent channel flow, and their transitions seamlessly. The simulated hydrological processes for several storm events fit well to those observed at the gage station. The capability would be useful for studies related to soil erosion and agro-pollutant transport. The model is currently used for watershed applications without considering interception, evapotranspiration, and infiltration. Additional work is needed to further extend the research in these areas.

### Acknowledgements

This work is supported in part by USDA Agriculture Research Service under the Research Project No. 6060-13000-025-00D (NCCHE) monitored by the USDA-ARS National Sedimentation Laboratory (NSL). Support is also in part by the Southeast Region Research Initiative (SERRI) project and the University of Mississippi (UM).

Author details

F. Douglas Shields Jr.<sup>3</sup>

Laboratory, MS, USA

\*, Tahmina Shirmeen<sup>1</sup>

University of Mississippi, MS, USA

3 Shields Engineering, LLC, MS, USA

provided the original work is properly cited.

\*Address all correspondence to: jia@ncche.olemiss.edu

, Martin A. Locke<sup>2</sup>

2 Water Quality and Ecology Research Unit, USDA-ARS National Sedimentation

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium,

1 National Center for Computational Hydroscience and Engineering, The

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model

DOI: http://dx.doi.org/10.5772/intechopen.80214

, Richard E. Lizotte Jr.<sup>2</sup> and

Yafei Jia<sup>1</sup>

51

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model DOI: http://dx.doi.org/10.5772/intechopen.80214

## Author details

experimental data due to steady and unsteady rainfall intensity, and applied to a real world watershed. Good agreement between the analytical solutions, experimental data, and numerical simulations were obtained. For the experimental cases involving complex watershed shapes, the numerical model has the ability to simu-

A numerical scheme was developed to correct the bilinear interpolation of the water surface elevation from its solutions at the staggered cell centers to the collocation nodes. The scheme was necessary and effective for obtaining good sheet runoff simulation results in watersheds with irregular topography. One would have to smooth the ground topography if a model requires the interpolation of water

The model was applied to an agricultural watershed in the Mississippi River alluvial plain. It was useful to identify the boundary of the monitored watershed and develop the rating curve at the gage station of the watershed. Several significant

hydrographs and the rating curves agreed well with those observed in the field. The sensitivity of the model to overland sheet flow friction was studied. An increase in the bed surface friction coefficient significantly diminishes the peak of runoff discharge, delaying its time of arrival. Values of n = 0.2–0.3 m1/3s for overland flow were found to be adequate to best fit the numerical simulations and the observed data in the studied watershed. With a high-resolution mesh, the model can predict the complex surface runoff pattern over the agricultural land. Ditch and stream channels in the domain are a connected channel network. The model is able to simulate sheet runoff, turbulent channel flow, and their transitions seamlessly. The simulated hydrological processes for several storm events fit well to those observed at the gage station. The capability would be useful for studies related to

runoff events were selected for simulation. Each of the simulated runoff

soil erosion and agro-pollutant transport. The model is currently used for watershed applications without considering interception, evapotranspiration, and infiltration. Additional work is needed to further extend the research in

This work is supported in part by USDA Agriculture Research Service under the Research Project No. 6060-13000-025-00D (NCCHE) monitored by the USDA-ARS National Sedimentation Laboratory (NSL). Support is also in part by the Southeast Region Research Initiative (SERRI) project and the University of

late runoff over the slope surfaces and the channel flows.

Soil Erosion - Rainfall Erosivity and Risk Assessment

surface solution under this condition.

these areas.

Acknowledgements

Mississippi (UM).

50

Yafei Jia<sup>1</sup> \*, Tahmina Shirmeen<sup>1</sup> , Martin A. Locke<sup>2</sup> , Richard E. Lizotte Jr.<sup>2</sup> and F. Douglas Shields Jr.<sup>3</sup>

1 National Center for Computational Hydroscience and Engineering, The University of Mississippi, MS, USA

2 Water Quality and Ecology Research Unit, USDA-ARS National Sedimentation Laboratory, MS, USA

3 Shields Engineering, LLC, MS, USA

\*Address all correspondence to: jia@ncche.olemiss.edu

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## References

[1] Kirkby MJ, editor. Hillslope Hydrology. Chichester: John Wiley; 1978. p. 389

[2] Schmidt J. Soil Erosion: Application of Physically Based Models. Springer Verlag Berlin and Heidelberg GmbH Co KG; 2000. p. 318

[3] Liggett JA, Woolhiser DA. Difference solutions of the shallow water equations. Journal of the Engineering Mechanics Division, ASCE. 1967;93(2):39-71

[4] Cea L, Puertas J. Experimental validation of two-dimensional depthaveraged models for forecasting rainfallrunoff from precipitation data in urban areas. Journal of Hydrology. 2010;382: 88-102

[5] Kivva SL, Zheleznyak MJ. Twodimensional modeling of rainfall-runoff and sediment transport in small catchments areas. International Journal of Fluid Mechanics Research. 2005; 32(6):703-716

[6] Chow VT, Maidment DR, Mays LW. Applied Hydrology. New York: McGraw-Hill Book Company; 1988. p. 572

[7] Singh VP. Kinematic wave modelling in water resources: A historical perspective. Hydrological Processes. 2001;15(4):671-706

[8] Singh VP. Is hydrology kinematic? Hydrological Processes. 2002;16(3): 667-716. 1

[9] Morris EM, Woolhiser DA. Unsteady one-dimensional flow-over a plane: Partial equilibrium and recession hydrographs. Water Resources Research. 1980;16(2):355-360

[10] Book DE, Labadie JW, Morrow DM. Dynamic vs. Kinematic Routing in Modeling Urban Storm Drainage. In: Proceedings of the Second International

Conference on Urban Storm Drainage. Yen BC editor. 14-19 June 1981. Urbana, IL, USA: University of Illinois at Urbana Champaign, 154-163

overland flow models. Advances in Water Resources. 2008;31:173-180

DOI: http://dx.doi.org/10.5772/intechopen.80214

Environmental and Water Resources Congress, 2015; 17-21 May 2015; Austin,

[27] Iwagaki Y. Fundamental Studies on the Runoff Analysis by Characteristics.

[28] Govindaraju RS, Kavvas ML, Jones SE. Approximate analytical solutions for overland flows. Water Resources Research. 1990;26(12):2903-2912

[30] Cea L, Puertas J, Pena L, Garrido M. Hydrologic forecasting of fast flood events in small catchments with a 2D-SWE model. Numerical model and experimental validation. In: World Water Congress 2008; 1-4 September

Vol. 10. Japan: Bulletin - Kyoto University Disaster Prevention Res.

[29] Singh VP. Kinematic Wave Modeling in Water Resources: Surface Water Hydrology. Chichester: John

Wiley and Sons Ltd.; 1996

2008; Montpellier, France

[31] Downer CW. Demonstration of GSSHA Hydrology and Sediment Transport at the Goodwin Creek Experimental Watershed. ERDC TN-SWWRP-08-x. Vicksburg, MS: U.S. Army Engineer Research and

Development Center; 2008 Available from: https://swwrp.usace.army.mil/

[32] Yeh GT, Shih DS, Cheng JRC. An integrated media, integrated processes watershed model. Computers & Fluids.

[33] Graham DN, Butts MB. Flexible, integrated watershed modelling with MIKE SHE. In: Singh VP, Frevert DK, editors. Watershed Models. CRC Press; 2005. pp. 245-272. ISBN: 0849336090

[34] Ewen J, Parkin G, O'Connell PE. SHETRAN: Distributed river basin flow and transport modeling system. Journal of Hydrologic Engineering, ASCE. July

2011;45(2011):2-13

2000;5(3):250-258

Inst.; 1955. pp. 1-25. 2

Texas. 2015

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model

[19] Esteves M, Faucher X, Galle S, Vauclin M. Overland flow and infiltration modeling for small plots unsteady rain: Numerical results versus observed values. Journal of Hydrology.

[20] Howes DA, Abrahams AD, Pitman EB. One- and two-dimensional modelling of overland flow in semiarid shrubland, Jornada basin, New Mexico. Hydrological

[21] Liu QQ, Chen L, Li JC, Singh VP. Two-dimensional kinematic wave model of overland-flow. Journal of Hydrology. 2004;291(1–2):28-41

[22] Costabile P, Costanzo C, Macchione F, Mercogliano P. Two-dimensional model for overland flow simulations: A case study. European Water. 2012;38:

[23] Nunoz-Carpena R, Miller CT, Parsons JE. A quadratic Petrov-Galerkin solution for kinematic wave overland flow. Water Resources Research. 1993;

[24] Venkata RK, Eldho TI, Rao EP. A diffusion wave based integrated FEM-GIS model for runoff simulation of small watersheds. Journal of Water Resource and Protection. 2009;1:391-399.

Published online December 2009. DOI:

generated overland flow using nonlinear shallow-water equations. Journal of Hydrologic Engineering. 2104:1943- 5584. https://doi.org/10.1061/(ASCE)

[25] Singh J, Altinakar MS, Ding Y. Numerical modeling of rainfall-

[26] Shirmeen T, Jia Y, Locke MA, Lizotte Jr. RE. Numerical modeling of rain induced overland flows. In: World

10.4236/jwarp.2009.16047

HE.1943-5584.0001124

53

Processes. 2006;20:1027-1046

2000;228:265-282

13-23

29(8):2615-2627

[11] Woolhiser DA, Liggett JA. Unsteady, one dimensional flow over a plane—The rising hydrograph. Water Resources Research. 1967;3(3):753-771

[12] Freeze RA. Mathematical models of hillslope hydrology. In: Kirkby MJ, editor. Hillslope Hydrology. New York: Wiley Interscience; 1978. pp. 177-225

[13] Cundy TW, Tento SW. Solution to the kinematic wave approach to overland flow routing with rainfall excess given by the Philip equation. Water Resources Research. 1985;21: 1132-1140

[14] de Lima JLMP, Singh VP. The influence of the pattern of moving rainstorms on overland flow. Advances in Water Resources. 2002;25:817-828

[15] Rousseau M, Cerdan O, Delestre O, Dupros F, James F, Cordier S. Overland flow modelling with the shallow water equation using a well-balanced numerical scheme: Adding efficiency or just more complexity? 2012. Available from: https://hal.archives-ouvertes.fr/ hall-00664535

[16] Du JK, Shunping X, Youpeng X, Xu C, Singh VP. Development and testing of a simple physically-based distributed rainfall-runoff model for storm runoff simulation in humid forested basins. Journal of Hydrology. 2007;336:334-346

[17] Zhang W, Cundy TW. Modeling of two-dimensional overland flow. Water Resources Research. 1989;25(9): 2019-2035

[18] Gottardi G, Vinutelli M. An accurate time integration method for simplified

Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model DOI: http://dx.doi.org/10.5772/intechopen.80214

overland flow models. Advances in Water Resources. 2008;31:173-180

References

1978. p. 389

KG; 2000. p. 318

88-102

32(6):703-716

[1] Kirkby MJ, editor. Hillslope Hydrology. Chichester: John Wiley;

[2] Schmidt J. Soil Erosion: Application of Physically Based Models. Springer Verlag Berlin and Heidelberg GmbH Co

Soil Erosion - Rainfall Erosivity and Risk Assessment

Conference on Urban Storm Drainage. Yen BC editor. 14-19 June 1981. Urbana, IL, USA: University of Illinois at Urbana

Unsteady, one dimensional flow over a plane—The rising hydrograph. Water Resources Research. 1967;3(3):753-771

[12] Freeze RA. Mathematical models of hillslope hydrology. In: Kirkby MJ, editor. Hillslope Hydrology. New York: Wiley Interscience; 1978. pp. 177-225

[13] Cundy TW, Tento SW. Solution to the kinematic wave approach to overland flow routing with rainfall excess given by the Philip equation. Water Resources Research. 1985;21:

[14] de Lima JLMP, Singh VP. The influence of the pattern of moving rainstorms on overland flow. Advances in Water Resources. 2002;25:817-828

[15] Rousseau M, Cerdan O, Delestre O, Dupros F, James F, Cordier S. Overland flow modelling with the shallow water

numerical scheme: Adding efficiency or just more complexity? 2012. Available from: https://hal.archives-ouvertes.fr/

[16] Du JK, Shunping X, Youpeng X, Xu C, Singh VP. Development and testing of a simple physically-based distributed rainfall-runoff model for storm runoff simulation in humid forested basins. Journal of Hydrology. 2007;336:334-346

[17] Zhang W, Cundy TW. Modeling of two-dimensional overland flow. Water Resources Research. 1989;25(9):

[18] Gottardi G, Vinutelli M. An accurate time integration method for simplified

equation using a well-balanced

Champaign, 154-163

1132-1140

hall-00664535

2019-2035

[11] Woolhiser DA, Liggett JA.

[3] Liggett JA, Woolhiser DA. Difference solutions of the shallow water equations. Journal of the Engineering Mechanics Division, ASCE. 1967;93(2):39-71

[4] Cea L, Puertas J. Experimental validation of two-dimensional depthaveraged models for forecasting rainfallrunoff from precipitation data in urban areas. Journal of Hydrology. 2010;382:

[5] Kivva SL, Zheleznyak MJ. Twodimensional modeling of rainfall-runoff

[6] Chow VT, Maidment DR, Mays LW. Applied Hydrology. New York: McGraw-

[7] Singh VP. Kinematic wave modelling

[8] Singh VP. Is hydrology kinematic? Hydrological Processes. 2002;16(3):

[9] Morris EM, Woolhiser DA. Unsteady one-dimensional flow-over a plane: Partial equilibrium and recession hydrographs. Water Resources Research. 1980;16(2):355-360

[10] Book DE, Labadie JW, Morrow DM. Dynamic vs. Kinematic Routing in Modeling Urban Storm Drainage. In: Proceedings of the Second International

and sediment transport in small catchments areas. International Journal of Fluid Mechanics Research. 2005;

Hill Book Company; 1988. p. 572

in water resources: A historical perspective. Hydrological Processes.

2001;15(4):671-706

667-716. 1

52

[19] Esteves M, Faucher X, Galle S, Vauclin M. Overland flow and infiltration modeling for small plots unsteady rain: Numerical results versus observed values. Journal of Hydrology. 2000;228:265-282

[20] Howes DA, Abrahams AD, Pitman EB. One- and two-dimensional modelling of overland flow in semiarid shrubland, Jornada basin, New Mexico. Hydrological Processes. 2006;20:1027-1046

[21] Liu QQ, Chen L, Li JC, Singh VP. Two-dimensional kinematic wave model of overland-flow. Journal of Hydrology. 2004;291(1–2):28-41

[22] Costabile P, Costanzo C, Macchione F, Mercogliano P. Two-dimensional model for overland flow simulations: A case study. European Water. 2012;38: 13-23

[23] Nunoz-Carpena R, Miller CT, Parsons JE. A quadratic Petrov-Galerkin solution for kinematic wave overland flow. Water Resources Research. 1993; 29(8):2615-2627

[24] Venkata RK, Eldho TI, Rao EP. A diffusion wave based integrated FEM-GIS model for runoff simulation of small watersheds. Journal of Water Resource and Protection. 2009;1:391-399. Published online December 2009. DOI: 10.4236/jwarp.2009.16047

[25] Singh J, Altinakar MS, Ding Y. Numerical modeling of rainfallgenerated overland flow using nonlinear shallow-water equations. Journal of Hydrologic Engineering. 2104:1943- 5584. https://doi.org/10.1061/(ASCE) HE.1943-5584.0001124

[26] Shirmeen T, Jia Y, Locke MA, Lizotte Jr. RE. Numerical modeling of rain induced overland flows. In: World Environmental and Water Resources Congress, 2015; 17-21 May 2015; Austin, Texas. 2015

[27] Iwagaki Y. Fundamental Studies on the Runoff Analysis by Characteristics. Vol. 10. Japan: Bulletin - Kyoto University Disaster Prevention Res. Inst.; 1955. pp. 1-25. 2

[28] Govindaraju RS, Kavvas ML, Jones SE. Approximate analytical solutions for overland flows. Water Resources Research. 1990;26(12):2903-2912

[29] Singh VP. Kinematic Wave Modeling in Water Resources: Surface Water Hydrology. Chichester: John Wiley and Sons Ltd.; 1996

[30] Cea L, Puertas J, Pena L, Garrido M. Hydrologic forecasting of fast flood events in small catchments with a 2D-SWE model. Numerical model and experimental validation. In: World Water Congress 2008; 1-4 September 2008; Montpellier, France

[31] Downer CW. Demonstration of GSSHA Hydrology and Sediment Transport at the Goodwin Creek Experimental Watershed. ERDC TN-SWWRP-08-x. Vicksburg, MS: U.S. Army Engineer Research and Development Center; 2008 Available from: https://swwrp.usace.army.mil/

[32] Yeh GT, Shih DS, Cheng JRC. An integrated media, integrated processes watershed model. Computers & Fluids. 2011;45(2011):2-13

[33] Graham DN, Butts MB. Flexible, integrated watershed modelling with MIKE SHE. In: Singh VP, Frevert DK, editors. Watershed Models. CRC Press; 2005. pp. 245-272. ISBN: 0849336090

[34] Ewen J, Parkin G, O'Connell PE. SHETRAN: Distributed river basin flow and transport modeling system. Journal of Hydrologic Engineering, ASCE. July 2000;5(3):250-258

[35] Vieira DA, Wu W. One-Dimensional Channel Network Model CCHE1D Version 3.0: User's Manual (Report NCCHE-TR-2002-2). Mississippi, USA: National Center for Computational Hydroscience and Engineering, The University of Mississippi; 2000

[36] Jia Y, Wang SSY, Xu YC. Validation and application of a 2D model to channel with complex geometry. International Journal of Computational Engineering Science. 2002;3(1):57-71

[37] Jia Y, Zhang Y, Wang SSY. Simulating bank erosion process using a depth averaged computational model. In: 7th Symposium on River, Coastal and Estuarine Morphodynamics (RCEM 2011), 2011; 5-8 September 2011; Beijing, China, CD-ROM

[38] Jia Y, Chao XB, Zhang YX, Zhu TT. Technical Manual of CCHE2D, V4.1, NCCHE-TR-02-2013. The University of Mississippi; 2013

[39] Jia Y, Altinakar M, Chao X, Zhang Y. Numerical simulations of spilled coal ash in the Dan River and the environmental impact of the incident. In: World Environmental and Water Resources Congress. 2016. pp. 114-125. DOI: 10.1061/9780784479865.012

[40] Ding Y, Wang SSY, Jia Y. Development and validation of a quasithree dimensional coastal area morphological model. ASCE, Journal of Waterway, Port, Coastal, and Ocean Engineering. ASCE. 2006;132(6):462-476

[41] Wang SSY, Roche PJ, Schmalz RA, Jia Y, Smith PE. Verification and Validation of 3D Free-Surface Flow Models. American Society of Civil Engineering; 2008. ISBN: 968-0-7844- 0957-2

[42] Jia Y, Wang SSY. Numerical model for channel flow and morphological change studies. Journal of Hydraulic

Engineering, ASCE. 1999;125(9): 924-933

[43] Jia Y, Shirmeen T. Simulation of rainfall-runoff process in watersheds using CCHE2D. In: 12th International Conference on Hydroscience & Engineering: Hydro-Science & Engineering for Environmental Resilience; 6-10 November 2016; Tainan, Taiwan

[44] Singh VP, Regl RR. Analytical solutions of kinematic equations for erosion on a plane. Ι. Rainfall of indefinite duration. Advances in Water Resources. 1981;6(1):2-10

[45] Singh VP. Analytical solutions of kinematic equations for erosion on a plane. ΙΙ. Rainfall of finite duration. Advances in Water Resources. 1983; 6(2):88-95

[46] Morgali JR, Linsley RK. Computer analysis of overland flow. Journal of Hydrological Division, Proc. American Society of Civil Engineers. 1965;HY3: 81-100

[47] Soil Survey Staff, Natural Resources Conservation Service, United States Department of Agriculture, Web Soil Survey. Available from: https:// websoilsurvey.nrcs.usda.gov/ [Accessed: 10-04-2017]

[48] Fisk HN. Geological Investigations of the Alluvial Valley of the Lower Mississippi River. Vicksburg, Mississippi: U.S. Army Corps of Engineers, Mississippi River Commission; 1944. p. 78

[49] Kalyanapu AJ, Burian SJ, McPherson TN. Effect of land use-based surface roughness on hydrologic model output. Journal of Spatial Hydrology. 2009;9(2):51-71. Fall 2009

**55**

Section 3

Assessment of Soil Erosion

Risk

## Section 3
