Influence of Rainfall Events and Surface Inclination on Overland and Subsurface Runoff Formation on Low-Permeable Soil

Citation: Gruchot, A.; Zydro´n, T.; Wał˛ega, A.; Paˇrílková, J.; Stanisz, J. Influence of Rainfall Events and Surface Inclination on Overland and Subsurface Runoff Formation on Low-Permeable Soil. Sustainability 2022,14, 4962. https://doi.org/ 10.3390/su14094962 Academic Editor: Li He Received: 17 February 2022 Accepted: 11 April 2022 Published: 20 April 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). sustainability Article Influence of Rainfall Events and Surface Inclination on Overland and Subsurface Runoff Formation on Low-Permeable Soil Andrzej Gruchot 1,* , Tymoteusz Zydro´n 1, Andrzej Wał˛ega 2, Jana Paˇrílková3and Jacek Stanisz 4 1Department of Hydraulic Engineering and Geotechnics, Faculty of Environmental Engineering and Land Surveying, University of Agriculture in Kraków, Mickiewicza 24/28, 30-059 Cracow, Poland; [email protected] 2Department of Sanitary Engineering and Water Management, Faculty of Environmental Engineering and Land Surveying, University of Agriculture in Kraków, Mickiewicza 24/28, 30-059 Cracow, Poland; [email protected] 3Institute of Water Structures, Faculty of Civil Engineering, Brno University of Technology, Veveˇrí331/95, 60200 Brno, Czech Republic; [email protected] 4Mineral and Energy Economy Research Institute, Polish Academy of Science, Józefa Wybickiego 7 A, 31-261 Cracow, Poland; [email protected].pl *Correspondence: andrzej.gr[email protected]; Tel.: +48-12-622-4136 Abstract: This paper presents the results of laboratory tests that allowed us to determine the effect of the soil surface inclination and its initial moisture content on the formation of overland and subsurface runoff. The experiments were carried out for the soil that is commonly present in the southern part of Poland, including the Outer Carpathians. The results of these measurements served as a reference for overland runoff calculations using the Richards model, simplified Green–Ampt model, and the empirical model (MSME). The results of the measurements showed that, for low-permeable soil, overland runoff is the dominant form. It was shown that a slope in the range of 2.5–5.0% does not have a significant effect on the amount of overland runoff, but affects its dynamics. The measurements also showed that the starting time and amount of overland runoff are strictly associated with the initial soil moisture content. High soil moisture content in the period preceding the onset of rainfall causes faster generation and an increase in overland runoff, which is caused by the saturation of the surface layer of the soil. This mechanism was confirmed by the results of calculations using the Richards model and measurements of the electrical resistance of the soil. Theoretical calculations showed that the results of the runoff calculations using the Richards and Green–Ampt models are strongly dependent on the hydraulic properties of the soil adopted for the analysis. It was also demonstrated that the modified MSME model satisfactorily estimates the amount of overland and subsurface runoff, but requires parameter calibration based on existing hydrological data. Keywords: overland flow; infiltration; subsurface flow; ponding of surface; surface inclination 1. Introduction Rainfall–runoff modeling has been developed over many decades. Despite the relatively simple theoretical basis, the models are still commonly used in hydrological applications worldwide. Values of the model parameters are estimated using recorded rainfall–runoff episodes. In many parts of the world, rainfall and runoff data are seldom adequate to determine a unit hydrograph for a catchment or a watershed. In the absence of rainfall–runoff data, unit hydrographs can be derived by synthetic means [1,2]. One of the crucial problems with the practical use of rainfall–runoff models is the uncertainty of the net or excess rainfall estimation. Many rainfall–runoff models are capable of assessing direct runoff, which can be interpreted as Hortonian overland flow or overland saturated flow and subsurface flow (throughflow or interflow) [3]. Sustainability 2022,14, 4962. https://doi.org/10.3390/su14094962 https://www.mdpi.com/journal/sustainability Sustainability 2022,14, 4962 2 of 27 Examples of rainfall–runoff models that can be used for direct runoff assessing include synthetic unit hydrographs (Snyder, SCS-UH, Clark-UH) [ 4 , 5 ] or conceptual models based on the cascade of Nash linear tanks, or double cascade of tanks (Wackermann model) [ 6 ]. An alternative to the rainfall–runoff models used so far may be the recently developed “Event-based Approach for Small and Ungauged Basins” (EBA4SUB) model. It allows for estimating the magnitude of the peak flow along with the characteristics of the design hydrograph [ 7 – 9 ] or new modification, named COSMO4SUB (“Continous Simulation Model for Small and Ungauged Basins”) for continuous simulation of runoff [ 10 ]. All of the mentioned models are common tools used in hydrological studies to design hydraulic structures like culverts, reservoirs, and flood areas. Depending on the catchment conditions, runoff can be dominated by overland runoff or subsurface drainage. In natural catchments, direct runoff can be dominated by subsurface runoff, which may also occur outside the immediate vicinity of watercourses [ 11 ]. Many commonly used empirical models, like the Soil Conservation Service Curve Number ( SSC-CN ) described by the Natural Resources Conservation Service [ 4 ], are only able to assess overland runoff. Yuan et al. [ 12 ] modified the SCS-CN method to calculate subsurface runoff by defining CN values for drainage flow for five poorly drained sites in east-central Illinois, USA. Wał˛ega et al. [ 13 ] proposed a modification of the SCS-CN method, named the MSME model, to predict overland and subsurface runoff from a small, forested watershed in North Carolina, USA. They modified the Sahu–Mishra–Eldo (SME) approach [ 14 ] by including a component to assess the depth of subsurface runoff. The key role here was played by the “a” parameter—that is, the factor determining the proportion of saturated area in the whole catchment. The method was successfully tested on a flat, coastal forested watershed in the South Carolina Atlantic Coastal Plain, USA. The method, however, has not yet been tested in catchments with higher slope inclination and less permeable soils. Based on the literature review, it is clear that there is a lack of studies that would describe simple rainfall–runoff models, including overland and subsurface runoff. The first stage of studies comprises experiments on runoff formation, including infiltration processes and soil properties. A novelty of the presented studies was that the MSME model was tested for its ability to estimate the overland and subsurface runoff in slopes with low-permeable soils and varied runoff surface inclinations. The MSME model was tested in laboratory conditions, which allowed us to show the relationship between the soil parameters and the model parameters; this method can be used to analyze in detail the formation of overland flow [15–18]. Deterministic models of water flow in unsaturated soil can also be used for surface and subsurface runoff calculations. The most popular physical infiltration model is that of Richards [ 19 ], which is an extension of Darcy’s law for the unsaturated medium. The model requires information on the soil retention and hydraulic conductivity functions, so its use in hydrology can be impractical or can make hydrological calculations complicated. Therefore, simplified physical models are a good alternative to Richards’. One of the most popular physical-based infiltration models used for calculation overland flow is the Green–Ampt model [ 20 ]. This model assumes the uniform distribution of water content with the soil profile, constant value of suction pressure in the unsaturated zone, and permeability coefficient of soil. These assumptions ensure that the model is not too sophisticated and does not have a significant impact on the accuracy of the calculation results. So, the model is often used not only in hydrology [ 21 – 24 ] but also in slope stability analysis [ 25 – 27 ]. The most reliable results of surface and subsurface runoff studies are provided by field and laboratory measurements. In most cases, research on experimental plots focuses on determining the impact of plot use or dimensions on the runoff [ 28 – 33 ], but the infiltration process is rarely monitored on them [ 34 – 36 ]. On the other hand, the best way to select soil conditions, control the intensity of rainfall, or choose the instrumentation of the measuring model is based on laboratory measurements. Poesen [ 37 ] described a laboratory experiment with the aid of simulated runoff, illustrating that the surface slope inclination is negatively correlated to the Hortonian overland runoff, while the infiltration rate is dependent not Sustainability 2022,14, 4962 3 of 27 only on the soil properties but also on relief parameters such as the surface slope inclination. Laboratory tests by Chu et al. [ 17 ] showed that surface microrelief, soil type, initial soil moisture content, and rainfall characteristics affect the overland flow generation. Nassif and Wilson [ 38 ] show that, for soils of low permeability, there is little effect of rainfall intensity on infiltration. Similar results were provided by Wang et al. [ 16 ], who revealed that an increase in the inclination of soil surface increases overland flow and decreases infiltration. The formation of overland and subsurface runoff significantly affects the dynamics of floods. In Poland and other countries, this problem mainly concerns sealed surfaces, i.e., urbanized areas [ 11 , 39 – 43 ], as well as mountain areas [ 44 – 47 ]. The sealing of the area causes an increase in the number of floods with a rapid course of a flood wave, but of short duration. The formation of such floods is heightened by climate change and other human activities, which causes an imbalance between precipitation and the processes of runoff, infiltration, and transpiration. As a result of this disturbance, the volume of water infiltrating the soil decreases and overland runoff develops. In southern Poland, including mountainous areas of the Outer Carpathians, the natural substrate is relatively often formed by silty (loesses) and clay soils of medium and low water permeability. These soils are the product of weathering of shales and sandstones that were formed in the Tetyda Ocean and lifted during Alpine orogenesis. The presence of this type of formation and the high inclination of slopes favor the formation of overland runoff. Another factor contributing to the increase in overland runoff in this region is human activity related to changes in land use, i.e., deforestation, sealing, and the development and reconstruction of transport networks. All kinds of earthworks accompanying human activities remove vegetation, permanently or temporarily exposing the soil, which becomes more susceptible to the formation of overland runoff [48–50]. Still common in practical calculations of runoff are simple empirical models based on a hydrological approach. In comparison to physical-based models like the Green–Ampt or Richards, the hydrological approach uses many simplification, mainly for the description of runoff formation. In this study we tested following hypothesis: the modified SCS-CN method, name MSME model, gives similar results of runoff to compare to the commonly using Green–Ampt and Richards models. The aim of this study was to assess the influence of the slope inclination, the initial moisture content of the soil on the formation of overland, and subsurface runoff. In addition, the usefulness of the MSME method for estimating overland and subsurface runoff formed on a slope with low-permeable soils was determined. The results from the MSME model were compared to commonly used models of water infiltration in the soil, e.g., the one-dimensional Green–Ampt model and the two-dimensional model described by the Richards equation. In the tests, the infiltration process was monitored using the electric impedance spectrometry (EIS) method. 2. Materials and Methods The tests were carried out in laboratory conditions on a runoff simulator for silty soil collected from the vicinity of Kraków, Poland. The Kraków region is located on the border of several large tectonic units, which gives it a complicated and varied geological structure [ 51 ]. The northern part of the city and the part south of the Vistula River are located in the Carpathian Foredeep zone. The southern part of the city is located in the Flysch Carpathians. The eastern part of Kraków is located in the Sandomierz Basin, where the river sediments of the Vistula and its tributaries predominate. In turn, the western part of Kraków is the area of the Krakowska Gate, which is made of a limestone framework separated by narrow tectonic ditches [ 52 ]. According to the Detailed Geological Map of Poland, Part: Kraków [ 53 ], in the northern part of Kraków, where the soil samples were taken, the surface formations are represented by loesses, and their presence is also noted in the southern zone, where they cover flysch formations. Loesses in the vicinity of Krakow are characterized by a varied content of the clay fraction (5–13%) [ 54 ], which allows them to be classified as silt and silty loams according to geotechnical standards [55]. Sustainability 2022,14, 4962 4 of 27 The total annual rainfall in Kraków (Poland) is approximately 700 mm [ 56 ]. Some of this rainfall is intense, causing rapid increases in water levels in watercourses and local flooding. The scale of this phenomenon can be proved by the data from the Institute of Meteorology and Water Management–National Research Institute [ 57 ]. In May–June 2021, the Hydrological Forecasting Office in Kraków issued almost 50 hydrological warnings indicating the occurrence of storms, including urbanized catchments located in southeastern Poland. 2.1. Research on Geotechnical Properties The scope of the research included the determination of the particle size distribution, bulk density, consistency limits, compaction parameters, and the filtration coefficient. The grain size composition was determined by the combined method, i.e., wet sieve analysis for grains larger than 0.063 mm and the hydrometric method for particles smaller than 0.063 mm [ 55 , 58 ]. The liquidity limit was determined by the Casagrande method [ 59 ] as the moisture content at which a fine-grained soil no longer flows like a liquid. The plasticity limit was determined by the rolling method [ 59 ], as the moisture content at which a finegrained soil can no longer be remolded without cracking. The compaction parameters (optimum moisture content and maximum dry density) were determined in the Proctor apparatus at a compaction energy of 0.59 J · cm −3 [ 60 ]. The compaction characteristics of a soil can be assessed by means of standard laboratory tests. In the Proctor test, the volume of the mold is 1 dm 3 and the soil is compacted by a rammer consisting of a 2.5 kg mass falling freely through 320 mm. The soil is compacted in three equal layers, with each layer receiving 25 blows with the rammer for at least five soil samples. After compaction, the bulk density and water content of the soil are determined and the dry density calculated. The filtration coefficient was determined in an oedometer on samples of 7 cm in diameter and 1.9 cm high, formed in the ring of the apparatus at the optimum moisture to obtain a compaction index of I S = 0.88, 0.92, 0.95, or 1.00, which corresponds to a soil porosity n= 0.40, 0.38, 0.35, or 0.32, respectively. The determination of the filtration coefficient was also carried out with the use of a Saturo infiltrometer (dual head) (METER Group, Pullman, DC, USA) (Figure 1). In this case, the tests were performed on samples formed in the cylinder of a medium-sized Proctor apparatus (Wille Geotechnik, Germany) with a diameter of 25 cm and a height of 20 cm at a moisture content of 10% until a compaction index of I S = 0.88 (n= 0.40) was obtained. The filtration coefficient from the infiltrometer tests was determined using a ring with a diameter of 14.4 cm stuck into the soil sample to a depth of 5 cm. The test was performed in duplicate with two and three measurement cycles. Oedometer tests were carried out in a laboratory under saturated conditions. The direction of water through the soil was from bottom to top of a soil sample and the hydraulic gradient during the test was variable. A Saturo infiltrometer was used to measure the infiltration rate of the water through soil pore spaces under field conditions. The measurement took place within a single ring that was pushed into the soil. The tests were carried out at the two values of ponding head, which enabled us to reduce the impact of lateral flow on the measured infiltration rate. 2.2. Surface and Subsurface Runoff Studies The physical model was prepared in a runoff simulator (Figure 2a), which enables the simulation of rainfall with a predetermined intensity for any soil or material used for the construction of communication pavements. The basic elements of the simulator include the water reservoir, the main bed in which the soil sample is formed and the reservoirs of feed water, and water flowing out of the bed in the form of surface and subsurface runoff. The active height of the main bed of the simulator is 0.216 m, with cross section dimensions of 0.64 m × 1.54 m. On one side of the bed, there are channels collecting water flowing down the ground surface and water draining from the soil (subsurface runoff) (Figure 2b). At 0.5 m above the main channel, there is a water tank with 372 holes (nozzles) from which Sustainability 2022,14, 4962 5 of 27 water falls onto the sample surface. The hydraulic flow of water supplied to the reservoir is regulated, and for the purposes of this research it was related to the bed (soil) surface and converted into the intensity of precipitation. Before starting the tests, the pump dosing the precipitation to the soil sample was calibrated. Sustainability 2022, 14, x FOR PEER REVIEW 5 of 31 (a) (b) Figure 1. General view of Saturo infiltrometer (dual head) (a) and model of soil permeability testing by the infiltration method (b). 2.2. Surface and Subsurface Runoff Studies The physical model was prepared in a runoff simulator (Figure 2a), which enables the simulation of rainfall with a predetermined intensity for any soil or material used for the construction of communication pavements. The basic elements of the simulator include the water reservoir, the main bed in which the soil sample is formed and the reservoirs of feed water, and water flowing out of the bed in the form of surface and subsurface runoff. The active height of the main bed of the simulator is 0.216 m, with cross section dimensions of 0.64 m × 1.54 m. On one side of the bed, there are channels collecting water flowing down the ground surface and water draining from the soil (subsurface runoff) (Figure 2b). At 0.5 m above the main channel, there is a water tank with 372 holes (nozzles) from which water falls onto the sample surface. The hydraulic flow of water supplied to the reservoir is regulated, and for the purposes of this research it was related to the bed (soil) surface and converted into the intensity of precipitation. Before starting the tests, the pump dosing the precipitation to the soil sample was calibrated. Two series of tests were carried out, in which the inclination of the sample (soil) surface was 2.5% and 5.0%. Both samples were formed at the same compaction, directly in the simulator’s main bed, in five layers. The soil was compacted at 10% moisture to obtain a volumetric density of 1.75 g·cm−3, which allowed us to obtain a soil compaction index of IS = 0.87. This index corresponds to a soil porosity of n = 0.40, which falls within the range of silty soil porosity given by Kaczyński [61] for upland and mountain areas of southern Poland. For each series, three precipitation simulations were performed, during which the simulator operated for 40 min. As part of the research, 30 mm rainfall was generated for a duration of 40 min, which corresponds to a rainfall intensity of 0.75 mm·min−1. Figure 1. General view of Saturo infiltrometer (dual head) ( a ) and model of soil permeability testing by the infiltration method (b). Two series of tests were carried out, in which the inclination of the sample (soil) surface was 2.5% and 5.0%. Both samples were formed at the same compaction, directly in the simulator’s main bed, in five layers. The soil was compacted at 10% moisture to obtain a volumetric density of 1.75 g · cm −3 , which allowed us to obtain a soil compaction index of I S = 0.87. This index corresponds to a soil porosity of n= 0.40, which falls within the range of silty soil porosity given by Kaczy´nski [ 61 ] for upland and mountain areas of southern Poland. For each series, three precipitation simulations were performed, during which the simulator operated for 40 min. As part of the research, 30 mm rainfall was generated for a duration of 40 min, which corresponds to a rainfall intensity of 0.75 mm · min −1 . According to Lambor [ 62 ], for an area with an annual precipitation of 700 mm, this intensity is characterized by a probability of about 5–10%, which corresponds to the probability for which road drainage devices are designed [ 63 ]. According to the Institute of Meteorology and Water Management–National Research Institute [ 64 ], rainfall of 30 mm is considered critical, leading to rising water in water courses and the start of surface runoff. Runoff measurements were conducted during each rainfall episode every 2 min, and continued for 10–40 min thereafter, depending on the dynamics of the recorded runoff. Before the research was started and after the end of each episode, the soil moisture content was controlled by taking soil samples from two extreme points of the model. Sustainability 2022,14, 4962 6 of 27 Sustainability 2022, 14, x FOR PEER REVIEW 6 of 31 According to Lambor [62], for an area with an annual precipitation of 700 mm, this intensity is characterized by a probability of about 5–10%, which corresponds to the probability for which road drainage devices are designed [63]. According to the Institute of Meteorology and Water Management–National Research Institute [64], rainfall of 30 mm is considered critical, leading to rising water in water courses and the start of surface runoff. Runoff measurements were conducted during each rainfall episode every 2 min, and continued for 10–40 min thereafter, depending on the dynamics of the recorded runoff. Before the research was started and after the end of each episode, the soil moisture content was controlled by taking soil samples from two extreme points of the model. (a) (b) Figure 2. General view (a) and scheme (b) of the runoff simulator. 2.3. Monitoring of the Infiltration Process Using the EIS Method The EIS method used in this study was applied to measure the electrical characteristics of the soil in a complex form, i.e., the resistance Rx (electrical resistance) and reactance X in the frequency range of the supply signal. Resistance is the actual (active) resistance in an electrical circuit and is the so-called electrical resistance, which depends on the soil type and the water content in its pores. On the other hand, reactance is the imaginary part of the impedance. A nonzero value of reactance causes a phase shift (phase difference) between the intensity and voltage of the electric current in the circuit. Reactance characterizes the changes in the soil caused by variations in its grain size or porosity [65]. The resistance measurements of the soil built into the model were carried out with two five-channel EIS measuring electrodes with a total length of 22.5 cm (Figure 3b) at a signal frequency of 2050 Hz. The electrodes were placed in the central part of the sample with a spacing of 20 cm. The measuring electrodes were paired up to a measuring probe with five vertical measuring segments (channels), 2.5 cm high each (current-conducting sections), spaced vertically every 2.5 cm. Due to the height of the soil sample, four measurement channels were used in the experiment. A Z-Meter (GEOtest, Brno, Czech Republic) was used to record the electrical resistance of the soil embedded in the model during water infiltration (Figure 3a). The device uses a measurement method that consists of comparing the measured impedance Z with normal resistance Rx with a known value of electrical resistance. The measurement results were recorded in a Microsoft Excel spreadsheet. Figure 2. General view (a) and scheme (b) of the runoff simulator. 2.3. Monitoring of the Infiltration Process Using the EIS Method The EIS method used in this study was applied to measure the electrical characteristics of the soil in a complex form, i.e., the resistance Rx (electrical resistance) and reactance X in the frequency range of the supply signal. Resistance is the actual (active) resistance in an electrical circuit and is the so-called electrical resistance, which depends on the soil type and the water content in its pores. On the other hand, reactance is the imaginary part of the impedance. A nonzero value of reactance causes a phase shift (phase difference) between the intensity and voltage of the electric current in the circuit. Reactance characterizes the changes in the soil caused by variations in its grain size or porosity [65]. The resistance measurements of the soil built into the model were carried out with two five-channel EIS measuring electrodes with a total length of 22.5 cm (Figure 3b) at a signal frequency of 2050 Hz. The electrodes were placed in the central part of the sample with a spacing of 20 cm. The measuring electrodes were paired up to a measuring probe with five vertical measuring segments (channels), 2.5 cm high each (current-conducting sections), spaced vertically every 2.5 cm. Due to the height of the soil sample, four measurement channels were used in the experiment. A Z-Meter (GEOtest, Brno, Czech Republic) was used to record the electrical resistance of the soil embedded in the model during water infiltration (Figure 3a). The device uses a measurement method that consists of comparing the measured impedance Z with normal resistance Rx with a known value of electrical resistance. The measurement results were recorded in a Microsoft Excel spreadsheet. 2.4. Calculations of Overland and Subsurface Runoff with the Use of Models Taking into Account the Infiltration Process The obtained results of overland and subsurface runoff measurements were compared with the results of calculations carried out using the one-dimensional Green–Ampt model and the two-dimensional soil model developed in the GeoStudio 2020 software (Seequent, a Bentley Company, Broomfield, CO, USA), in which the Richards [ 19 ] equation was adopted to calculate the flow of water in the soil: ∂θ ∂t=∂ ∂xkx∂H ∂x+∂ ∂zkz∂H ∂z+Q, (1) Sustainability 2022,14, 4962 7 of 27 where θ (-) is the volumetric water content, t(s) is time, k x (m · s −1 ) is the hydraulic conductivity of the soil in the x-direction, k z (m · s −1 ) is the hydraulic conductivity of the soil in the z-direction, H(m) is the water total head, and Q(m·s−1) is the boundary flux. Sustainability 2022, 14, x FOR PEER REVIEW 7 of 31 (a) (b) Figure 3. View of the Z-Meter IV apparatus (a) and measuring electrodes (b). 2.4. Calculations of Overland and Subsurface Runoff with the Use of Models Taking into Account the Infiltration Process The obtained results of overland and subsurface runoff measurements were compared with the results of calculations carried out using the one-dimensional Green–Ampt model and the two-dimensional soil model developed in the GeoStudio 2020 software (Seequent, a Bentley Company, Broomfield, CO, USA), in which the Richards [19] equation was adopted to calculate the flow of water in the soil:  =  +  + , (1) where θ (-) is the volumetric water content, t (s) is time, kx (m·s−1) is the hydraulic conductivity of the soil in the x-direction, kz (m·s−1) is the hydraulic conductivity of the soil in the z-direction, H (m) is the water total head, and Q (m·s−1) is the boundary flux. The calculation of the infiltration by the Green–Ampt method [20] consists of determining the soil capacity rate (potential infiltration rate): =∙=∙++ [·] (2) where ks (m·s−1) is the hydraulic conductivity of a soil, given as half of the value of the soil coefficient of permeability [66]; i (-) is the hydraulic gradient; Hp (m) is the height of ponding, often assumed to be equal to 0, because surface ponding triggers overland runoff; zf (m) is the depth of the wetting front location; and ψf (m) is the weight of soil suction pressure at the base of the wetting front. The height of the soil suction pressure at the base of the wetting front can be determined in various ways [3,67]. In this study, the formula given by Maidmend [3] was used, as shown in Equation (3): Figure 3. View of the Z-Meter IV apparatus (a) and measuring electrodes (b). The calculation of the infiltration by the Green–Ampt method [ 20 ] consists of determining the soil capacity rate (potential infiltration rate): f=ks·i=ks·zf+ψf+Hp zfhm·s−1i(2) where k s (m · s −1 ) is the hydraulic conductivity of a soil, given as half of the value of the soil coefficient of permeability [ 66 ]; i(-) is the hydraulic gradient; H p (m) is the height of ponding, often assumed to be equal to 0, because surface ponding triggers overland runoff; z f (m) is the depth of the wetting front location; and ψf (m) is the weight of soil suction pressure at the base of the wetting front. The height of the soil suction pressure at the base of the wetting front can be determined in various ways [ 3 , 67 ]. In this study, the formula given by Maidmend [ 3 ] was used, as shown in Equation (3): ψf=0.01·6.53 −7.326·n+0.00158Cl2+3.809·n2+0.000344·Sa·Cl −0.04989·Sa·n+0.0016 ·Sa2·n2+0.0016·Cl2·n2−0.0000136·Sa2·Cl −0.00348·Cl2·n−0.000799·Sa2 ·n] [cm] (3) where n(-) is the soil porosity, S a (%) is the sand fraction content, and Cl (%) is the clay fraction content. The equation is easy to calculate and there is no need for information on the soil water characteristic curve (retention curve). The above equation assumes that the value of the suction pressure at the base of the wetting front is independent of the soil moisture content. Therefore, for the purposes of the analyses, Equation (3) was modified to take into account the soil moisture content: ψf m =(1−Sr)·ψf[cm], (4) where Sr(-) is the degree of soil saturation. In the calculations, the infiltration capacity is compared with the intensity of precipitation, assuming that, until the values of both parameters are equal, the total amount of Sustainability 2022,14, 4962 8 of 27 rainwater infiltrates the soil profile. This is referred to as the surface flooding time (t p ). After this, the overland runoff begins, and the amount of water accumulated (F) in the soil profile can be calculated from the following relationship: F=Fp+ks·t−tp+(θs−θi)·ψf m·ln (θs−θi)·ψf m +F (θs−θi)·ψf m +Fp , (5) where F p (mm) is the accumulation of rainwater in the soil at the moment of ponding of the soil surface, θs (-) is the soil moisture content at full saturation, and θi (-) is the soil moisture content before rainfall. The solution to Equation (5) was found by using the iterative method in Microsoft Excel. The values of initial soil moisture before the first rainfall episode were determined on the basis of the soil porosity and mass moisture, while for the next two episodes the balance of water infiltrating the soil and the amount of overland and subsurface runoff were taken into consideration. The second method of calculating overland runoff was based on numerical calculations using the physical Richards equation [ 19 ] of water flow in the zone of incomplete saturation. The calculations were conducted in the SEEP/W module (Seequent, a Bentley Company, Broomfield, CO, USA) for the two-dimensional model of the plot (Figure 4), assuming that the outflow of water from the ground occurs in the lower part of one of the walls. On the other hand, the amount of rainfall, its duration, and the periods between individual rainfall episodes were the same as in the actual studies. As a result of the calculations, the volume of water stored and drained from the soil was obtained, and the amount of overland runoff was calculated as the difference between the rainfall and the amount of water infiltrating the soil: OF =R−S, (6) where OF (mm) is the overland runoff, R(mm) is the rainfall, and S(mm) is the height of the water column stored in the ground. Sustainability 2022, 14, x FOR PEER REVIEW 9 of 31 Figure 4. Cross section of the calculation model scheme (5.0% soil surface inclination). For the purpose of calculating the coefficient of permeability of soil, values of (2 ÷ 5) × 10−6 m·s−1 were adopted from the Saturo infiltrometer tests, and the retention characteristics were determined using the retention parameters proposed by the SEEP/W program for silty soils. 2.5. Estimation of Overland and Subsurface Runoff Using the MSME Model The concept of the model was introduced by Wałęga and Amatya [68] for a flat coastal watershed. The method is based on a common SCS-CN method used to calculate direct runoff. In the present study, the concept was slightly modified. In the MSME model, the direct runoff was estimated for every event using the following equations: =()∙()  > ; (7) Supposing =0 then: =∙(−) , (8) where S is the maximum retention capacity, and P (mm) is the sum of the precipitation during the event. In Equation (7), the higher the antecedent moisture (M), the lower the initial abstraction (Ia), and vice versa. The antecedent moisture content is given as follows: =(5−) (5−)+5> (9) =05< , (10) where β (-) and λ (-) are parameters that are optimized, and P5 is the amount of water that infiltrated during the previous five-day period. Subsurface runoff (MSMEQsubs) was calculated using the following equations:  =(−)∙(−+) (−+)> (11) 1=∙5+(−) (12)  =0< (13) Figure 4. Cross section of the calculation model scheme (5.0% soil surface inclination). For the purpose of calculating the coefficient of permeability of soil, values of (2 ÷5) ×10−6m·s−1 were adopted from the Saturo infiltrometer tests, and the retention characteristics were determined using the retention parameters proposed by the SEEP/W program for silty soils. 2.5. Estimation of Overland and Subsurface Runoff Using the MSME Model The concept of the model was introduced by Wał˛ega and Amatya [ 68 ] for a flat coastal watershed. The method is based on a common SCS-CN method used to calculate direct runoff. In the present study, the concept was slightly modified. In the MSME model, the direct runoff was estimated for every event using the following equations: Q=(P−Ia)·(P−Ia+M) P−Ia+Si f P >Ia; (7) Sustainability 2022,14, 4962 9 of 27 Supposing Q=0 then: Ia=λ·(S−M), (8) where Sis the maximum retention capacity, and P(mm) is the sum of the precipitation during the event. In Equation (7), the higher the antecedent moisture (M), the lower the initial abstraction (Ia), and vice versa. The antecedent moisture content is given as follows: M=β(P5−λS)λS (P5−λS)+Sf or P5>S(9) M=0f or P5<S, (10) where β (-) and λ (-) are parameters that are optimized, and P5 is the amount of water that infiltrated during the previous five-day period. Subsurface runoff (MSMEQsubs) was calculated using the following equations: MSMEQsubs =(P−Ia1)·(P−Ia1+M) (P−Ia1+Sa)i f P >Ia1(11) Ia1=a·(PET5+(Sa−M)) (12) MSMEQsubs =0i f P <Ia1(13) M=β·(P5−PET5)·Sa (P5−PET5)+Sa(14) Sa=a25400 CN −254, (15) where I a1 (mm) is the initial abstraction for subsurface runoff, PET5 (mm) is the sum of five days’ potential evapotranspiration computed using the Penman–Monteith method, S a (mm) is the maximum potential retention for the area where subsurface runoff occurs, and a(-) is the coefficient for the proportion of area with saturated soil. Overland runoff (MSMEQsurf) was calculated using the following equations: MSMEQsur f =(P−Ia2)·(P−Ia2) P−Ia2+Sb i f P >Ia2(16) Ia2=(1−a)·(PET5+Sb)(17) MSMEQsur f =0i f P <Ia2(18) Sb=(1−a)·25400 100 −CN −254, (19) where I a2 (mm) is the initial abstraction for overland runoff, S b (mm) is the maximum potential retention for the area where overland runoff occurs, and CN (-) is the curve number (CN) parameter. Total runoff was calculated as the sum of overland and subsurface runoff: MSMEtot =MSMEQsur f +MSMEQsubs. (20) In the initial phase of rainfall, when the soil is dry with a large soil water deficit (storage), it is assumed that runoff does not occur. In the case of highly impermeable soil, after the first soil saturation threshold is exceeded, overland runoff begins. The factor determining the proportion of saturated area is parameter “a”, which is calibrated based on the observed rainfall and runoff from the laboratory model. If the soil is almost fully saturated ( Sr∼ = 1.0), subsurface runoff may occur after overland runoff because the wetting front is moving in the soil profile and saturation increases. If the rainfall continues, the soil will continue to be saturated, achieving the second threshold and resulting in subsurface Sustainability 2022,14, 4962 16 of 27 were carried out using two calculation methods, taking into account the physical properties of the soil (Table 2). The results of the calculations show a significant influence of the coefficient of permeability on the results of the overland and subsurface runoff calculations. Generally, for lower values of the coefficient of permeability, higher values of overland runoff and lower values of subsurface runoff were obtained. The exemplary results of overland runoff calculations using the Green–Ampt method (Figure 12) indicate that the greatest discrepancies between the calculation and observation results were obtained in the first rainfall episode, where the calculated overland runoff value was significantly greater than the calculated value. The root mean square error and the Nash–Sutcliffe coefficient showed that, in the case of overland runoff calculations, the most accurate values of both parameters were obtained with a permeability coefficient of 2.10 × 10 −6 m · s −1 . The value of the EF parameter (Table 3) for the model with a 2.5% inclination was 0.77–0.80 in the Green–Ampt and Richards method, respectively, and 0.95–0.96 for the model with a 5.0% inclination of the ground surface. On the other hand, acceptable EF values (0.64–0.79) were obtained for the model with a surface inclination of 5% for the two highest values of the coefficient of permeability. This is taking into account that, in the analyzed cases, overland runoff dominated, so the results of the calculations were used for further analyses, with a coefficient of permeability of 2.10 ×10−6m·s−1. Sustainability 2022, 14, x FOR PEER REVIEW 17 of 31 Figure 10. Changes in soil resistance during three rainfall episodes at a 2.5% inclination of the soil surface. Figure 11. Changes in soil resistance during three rainfall episodes at a 5.0% inclination of the soil surface. 3.4. Verification of Calculation Models 3.4.1. Green–Ampt and Richards models The research results confirmed the effect of individual rainfall episodes on the amount of overland runoff generated. For comparison, theoretical calculations of water infiltration were carried out using two calculation methods, taking into account the physical properties of the soil (Table 2). The results of the calculations show a significant Figure 11. Changes in soil resistance during three rainfall episodes at a 5.0% inclination of the soil surface. Figures 13a and 14a compare the results of numerical calculations of overland and subsurface runoff from the Richards model with the results of laboratory tests. The calculations were carried out as a single time series, which started at least ten minutes before the beginning of the first rainfall episode and ended several dozen minutes after the end of the last rainfall episode. The time intervals between the individual rainfall episodes were identical to those used in the laboratory experiments. The calculation results for both inclinations of the runoff soil surface showed that the total values of overland and subsurface runoff were slightly lower than those obtained from the tests. It could be noticed that the overland runoff values calculated for the first rainfall episode were overestimated compared to those measured in laboratory conditions. On the other hand, the values calculated for the second rainfall episode were similar, while those for the third episode were underestimated compared to the results obtained from the measurements. In the case of the last rainfall episode, the measured overland runoff was 27 mm, representing almost 90% of the precipitation. Sustainability 2022,14, 4962 17 of 27 Table 2. Summary of measurement and calculation results for overland and subsurface runoff. Number of the Rainfall Episode Type of Runoff Soil Coefficient of Permeability Used in the Calculations (m·s−1) Inclination of Soil Surface 2.5% 5.0% Observations Model Observations Model Green– Ampt Richards Green– Ampt Richards Runoff Value (mm) 1overland 4.0 ×10−6 2.0 ×10−6 1.0 ×10−6 5.0 ×10−7 0.61 0.00 3.02 9.76 15.49 0.96 1.32 6.85 12.96 2.03 0.00 3.02 9.76 15.49 0.20 0.88 6.76 12.91 subsurface 4.0 ×10−6 2.0 ×10−6 1.0 ×10−6 5.0 ×10-7 0.00 - 0.00 0.00 0.00 0.00 0.0 - 0.00 0.00 0.00 0.00 2overland 4.0 ×10−6 2.0 ×10−6 1.0 ×10−6 5.0 ×10−7 16.11 6.83 14.24 19.37 22.78 7.38 8.33 14.46 19.32 14.15 6.10 13.59 18.87 22.41 8.25 8.95 15.27 19.91 subsurface 4.0 ×10−6 2.0 ×10−6 1.0 ×10−6 5.0 ×10−7 0.39 - 0.19 0.00 0.00 0.00 0.0 - 0.56 0.19 0.00 0.00 3overland 4.0 ×10−6 2.0 ×10−6 1.0 ×10−6 5.0 ×10−7 25.17 11.23 17.94 22.14 24.78 28.25 28.36 18.14 18.01 26.78 14.61 20.60 24.10 26.20 28.11 28.31 17.88 18.00 subsurface 4.0 ×10−6 2.0 ×10−6 1.0 ×10−6 5.0 ×10−7 0.75 - 0.90 0.77 0.00 0.00 1.53 - 1.42 0.77 0.00 0.00 Total overland 4.0 ×10−6 2.0 ×10−6 1.0 ×10−6 5.0 ×10−7 43.72 18.1 35.2 51.3 63.1 36.59 38.00 39.46 50.29 42.96 20.71 37.21 52.73 64.10 36.56 38.14 39.91 50.82 subsurface 4.0 ×10−6 2.0 ×10−6 1.0 ×10−6 5.0 ×10−7 5.20 - 1.49 0.91 0.00 0.00 2.66 - 2.42 0.78 0.00 0.00 3.4.2. MSME Model (Verification of Suitability for the Total Runoff Estimation) Table 4shows the mean values of the overland and subsurface runoff obtained based on the results of laboratory tests and calculations using the MSME model. The test results showed that the subsurface runoff in the tested soil was small and constituted just over 3.3% and 3.7% of the total runoff for the 2.5% and 5.0% soil surface inclination, respectively. The amount of overland and subsurface runoff from both the observation and the MSME model increased with the inclination of the runoff surface, which is consistent with the study of Chen et al. [ 26 ]. The results obtained from the MSME model were similar to the observations, regardless of the slope of the runoff surface. The RMSE error value was 0.52 and 0.05 mm for the 2.5% and 5.0% slope, respectively, while the EF coefficient was 0.99 in both cases. Therefore, it can be concluded that the MSME model allowed for a very good estimation of the total runoff, which proves the correct concept of runoff formation described by this model. The inclination of the runoff surface affects the soil degree of saturation, expressed by the parameter “a”, which for the 2.5% inclination of the runoff surface was 0.84, while for 5.0% it was 0.90. Thus, in the case of a surface with higher inclination, the subsurface runoff may be formed at a higher soil saturation. However, the differences in soil moisture content are small due to the fact that, during heavy rainfall, Sustainability 2022,14, 4962 18 of 27 the infiltration rate depends on the runoff surface slope under conditions where it exceeds 10 ◦ [ 26 ]. The value of the CN parameter, understood as the catchment’s potential for runoff formation, is as high as 89.7 for a 2.5% surface slope and 93.0 for a 5% slope, which is understandable as the investigated land is characterized by a low infiltration capacity. At the same time, it should be noted that the inclination of the runoff surface affects the value of the CN parameter. The value of the M parameter related to the initial water content in the soil profile, coming from the infiltration of rainwater before the runoff occurred, is 0.00 for both variants of the runoff surface slope. This is caused by the methodology used to prepare the soil for testing, which is characterized by a relatively low moisture content. The values of S a and S b retention parameters were significantly lower for soil with a 5.0% slope in comparison to a 2.5% slope. It is related, as in the case of the CN parameter, to a lower water retention capacity in soils with a higher slope of the runoff surface, which was confirmed in this study. As a consequence, the initial losses of I a1 parameter for subsurface runoff and Ia2parameter for overland runoff were also reduced. Sustainability 2022, 14, x FOR PEER REVIEW 19 of 31 Figure 12. Dependence of overland runoff on the duration of a rainfall episode: comparison of the results of measurements and calculations using the Green–Ampt method. Figure 12. Dependence of overland runoff on the duration of a rainfall episode: comparison of the results of measurements and calculations using the Green–Ampt method. Sustainability 2022,14, 4962 19 of 27 Table 3. Values of metrics—RMSE and Nash–Sutcliffe EF [70]. Type of Runoff Soil Coefficient of Permeability Used in the Calculations (m.s−1) Inclination of Soil Surface 2.5% 5.0% Model Model Green–Ampt Richards Green–Ampt Richards Root Mean Square Error, RMSE (mm) (Equation (21)) overland 4.0 ×10−6 2.0 ×10−6 1.0 ×10−6 5.0 ×10−7 162.15 35.38 59.77 153.62 49.50 41.10 52.57 123.60 125.29 22.79 51.50 144.27 23.07 17.70 59.35 132.08 subsurface 4.0 ×10−6 2.0 ×10−6 1.0 ×10−6 5.0 ×10−7 - 0.25 0.30 0.63 0.63 - 0.19 0.33 1.35 1.35 Modeling efficiency, EF [-] (Equation (22)) overland 4.0 ×10−6 2.0 ×10−6 1.0 ×10−6 5.0 ×10−7 0.09 0.80 0.66 0.14 0.72 0.77 0.70 0.31 0.74 0.95 0.89 0.70 0.95 0.96 0.87 0.72 subsurface 4.0 ×10−6 2.0 ×10−6 1.0 ×10−6 5.0 ×10−7 - −0.54 −0.86 −2.86 −2.86 - 0.79 0.64 −0.50 −0.50 3.5. Discussion When analyzing the influence of the surface inclination on the measured values of overland runoff, it can be concluded that it was not very significant. These results are consistent with the results of the analyses by Chen and Young [ 26 ], who showed that, at soil inclinations smaller than 10 ◦ , the differences in the amount of water infiltrating the soil are small, so the values of overland runoff will be similar. Similar results were generated in the field tests conducted by Smolska [ 32 ]. Laboratory tests carried out by Wang et al. [ 16 ] provide slightly different relationships, as they indicate that there are differences between the measured values of effective precipitation and the amount of water infiltrating the soil for surface inclinations of 3 ◦ and 5 ◦ . In general, the amount of overland runoff increases with increasing surface inclination, while the amount of water infiltrating the soil decreases. Unfortunately, Wang et al. [ 16 ] do not provide the physical properties of the soil they examined. Interesting data were generated by numerical calculation. The results of soil degree of saturation calculations (Figures 13b and 14b) indicate that, during the first rainfall episode, there were significant changes in the surface part of the soil sample, regardless of the soil slope. In the lower part of the sample, at a depth of 15.5 cm, these changes were very small and appeared at the end of the precipitation event. This seems to be consistent with the results of measurements of the electrical resistance of the soil (Figures 10 and 11). The results of calculations of soil moisture content in the upper part of the soil sample indicate that, in the final phase of the first episode, the soil becomes saturated, which coincides with the beginning of the overland runoff. This relationship indicates that the water flow in the soil is consistent with the assumptions of infiltration calculations using piston models, in which the wetting front formed in the upper part of the profile moves downwards with the existing water supply (precipitation) to the soil. In the second rainfall episode, the saturation of the surface layer of soil occurred in the initial period of the rainfall, similar to when the overland runoff began. On the other hand, the soil moisture content in the lower part of the soil sample systematically increases until Sustainability 2022,14, 4962 20 of 27 soil saturation was achieved, at the end of the rainfall episode. The calculations show that, after the episode ends, the subsurface runoff starts, which is consistent with the test results obtained for the model with a runoff slope of 2.5%. The values resulting from the electrical resistance measurements in the lower part of the sample with a slope of 2.5% decreased noticeably at the end of this episode, to a value that was maintained until the end of the study. This dependence may indicate the saturation of the lower part of the sample. In the case of the soil sample with a 5.0% slope, the resistance value systematically decreased, but the minimum value was achieved only during the last rainfall episode. Sustainability 2022, 14, x FOR PEER REVIEW 22 of 31 Figure 13. The results of overland and subsurface runoff calculations for the model with a soil surface inclination of 2.5%. Figure 13. The results of overland and subsurface runoff calculations for the model with a soil surface inclination of 2.5%. The numerical calculations indicate that, after the second rainfall episode, the soil was fully or nearly saturated in virtually the entire sample profile. Therefore, during the third episode, most of the rainfall was transformed into overland runoff. When analyzing the results of the numerical calculations, it can be noticed that, before the beginning of the rainfall, the lower part of the sample was characterized by a higher moisture content than its upper part. This observation seems to be consistent with the results of soil resistance measurements, in which higher resistance values as well as lower soil moisture content were recorded closer to the surface. The numerical calculations indicate that, in the third Sustainability 2022,14, 4962 21 of 27 rainfall episode, the infiltration of rainwater does not have the nature of a piston movement, but causes the supply of water to the lower, saturated zone of the soil. Sustainability 2022, 14, x FOR PEER REVIEW 23 of 31 Figure 14. The results of overland and subsurface runoff calculations for the model with soil surface inclination of 5.0%. The calculation results for both inclinations of the runoff soil surface showed that the total values of overland and subsurface runoff were slightly lower than those obtained from the tests. It could be noticed that the overland runoff values calculated for the first rainfall episode were overestimated compared to those measured in laboratory conditions. On the other hand, the values calculated for the second rainfall episode were similar, while those for the third episode were underestimated compared to the results obtained from the measurements. In the case of the last rainfall episode, the measured overland runoff was 27 mm, representing almost 90% of the precipitation. 3.4.2. MSME model (verification of suitability for the total runoff estimation) Table 4 shows the mean values of the overland and subsurface runoff obtained based on the results of laboratory tests and calculations using the MSME model. The test results showed that the subsurface runoff in the tested soil was small and constituted just over 3.3% and 3.7% of the total runoff for the 2.5% and 5.0% soil surface inclination, respectively. The amount of overland and subsurface runoff from both the observation and the MSME model increased with the inclination of the runoff surface, which is consistent with the study of Chen et al. [26]. The results obtained from the MSME model were similar to the observations, regardless of the slope of the runoff surface. The RMSE error value was 0.52 and 0.05 mm for the 2.5% and 5.0% slope, respectively, while the EF coefficient was 0.99 in both cases. Therefore, it can be concluded that the MSME model allowed for a very Figure 14. The results of overland and subsurface runoff calculations for the model with soil surface inclination of 5.0%. Table 4. Average runoff characteristics from the experimental laboratory model and MSME model. Slope Inclination Qsurfobs Qsubsurfobs Qtotobs Qsurcalc Qsubsurfcalc Qtotcalc CN M SaSbIa1Ia2 - mm - mm 2.5 10.8 0.40 11.16 11.19 0.21 11.40 89.7 0.00 27.1 787.1 22.9 122.6 5.0 12.1 0.47 12.53 11.91 0.59 12.51 93.0 0.00 18.2 654.7 16.5 62.2 Notes: Q surfobs is the observed overland runoff, Q subsurfobs is the observed subsurface runoff, Q totobs is the observed total runoff (sum of overland and subsurface runoff), Q surfcalc is the calculated observed overland runoff, Q subsurfcal is the calculated observed subsurface runoff, Q totcalc is the calculated total runoff (sum of overland and subsurface runoff), I a1 is the initial abstraction for subsurface runoff, Mis the antecedent moisture content, S a is the maximum potential retention for the area where subsurface runoff occurs, I a2 is the initial abstraction for overland runoff (in mm), S b is the maximum potential retention for the area where overland runoff occurs (mm), and CN is the calculated curve number. It should also be noted that, in the third episode, both the measured and calculated values of the subsurface runoff were higher than in previous episodes. On the other hand, the calculated values of subsurface runoff were clearly smaller than the measured results. It seems that the differences may be significantly influenced by the accuracy of the preparation of the soil sample in the critical places, i.e., where the soil is exposed by the grooves capturing the water flowing out of the soil. Sustainability 2022,14, 4962 22 of 27 The numerical calculation of infiltration and overland runoff, which were simulated for laboratory runoff apparatus, was presented by Mendes et al. [ 34 ]. Analysis results indicate that the initiation of runoff generation is significantly related to the saturation of the soil sample in its upper part. They also stated that the simulated runoff results are insensitive to the length of the sample length, which differs from some site measurement results (e.g., [ 29 , 30 , 82 ]), which found that higher runoff coefficients were found on small plots than on large plots. The results of the experiments conducted indicate that, during the rainfall, the amount of overland runoff for the model with a soil surface inclination of 5.0% was slightly higher. On the other hand, the total amount of rainfall was similar (Table 2). These results may prove the influence of the soil surface inclination on the delay in the outflow of water from the catchment, as confirmed by Wang et al. [16]. The obtained results also prove the significant influence of the soil moisture content on the conditions of formation and the amount of overland and subsurface runoff. The occurrence of a series of short-term and intense rainfall episodes favors an increase in the soil moisture content and a reduction in its retention capacity. Similar relationships are observed in natural conditions (e.g., [ 36 , 48 ]) and when taking into account the climate change observed over many years, e.g., short and intense rainfall in Central Europe [ 83 ] and other regions [ 84 , 85 ]. They generate overland runoff, often causing flooding, usually of a local nature. The conducted experiments showed that in extreme cases more than two-thirds of rainfall is transformed into overland runoff. Similar dependencies in natural conditions are presented, among others, by Kijowska-Strugała and Kiska [ 79 ], Scherrer et al. [ 36 ], and Zydro´n et al. [ 86 ]. It should be emphasized that the formation of this process is greatly influenced by the soil permeability coefficient, and, in built-up areas, by the degree of surface sealing. In the presented studies, the precipitation intensity was, on average, 0.75 mm · min −1 , and the value of the coefficient of permeability was estimated at 0.12 mm · min −1 . In the case of less permeable soil, one should expect a higher value of overland runoff. It should be noted that the results of laboratory tests concerned soil without vegetation cover, which may correspond to areas covered by construction work or ground communication routes. Similar conditions may also apply to bare soil. Measurements on bare soil [ 31 , 32 ] indicate that overland runoff on this kind of surface is higher than that on cultivated land or soil with vegetation cover. The research carried out by Mounirou et al. [ 43 ] showed the effect of land use and land cover changes on the water cycle in a small Sahelian watershed. Simulation results provided an explanation of the Sahelian area paradox. The authors proved that a change in the sealing of the catchment area could lead to an increase in surface runoff despite the drop in rainfall. Results of runoff from the Green–Ampt and Richards equations were compared with the hydrological model MSME. The quality of the MSME model, despite its considerable simplicity compared to the physical models, was better than that of the Richards and Green–Ampt models. Therefore, the MSME model correctly describes the total runoff and its individual parts in the case of low-permeable soils. In the MSME model, the soil degree of saturation parameter “a” plays an essential role in separating subsurface and overland runoff. The value of the “a” parameter depends on the soil inclination. The values for the “a” parameter obtained in this study were much lower compared to in forest catchments with gentle slopes, as described by Wał˛ega et al. [ 13 ]. However, it should be taken into consideration that, in the case of a forest catchment, midsoil runoff is predominant, while overland runoff is formed only after very intense rainfall. In the original methodology [ 4 ], the runoff surface slope is not taken into account in the calculation of the CN parameter; however, as the slope of the runoff surface increases, the runoff rate increases, and therefore the infiltration capacity decreases. A similar tendency was observed in the present study, i.e., that the water infiltration capacity was reduced with an increase in the runoff surface inclination. This was manifested by an increase in the value of the CN parameter. The impact of an increase in the inclination of the runoff surface Sustainability 2022,14, 4962 23 of 27 on the value of the CN parameter has been the subject of many studies. For example, Mły´nski [ 87 ] showed that taking into account the slope of a catchment area improves the quality of adjusting the CN to the observed values. The obtained results may be a good indicator of areas subject to high human pressure. The obtained results indicate that high-intensity overland runoff is accompanied by soil denudation processes, and therefore, in order to reduce the intensity of these processes, it seems reasonable to use appropriate biogeotechnical treatments. 4. Conclusions Based on the results of tests conducted on the coarse clay silt, it can be concluded that overland runoff is the dominant type of runoff in the analyzed soil. The impact of subsurface runoff on the water balance will be small due to the low water permeability of the examined soil. This may have a significant impact on the formation of a flood wave in a watercourse in areas where there are soils with similar geotechnical characteristics, in particular in terms of filtration. Our tests proved that an increase in the soil surface inclination from 2.5% to 5.0% did not cause a significant change in the amount of overland runoff. However, increased dynamics of its formation was found, which was certainly influenced by the lack of a vegetation cover that would reduce and delay its formation. The obtained results indicate that, at a low initial soil moisture content, the runoff is of Horton overland flow character, in which the subsurface layer of soil is mostly saturated, and in the case of high moisture content it is caused by the saturation of the entire soil profile. This is confirmed by the results of electrical resistance measurements and numerical calculations. It was demonstrated that intensive rainfall repeating in short intervals will contribute to the generation of significant overland runoff, and in extreme cases more than two-thirds of rainfall can be transposed into overland runoff. It was also demonstrated that the Green–Ampt and Richards models correctly describe the course of the infiltration process and overland runoff formation. The Richards model, however, tended to underestimate the amount of subsurface runoff. The filtration coefficient adopted for the calculations has a very significant influence on the results of the calculations, especially in the case of the Green–Ampt model. The MSME model can be used to estimate the overland runoff (Qsurf) and subsurface runoff (Q subsr ) of impermeable soils, and its quality is comparable to that of the Green– Ampt and Richards models. The Green–Ampt and Richards method are based on physical assumptions and soil properties, but the simple empirical model MSME has parameters that are optimized based on observed rainfall–runoff events. This may be the main reason why the results of the MSME model are slightly better than those achieved by physical approaches. A limitation of the MMSE model in relation to the physical models is the inability to determine the course of infiltration over time and the need to calibrate the parameters based on the observed rainfall–runoff episodes. However, since it is a simple empirical model, it can be used by hydrologists to estimate the runoff in catchments. So, it is necessary to conduct further research on the application of this model on a real scale in various catchments to prove its usefulness in engineering calculations. Also, the authors will consider performing a study on the influence of artificial drainage on runoff formation and would include these results in the analyzed models. In a real catchment, the use of modern technique, like satellite imagines will be considered to indirectly assess soil moisture and link them with the parameters of the MSME model. Author Contributions: Conceptualization, A.G., T.Z. and A.W.; methodology, A.G., T.Z. and A.W.; software, A.G., T.Z., A.W. and J.P.; validation, A.G., T.Z., A.W. and J.S.; formal analysis, A.G., T.Z., A.W. and J.P.; investigation, A.G., T.Z., A.W. and J.P.; resources, A.G., T.Z. and A.W.; data curation, A.G., T.Z. and J.P.; writing—original draft preparation, A.G., T.Z. and A.W.; writing—review and editing, A.G., T.Z. and A.W.; visualization, T.Z.; supervision, T.Z. and J.S. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Sustainability 2022,14, 4962 24 of 27 Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Not applicable. Acknowledgments: The authors would like to thank Adam Kapała for his help carrying out the laboratory tests presented in this paper. The authors would like to thank the anonymous reviewers for their efforts towards improving the manuscript. Conflicts of Interest: The authors declare no conflict of interest. References 1. Egiazarova, D.; Kordzakhia, M.; Wał˛ega, A.; Dro˙ zd˙ zal, E.; Milczarek, M.; Radecka, A. Application of Polish experience in the implementation of the flood directive in Georgia—Hydrological calculations. Acta Sci. Pol. Form. Circumiectus 2017 ,16, 89–110. [CrossRef] 2. G ˛adek, W.; Bodziony, M. The hydrological model and formula for determining the hypothetical flood wave volume in non-gauged basins. Meteorol. Hydrol. Water Manag. 2015,3, 3–10. [CrossRef] 3. Maidmend, D.R. Handbook of Hydrology; CRC Press: Boca Raton, FL, USA, 1993. 4. USDA Natural Resources Conservation Service. Hydrology. In National Engineering Handbook; Chapter 10; USDA Soil Conservation Service: Washington, DC, USA, 2004. 5. Karabová, B.; Sikorska, A.; Banasik, K.; Kohnová, S. Parameters determination of a conceptual rainfall-runoff model for a small catchment in Carpathians. Ann. Wars. Univ. Life Sci.-SGGW. Land Reclam. 2012,44, 155–162. [CrossRef] 6. Mły´nski, D.; Wałega, A.; Ksi ˛a˙ zek, L.; Florek, J.; Petroselli, A. Possibility of using selected rainfall-runoff models for determining the design hydrograph in mountainous catchments: A case study in poland. Water 2020,12, 1450. [CrossRef] 7. Grimaldi, S.; Petroselli, A. Do we still need the rational formula? An alternative empirical procedure for peak discharge estimation in small and ungauged basins. Hydrol. Sci. J. 2015,60, 67–77. [CrossRef] 8. Piscopia, R.; Petroselli, A.; Grimaldi, S. A software package for the prediction of design flood hydrograph in small and ungauged basins. J. Agric. Eng. 2015,432, 74–84. [CrossRef] 9. Petroselli, A.; Grimaldi, S. Design hydrograph estimation in small and fully ungauged basin: A preliminary assessment of the EBA4SUB framework. J. Flood Risk Manag. 2018,11, 197–201. [CrossRef] 10. Petroselli, A.; Grimaldi, S.; Piscopia, R.; Tauro, F. Design hydrograph estimation in small and ungauged basins: A comparative assessment of event based (EBA4SUB) and continuous (cosmo4sub) modelling approaches. Acta Sci. Pol. Form. Circumiectus 2019 , 18, 113–124. [CrossRef] 11. Szymczak, T.; Kr˛e˙ załek, K. Prognostic model of total runoff and its components from a partially urbanized small lowland catchment. Acta Sci. Pol. Form. Circumiectus 2018,18, 185–203. [CrossRef] 12. Yuan, Y.; Mitchell, J.K.; Hirschi, M.C.; Cooke, R.A. Modified SCS curve number method for predicting subsurface drainage flow. Trans. ASAE 2001,44, 1673–1682. [CrossRef] 13. Wał˛ega, A.; Amatya, D.M.; Caldwell, P.; Marion, D.; Panda, S. Assessment of storm direct runoff and peak flow rates using improved SCS-CN models for selected forested watersheds in the Southeastern United States. J. Hydrol. Reg. Stud. 2020 ,27, 100645. [CrossRef] 14. Sahu, R.K.; Mishra, S.K.; Eldho, T.I. Performance evaluation of modified versions of SCS curve number method for two watersheds of Maharashtra, India. ISH J. Hydraul. Eng. 2012,18, 27–36. [CrossRef] 15. De Lima, J.L.M.P.; Singh, V.P. Laboratory experiments on the influence of storm movement on overland flow. Phys. Chem. Earth 2003,28, 277–282. [CrossRef] 16. Wang, A.; Jin, C.; Pei, J. A modified hortonian overland flow model based on laboratory experiments. Water Resour. Manag. 2006 , 20, 181–192. [CrossRef] 17. Chu, X.; Padmanabhan, G.; Bogart, D. Microrelief-controlled overland flow generation: Laboratory and field experiments. Hindawi Publ. Corp. Appl. Environ. Soil Sci. 2015,2015, 642952. [CrossRef] 18. Danino, D.; Svoray, T.; Thompson, S.; Cohen, A.; Crompton, O.; Volk, E.; Argaman, E.; Levi, A.; Cohen, Y.; Narkis, K.; et al. Quantifying shallow overland flow patterns under laboratory simulations using thermal and LiDAR imagery. Water Resour. Res. 2021,57, e2020WR028857. [CrossRef] 19. Richards, L.A. Capillary conduction of liquids in porous mediums. Physics 1931,1, 318–333. [CrossRef] 20. Green, W.H.; Ampt, G.A. Studies of soils physics I. The flow of air and water through soils. J. Agric. Sci. 1911,4, 1–24. 21. Mein, R.G.; Larson, C.L. Modeling infiltration during a steady rain. Water Resour. Res. 1973,9, 2, 384–394. [CrossRef] 22. Chow, V.T.; Maidment, D.R.; Mays, L.W. Applied Hydrology; McGraw-Hill Book Company: New York, NY, USA, 1988. 23. Chen, L.; Young, M.H. Green-Ampt infiltration model for sloping surfaces. Water Resour. Res. 2006,42, W07420. [CrossRef] 24. Chorma´nski, J.; Ignar, S.; Caba´nski, P. Zastosowanie modelu infiltracyjnego Green’a i Ampt’a oraz metod GIS do okre´slania opadu efektywnego w modelowaniu opad-odpływ na przykładzie zlewni górnej wilgi. Zesz. Probl. Post˛epów Nauk. Rol. 2008,532, 77–89. 25. Cho, S.E.; Lee, S.R. Evaluation of surficial stability for homogeneous slopes considering rainfall characteristics. J. Geotech. Geoenviron. Eng. 2002,128, 756–763. [CrossRef] Sustainability 2022,14, 4962 25 of 27 26. Cho, S.E. Infiltration analysis to evaluate the surficial stability of two-layered slopes considering rainfall characteristics. Eng. Geol. 2009,105, 32–43. [CrossRef] 27. Muntohar, A.S.; Liao, H.J. Analysis of rainfall-induced infinite slope failure during typhoon using a hydrological—Geotechnical model. Environ. Geol. 2009,56, 1145–1159. [CrossRef] 28. Bochenek, W.; Gil, E. Water circulation, soil erosion and chemical denudation in flysh catchment area. (In Polish: Procesy obiegu wody, erozji gleb i denudacji chemicznej w zlewni Bystrzanki). Przegl ˛ad Nauk. In˙ z. Kształt. Sr. 2007,16, 28–42. 29. Bochenek, W.; Gil, E. The diversity of overland flow and soil wash on experimental plots of different lengths (Szymbark, Low Beskidy Mts.) (In Polish: Zró˙ znicowanie spływu powierzchniowego i spłukiwania gleby na poletkach do´swiadczalnych o ró˙ znej długo´sci (Szymbark, Beskid Niski)). Pr. Studia Geogr. 2010,45, 265–278. 30. Martinez, G.; Weltz, M.; Pierson, F.B.; Spaeth, K.E.; Pachepsky, Y. Scale effects on runoff and soil erosion in rangelands: Observations and estimations with predictors of different availability. Catena 2017,151, 161–173. [CrossRef] 31. Mounirou, L.A.; Zoure, C.O.; Yonaba, R.; Paturel, J.E.; Mahe, G.; Niang, D.; Yacouba, H.; Karambiri, H. Multi-scale analysis of runoff from a statistical perspective in a small Sahelian catchment under semi-arid climate. Arab. J. Geosci. 2020 ,13, 154. [CrossRef] 32. Smolska, E. Runoff and soil erosion on sandy slope in the last-glacial area—Plots measurements (Suwałki Lake land, NE Poland. (in polish: Spływ wody i erozja gleby na piaszczystym stoku w obszarze młodoglacjalnym—Pomiary poletkowe (Pojezierze Suwalskie, Polska NE)). Pr. Studia Geogr. 2010,45, 197–214. 33. ´ Swi˛echowicz, J. Slopewash on agricultural foothill slopes in hydrological years 2007–2008 in Łazy (Wi´snicz Foothills). (In Polish: Spłukiwanie gleby na u˙ zytkowanych rolniczo stokach pogórskich w latach hydrologicznych 2007–2008 w Łazach (Pogórze Wi´snickie)). Pr. Studia Geogr. 2010,45, 243–263. 34. Mendes, T.A.; Gitirana, G.F.N.J.; Rebolledo, J.F.R.; Vaz, E.F.; da Luz, M.P. Numerical evaluation of laboratory apparatuses for the study of infiltration and runoff. Braz. J. Water Resour. 2020,25, e37. [CrossRef] 35. Rahardjo, H.; Lee, T.T.; Leong, E.C.; Rezaur, R.B. Response of a residual soil slope to rainfall. Can. Geotech. J. 2005 ,42, 340–351. [CrossRef] 36. Scherrer, S.; Naef, F.; Faeh, A.O.; Cordery, I. Formation of runoff at the hillslope scale during intense precipitation. Hydrol. Earth Syst. Sci. 2007,11, 907–922. [CrossRef] 37. Poesen, J. The influence of slope angle on infiltration rate and Hortonian overland flow. Geomorphology 1984,49, 117–131. 38. Nassif, S.H.; Wilson, E.M. The influence of slope and rain intensity on runoff and infiltration (L’influence de l’inclinaison de terrain et de l’intensitéde pluie sur l’écoulement et l’infiltration). Hydrol. Sci. J. 1975,20, 539–553. [CrossRef] 39. Kowalczak, P.; Kundzewicz, Z.W. Water-related conflicts in urban areas in Poland. Hydrol. Sci. J. 2011,56, 588–596. [CrossRef] 40. Skotnicki, M.; Sowi´nski, M. The influence of depression storage on runoff from impervious surface of urban catchment. Urban Water J. 2015,12, 207–218. [CrossRef] 41. Jarosi´nska, E. Local flooding in the USA, Europe, and Poland—An overview of strategies and actions in face of climate change and urbanisation. Infrastruct. Ecol. Rural. Areas 2016,3, 801–821. [CrossRef] 42. Walczykiewicz, T.; Skonieczna, M. Rainfall flooding in Urban areas in the context of geomorphological aspects. Geoseinces 2020 ,10, 457. [CrossRef] 43. Yonaba, R.; Biaou, A.C.; Koïta, M.; Tazen, F.; Mounirou, L.A.; Zouré, C.O.; Queloz, P.; Karambiri, H.; Yacouba, H. A dynamic land use/land cover input helps in picturing the Sahelian paradox: Assessing variability and attribution of changes in surface runoff in a Sahelian watershed. Sci. Total Environ. 2021, 757. [CrossRef] 44. Starkel, L. Geomorphic hazards in the Polish Flysch Carpathians. Studia Geomorphol. Carpatho-Balc. 2006,11, 7–19. 45. Bodziony, M.; Baziak, B. Heavy rain effects on the example of 1997 and 2005 floods in the Wielka Puszcza basin. (In Polish: Skutki deszczy nawalnych na przykładzie powodzi w zlewni rzeki Wielkiej Puszczy w latach 1997 i 2005). Czas. Tech. Sr. 2007 ,2, 13–28. 46. Wałega, A.; Cupak, A.; Amatya, D.M.; Dro˙ zd˙ zal, E. Comparison of direct outflow calculated by modified SCS-CN methods for mountainous and high land catchments in upper vistula basin, Poland and lowland catchment in South Carlina, USA. Acta Sci. Pol. Form. Circumiectus 2017,16, 187–207. [CrossRef] 47. Szwagrzyk, M.; Kaim, D.; Price, B.; Wypych, A.; Grabska, E.; Kozak, J. Impact of forecasted land use changes on flood risk in the Polish Carpathians. Nat. Hazards 2018,94, 227–240. [CrossRef] 48. Gil, E. Water Circulation and Wash Down on the Flysch Slopes Used fo r Fanningpurposes in 1980–1990 Years (Results of Investigation on Experimental Plots at ResearchStation of Institute of Geography and Spatial Organization Polish Academy of Sciences in Szymbark) (In Polish: Obieg Wody i Spłukiwanie na Fliszowych Stokach U˙ zytkowanych Rolniczo w Latach 1980–1990 (Wyniki Bada´n Przeprowadzonych na Poletkach Do´swiadczalnych na Stacji Naukowej IGiPZ PAN w Szymbarku)); Zeszyty IGiPZ PAN: Warszawa, Poland, 1999. 49. Moeyersons, J.; Makanzu, F.M.; Imwangana; Dewitte, O. Siteand rainfall-specific runoff coefficients and critical rainfall for mega-gully development in Kinshasa (DR Congo). Nat. Hazards 2015,79, S203–S233. [CrossRef] 50. Mahmoud, S.H.; Mohammad, F.S.; Alazba, A.A. Determination of potential runoff coefficient for Al-Baha Region, Saudi Arabia using GIS. Arab. J. Geosci. 2014,7, 2041–2057. [CrossRef] 51. Rutkowski, J. Budowa geologiczna regionu Krakowa. Prz. Geol. 1989,37, 302–308. 52. Wójcik, A.; Kamieniarz, S.; Wódka, M.; Biajgo, A.; Janeczek, A.; Walatek, M. Atlas Osuwisk Miasta Krakowa; Urz ˛ad Miasta: Kraków, Poland, 2019. 53. Rutkowski, J. Szczegółowa Mapa Geologiczna Polski, Arkusz 973—Kraków; Pa´nstwowy Instytut Geologiczny: Warszawa, Poland, 1989.