Improved calibration of Green-Ampt infiltration in the EROSION-2D/3D model using a rainfall-runoff experiment database

EROSION-2D/3D model using a rainfall-runoff experiment database Hana Beitlerová1, Jonas Lenz2, Jan Devátý3, Martin Mistr1, Jiří Kapička1, Arno Buchholz2, Ilona Gerndtová4, and Anne Routschek2 1Research Institute for Soil and Water Conservation, Prague, Czech Republic 2Soil and Water Conservation Unit, TU Bergakademie Freiberg, Freiberg, Germany 3Czech Technical University in Prague, Prague, Czech Republic 4Research Institute of Agricultural Engineering, Prague, Czech Republic Correspondence: Hana Beitlerová (beitlerova.hana@seznam.cz)

in EROSION-2D/3D requires input parameters that can be measured or predicted with common methods (i.e., bulk density, initial soil moisture, grain size distribution, and organic bound carbon) and the skinfactor calibration parameter. The skinfactor can be determined from rainfall-runoff or infiltration experiments with the hillslope simulation tool EROSION-2D (Michael et al., 1996). This process requires extended time and demands manual labour, limiting the skinfactor determination to a relatively small number of combinations of soil and vegetation conditions. 15 Previous studies have focused on estimating skinfactors for those other than measured conditions. The studies are based on 116 rainfall experiments conducted in Saxony (Germany) between 1992 and 1995, which are published in the EROSION-3D Catalogue of Input Parameters (Parameter Catalogue) (Michael et al., 1996). Michael et al. (1996) andvon Werner (2009) estimated the skinfactors using information on German KA5 soil textural classes (Sponagel and Ad-hoc-Arbeitsgruppe Boden, 2005), initial soil saturation (dry or wet conditions), plant development stages, management practices, and field conditions. 20 All of the predictors were factorial variables. The resulting matrix of skinfactor values provides guidance for a limited number of vegetation and soil condition combinations, which is available in the Parameter Catalogue for model users. However, the statistical background of the matrix and the selection of the predictors were not published and are not traceable. For other conditions, users must estimate values by themselves from the limited and incomplete matrix. Another approach (Michael, 2000;Schlegel, 2012) was to predict skinfactors from the numeric soil input parameters of the infiltration module (i.e., clay, 25 silt, sand, organic carbon, bulk density, and soil moisture). Both studies used regression models to analyse the strongest predictors for different groups of experiments according to the soil types, management practices, and moisture conditions. The entire dataset shows the strongest correlation between the skinfactor and the bulk density, soil moisture, and silt content, but with a low statistical significance and small correlation coefficient. Analysis of specific groups of experiments (e.g., sandy soils and conservational management practices) exhibits better results, but are based on an insufficient number of 30 experiments. For this study, an R package toolbox.e3d was developed to enable automatic and batch determination of the skinfactors for multiple rainfall-runoff infiltration experiments. An extensive rainfall-runoff experiment database was processed by the package, creating a sufficient amount of data to statistically analyse the relationships between the skinfactor and other parameters describing the soil and vegetation conditions of the experiments. The aim of this study is to improve the performance of EROSION-2D/3D by providing easy to use transfer functions to calibrate the infiltration module of the model. This paper reports the skinfactor transfer functions derived from currently available data; however, this process is fully reproducible using the R code provided in the supplementary material of this paper, such that the functions can be improved and more robustly validated using the growing dataset of rainfall simulations.
2 Data and methods 5

Skinfactor
The infiltration submodule used in EROSION-2D/3D was developed by Schmidt (1996) based on the Green-Ampt infiltration approach (Green and Ampt, 1911), which includes a simplification of the infiltration process by assuming that rainwater penetrates the soil in a piston-like flow and completely saturates the available pore space. Empirical functions are used to estimate the matrix potential Van Genuchten, 1980) and saturated hydraulic conductivity (Campbell, 10 1985). A full description of the infiltration submodule is given in Schmidt (1996). As the theoretical concept of infiltration assumes a rigid soil matrix, time variable structural processes, such as soil compaction, slaking, and crusting or macropores due to shrinking and biological activities, should be considered using an empirical factor, known as the skinfactor. This factor is used to adjust the saturated hydraulic conductivity k s according to Schindewolf and Schmidt (2012) as follows: Values of the skinfactor <1 reduce the infiltration rate to consider the effects of soil slaking and crusting, as well as anthropogenic compaction. Values of the skinfactor >1 cause a positive correction of the infiltration rate, e.g., to consider increased infiltration in macropores due to soil shrinking, biological activity, or tillage impact. Two methods of deriving 20 the skinfactors from rainfall-runoff experiments were established in previous studies, both yielding slightly different values, resulting in different surface runoff rates. The first established method uses the skinfactor to adjust the amount of cumulative runoff from the plot area (skinfactor runoff ) (Michael, 2000). The second established method uses the skinfactor to adjust a certain infiltration rate, usually the final infiltration rate at the end of the experiment (skinfactor inf ) (Schindewolf and Schmidt, 2012). The best method remains a topic of debate among model developers. In this study, we used both methods to derive 25 the skinfactors for the analysis. Transfer functions for the skinfactor inf showed a better fit to the validation datasets and are therefore presented in this study.

Rainfall-runoff data
An open database for storing, maintaining, and sharing protocols from rainfall-runoff experiments is being developed in parallel to this study (Devátý et al., 2020). Currently, the database contains protocols from three working groups: The Technical The CTU data do not contain organic carbon content and bulk density and were thus not used in this study. Another 44 RISWC and TUBAF experiments were excluded from further analysis due to missing input parameters, no generated runoff, or the use of non-standard experimental conditions. Factorial predictors of crops and management practices were fractionated into many levels represented by a few to tens of cases. For better statistical representation, the predictors were categorized into subgroups based on their similar behaviour during the erosion process (Table 1). The complete and consolidated dataset for statistical 10 analysis contains 273 RISWC and TUBAF experiments. Parameters included in the statistical analysis and respective data acquisition methods used by the working groups are listed in Table 2.

Skinfactor prediction
The skinfactor has a nearly logarithmic distribution, with values ranging from 0.001 to 100 in the dataset. The assumption of normally distributed residuals in the linear mixed effects models used in this study is violated when using untransformed 15 skinfactors. Logarithmic transformation of skinfactors produces a near normal distribution for the residuals. Therefore, this transformation was used for all skinfactor values in the statistical analysis.
The dependency of the skinfactor on single predictors was tested in the correlation matrix for the numerical predictors and via an ANOVA analysis for the factorial predictors to obtain the first insight into the relationships. Numerical variables of the initial soil moisture, bulk density, and soil texture were correlated with the skinfactor. Multicollinearity was observed between the sand and silt content and between the vegetation cover and time of consolidation. The sand content was removed as the soil texture predictor has less of a correlation with the skinfactor than the silt content. The time of consolidation was removed as a parameter because it is harder to obtain for model users than the vegetation cover. Among the factorial variables, the significant impact that soil saturation (dry/wet experiments) has on the skinfactor was detected, which corresponds to the correlation 5 between the skinfactor and the initial soil moisture. Dry soil leads to lower skinfactors than saturated soils. However, it is important to consider the soil saturation not only in the context of the soil moisture (low x high), but also in the context of the state of the topsoil. While dry experiments represent the natural conditions of the soil cover, wet experiments represent the soil cover after rainfall and impacts from the destruction of soil aggregates and soil crust, loss of trapped air, or water repellence. The crop type and soil texture group also have an impact on the skinfactor, but only on the inter-level stage. For 10 the crop predictor, unlikely relations were observed. Differences between similar crop groups (e.g., catch crops versus cereals) were more significant than the differences between highly diverse crop groups (e.g., catch crops versus seedbed). Significantly different skinfactor values were also observed between working groups.
To determine the transfer functions for the skinfactor, linear mixed-effect models (Galecky and Burzykowski, 2013) were applied. All numerical soil input parameters and categorical variables used in previous studies were included in the analysis 15 as fixed effects. Furthermore, two nested random effects were included in the model to account for the interdependency and hierarchy of the data. The first random effect is the working group. Results of the experiments can be affected by the use of a specific rainfall-runoff simulator. The rainfall parameters and methodology for data acquisition differ between the working groups ( Table 2). The second random effect is the plot ID, which is nested in the working group. Both working groups usually repeat their measurements twice on an identical plot to obtain data under the dry and wet conditions. Measurements with the 20 same plot ID are thus interdependent.

Model selection
The experimental dataset was divided into the training subset, containing 75% of the randomly selected experiments, and validation subset, containing the remaining 25% of the experiments. Various models were fitted using the training subset.
Model ORIG, with factorial predictors originally used in the Parameter Catalogue, was fitted to statistically evaluate the current 25 skinfactor prediction method available for model users (Michael et al., 1996). The dataset structures used in the Parameter Catalogue and presented in this study are not identical; therefore, the equivalents of the predictors were used to remain as close to the Parameter Catalogue approach as possible (e.g., factorial predictor plant development is not available for RISWC data; therefore, it was substituted by the numerical variable, vegetation cover). STEP1-STEP3 represent the group of models manually selected using the stepwise method from the initial model containing all factorial predictors in the interactions with all 30 numerical predictors. The manually controlled backward elimination approach was followed. Single predictors with the lowest significance were continuously removed from the model while controlling for the significance of the remaining predictors and interactions, the Akaike Information Criterion (AIC) (Akaike, 1987) and the environmental sensitivity of the selected predictors. STEP1 was the most complex model, whereas STEP2 and STEP3 were selected by simplifying model STEP1 to provide a suitable model for EROSION-2D/3D users according to information on the study area and available predictors. The simplest model, i.e., STRONG, contains only the two most significant predictors. season, six rainfall events produced runoff. However, runoff was never recorded in all three plots, which shows high variability in the rainfall-runoff processes even within a very small area. The parameters of the events are presented in Table 3. Each rainfall event was modelled by Erosion-3D with the skinfactor predicted by transfer functions STEP1-3 and STRONG; for each function, the skinfactor was corrected by the positive and negative MAPE error to account for the uncertainties in the 5 predictions. 3 Results

Skinfactor prediction
Five models were fitted to evaluate the skinfactor estimation method given in the Parameter Catalogue and determine new transfer functions for predicting skinfactors using the most significant predictors (Fig. 1). Table 4 lists the evaluation of 10 the model performance based on the validation dataset and the model predictors with the coefficients for transfer function construction. The ORIG model, fitted to the predictor equivalents from the Parameter Catalogue, has low explanatory significance (variance explained by fixed effects R 2 = 0.12). Only soil saturation (dry or wet experiment) is a highly significant predictor.
The new transfer functions provide significant improvement to the accuracy of skinfactor prediction. Soil moisture and bulk density were determined as by far the most significant predictors, explaining together 79% of the skinfactor variability. The

15
skinfactor increased with an increase in both of the predictors (Fig. 2). Other significant predictors, e.g., silt content, soil texture group, and soil saturation, slightly improved the model fit. The most complex STEP1 model containing all of the significant predictors, including the interactions, explains only an additional 4% of the skinfactor variability. All four transfer functions performed well according to the interpretation of the RSR indicator by (Moriasi et al., 2007). The mean absolute percent error was between 66% and 72% for the new transfer functions while it was 192% for the ORIG function.

Error propagation
Error propagation of the predicted skinfactor for the surface runoff and sediment volume simulated by EROSION-3D was evaluated on the hypothetical 400 m long slope. The skinfactors input into ORIG model produced no runoff for 24 out of the 64 validation datasets while the skinfactors input into the new transfer functions produced no runoff only for 1-3 datasets (Fig.   3). There is not a large difference in the error propagation between models STEP1-3 and STRONG (Fig. 5). This indicates the 25 major impact of the two strongest predictors, i.e., initial soil moisture and bulk density. The median error of the surface runoff was 44-46% while that of the sediment volume was 52-56% (for the ORIG model these were 93 and 100%, respectively).
Errors below 100% characterised 78% of the datasets for surface runoff and 70% of the datasets for sediment volume, whereas, for the ORIG model, these values were 50 and 42%, respectively. and 4). The simplest model, i.e., STRONG, produced better results for certain metrics than more complex models. In general, all of the new transfer functions showed similar error propagation values, such that they can be used to predict the skinfactor.
The results suggest that the simplest function does not necessarily lead to the poorest result.

Validation with real events
Real rainfall-runoff events were modelled using the new transfer functions. To account for the potential error in the functions, each event was simulated with the predicted skinfactor and the skinfactor corrected by +MAPE error and -MAPE error.      skinfactor corrected by -MAPE to increase the infiltration rate, produced no runoff for all events. Only events 14.5. and 15.7. produced runoff (Table 3). For all of the transfer functions, the modelled runoff was within the range or close to the runoff value recorded by the trap devices. The STRONG model simulated less runoff than the other models and only the simulations with +MAPE correction produced runoff. The recorded runoff values for events 5.5., 2.7., and 11.7. are questionable, because 10 the rainfall data had very low volume and intensity, significantly lower than the erosion causing rainfall, as defined by (Janeček et al., 2012) (12.5 mm volume or 6 mm/15 min intensity). Event 29.6. had one of the highest volumes, but had a relatively long duration and low intensity. While this event produced the largest runoff, as recorded by a trap device, EROSION-3D simulated no runoff. Crust on the topsoil was recorded by field workers for the last four events, which likely initiated runoff from the low-volume and low-intensity rainfall events. The fact that runoff was never recorded in three trap devices during the same event shows the high natural variability of the rainfall-runoff process within a small area. More validation datasets  Validation at the field or the catchment scale is appropriate because the measured runoff data represent average conditions, where site-to-site changes, as recorded using the trap device, are blurred.

Discussion
The joint rainfall simulation dataset of TUBAF and RISWC provides a sufficient amount of data to statistically analyse the 10 relationships between the skinfactor calibration parameter and commonly measured soil and vegetation parameters, as well as to derive the transfer functions for the skinfactor.
The current skinfactor prediction method published in the Parameter Catalogue is based on easily and accurately measurable factorial variables, i.e., crop, management practice, soil saturation, development stage of vegetation, and soil texture class.
The results of this study show that the variables, except for the soil saturation have statistically negligible evidential influence on the skinfactor. The most significant predictors identified in this study, i.e., the initial soil moisture and bulk density, are highly variable in time and space and cannot be easily obtained. The initial soil moisture can be calculated from antecedent precipitation (Heggen, 2001) and other soil, vegetation, and relief properties (Pan et al., 2003;Zhao et al., 2011). Tramblay et al. (2011 used external software to derive initial soil moisture as an input parameter for the runoff model. The number of 5 projects developing methods and producing open data for soil moisture based on remote sensing techniques is increasing (e.g.,   in Das (2013). Ballabio et al. (2016) presented the European map of bulk density.
The relationship between the skinfactor and the soil moisture and bulk density indicates that infiltration rates are overestimated 5 at low soil moisture and low bulk density values and underestimated at high bulk density values by the infiltration module used in EROSION-2D/3D. Previous studies have also discussed the dependency on both the soil moisture and bulk density. Soil moisture has been explained by the stability of aggregates (Michael, 2000). Dry aggregates are prone to destruction by enclosed air, which becomes compressed by water infiltrating into the aggregates. The smaller particles from the destroyed aggregates then cause surface sealing and smaller skinfactors. Wet aggregates are more stable because their matrix potential is lower and 10 the infiltrating water does not produce such high, destructive pressure in the aggregates. Schindewolf and Schmidt (2012) used air trapping on a larger scale. Air trapping occurs when the wetting front enters the soil. The enclosed soil air then hinders, to a certain extent, the infiltration. A further theoretical explanation is hydrophobicity, which results from hydrophobic particles (mainly organic matter) in the soil matrix. Once dried, particles are harder to rewet than hydrophilic particles (Hallett, 2007;Seidel, 2008;Kuhnert, 2008;Schmidt, 2009). All of these effects would decrease the infiltration rates for dry 15 soils, but are not considered in model algorithms. Therefore, all of these theories can be considered reasonable explanations for the dependency of the skinfactor on soil moisture, but none of them are validated in the rainfall experiments. An alternative explanation is the misfit of the empirical estimation functions for the saturated hydraulic conductivity and matrix potential. The experimental basis behind Campbell's model is unknown (Campbell, 1985). The equations for the matrix potential estimation are based on the measurements of 40 important Belgian soil series. They represent a local dataset, which may be comparable 20 to other regions, but validation is required . Schmidt (1996) showed that these equations lack accuracy for very dry conditions (pF >4).
The existing prediction methods for the skinfactor fail to include this dependency on soil moisture. They distinguish only between dry and wet run conditions (Michael et al., 1996;von Werner, 2009), which can rather correspond with the impact of the rainfall on soil cover (e.g., soil sealing and broken aggregates) as opposed to moisture (Fiener et al., 2011).

25
The dependency of the skinfactor on bulk density is associated with macropores and surface sealing (Michael, 2000). Soils with high bulk density values are likely treated with reduced tillage practice and therefore are rich on macropores, which enhance infiltration and lead to greater skinfactors. Soils with low bulk density are likely freshly ploughed and therefore are prone to surface sealing, which hinders infiltration and leads to lower skinfactors.
Previous studies associated skinfactor values greater than one with macropores and values smaller than one with surface 30 sealing (Michael, 2000;Seidel, 2008;Schindewolf and Schmidt, 2012). This study indicates that skinfactor values do not systematically relate to these conditions; all of the experiments at dry conditions had skinfactors smaller than one, including those with reduced tillage, which tend to develop more macropores. However, this cannot be proven because surface sealing or macropore conditions were not recorded. Previous studies have attempted to determine the empirical equations for skinfactor prediction (Michael, 2000;Schlegel, 2012). Although these authors do not recommend the use of these equations (Schlegel, 2012), as well as the fact that certain predicted values are unreasonable (Lenz et al., 2018), the initial soil moisture and bulk density were identified as the most important predictors, which consistent with this study. In these attempts, experiments were 5 grouped into subsets based on texture, management practices, and the type of run to derive regression models for the subsets.
This method reduces the number of experiments and achieves higher R 2 values for each single subset, as compared with the method applied in this study, which uses categorical variables as covariables in linear models. Previous studies determined different dependencies for the prediction parameters (e.g., the intercept of soil moisture) on each single subset, whereas this study assumed an equal dependency on each parameter for the entire dataset.

Conclusion
This study aimed to increase the accuracy of the infiltration module of the EROSION-2D/3D soil erosion simulation tool by introducing new transfer functions to estimate the skinfactor calibration parameter. The relationship of the skinfactor with soil, vegetation, and farm management parameters was analysed using the linear mixed effect models based on 273 rainfall-runoff experiments. The initial soil moisture and bulk density were found to be the most important predictors, together explaining 79% 15 of the skinfactor variability. These parameters are not considered in currently available prediction methods provided in (Michael et al., 1996). Other significant predictors of soil texture (i.e., the silt content and KA5 soil texture group) and the impact of previous rain events only slightly improved the skinfactor prediction. Four transfer functions with different complexities and number of predictors were presented, such that the users can make a selection according to the available data in their study area. The proposed transfer functions present significant increases in the skinfactor prediction accuracy, as compared with 20 currently available methods (decrease in the MAPE error from 192 to 66-72%). Error propagation of the estimated skinfactors indicates substantial improvements to surface runoff and soil loss simulations. Real rainfall-runoff events were modelled by EROSION-3D with the skinfactors predicted by the proposed functions, exhibiting good model performance for events with higher total precipitation and intensity.
. This paper was compiled using the RMD-template by Nuest (Allaire et al., 2020) in RStudio (RStudio Team, 2020). The source file with all calculations performed in R (R Core Team, 2020) and not open accessible input data are available in the supplementary materials.
. AR, JD, MM, and AB made rainfall experiments, HB, JL and JD processed rainfall experiments data, JL automatized skinfactor determination,

5
HB, JL and JD made the statistical analysis, IG provided data for validation on real events, HB and JL wrote the code and prepared manuscript, AR and JK consulted the whole process.
. The authors declare no competing interests.
. This study was supported by the Ministry of Agriculture of the Czech Republic (QK1810341, MZE-RO0218) and by the European Social Fund in the Free State of Saxony (Förderbaustein: Promotionen)