Splash erosion affected by initial soil moisture and surface conditions under simulated rainfall

Soil erosion by water is one of the most severe soil degradation processes. Splash erosion is the initial stage of soil erosion by water, resulting from the destructive force of rain drops acting on soil surface aggregates. Apart from rainfall properties, constant soil physical properties (texture and soil organic matter) are crucial in understanding the splash erosion. However, there is lack of information about the effect of variable soil properties such as soil initial water content and surface condition (seal formation) on splash erosion. The objective of the present study was to determine how initial water content and surface condition affected soil splash erosion under simulated rainfall. The changes in soil surface condition were characterized by hydraulic variability (saturated hydraulic conductivity) due to surface seal formation. Slit loam and loamy sand soil textures were used in the experiment. The soil samples were collected from the top layer; air dried, and filled into modified Morgan splash cups for splash erosion measurements. Rainfall was created in the laboratory using two types of rainfall simulators covering intensity range from 28 to 54 mm h and from 35 to 81 mm h. The soil samples were exposed to three consecutive rainfall simulations with different time intervals between simulations and different initial water content and surface conditions (air-dried, wet-sealed, and dry-crusted). Wet-sealed soil samples had up to 70% lower splash erosion rate compared to air-dried samples, due to surface ponding followed by seal formation. A significant decrease in soil saturated hydraulic conductivity indicated the formation of surface seal for silt loam soils. A non-significant decrease in saturated hydraulic conductivity for loamy sand soil was attributed to earlier formation of stable seals. Two different rainfall simulators produced different amount of splash erosion rates; however, the splash erosion development for increasing rainfall intensity was almost equal considering same initial surface condition. These results provide insight into dynamic changes of individual soil parameters affected by rainfall, and could find wider application for more complex soil erosion prediction models.


Introduction
Detachment of soil by rain drop impact is the first stage of the soil erosion process by water (Quansah, 1981). According to Rose (1960) and Hairsine and Rose (1991), splash has more influence on detached soil particles than surface runoff, before the stage of rill and gully erosion is reached. Bare soil surface exposed to rain drops changes its structural and hydraulic properties, which remarkably influences soil infiltration, soil water repellency, overland flow, and final soil erosion rates (Fernández-Raga et al., 2019). The main driver for the splash detachment process is the kinetic energy (KE) of rainfall, which depends on the amount, size, and fall velocities of the drops according to Wischmeier et al. (1971), Ghadiri and Payne (1977) and Morgan (2005). Together with rainfall characteristics, soil physical parameters are crucial in defining the soil erosion process. Wischmeier and Smith (1978) concluded that particle size distribution and organic matter content were the most dominant indicators of soil erodibility. Le Bissonnais (2016) also reported that soil mineralogy, soil texture, organic matter content and initial water content (θ a ) influence the formation of aggregates, where higher θ a increases the resistance of aggregates against the rain drop impact.
At the beginning of rainfall, KE of rain drops has to exceed a threshold, in order to induce destruction of soil aggregates (Kinnell, 2005). If detached soil particles are deposited on the soil surface, they can cause pore filling and clogging by wash-in of fine soil particles (Assouline, 2004;McIntyre, 1958). Soil infiltration rate is then reduced with incipient surface ponding. During ponding the soil surface is rapidly sealed and the crust formation is dependent on cumulative rainfall KE (Baumhardt et al., 1990). Development of seal on soil surface depends on the soil characteristics, being mostly common for soils with high silt content and low organic matter (Armenise et al., 2018). When observing the interaction of splash erosion and seal development, Cheng et al. (2008) discovered that the splash erosion fluctuations are related to surface crust development, where the final splash rates were lower for the soils suspected to surface sealing. Furthermore, several researches have investigated the effect of ponding on soil erosion. For instance, Gao et al. (2003) proved that shallow water layer can accelerate the rain drop impact force on soil erosion, until it reaches a critical depth when the soil detachability decreases. However, they did not consider surface seal formation during ponding phase. Proposed by Guy et al. (1987), the defined critical ponding depth on soil surface is rainfall dependent and equal to three drop diameter. During the rainfall, soil hydraulic and structural properties define the soil erosion and runoff production. Different approaches have been developed to describe the water movement through the soil affected by surface sealing. Assouline and Mualem (1997) related the seal formation to increase in soil bulk density. Furthermore, in later study (Assouline and Mualem, 2002), they proved that the use of saturated hydraulic conductivity (K s ) for an infiltration model can sufficiently describe the heterogeneity in soil hydraulic properties.
Many experimental studies (e.g. Kinnell, 1982;Salles et al., 2001;Salles and Poesen, 2000;Sharma et al., 1991) have been dealing with rainfall properties controlling splash erosion. However, there is still scarce information in most of these studies about the influence of soil physical properties (moisture, texture, structure, infiltration capacity etc.) on splash erosion. Few authors (Beczek et al., 2019;Truman and Bradford, 1990;Vermang et al., 2009) investigated the effect of different soil moisture on splash erosion, however, different results reported could be contributed to particular conditions in the experiment. Variable results were mainly concerning the different preventing mode, soil organic carbon and clay content of soil samples (Vermang et al., 2009), which can affect the result interpretation. Furthermore, there is lack of studies relating splash erosion with sealing and crust formation and its influence on infiltration (Fernández-Raga et al., 2017).
The main aim of this study is to investigate the effect of different soil moisture content and surface condition (seal formation) on splash erosion for three soils under simulated rainfall. The changes on soil surface were related to changes in soil K s . Consequently, second objective is to quantify the differences in soil K s affected by different rainfall intensities and corresponding KE. According to the recent knowledge of authors, the influence of different rainfall characteristics produced by two rainfall simulators on splash erosion development has never been investigated within one study. Therefore, a third objective is to quantify the differences between rainfall characteristics produced by two rainfall simulators and their impact on soil splash erosion for different soil initial conditions.
Samples for chemical analyses were air-dried, crushed, mixed and passed through a 2 mm sieve. Total carbon content was analysed by dry combustion (ÖNORM L 1080(ÖNORM L , 2013(ÖNORM L , 2013) using a C/N Analyzer (Vario Max CN, Elementar). Soil organic carbon content (OC) was obtained by subtracting inorganic carbon content measured volumetrically by the Scheibler method with a Calcimeter (ÖNORM L 1084, 2016) from total carbon content measured by C/N Analyzer. Further, the calcium carbonate was measured with the calcimeter using the method by Scheibler (ÖNORM L 1084. Soil pH was determined in a 1:2.5 soil to water solution, using a glass electrode by Metrohm. Cation exchange capacity (CEC) was estimated by extraction of the effective exchangeable cations by barium chloride solution following ÖNORM L 1086-1, 2014. Aggregate stability (AS) of soils was determined with the wet sieving method according to Kemper and Koch (1966) in modified form. The physical and chemical properties of soils obtained from laboratory analysis are listed in Table 1.

Experimental design
The experiments took place at the Institute for Soil Physics and Rural Water Management at University of Natural Resources and Life Sciences in Vienna, Austria (BOKU) and the Institute of Land and Water Management Research in Petzenkirchen, Austria (BAW). Modified Morgan splash cups (Morgan, 1981) were used for the splash erosion measurements. The splash cups were constructed from PVC drainage pipe tops with the inner diameter of 10.5 cm ( Fig. 1-a) and height of 5 cm ( Fig. 1-b). In the bottom of the splash cup, holes were drilled, to ensure water drainage through the soil layer during the rainfall. First, two nylon meshes (500 and 1000 µm) were placed on the bottom of the splash cup, to prevent soil loss trough the holes ( Fig. 1-a). Air-dried and sieved (aggregates < 10 mm) soil material was filled into the splash cups. The soil was loosely packed in three layers to reach similar bulk density as in seedbed condition (1-1.2 g cm −3 ) for each corresponding soil. Soil layers were levelled with a long metal needle. Top layer was filled up to 1 cm below the splash cup edge to prevent surface overflow during high intensities ( Fig. 1-c). At the top layer soil aggregates were randomly distributed to achieve heterogeneous arrangement of the all fractions (from fine to coarse).
Splash cups were placed in the middle of a splash collector, with standing height of 30 cm and diameter of 45.5 cm ( Fig. 1-d). Splash collector had an outlet, ensuring the drainage of rainfall water with detached soil into buckets placed underneath. After the exposure to rainfall, eroded soil was rinsed from the collector rim, filtered and oven dried at 105°C. The splash weighted from the filters corresponds to splashed soil from the splash cup surface of 86.6 cm 2 . Fig. 2 illustrates schematic overview of experiments made with two rainfall simulators. The three soils (MI, ZW, BK) were subjected to three rainfall simulations of 30 min, with various rainfall intensities across the positions under the rainfall simulator and different initial soil surface conditions. The first rainfall simulation was performed on air-dried soil surface (AD, θ < 5%). Second rainfall simulation followed after 24 h on wet-sealed (WS, θ a~3 0%) and final (third) simulation on dry-1 Primary classifications of particle size distribution following ÖNORM L 1050, 2016 (1994): sand (particle sizes from 0.2 mm − 0.063 mm diameter), silt (particle sizes from 0.063 mm − 0.002 mm diameter), clay (particle sizes smaller than 0.002 mm). N. Zambon, et al. Catena 196 (2021) 104827 crusted soil surface (DC, 5% < θ a < 10%) after approximately ten days of drying. A total of 9, 6 and 3 replications of splash erosion measurements were made for the AD, WS and DC surface condition, respectively. At the end of each rainfall simulation, K s of soil samples was measured using constant head method, according to description by Klute and Dirksen (1987) at the BAW. Each measurement of soil K s was replicated three times. Furthermore, the changes in soil infiltration rate were observed by measuring the time from the beginning of the rainfall until the appearance of surface ponding (accumulated water layer). This was defined as ponding time (t p ). Ponding was recorded by two cameras installed on the corners of experimental area under the BAW rainfall simulator. The cameras recorded a photo in one-minute interval during the rainfall simulation. Photos from the camera were analysed and the t p was registered. Additionally, ponding was visually observed during the rainfall simulation where t p for each soil sample was noted.
The final t p form camera and visual observation varied in ± 1 min, therefore t p form the camera observations were taken as a reference.
Ponding was not temporary observed under the rainfall simulator at BOKU due to technical reasons; it was only noted if the samples had surface ponding or not after each experiment.

Rainfall simulators
For generating artificial rainfall two types of rainfall simulators were used. Norton Ladder type of rainfall simulator located at the BOKU consisted of four oscillating VeeJet 80100 spray nozzles, arranged in two rows and elevated 2.3 m above the splash cups. The simulated rainfall was operated with a pressure of 0.45 bar at the nozzles and the water was distributed from deionized water supply. Totally, nine different positions were arranged under the rainfall simulator for splash erosion measurements, as presented on the Fig. 3-a.
The rainfall simulator at BAW was equipped with one FullJet nozzle (½ HH-30WSQ), where intensity was controlled electronically by discontinuous spraying (Strauss et al., 2000). This design was deviating from normal use (three nozzles) to produce as much heterogeneity in rainfall intensity as possible. The nozzle was elevated 2.3 m above the splash cups and six positions were selected for splash erosion measurements ( Fig. 3-b). The BAW rainfall simulator used deionized water with a constant water pressure of 0.25 bar at nozzles.
Rainfall simulations were performed for 30 min with a constant rainfall intensity rate. The maximum rainfall intensity tended to be between 40 and 80 mm h −1 , since it was attempted to keep the average intensity equal to those measured (data not shown) under natural rainfall at the locations from which the soil samples were collected. Real distribution of rainfall intensities below the simulators was measured for each position defined for splash erosion measurements, as shown on Fig. 3. Intensities were calculated from the accumulated water volume captured by the splash collectors and drained into the buckets beneath during 30 min rainfall simulation. Total volume of accumulated water in the buckets below the splash collectors was divided by the splash collector area (1625 cm 2 ). The intensities under the rainfall simulator at BOKU were measured after each experiment with soil samples. The mean of totally three replicates per soil type and initial condition for defined position was used in the results. Under the Table 1 Physical and chemical properties of soil material (0-10 cm soil depth) including: particle size distribution, aggregate stability (AS), alkalescence (pH), calcium carbonate content (CaCO 3 ), organic carbon (OC) and cation exchange capacity (CEC).  N. Zambon, et al. Catena 196 (2021) 104827 rainfall simulator at BAW, the intensity for defined positions of splash cups was calculated as the mean value obtained from totally six rainfall simulations (independent from soil or initial condition). Rainfall KE, mean drop diameter, median drop diameter and mean drop fall velocity were derived from drop size distribution (DSD) measured with an optical laser disdrometer Weather Sensor OTT Parsivel Version 1 (Parsivel) by OTT Messtechnik. Parsivel disdrometer has a measuring area of 54 cm 2 and categorises the drops into 32 drop size and velocity classes (OTT, 2005). The KE under both rainfall simulators was obtained for each position defined for splash erosion measurements, as shown in Fig. 3. Disdrometer was centred on the positions of corresponding splash cups, which were previously marked on the ground. It was ensured that the height of the disdrometer laser beam is equal to the height of soil surface in the splash cup. For each position 15 min of rainfall was measured with the disdrometer.
The kinetic energy KE i,j (J m −2 ) of rainfall per minute was computed for the diameter class i and velocity class j, that are provided by disdrometer, as follows: where m i is the mean mass [g] corresponding to the drop diameter class i; N is number of detected raindrops of a certain size class i and velocity class j; A is the is the sampling area of the disdrometer (m 2 ); ρ w is density of water (g cm −3 ); D i is mean drop diameter (mm) of size class i; and v j is mean fall velocity (m s −1 ) of velocity class j. The mass of raindrop was calculated assuming a spherical drop shape. Total KE is the sum of KEs for each drop size and velocity, multiplied with number of drops in the corresponding classes.

Statistical analysis
Statistical analyses were performed in R Studio (R Development Core Team, 2015). The Kruskal and Wallis (1952) test (one-way ANOVA on ranks) was used to identify statistical differences in K s for different intensity rates (positions). Each position under the rainfall simulator represented one group of data for which the K s was obtained. The positions (groups) differ in intensity rates which they are exposed to. The comparison is based on three replicates (n = 3) of K s values obtained for each group (intensities). Multiple comparison between groups (post hoc) was further conducted with Dunn's test (Dunn, 1964) with p-adjustments of Benjamini-Hochberg (Benjamini and Hochberg, 1995). The differences in splash erosion rates between the three soils were analysed using the Student's t-test. N. Zambon, et al. Catena 196 (2021) 104827 3. Results

Intensity and kinetic energy
According to intensity measurements, the rainfall simulators at BOKU and BAW covered intensity range from 28.2 to 54.2 mm h −1 ( Table 2) and 35.3 to 81.2 mm h −1 (Table 3), respectively. Drop size distribution indicated that the rainfall simulator at BAW produced higher percentage of drops with drop diameter larger than 0.6 mm ( Fig. 4-a), however, rainfall simulator at the BOKU had higher percentage of drops with velocity greater than 5.2 mm s −1 (Fig. 4-b). This resulted in larger mean and median drop size obtained for positions under the BAW rainfall simulator (Table 2) and higher mean raindrop velocity for positions under the BOKU rainfall simulator (Table 3). Distribution of mean drop velocity per drop size class (Fig. 5), indicated that 14% of total drops with diameter larger than 1.3 mm under the BAW rainfall simulator did not reach their terminal velocity, as specified by Atlas et al. (1973). Since the rainfall simulator at the BOKU produced higher velocity of larger drops, the resulting KE per mm of rainfall was on average 62% higher (Table 2) compared to KE per mm of rainfall for rainfall simulator at the BAW (Table 3). Total KE ranged between 504.4 and 923.1 J m −2 h −1 for the BOKU rainfall simulator (Table 2), and between 375.8 and 961.8 J m −2 h −1 for the BAW rainfall simulator (Table 3).

Splash erosion
Due to higher KE per mm of rainfall measured under the rainfall simulator at the BOKU, overall splash erosion rates were almost twice as high compared to results obtained with BAW simulator as shown on Fig. 6 (excluding the ZW soil for AD condition on Fig. 6-d).
For the rainfall simulator at BOKU, the MI soil had widest range of splash rates (0.06-0.58 kg m −2 h −1 ) compared to ZW (0.02-0.31 kg m −2 h −1 ) and BK (0.12-0.43 kg m −2 h −1 ) soil. Splash erosion under the BAW rainfall simulator ranged similarly between three soils (0.02-0.23 kg m −2 h −1 ) for AD initial condition ( Fig. 6a,d,g). The highest splash erosion was measured for samples exposed to KE of 667 J m −2 h −1 with largest drop diameter (Table 3). Relationship between splash erosion for three soils and KE was linear for the BOKU rainfall simulator, and non-linear for the BAW rainfall simulator (Table 4). Highest coefficient of determination (R 2 ) was calculated for ZW soil, considering both rainfall simulators. Between 60 and 72% less splash erosion, compared to simulation with AD initial condition, was measured for three soils with WS initial condition under both rainfall simulators. This reduction was evident for samples exposed to KE > 660 J m −2 h −1 , which were affected by surface ponding. However,  N. Zambon, et al. Catena 196 (2021) 104827 splash erosion increased for samples exposed to lower KE. This was particularly noted for MI and ZW soil under the BOKU rainfall simulator ( Fig. 6-b,e). Reduction in splash erosion for WS initial condition resulted in negative linear function between splash and KE with low R 2 (mostly < 0.10) considering both rainfall simulators ( Table 4). The splash erosion for the samples with DC surface condition increased with KE, except for BK soil under the BOKU rainfall simulator (Fig. 6-i). Similarly to splash erosion results with WS condition, MI and ZW soil showed increase of roughly 50% in splash erosion rates for the samples exposed KE < 660 J m −2 h −1 (Fig. 6-c,f). Highest R 2 between splash erosion for three soils and KE was obtained for samples under the BAW rainfall simulator (Table 4).

Saturated hydraulic conductivity
Splash erosion affected by surface ponding during the rainfall simulations was contributed to different θ a and changes in soil hydraulic properties. To quantify this changes soil K s was measured after each rainfall simulation with different soil initial condition. Fig. 7 represents average K s (n = 3) for the three soils and three initial conditions. Highest K s between three soils were measured for ZW soil with the maximum of 1,102.9 mm h −1 obtained for initial AD surface condition ( Fig. 7-d). Lowest values were measured for DC surface condition ( Fig. 7-f). Similar trend was observed for MI soil, where highest K s of 461.7 mm h −1 decreased to 56.3 mm h −1 for DC surface condition ( Fig. 7-a,c). Generally, BK soil exhibited lowest K s between three soils with maximum K s of 193 mm h −1 for WS initial surface condition ( Fig. 7-g-i). When comparing the K s values measured for individual surface conditions it was observed that the samples exposed to lower (< 56.7 mm h −1 ) intensities had higher K s than samples exposed to high intensities. These differences were significant (P < 0.05) for MI soil with AD samples (Fig. 7-a), including ZW soil with WS and DC samples (Fig. 7-e,f). Furthermore, high variations between the replicates, typical for the samples exposed to lower intensities (< 56.7 mm h −1 ), could result in no significant (P > 0.05) difference between K s obtained for AD and WS surface condition considering three soils ( Fig. 7-b,d,h).

Relationship between soil saturated hydraulic conductivity and rainfall kinetic energy
As shown in previous results, the K s of soil samples tend to decrease with the increasing intensity and after subsequent exposure to rainfall. To investigate the impact of KE on K s , total KE applied on soil surface during the three rainfall simulations was plotted against K s measured for each soil initial condition (Fig. 8). The relationship obtained for MI, ZW and BK soil can be described as a power function of accumulated KE, which agrees best for MI and ZW soil with R 2 of 0.54 and 0.74, respectively, where for BK soil low R 2 of 0.22 was obtained.
From the KE-K s functions obtained, the decrease in K s with increasing KE was observed until reaching a steady state, where KE had no impact on further K s reduction. Assuming that the decrease in K s was accompanied with a surface sealing, constant K s values would indicate a final stage in the surface seal formation. The amount of KE needed for fully developed surface sealing could be calculated by considering the impact of KE on soil K s . After a certain amount of KE applied on the soil    Atlas et al. (1973). The velocities are shown for drop size classes within the total drops amount was larger than 100.
surface, K s significantly (P < 0.05) varied for increasing KE until reaching a constant value (Fig. 8). This amount of KE was addressed as the threshold value required for fully developed seal formation and it represented the transformation point between variable and constant soil K s . Accordingly, the calculated threshold value of accumulated kinetic energy (KE 0 ) was found to be 0.39 kJ m −2 for MI soil and 0.59 kJ m −2 for ZW soil ( Fig. 8-a,b). The BK soil showed no significant impact of increasing KE on K s ; therefore, the threshold value could not be obtained (Fig. 8-c).

Surface ponding
Present surface ponding was mostly concentrated within samples in central positons under the rainfall simulator at BOKU (e.g. positon 2, 4, 5, 6, and 8), where the drops overlapped (Table 2). For these positons higher KE and lower final splash erosion rates were measured, compared to positons with no surface ponding. During the rainfall simulation at the BAW, t p was recorded for each soil sample in splash cup (Table 5). Generally, highest t p was measured for the AD and DC surface condition, and lowest for WS surface condition. This was expected considering that samples with AD and DC initial condition had higher infiltration capacity due to lower θ a and samples with WS initial condition had lower infiltration capacity due to higher θ a . Furthermore, this results correlated to the K s values obtained for the three soils. The shortest t p was measured for BK soil, which had the lowest K s and longest t p for the ZW soil with highest K s compared to other soils.
Measured t p for the samples with WS initial condition was related to Whiskers indicate ± standard deviation of mean splash erosion within replicates (n = 3, 6, 9).

Table 4
Linear and nonlinear regression with associated determination coefficient (R 2 ) between kinetic energy (KE) (kJ m −2 h −1 ) and splash erosion (S) (kg m −2 h −1 ) obtained for two rainfall simulators, at BOKU and BAW institute, and Mistelbach (MI), Zwerbach (ZW) and Býkovice (BK) soil with air-dried (AD), wet-sealed (WS) and dry-crusted (DC) soil initial condition (IC).   N. Zambon, et al. Catena 196 (2021) 104827 the intensity rates (Fig. 9). The results showed high agreement between t p and rainfall intensity with R 2 of 0.90 and 0.91 obtained for MI and BK soil, respectively. The R 2 for ZW soil was remarkably lower, due to high difference in t p between the positons 3 and 6, with same intensity rate of 56 mm h −1 (P3, P6 on Fig. 9). Larger drop diameter measured for position 3 (Table 3) could contribute to greater surface compaction during the previous simulation with AD initial condition. This eventually resulted in lower K s for positon 3 and therefore, shorter t p during the simulation with WS samples (Fig. 7-e). Similar results of t p for position 3 and 6 were also observed for MI soil during the simulation with AD initial condition (Table 5).

Effect of rainfall kinetic energy produced by two rainfall simulators on splash erosion
The results of splash erosion obtained from both rainfall simulators revealed comparable trend of splash development in relation with KE (Fig. 6). However, twice as high splash erosion rates were measured for the rainfall simulator at BOKU. The velocities of raindrops were crucial in defining rainfall kinetic energy, where higher velocities contributed to higher KE for BOKU rainfall simulator (Fig. 4-b). The differences in drop size-velocity distribution between different rainfall simulators can result in the significantly different splash erosion rates. This should be considered when comparing the splash erosion results from different studies.
The rainfall KE produced by rainfall simulator is often defined by (uniform) drops of different size and fall height, while the KE obtained in the field conditions depends primarily on drop size distribution (Van Dijk et al., 2003). Similar was observed with the two rainfall simulators used in present study. The drop distribution of rainfall simulator at the BAW was too uniform to describe the natural rainfall and the velocities of larger drops (> 1.3 mm) were far from the terminal velocity, defined by Atlas et al. (1973). On the contrary, the BOKU rainfall simulator produced higher velocities for most of drops in the same diameter class. However, large drops were still under the terminal velocity line (Fig. 5). Related to drop terminal velocity under the simulated rainfall, Iserloh et al. (2013) concluded that larger drops produced by different rainfall simulators are usually not able to reach their terminal velocity (mostly due to low fall height). For this reason, the rainfall KE produced in the laboratory cannot completely represent KE of natural rainfall.

Effect of soil water content and surface condition on splash erosion
Fluctuations in splash rates obtained from individual rainfall simulator were combined effect of θ a and seal development affected by KE. Lower θ a (AD initial condition) contributed to greater splash erosion for samples exposed to high intensities and KE. Therefore, splash erosion could be described as a linear or power function of increasing KE for AD initial condition, which was also obtained in previous studies (Fernández-Raga et al., 2010;Sharma et al., 1991;Zumr et al., 2019). However, appearance of surface ponding in the last stage of rainfall simulation indicated the reduction in soil infiltration rate, characterized by lower soil permeability due to surface seal formation (Assouline, 2004;Liu et al., 2011;Sharma et al., 1995). Surface seal was easily observable on the soil samples exposed to high intensities (> 50 mm h −1 ) with completely smooth soil surface after the rainfall simulation. In addition, significantly lower K s (Fig. 7) for the samples exposed to high intensities indicated the beginning of surface seal formation.
In the following rainfall simulation with WS soil samples splash erosion decreased with increased rainfall KE (Fig. 6-b,e,h). The final splash erosion rates for samples exposed to highest KE were evidently lower compared to rates for AD surface condition. According to results in Table 5, surface ponding was initiated earlier than for AD and DC initial condition. This was contributed to high θ a , where early surface ponding at higher θ a results from rapid decline in the hydraulic gradient for intensive rainfall in combination with a lower storage capacity (Liu et al., 2011;Vermang et al., 2009). Furthermore, partly formed seals (for MI and ZW soil) from previous simulation with AD soil samples probably resulted in higher resistance of soil surface to splash erosion prior to surface ponding. Following surface ponding eventually lowered the raindrop impact on soil surface preventing the further splash erosion (Poesen, 1981). Similar results were reported by Vermang et al. (2009), where soil detachment decreased along with surface seal development and formation of shallow water layer by surface runoff.
During the simulation with DC initial condition, the influence of the  θ a on reduction in t p could be excluded, since the initial θ a was approximately three times lower. Therefore, the reduction in infiltration rates followed by surface ponding (Table 5) may be due to already developed seals, after previous exposure of the samples to rainfall (with AD and WS initial condition). Unlikely to splash erosion for WS initial condition, the reduction in splash erosion rates was noticed only for BK soil under the BOKU rainfall simulator (Fig. 6-i). Results for other soils, however, suggested the increase in splash erosion compared to samples with AD initial condition under lower rainfall intensities. Opposite results were reported by Le Bissonnais et al. (1995), where drying of aggregates increased the stability against aggregate breakdown. According to Lado et al. (2004), the wetting rate in combination with initial water content determined the magnitude of slaking forces causing aggregate break down. These forces are greater if the aggregates are drier. Therefore, it could be hypothesized that the higher splash erosion rates during the rainfall simulation with DC initial condition for the MI and ZW samples could be contributed to aggregate destruction by slaking. Generally splash erosion increased or did not vary between the simulations on samples with different initial conditions which were not affected by surface ponding. This was related mostly to samples exposed to lower intensity and KE (< 660 J m −2 h −1 ). High θ a for WS samples may result in higher splash erosion (compared to AD initial condition) for the slow ponding soil samples considering both ZW and MI soil ( Fig. 6-a-f). The findings of Beczek et al. (2019) and Gao et al. (2003) confirmed that the mass of monitored eroded material increased with higher θ a due to lower cohesion between soil particles, considering that no surface seal was present.
Following the above lines discussed, it is difficult to select one factor that describes the differences in splash erosion rates for different initial conditions and soils. Soil conditions before rainfall and their variability during the rainfall will contribute to variation in splash characteristics. Depending on rainfall, soil structure, physico-chemical properties and θ a , response of the soil to raindrop impact will vary. Nevertheless, for the three soils studied under rainfall conditions in this experiment can be concluded that reduction in splash erosion was primarily attributed to surface ponding initiated by high rainfall intensities and subsequent sealing formation. High initial θ a reduced soil infiltration capacity and induced faster surface ponding (for WS initial condition) with decrease in splash erosion rates (Walker et al., 2007). In the case of the three soils studied here, initially lower K s for MI and BK soil could be attributed to lower porosity considering lower clay content (12-18%) than for ZW soil, as indicated by Wei et al. (2015). Furthermore, it should be noted that beside the soil texture other soil properties, such as OC and calcium cations (Ca 2+ ) could affect the soil hydraulic properties (Wuddivira and Camps-Roach, 2007) and eventually splash erosion rates. However, we could not establish correlation between K s and OC or CaCO 3 (contains Ca 2+ ) for the three soils presented here (Table 1, Fig. 7).

Effect of soil properties on surface seal formation
Soil detachment is depended on aggregate strength and after seal development on seal strength (Vermang et al., 2009). The both are highly related to soil physico-chemical properties (Mualem et al., 1990). Stability of aggregates obtained for three soils could not explain differences in the splash rates or the seal development in the present study. For example, BK soil had highest ratio of stable aggregates (Table 1), however, the splash erosion rates for AD surface condition did not significantly (P < 0.05) differ compared to MI soil with lowest AS. High OC and CaCO 3 content could be favourable for higher soil AS (Le Bissonnais et al., 1993). Similar OC content between the soils (Table 1) and low CaCO 3 soil could not explain higher AS for BK. However, it might be considered that higher clay content for MI and ZW soil could promote the slaking forces during the wet-sieving and decrease of the AS. Furthermore, the K s obtained for BK was lowest compared to other two soils and did not vary between the simulations (Fig. 8-c, Fig. 7-g-i). On the one hand, this could indicate that the KE did not have major effect on the surface seal formation. On the other hand, stable surface seal formation could be initiated in early stage of rainfall simulation with AD soil samples. Therefore, further development (decrease) was not detected. In addition to that, the BK soil had visibly smaller aggregates compared to other two soils. For this reason, the lower KE was needed to destroy aggregates and to form crust. In the study by Fox et al. (2007), smaller fractions were more susceptible to surface crusting and splash erosion than the coarser fractions. Furthermore, lower surface roughness due to smaller aggregates might indicate smaller depressional storage, which induced earlier surface runoff (Truman and Bradford, 1990).
ZW soil showed lowest splash erosion rates and high K s for AD surface condition compared to three soils. Le Bissonnais et al. (1995) reported that soils with high OC content have lowest erosion rate in AD conditions. Similar values of OC content were obtained for three soils (Table 1). Considering this, low splash erosion rates for ZW soil in AD condition were not contributed to OC content. However, according to Vermang et al. (2009) andXiao et al. (2018) during the fast wetting period on dry surface, formation stable aggregates is affected by cementing effect of clay particles trough the surface water tension. This may be the reason for higher resistance of the ZW soil to soil erosion and seal development.
Manny studies provided the evidence that formation of sealing was characteristics of soils with high silt content (Cheng et al., 2008;Rodrigo Comino et al., 2017;Truman and Bradford, 1990). This could be also applied on results obtained in this study for MI soil, characterized by high silt content and low AS ( Table 1). Assuming that the decrease in splash erosion and K s was affected by surface seal formation, the MI soil showed highest reduction in K s between the rainfall simulations among the three soils.

Effect of rainfall kinetic energy on surface seal formation
The obtained threshold values of KE required to form stable seal formation, confirmed the assumption of greater seal development under the soil surfaces exposed to higher KE. This was also stated by Bedaiwy (2008), where the important influence of KE on surface formation was highlighted. Constant or increasing splash erosion rates obtained for all three soils exposed to low intensity rates (< 35.3 mm h −1 ) implied that accumulated KE throughout the simulations was lower than critical KE to initiate surface seal formation. However, high deviations between the replicates for K s require more measurements in order to more precisely describe the process of surface crust development for certain scenarios.

Conclusion
In this study different scenarios of splash erosion development were obtained by applying simulated rainfall produced from two different rainfall simulators on soil samples with air-dried, wet-sealed and drycrusted surface condition. Both rainfall simulators exhibited same trend of splash erosion development by applying similar range of intensities, though the total amount of splash rates where significantly different. This was contributed to differences in drop and velocity spectrum resulting in different kinetic energy produced for the same rainfall amount. Since splash erosion in primarily affected by raindrop impact, a special attention should be given when comparing the results obtained with different rainfall simulators due to variabilities in drop and velocity distribution.
The splash erosion rates increased with increasing kinetic energy for air-dried and dry-crusted soil samples with lower initial water content. Higher initial water content contributed to decrease (up to 70%) in splash erosion rates for the wet-sealed initial condition due to early surface ponding and sealing. Time to ponding measurements verified that decrease in soil infiltration rate for wet-sealed condition is the function of increasing intensity with R 2 of 0.90 and 0.91 for silt loam and loamy sand soil, respectively. The formed seal layer can be reflected on decrease in saturated hydraulic conductivity, where the predominant factor for its reduction was kinetic energy. We identified a threshold rainfall kinetic energy from which stable surface seals are formed, which equals to 0.39 and 0.59 kJ m −2 h −1 for two silt loam soils used in this experiment. Considering high variability among replicates for saturated hydraulic conductivity measurements, it is difficult to identify clear relationship for some scenarios. Further experiments comparing more soil types should be conducted in order to precisely define soil specific properties controlling the splash erosion. In natural conditions soil is exposed to frequent changes in soil moisture and surface structure due to variable weather impacts. This research presents that the conditions before rainfall such as initial water content and surface condition highly affect soil erodibility, infiltration and final splash erosion. Therefore, including these parameters into soil erosion prediction models could significantly improve their accuracy.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.