One-Dimensional Numerical Cascade Model of Runoff and Soil Loss on Convergent and Divergent Plane Soil Surfaces: Labora-Tory Assessment and Numerical Simulations

1. Introduction

Rainfall-induced soil erosion poses a serious environmental threat, leading to reduced fertility of arable lands, lower crop yields, intensified surface runoff, and the transportation of harmful pollutants to freshwater bodies. Consequently, soil conservation is crucial for sustainable agriculture, necessitating accurate estimation of soil loss and identification of critical areas for implementing soil conservation measures, especially on complex soil surface geometries.

Mathematical models offer cost-effective tools to understand rainfall-induced soil erosion processes, evaluate soil and water conservation practices under different environmental conditions, and assist in policy development for future climate scenarios. However, no universal model can accurately predict soil erosion across all temporal and spatial scales.

Characterizing rainfall-induced runoff response is essential for estimating soil erosion. The kinematic wave (KW) theory is widely used for modelling overland flow, with one-dimensional (1D) KW approximations proving effective on small slopes [1,2,3].

Sediment transport data, which are crucial for testing soil erosion models, can be obtained from field plots and laboratory soil flumes, taking into account factors such as rainfall characteristics, slope geometry, soil properties and land cover. While numerous studies have evaluated, under controlled laboratory conditions, physically based soil erosion and transport models on simple geometries such as square or rectangular plane soil surface geometries [4,5,6,7,8], complex soil surfaces, common in nature, are less studied.

Only a few researchers [9,10,11,12,13] have examined erosion on curved rectangular surfaces (convex, plane and concave). These types of geometries are similar to those often found in rocky catchments where landslides are common [14,15]. To the best of our knowledge, to date, no experimental or numerical studies have assessed the impact of varying converging and diverging angles of plane soil surfaces on rainfall-induced runoff and soil loss. Additionally, no physically based 1D model exists for estimating erosion on these surfaces.

This paper has two objectives: first, to develop a 1D cascade model for estimating rainfall-induced runoff and soil erosion on converging and diverging plane surfaces, and to test it using controlled laboratory experiments. The physically based mathematical model, grounded in mass conservation principles, simulates sediment detachment, transport and deposition processes [16,17]. Laboratory experiments were conducted using a rainfall simulator and a soil flume with a planar circular sector geometry to create converging and diverging flow conditions. Second, the paper evaluates the influence of different plane surfaces’ converging and diverging angles on runoff and soil loss through numerical simulations on seven plane surfaces, including one rectangular/linear configuration (angle along the flow direction θ=0°), and surfaces with three unique converge and divergence angles (θ=8.5°, 16.7°, and 24.2°). Soil surface area, rainfall intensity and slope gradient were kept constant across all numerical simulation.

2. Proposed Cascade Model

Soil erosion is a complex phenomenon influenced by the geometry and geomorphology of soil surfaces. The 1D cascade model used in this work to estimate rainfall-induced runoff and soil erosion represents the soil surface geometry as a series of cascading rectangular planes. The model consists of three components: (i) Horton’s equation for soil infiltration; (ii) the KW approximation of the 1D Saint-Venant shallow water equations for overland flow generation; and (iii) the 1D sediment transport equation for soil loss calculation.

2.1. Infiltration Model

Horton [18,19,20] introduced an exponential decay equation to describe how the infiltration rate varies over time after a soil surface is exposed to constant rainfall:

$$f_c\left(t\right)=f_1+\left(f_0-f_1\right)e^{-kt}$$

(1)

where $f_c\left(t\right)$ is the infiltration rate at any time t, $f_0$is the initial infiltration rate when rainfall starts (at t=0), $f_1$ is the final or equilibrium infiltration rate, and $k$ is the exponential time decay coefficient. Cumulative infiltration $F$ a at any time t is:

$$F={\int }_0^t{f}_c\left(t\right)dt=f_{1}t+{\displaystyle \frac{\left(f_0-f_1\right)}{k}} \left(1-e^{-kt}\right)$$

(2)

Using Equation (1) to eliminate t in Equation (2), the resultant equation becomes:

$$F={\displaystyle \frac{f_0-f_c}{k}}-{\displaystyle \frac{f_1}{k}}\ln \left({\displaystyle \frac{f_c-f_1}{f_0-f_1}}\right)$$

(3)

Before ponding occurs, cumulative infiltration is given by $F=It$, where $I$ is the mean rainfall intensity. Ponding occurs when the infiltration rate equals the rainfall intensity $I$. Setting $f_c=I$ in Equation (3), one obtains:

$$F_p={\displaystyle \frac{f_0-I}{k}}-{\displaystyle \frac{f_1}{k}}\ln \left({\displaystyle \frac{I-f_1}{f_0-f_1}}\right)$$

(4)

where $F_p$ is the cumulative infiltration at ponding time $t_p$ [21,22], which is:

$$t_p={\displaystyle \frac{F_p}{I}}={\displaystyle \frac{f_0-I}{kI}}-{\displaystyle \frac{f_1}{kI}}\ln \left({\displaystyle \frac{I-f_1}{f_0-f_1}}\right)$$

(5)

2.2. Overland Flow Model

Overland flow on any given plane is described by the KW approximation of the 1D Saint-Venant equation [1,2,23]:

$${\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{1}{W}}{\displaystyle \frac{\partial Q}{\partial x}}=r$$

(6)

where $h$ is the water depth, $W$ is the runoff width, $Q$ is the runoff discharge, $r$ is the rainfall excess intensity, $t$ is time, and $x$ is the distance in the flow direction.

The KW equation assumes that the friction slope equals the bed slope. By using existing open channel flow equations, a relationship between discharge $Q$ and water depth h at any point on a plane can be expressed as:

$$Q=\alpha W{h}^{\beta }$$

(7)

where α is the KW resistance parameter and β is the Bakhmeteff dimensionless coefficient, assumed to be equal to 5/3. The parameter α is related to Manning’s roughness coefficient $n$ and can be calculated as $\alpha =n S^{1/2}$, where $S$ is the bed slope of the plane.

Substituting Equation (7) in Equation (6), the KW equation can be rewritten as:

$${\displaystyle \frac{\partial h}{\partial t}}+\beta \alpha h^{\beta -1}{\displaystyle \frac{\partial h}{\partial x}}=r$$

(8)

A cascade of rectangular n-planes representing converging or diverging surfaces exposed to rainfall is shown in Figure 1. Conceptually, these n-planes could vary in their lengths, widths, and slopes.

Here, the numerical scheme for overland flow computation is described for the converging cascade planes as an example. The governing equation for overland flow (Equation (8)) is derived using the second-order single-step Lax–Wendroff numerical scheme, expressed in finite difference form as [23]:

$$\begin{gathered} h_j^{i+1}=h_j^i+∆t\left(r_j^i-\beta \alpha {\displaystyle \frac{h_{j+1}^{i^{\beta -1}}+h_{j-1}^{i^{\beta -1}}}{2}} {\displaystyle \frac{h_{j+1}^i-h_{j-1}^i}{2∆x}}\right)+{\displaystyle \frac{{\left(∆t\right)}^2}{2}} {\displaystyle \frac{r_j^{i+1}-r_j^i}{∆t}} \\ -\beta \alpha {\displaystyle \frac{{\left(∆t\right)}^2}{2∆x}} \left[{\displaystyle \frac{h_{j+1}^{i^{\beta -1}}+h_j^{i^{\beta -1}}}{2}} \left({\displaystyle \frac{r_{j+1}^i+r_j^i}{2}}-\beta \alpha {\displaystyle \frac{h_{j+1}^{i^{\beta -1}}+h_j^{i^{\beta -1}}}{2}}{\displaystyle \frac{h_{j+1}^i-h_j^i}{∆x}}\right) \\ -{\displaystyle \frac{h_j^{i^{\beta -1}}+h_{j-1}^{i^{\beta -1}}}{2}}\left({\displaystyle \frac{r_j^i+r_{j-1}^i}{2}}-\beta \alpha {\displaystyle \frac{h_j^{i^{\beta -1}}+h_{j-1}^{i^{\beta -1}}}{2}} {\displaystyle \frac{h_j^i-h_{j-1}^i}{∆x}}\right)\right] \end{gathered}$$

(9)

where i and j represent the grid points in time and space, respectively. For the downstream boundary of each plane, the following first order scheme is used [23]:

$$h_{jL}^{i+1}=h_{jL}^i+∆t\left(r_{jL}^i-\beta \alpha {\displaystyle \frac{h_{jL}^{i^{\beta -1}}+h_{jL-1}^{i^{\beta -1}}}{2}} {\displaystyle \frac{h_{jL}^i-h_{jL-1}^i}{∆x}}\right) \mathrm{for} x=L$$

(10)

where $h_{jL}$ refers to the last point in space $x$, and $L$ is the length of each plane along the flow direction. The initial condition for every plane is given as:

$$h\left(x,0\right)=0, \mathrm{for} 0\le x\le L \mathrm{and} t=0$$

(11)

The upstream boundary condition for the first plane k is given as:

$$h_k\left(0,t\right)=0, \mathrm{for} x=0 \mathrm{and} 0\le t\le \infty$$

(12)

For the remaining cascade planes, the discharge leaving the downstream boundary of each plane enters the upstream boundary of the next plane, serving to establish the boundary condition required for the computation of water depth of that plane.

As an example, the computation of the upstream boundary condition for plane n is given as follows, respectively, for the upstream boundary entering discharge and water depth:

$$Q_{n-l}\left(x,t\right)={\alpha }_{n-l}{W}_{n-l}{h}_{n-l}^{\beta }\left(x,t\right), \mathrm{for} x=L_{n-l} (\mathrm{plane} nl) \mathrm{and} 0\le t\le \infty$$

(13)

$$h_n\left(0,t\right)={\left({\displaystyle \frac{Q_{n-l}}{{\alpha }_n{W}_n}}\right)}^{{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\beta $}\right.}}, \mathrm{for} x=0 (\mathrm{plane} n) \mathrm{and} 0\le t\le \infty$$

(14)

The numerical scheme needs to maintain the Courant number for stability, i.e.:

$${\displaystyle \frac{\mathsf{\Delta }t}{\mathsf{\Delta }x}}<{\displaystyle \frac{1}{\alpha \beta h^{\beta -1}(x,t)}}$$

(15)

2.3. Soil Erosion Model

Rainfall-induced soil erosion involves the processes of sediment detachment, transportation, and deposition [24]. Sediment detachment and transportation occur due to the impact of raindrops and the shear stress of overland flow, known as interrill and rill erosion. Deposition refers to the settling of eroded sediments. These processes are described by the 1D continuity equation of sediment transport [24,25]:

$${\displaystyle \frac{\partial \left(ch\right)}{\partial t}}+{\displaystyle \frac{\partial \left(cq\right)}{\partial x}}=e_R-d+e_I$$

(16)

where $h$ is water depth, $q$ is local discharge per unit width, $c$ is sediment concentration in overland flow, $e_R$ is the rill erosion rate, $d$ is the sediment deposition rate, $e_I$ is the interrill erosion rate, $t$ is time, and $x$ is the distance in the flow direction. The interrill erosion rate $e_I$ is given by [26]:

$$e_I=K_{I }I r$$

(17)

where $K_{I }$ is the interrill coefficient, $I$ is rainfall intensity, and $r$ is rainfall excess intensity. Rill erosion occurs when sediments cannot resist the shear stress of overland flow, and its rate can be expressed as [24,27]:

$$e_R=K_R {\tau }_e^b$$

(18)

where $K_R$ is a soil detachability factor for shear stress, ${\tau }_e$ the average "effective" shear stress assuming broad shallow flow, and $b$ is an exponent typically ranging from 1.0 to 2.0.

The sediment deposition rate $d$ is expressed as [28]:

$$d=\in c{ V}_s$$

(19)

where $\in$ is a dimensionless coefficient depending on soil and fluid properties, $c$ is the sediment concentration in overland flow, and${ V}_s$ is the settling velocity of sediment particle, which can be expressed as:

$$V_s=W_o\sqrt{{\displaystyle \frac{\left({\gamma }_s-\gamma \right)}{\gamma }}g{d}_s}$$

(20)

with

$$W_o=\sqrt{{\displaystyle \frac{2}{3}}+{\displaystyle \frac{36{\nu }^2}{g d_s^3\left(\frac{{\gamma }_s}{\gamma }-1\right)}}}-\sqrt{{\displaystyle \frac{36{\nu }^2}{g d_s^3\left(\frac{{\gamma }_s}{\gamma }-1\right)}}}$$

(21)

where ${\gamma }_s$ and $\gamma$ are the specific weights of sediment and water, respectively, $\nu$ is the kinematic viscosity of water, $g$ is the acceleration due to gravity, and $d_s$ is the size of settling sediment.

The numerical procedure proposed to solve the sediment continuity equation (Equation (16)) uses the four point-implicit scheme [16,29]:

$$\begin{gathered} {\displaystyle \frac{\varphi }{∆t}}\left[{\left(ch\right)}_{j+1}^{i+1}-{\left(ch\right)}_{j+1}^i\right]+{\displaystyle \frac{1-\varphi }{∆t}} \left[{\left(ch\right)}_j^{i+1}-{\left(ch\right)}_j^i\right]+{\displaystyle \frac{\omega }{∆x}}\left[{\left(cq\right)}_{j+1}^{i+1}-{\left(cq\right)}_j^i\right]+{\displaystyle \frac{1-\omega }{∆x}}\left[{\left(ch\right)}_{j+1}^i-{\left(cq\right)}_j^i\right]= \\ \omega \left[\varphi e_{R_{j+1}}^{i+1}+\left(1-\varphi \right)e_{R_j}^{i+1}\right]+\left(1-\omega \right) \left[\varphi e_{R_{j+1}}^i+\left(1-\varphi \right)e_{R_j}^i\right]-\omega \left[\varphi d_{j+1}^{i+1}+\left(1-\varphi \right) d_j^{i+1}\right] \\ -\left(1-\omega \right) \left[\varphi d_{j+1}^i+\left(1-\varphi \right) d_j^i\right]+\omega e_I^{i+1}+\left(1-\omega \right) e_I^i \end{gathered}$$

(22)

where $\varphi$ is the weighting factor for space, and ω is the weighting factor for time. If $\varphi$=1 and ω=0, the numerical scheme becomes explicit and needs Courant number (Equation (15)) to be maintained for stability. If $\varphi$=0.5 and ω > 0.5, the scheme becomes implicit and is unconditionally stable; however, for accuracy, Courant number must still be maintained. For our simulations, we used $\varphi$=0.5 and ω=0.6, which made the numerical scheme stable for the computation of small volumes of runoff and sediment, typically associated with laboratory conditions, which were used in this study to test the performance of the model.

The initial and upstream boundary conditions for each cascade plane are [16]:

$$c\left(x,t_p\right)={\displaystyle \frac{K_{I }I \left(t_p\right) r\left(t_p\right)}{\in V_s+r\left(t_p\right)}}, \mathrm{for} x\ge 0\mathrm{a}\mathrm{n}\mathrm{d}t=t_p$$

(23)

$$c\left(0,t\right)={\displaystyle \frac{K_{I }I \left(t\right) r\left(t\right)}{\in V_s+r\left(t\right)}},\mathrm{f}\mathrm{o}\mathrm{r}x=0\mathrm{a}\mathrm{n}\mathrm{d}t\ge t_p$$

(24)

These boundary conditions imply that the numerical scheme will start computing sediment concentration from the ponding time $t_p$. The downstream boundary condition was considered open for each cascade plane.

2.4. Calibration

To achieve a more accurate representation of the real-world system being studied and ensure that the model outputs align closely with the observed data, several model parameters can be adjusted, in particular infiltration parameters ($f_0$, $f_1$, and $k$, in Equation (1)) and two soil erosion parameters ($K_{I }$ and $K_R$, in Equations (17) and (18), respectively). This process enhances the reliability and predictive capability of the model for the specific conditions under which it was calibrated.

3. Materials and Methods

The performance of the model proposed in section 2 was evaluated by comparing numerical results with laboratory data from experiments using a rainfall simulator and a soil flume. This chapter presents the model testing procedure (section 3.1), including a description of the laboratory experiments and the corresponding scheme used in the numerical model. Additionally, it describes the numerical simulations used to explore the effect of several parameters on the model’s response, along with the scenarios adopted for this purpose (section 3.2).

3.1. Model Testing Procedure

3.1.1. Experimental Data

Data for testing the model was obtained from laboratory experiments conducted on a soil flume under controlled conditions, using simulated rainfall (Figure 2). The rainfall simulator’s sprinkler was equipped with a 3.58 mm orifice diameter HH-22 FullJet nozzle, installed 2.2 m above the soil flume surface, producing a full cone spray. This simulator has been used in several other studies [11,30,31,32]. [32] calculated the raindrop mean diameter and mean fall velocity as 0.74 mm and 1.92 m s⁻¹, respectively, at a constant nozzle operational pressure of 50 kPa for the same nozzle used in this study. To maintain this constant pressure, a hydraulic system was employed, pumping water from a constant-head reservoir with a submerged pump. The rainfall spatial distribution was measured by placing rain gauges 0.15 m apart on the soil flume and measuring the collected rainfall volume during a 2-minute rainfall event. The uneven rainfall spatial distribution on the flume was a direct consequence of using a single full-cone nozzle, positioned in the middle of the flume surface approximately 1 m from the outlet, along the slope.

The soil flume had a planar circular sector geometry with a length of 2.02 m, limited by arc lengths of 0.10 m and 2.10 m, leading to the formation of converging and diverging flow conditions. The flume’s depth was 0.12 m, and it was constructed from zinc-coated iron metal sheets. A geotextile rug was placed over the iron mesh at the bottom of the soil flume to allow for free drainage. However, drainage flux was not measured, as the focus was solely on surface runoff measurements needed to validate the numerical model’s overland flow generation component. The soil flume had runoff collection outlets at both ends, and by tilting the flume in opposite directions, experiments could be performed on both converging and diverging plane surfaces.

After sieving through a 5 mm mesh, sandy loam soil was spread manually in the flume in layers and gently compacted with a steel plate, The soil surface was levelled using a wooden blade, with a final soil depth of 0.1 m. Before each experiment, the soil was gently saturated with water (using a hose connected to tap water) until ponding occurred, followed by 1.5 hours of drainage. After each experiment, the soil was removed, replaced with air-dried soil, and the saturation procedure was repeated to ensure consistent soil moisture across all experiments.

A total of four experiments were conducted: two on converging and two on diverging plane surfaces. Each experiment was performed once without repetition. The data on applied mean rainfall intensity and the slope of the converging and diverging plane surfaces are given in Table 1. Rainfall duration for each experiment was 7 minutes. High rainfall intensity (≈ 66 mm/h) and a steep gradient (20%) created conditions prone to high soil erosion, while lower rainfall intensity (≈ 37 mm/h) with a gentle gradient (5%) ensured low soil loss conditions on both surface geometries. This approach was used to test the model’s ability to simulate runoff and erosion under these contrasting scenarios, for converging and diverging plane surfaces.

Runoff samples were collected at the flume outlet every minute for 10 seconds, from the initiation of runoff until the end of rainfall. Two additional samples were collected at 10-second intervals after the rainfall stopped. The runoff samples were weighted, oven-dried at 80°C for 24 hours and weighted again to estimate runoff and soil loss.

3.1.2. Numerical Simulations Data

For the numerical simulations aiming at calibrating and testing the model (section 2) based on the laboratory data, the surface was represented by a cascade of four planes (Figure 3), with a total area of 2.2 m², which was approximately equal to the surface area of the soil flume (2.24 m²). Each plane had a constant length ($L$) of 0.5 m along the flow direction, while the widths ($W$) of the planes varied. For a given cascade (either converging or diverging), all planes had a similar slope. Since the rainfall spatial distribution over the flume was not uniform, the simulated rainfall distribution on the cascade planes was also uneven (Figure 3).

The infiltration parameters ($f_0$, $f_1$, and $k$) and two soil erosion parameters ($K_{I }$ and $K_R$) were optimized by trial and error to fit the observed hydrographs and sediment graphs (Table 1). The remaining soil erosion parameters ($b$, $\in$, ${\gamma }_s$, $\gamma$, $d_s$ and $\nu$) were taken constants for all simulations (Table 2), which is justified as some of these parameters are physical constants, and some are expected to have minimal impact on the model results.

3.1.3. Model Evaluation

The model’s skill in simulating runoff peak, total runoff volume, soil loss peak, and total sediment mass was assessed using the Percentage of Deviation, $Dev(\%)$:

$$Dev(\%)={\displaystyle \frac{Sim-Obs}{Obs}}\times 100$$

(25)

where $Obs$ and $Sim$ represent the observed and simulated data, respectively.

The goodness of fit for the shape of the hydrographs and sediment graphs was evaluated using the Nash-Sutcliffe coefficient of efficiency ($NS$):

$$NS=1-{\displaystyle \frac{{\sum }_{t=1}^n{\left({Obs}_t-{Sim}_t\right)}^2}{{\sum }_{t=1}^n{\left({Obs}_t-{Obs}^m\right)}^2}}$$

(26)

where ${Obs}_t$ and ${Sim}_t$ represent the observed and simulated data at time t, ${Obs}^m$ is the mean of the observed data, and n is the total number of data points (n=9 in the experiments). The $NS$ values range from −∞ to 1. As the $NS$ value approaches 1, the model’s predictive efficiency improves. A value of 1 indicates a perfect fit between the model’s prediction and the observed data.

3.2. Numerical Simulation Scenarios for Varying Surface Convergence and Divergence Angles

A total of seven soil geometries with plane surface profile curvatures were numerically investigated, as shown in Figure 4. These included one rectangular/linear geometry, three diverging geometries, and three converging geometries. Each soil surface geometry had a constant area of 8 m². Three distinct angles of convergence or divergence (θ along the flow direction), specifically 8.5°, 16.7°, and 24.2°, resulted in three different geometries for both converging and diverging configurations. For the numerical analysis, the converging and diverging soil surface geometries were modelled using a cascade arrangement composed of four planes. Each individual plane maintained a consistent length of 1 m along the flow direction, with variations introduced in the widths of these planes. The rectangular surface was modelled as a single plane measuring 4 m in length and 2 m in width. For any given cascade arrangement, all planes had an identical slope of 20%. The total simulation time for each soil surface geometry was 9 minutes, with a uniform rainfall intensity of 45 mm/h simulated for 5 minutes on all cascade planes. The infiltration ($f_0$, $f_1$, and $k$) and soil erosion ($K_{I }$and $K_R$) parameters used were kept constant for all numerical simulations and are shown in Table 2. The remaining soil erosion parameters ($b$, $\in$, ${\gamma }_s$, $\gamma$, $d_s$ and $\nu$) used are also listed in Table 2, and coincide with the values used to test and calibrate the model (section 3.1.2). Identical simulation conditions ensured that only the variation in soil surface geometry influenced the runoff generation and soil loss estimation.

4. Results

4.1. Calibration and Testing of the Numerical Model

The infiltration parameters $f_0$, $f_1$, and $k$, as well as the two soil erosion parameters $K_{I }$ and $K_R$, determined during the model calibration (sections 2.4), are listed in Table 1. The remaining soil erosion parameters ($b$, $\in$, ${\gamma }_s$, $\gamma$, $d_s$ and $\nu$; see sections 2.3 and 3.1.2) used are listed in Table 2.

Simulated and observed hydrographs and sediment graphs for different rainfall intensity and slope combinations on converging and diverging plane surfaces are presented in Figure 5 and Figure 6, respectively. The numerical model’s simulated hydrographs aligned well with the observed ones, as indicated by the $NS$ values close to 1 (Figure 5). The rising limbs of the simulated hydrographs matched the observed values closely, although there was a delay in the recession limbs.

Overall, the model effectively reproduced the main characteristics of sediment transport. All simulated sediment graphs showed good agreement with the observed data, supported by $NS$ values consistently higher than 0.7 (Figure 6). However, the cascade model revealed a better capability to simulate runoff rates compared to soil loss rates.

Table 3 and Table 4 present the simulated and observed runoff peaks, runoff volumes, soil loss peaks, and total soil losses, along with the corresponding $Dev(\%)$ values, indicating that the model performed reasonably well. There were two exceptions: the total runoff in experiment 4 and the total soil loss in experiment 2, where $Dev(\%)$ values were higher. These cases involved lower rainfall intensities (≈ 37 mm/h) and gentle slopes (5%), which resulted in lower runoff volumes and sediment transport yields. Further adjustment of the model’s parameters could reduce the $Dev(\%)$ values, but this might decrease the $NS$ values for the corresponding hydrographs and sediment graphs. Although the $NS$ value, used as a model performance indicator, suggests that the model may face challenges in accurately estimating runoff volumes and total soil loss for smaller quantities, further testing of the model should be conducted under a broader range of conditions and parameters.

4.2. Numerical Simulations for Plane Soil Surfaces with Varying Converging and Diverging Angles

The infiltration and soil erosion parameters used for each numerical simulation were constant and are listed in Table 2. The simulated values for peak runoff, and peak and total soil loss across all simulated plane surface geometries, are listed in Table 5.

Among the three converging geometries, the soil surface with a 24.2° angle exhibited the lowest runoff peak while displaying the highest peak and total soil loss. Conversely, within the diverging geometries, the soil surface with an 8.5° angle recorded the lowest runoff peak, and the highest peak and total soil loss.

Figure 7 displays simulated hydrographs for seven different plane soil surface geometries. The rising and recession limbs of the hydrograph for the rectangular surface were steeper compared to those of the converging surfaces. Conversely, the rising and recession limbs of the hydrographs for the diverging surfaces were steeper than those of the rectangular surface.

The trend in the runoff peaks for all plane soil surface geometries are better depicted in Figure 8. Runoff peaks for converging surfaces were consistently lower than those for the rectangular surface; for the experimental data tested, results were 0.5–1.4% lower. Notably, an increase in the converging angle of the plane soil surface geometry led to a linear decrease in the runoff peak (R²=0.99). Conversely, for diverging soil surface geometries, the runoff peaks were higher than those for the rectangular one; for the data tested, results were 0.5–1.5% higher. An increase in the diverging angle of the soil surface geometry resulted in a linear increase in the runoff peak (R²=0.99).

Figure 9 and Figure 10 present the simulated sediment graphs and peak soil loss for the various plane soil surface geometries, respectively. These outputs also indicate a behavioral trend. For the varied tested scenarios, which led to a wide range of results, peak soil loss for converging surfaces were remarkably higher by 36–657% compared to the rectangular surface. The peak soil loss exhibited an exponential increase with the increasing converging angle of the soil surface geometry (R²=0.91). In contrast, diverging surfaces had peak soil loss 20–41% lower than the rectangular surface, and an increase in the diverging angle of the soil surface geometry led to an exponential decrease in the peak soil loss (R²=0.99). Soil loss peaks of converging surfaces exceeded those of diverging surfaces by a substantial 168–1250% for the converging/diverging angles studied.

5. Discussion

In this study, the concept of representing soil surfaces with convergent and divergent geometries using a cascade of n-planes accommodates spatial changes in these surfaces. Although the model is used here to simulate processes on plane surfaces, it can be adjusted for curved surfaces by varying the slope of the model planes, which broadens its applicability.

For the model version explored, the results reveal that when flow leaves the downstream boundary of one plane and enters the upstream boundary of the next plane, it experiences an abruptly change (shock wave phenomenon [33]) due to the change in cascade geometry. This issue, observed with a model consisting of four cascading planes, can also affect the shape of the hydrographs. Therefore, the number of cascade planes used for any specific study area needs to be increased until such abrupt oscillations are minimized.

Furthermore, the numerical model was tested with observed laboratory data from single events (without replicates). Therefore, experimental errors (such as those from soil preparation, rainfall simulator handling, and runoff sampling) cannot be ruled out. Future work should consider varying conditions (such as rainfall intensities and duration, slopes, surface shapes and dimensions, and soil types) to further evaluate the model performance. Additionally, a larger number of observed simulations should be tested.

6. Conclusions

A 1D cascade model for estimating rainfall-induced runoff and soil erosion on various soil surface geometries (rectangular/plane, converging planes, and diverging planes) was proposed and tested using a soil flume and a rainfall simulator. Overall, for the tested slopes and rainfall intensities, the simulated hydrographs and sediment graphs compared well with the observed ones, as suggested by the Nash-Sutcliffe coefficients of efficiency ($NS$). However, the model’s ability to simulate the observed hydrographs ($NS$ > 0.87) was better than to simulate the sediment graphs ($NS$ > 0.70).

To examine how varying the converging or diverging angle of plane soil surfaces affects runoff and soil loss, a series of simulations using the proposed model were conducted with seven different geometric configurations. These configurations included a rectangular/linear surface (angle along the flow direction θ=0°), three surfaces each with distinct converging angles, and three surfaces each with unique diverging angles. The studied converging and diverging angles were 8.5°, 16.7°, and 24.2°. The converging and diverging plane surfaces were modelled with a 4-plane cascade. Throughout all simulation runs, the soil surface area, rainfall intensity, slope gradient, and infiltration and soil parameters remained constant. As the convergence angle of the plane soil surface increased, there was a reduction in the runoff peak but an increase in both the peak and total soil loss. Conversely, increasing the divergence angle of the plane soil surface led to an increase in the runoff peak while causing a decline in both the peak and total soil loss. Compared to the rectangular soil surface, all converging geometries showed higher peak and total soil loss, while the diverging geometries displayed lower values. This indicates that converging plane soil surfaces are the most hazardous in terms of soil loss generation.

The presented model is thought to be a useful step towards simulating overland flow generation and sediment transport on convergent and divergent plane surfaces, avoiding the complexity of more advanced 2D models. Future studies should evaluate the ability of this model to simulate runoff and sediment loss on various curvatures (convex and concave) of different soil surface shapes (converging, linear, and diverging) under varying rainfall intensities, soil types, and slopes. The application of the model to larger areas, while also considering a larger number of planes in the cascade model, should be explored, as this will allow for the modelling of runoff and soil loss at the hillslope scale.