Numerical evaluation of laboratory apparatuses for the study of infiltration and runoff

Runoff apparatuses (RA) are developed to study infiltration, runoff generation, and erosion processes. Several RA designs are available, but limited attention has been given to the effects of the equipment scale, initial, and boundary conditions on measured runoff. This paper presents a model-based evaluation of RAs using a finite element solution for Richard’s equation and a novel ground surface boundary condition designed to accommodate unsaturated soil behavior. The hydraulic properties of two tropical soils were considered, with multiple combinations of initial water contents, specimen dimensions, and sloping angle. The numerical exercises indicate that soils with lower air-entry values require an equilibrium stage for the establishment of initial conditions. Testing protocols with equilibrium times of 48 hours are recommended. Moisture flow produced by gravity when sloping the specimen was shown to potentially affect surface conditions and, consequently, runoff. Testing specifications to minimize the effects of specimen sloping are presented. The runoff mechanism in an RA was shown to have up to three stages, all with clear physical meaning. The third stage is an undesirable consequence of the influence of the RA’s impervious bottom. The establishment of the minimum specimen thickness that prevents boundary effects was shown to have major importance to testing results.


INTRODUCTION
Runoff is a primary variable associate with landscape evolution. The study of rainfall-runoff processes involves numerous variables and complex processes associated with infiltration, ponding, and vadoze zone behavior. Runoff depends on the soil properties, soil surface and subsurface conditions, and meteorological characteristics (Smets et al., 2008a(Smets et al., , 2008bSilva Júnior, 2015). The interaction between the soil and the weather conditions is dynamic and varies in time and space. Nevertheless, runoff assessments used in practice are predominantly based on simplified approaches (Young et al., 2009).
The need to replaced crude estimates of runoff by more accurate field and laboratory experimental data is dictated by the problem type and scale (Martinez et al., 2017;Liu et al., 2019).
Poorly representative values of the runoff coefficient (C) may result, for example, in inaccurate hydrographs, errors in the computation of subsurface moisture conditions required for ground behavior evaluations, and poor estimation of erosivity (Sen, 2008).
A direct measurement of runoff may be obtained by means of a runoff apparatus (RA), also referred to in the literature as physical model or rainfall simulator (Montebeller et al., 2001;Alves Sobrinho et al., 2008;Sousa Júnior et al., 2017;Chouksey et al., 2017). An RA must reproduce, under a controlled laboratory environment, rainfall and geotechnical conditions with characteristics that are similar to those found in the field. Such accurate representation of in-situ conditions must consider the numerous variables involved, such as initial antecedent moisture, soil type and thickness, bare surface or vegetation cover characteristics, and surface slope, among others.
There are numerous studies with respect to the characteristics of the simulated rains. These studies have dealt with the limits and tolerances allowed regarding the uniformity and diameters of water drops, and the speed and kinetic energy of the rain (Abudi et al., 2012;Spohr et al., 2015;Sousa Júnior et al., 2017;Kim et al., 2018). Meanwhile, research about the evaluation of scale effects in RA's has been presented by Martinez et al. (2017), Liu et al. (2019), and Langhans et al. (2019), among others. It is also important to note that most studies employ field apparatuses instead of laboratory equipment. A variety of conditions has been considered by previous in-situ studies, such as land covered by gravel (Poesen et al., 1994), tillage (Sadeghi et al., 2015), vegetation (Cerdan et al., 2002(Cerdan et al., , 2004Mingguo et al., 2007), multiple cover types (Smets et al., 2008a,b), and bare soil surface (Fu et al., 2011).
The existing RA's are designed based on empirical observations (Egeli & Pulat, 2011;Cecconi et al., 2015;Montrasio et al., 2016). There is an apparent lack of details regarding the design process and little effort has been directed towards the optimization of RA's geometry. None of the studies cited herein have established ideal testing conditions or made use of computer-aided design, based on numerical modeling of the rainfall-runoff interaction. Some of the main variables of interest are the characteristics of the simulated rainfall, and the scale and boundary effects. The adequacy of the apparatus should also consider testing parameters at varying initial moisture conditions. Runoff measurements may depend on the soil sample thickness, on the distance of specimen boundaries, on the uniformity of the initial water content, and finally on the uniformity of the distribution of runoff generated along the sample surface. These variables need to be carefully evaluated so that test results are adequately interpreted, assuring the representativeness of obtained runoff coefficients.
Within the context of model-based experimental design, the objective of this paper is to evaluate the theoretical behavior and performance of laboratory runoff apparatuses under various conditions, using a comprehensive and rigorous numerical model for the analysis of runoff and infiltration. The sensitivity of the runoff process to the physical characteristics of the apparatuses and to the soils tested are numerically investigated. Special attention is given to scale effects and to the determination of the minimum specimen dimensions (i.e., length and height) required to represent field conditions. Finally, the manner how testing results are interpreted is evaluated, considering the design characteristics of the devices.

MODELING METHODOLOGY FOR RUNOFF TESTS
The evaluation of the behavior of runoff apparatuses presented in this paper is based on a mechanistic, phenomenological model of water flow under unsaturated/saturated conditions. The equations governing water flow are combined with boundary conditions specially developed to represent precipitation, infiltration, and runoff rates. The infiltration and runoff rates are based on a mechanistic nonlinear relationship that considers the pore-water pressure at the ground surface which in turn depends on the precipitation rate and on the ability of the soil to allow infiltration. The following sections present the approach adopted in detail. Figure 1 presents a general view of the flow components on an RA. A time and spatially uniform precipitation rate, P P [L 3 L -2 T -1 ], is generally applied to a specimen whose characteristics (i.e., initial conditions and internal and surface properties) are representative of in-situ conditions. The specimen surface slope is determined by α D . During the precipitation event, infiltration and runoff may be produced. The point infiltration rate, I P [L 3 L -2 T -1 ] and point runoff rate, R P [L 3 L -2 T -1 ] are usually time and spatially non-uniform variables that correspond to flux rates at any given point. Note that RA's are capable of measuring only accumulated Mendes et al. runoff generated over the entire specimen surface. Therefore, the interpretation of results depends on the spatial uniformity of generated runoff and scale effects become important issues.

Water flow in saturated/unsaturated media
The water flow through the soil specimen during a runoff test may be represented by a partial differential equation (PDE) that is obtained combining the differential equation of conservation of water mass, Darcy-Buckingham flow law, and a constitutive equation describing the volume of water stored in the soil. Considering two-dimensional conditions, the following equation is obtained (Fredlund & Rahardjo, 1993): where k(θ) is the hydraulic conductivity function [L 3 L -2 T -1 ], θ is the volumetric water content [L 3 L -3 ], h is the total hydraulic head [L], h = u w /γ w + z, u w is the pore-water pressure [M L -1 T -2 ], z is the elevation [L], and t is the time [T]. The pore-air pressure has been assumed constant and the soil is considered isotropic. This equation is known as Richard's equation. The governing equation is presented using total hydraulic head as the primary variable. Two nonlinear soil properties are required, namely, the hydraulic conductivity function and the water retention or soil-water characteristic curve (SWCC). A variety of fitting or prediction functions may be employed to represent these properties, such as those proposed by Burdine (1953), Brooks & Corey (1964), Mualem (1976), van Genuchten (1980, ,  to name only a few. The adequacy of any selected property function depends on the type of soil. For this study, the SWCC equation proposed by Gitirana Jr. & Fredlund (2004) was adopted, along with Brooks & Corey (1964) equation for the prediction of the hydraulic conductivity function. The soil types considered in this study will be presented later on, along with a comparison between the proposed hydraulic conductivity function and the prediction using the van Genuchten (1980) model.

Ground surface boundary condition
The equation governing the saturated/unsaturated internal flow must be solved considering the ground surface boundary condition. The boundary condition must consider the water flow balance across the ground surface, as follows: where P P is the rate of precipitation that reaches the surface at a given point [L 3 L -2 T -1 ]; α D is the surface slope [deg], I P is the rate of infiltration at any given point [L 3 L -2 T -1 ], and R P is the rate of runoff at any given point [L 3 L -2 T -1 ].
In the context of the governing PDE, the dimensions of P P , I P , and R P correspond to flux rate values averaged over an infinitesimal area taken around a point "P". These flux rates are referred to here as "point values" and may vary with the position on the ground surface.
It is assumed that the effect of evapotranspiration can be neglected during a precipitation event. It is also important to note that the water balance expressed by Equation 2 is not a reservoir balance as is usually considered but a surface rate local balance.
The ground surface water balance requires the knowledge of two of its three flow components. Assuming that P P and α D are known, either I P or R P must be also established. The general modeling approach adopted in this paper considers that the water balance may be rendered determined by computing I P using the following conditional equation: P n 1 P n 1 D w,gs n w p w n 1 w p w,gs n w p where n and n + 1 refer to consecutive time-steps, u w,gs is the pore-water pressure at the ground surface, γ w is the unit weight of water [M L -2 T -2 ], and d P is the allowed depth of ponding [L]. A ponding depth equal to zero may be assumed, depending on the surface drainage conditions. A natural boundary condition (i.e., rate of infiltration) is applied to the ground surface as long as the pore-water pressure at the ground surface is equal or below that corresponding to the specified maximum ponding depth. During the numerical integration of the governing equations, as soon as the ground surface pore-water pressure reaches the limiting values for a given time step the boundary condition is modified to a Dirichlet type boundary condition, imposing u w = γ w d P . The flux rate is obtained using Darcy's equation and the gradients produced at the ground surface. Finally, having obtained P P and I P , the rate of runoff if computed using Equation 2.
The physical basis of the conditional boundary condition presented by Equation 3 is sound. Unfortunately, it presents at least two limitations. First, once the ground surface pore-water pressure reaches zero, such boundary condition would never allow the return to an imposed flux rate. Once the Dirichlet boundary condition is applied, the system would never return to the natural boundary condition because it would be impossible for the pore-water pressures to decreased at the ground surface. Note that actual precipitation rates that reduce over time could result in a return to negative pore-water pressures. Secondly, such boundary condition imposes a time discontinuity that may produce water balance errors.
An equivalent boundary condition that eliminates these two limitations has been proposed by Gitirana Jr. (2005) and is adopted in this study. The boundary condition is represented by the following equation: where λ controls the smooth transition between positive and negative infiltration rates, as shown in Figure 2. The larger the value of λ, the more abrupt is the transition. During a precipitation event the observed rate of precipitation is applied to the surface as long as positive values of (d P − u w,gs ) are taking place. As soon as negative values of (d P − u w,gs ) are obtained for any given time step, a negative flux is applied, removing the excess water volume and pore-water pressure. This results in a controlled oscillatory behavior that follows the physical basis originally represented by Equation 3 while producing better stability and water balance, as shown by Gitirana Jr. (2005). This boundary condition may be used with explicit or implicit time integration schemes. The usual benefits of using implicit schemes associated with a Newton-Raphson solution for the nonlinear system of equations are particularly important when this type of conditional boundary condition is used. Therefore, the implicit scheme was adopted in this study.

Ground surface total flow volumes and runoff coefficient
The total volumes of precipitation, infiltration, and runoff may be obtained by integrating the point flow rate values along the ground surface (s), and time interval (t): where Q S [L 3 T -1 ] is the instantaneous flow for a given surface, Q P [L 3 L -2 T -1 ] is the flow rate point value, and Q T [L 3 ] is the total accumulated flow over a surface and time interval.
Applying the same integration procedures, point, surface, and total runoff values (R P , R S , and R T respectively) may be computed as follows: Runoff coefficients corresponding to point, surface, and total flow values (C P , C S , and C T respectively) are defined as follows: The various definitions of runoff coefficient are useful when interpreting the results obtained using runoff apparatuses. The adequacy of testing measurements taken over a surface and during a time interval may be investigated comparing simulated values of point, surface, and total runoff.

Solution of governing equations
The solution for the unsaturated soil flow model was implemented using FlexPDE version 6, a general-purpose partial differential equation solver (PDE Solutions, 2016). Gitirana Jr. (2005) and Silva Júnior (2015) verified similar implementations of ground surface flow models using FlexPDE. The user must Figure 2. Normalized infiltration rate as a function of the ground surface pore-water pressure.

5/16
develop a script presenting the governing equations, domain, boundary conditions, and initial conditions of the problem. FlexPDE builds the finite element model, solves the system of algebraic equations, and enables graphical output. First and second order PDEs in steady-state or transient conditions and for nonlinear problems can be solved. Cartesian, cylindrical, and spherical coordinate systems may be used, for one-, two-, or three-dimensional problems.
The general mass-conservative mixed form of the governing PDE, proposed by Celia et al. (1990) was adopted. The PDE presented in Equation 1 was coupled with the boundary condition presented in Equation 4. The Galerking weighted-residuals approach was adopted for the finite element solution. Isoparametric six-noded triangular elements were used for a quadratic polynomial representation of the primary variable. A three-step Green method was used for the time integration of the finite element equations. Newton-Raphson's method was used for the solution of the nonlinear system. Automatic mesh and time-step refinement procedures available in FlexPDE were used.
The discretized equations at the mesh nodes are formed by taking the integral of the weak form of the governing PDE, G(h). The individual cell integral is used as a measure of the mesh quality. This integral must be satisfied for each cell, according to a convergence criterium corresponding to the ratio dh/h, where dh is the error in h. The adopted tolerance for dh/h was 0.0001. If the convergence criterium is not met in a given cell, the mesh is refined by splitting the element. Note that only the diagonal components of the Jacobian matrix of derivatives of the Galerkin integral are used.
Automatic time-step refinement was used to meet the temporal error tolerance. An additional timestep of values of h, previous to the three points of the quadratic Green solution, is used. The size of the cubic term of a cubic fit implies the error in the quadratic solution. This error value is used to change the timestep, following a tolerance of 0.0001. More details regarding the refinement criteria may be found in PDE Solutions (2016). The adopted number of elements ranged from 12,210 to 22,130. The initial time step varied between 0.00001 and 0.001 minutes and changed over the course of the solution by a factor of up to 100, typically. The automated mesh and time refinement are related to the hydraulic head gradients, leading to a varying number of elements and varying time-steps. The observed hydraulic gradients and the resulting number of elements depended on the soil properties and initial condition, among other variables. For example, steeper SWCCs lead to higher gradients and result in denser meshes along the infiltration front. Mass balance was verified in all analyses and the obtained errors were always lower than 0.1%.

Proposed analyses
A general RA must be designed considering interdisciplinary concepts and requirements, mainly those from the hydrology and geotechnical disciplines. From a hydrologist point of view, the runoff coefficient is more commonly acknowledged as an average value along a relatively large surface. However, from a geotechnical modeler's perspective, runoff is often viewed as a local value at any point of a given surface. Both disciplines may benefit from small-scale studies and modeling as it allows an understanding of the physical mechanisms controlling runoff.
In order to adequately interpret runoff in a local scale, surface and subsurface measurements are required and numerical interpretation of test results is desirable. It will be demonstrated in this paper that the combination of surface and subsurface observations provides a better understanding of soil-weather interactions. The basic design of an RA is shown in Figure 3 and follows the recommendations presented by ASTM (American Society for Testing and Materials, 2015). The RA must be composed of a structure capable of supporting the specimen load and imposing varying surface slopes that are conveniently adjustable. The specimen container must be capable of withstanding the specimen load and compaction effort in the case of remolded materials.
The numerical analysis exercises presented herein are aimed at establishing the behavior of a specimen in terms of its pore-water pressure, degree of saturation and runoff coefficient. The specimen behavior will be evaluated considering variations in specimen length, L D , thickness, H D , and inclinations, α D . Table 1 presents the variations in RA geometry that are used in the parametric studies presented herein. The selected testing specifications are considered feasible in typical laboratory conditions.
A rainfall intensity of P P = 200 mm h -1 was adopted, corresponding to an extreme precipitation with a typical duration of 20 minutes and return period of 100 years for the city of Goiânia, GO, Brazil (Costa & Prado, 2003). Equation 4 was used as the boundary condition for the top of the specimen. Figure 3 shows that the boundary conditions adopted for the bottom and sides of the specimen correspond to impermeable surfaces (q = 0).
The use of an RA requires testing stages for specimen setup and rainfall simulation. It is desirable to interpret all moisture flow processes that take place during sample manipulation and evaluation. Three common testing stages were simulated, as depicted in Figure 3: • Stage 1: specimen setup and equilibrium (Figure 3a). This stage is only relevant to remolded soils that have uniform water contents right after specimen preparation and setup. The main objective is to evaluate the movement of moisture due to gravity and evaluate the time required for equilibrium. The equilibrium time was defined as the time when the difference between the hydraulic heads at the bottom and at the top of a specimen is lower than 0.01%. The selected initial conditions are presented in Table 2 and correspond to a constant water content and pore-water pressure, a horizontal sample (α D = 0), and no rainfall (P P = 0); • Stage 2: specimen sloping prior to rainfall (Figure 3b).
Evaluation of the changes in the equilibrium conditions generated by changing the slope of the sample. The aim is to evaluate now much time it takes until the equilibrium generated in the first stage is significantly altered after specimen sloping. The disturbance time is defined as the time required for the development of a variation of 10% in the pore-water pressures at the specimen top surface at its elevated top corner. It is important to emphasize that moisture movement due to specimen sloping could interfere with runoff measurements during the next testing stage. The analysis of stage 2 will demonstrate in which scenarios undesirable overlapping effects of precipitation and moisture redistribution due to specimen sloping would take place. Multiple specimen dimensions and slopes were considered, as presented in Table 1. The initial pore-water pressure conditions were assumed to be hydrostatic in the horizontal specimen position; • Stage 3: rainfall simulation (Figure 3c). Once hydraulic equilibrium is achieved in the first stage, the rainfall (P P = 200 mm h -1 ) is applied considering the multiple initial conditions and specimen geometries presented in Tables 1 and 2. The coefficient of runoff (C) is evaluated at the microplot scale. The rainfall simulation is aimed at allowing the interpretation of the water budget and runoff phenomena along the specimen surface using Equations 10 to 12.
To perform the numerical simulations, it is necessary to establish the soil-water characteristic curve (SWCC) and the hydraulic function of the soil. Two different tropical soils from Goiânia, GO, Brazil were considered, both with bimodal SWCCs. The first material, denominated Soil 1 (S1), is a highly weathered and permeable soil that was studied by Araújo (2013). The second material, denominated Soil 2 (S2), is a compacted residual soil that was studied by Laguna (2015). The properties of the two soils used in the numerical simulations are presented in Table 2 and Figures 4 and 5. The selected properties offer a wide range of possible soil characteristics, making it possible to evaluate the behavior of RA's under different conditions. The SWCCs where represented using a best-fit analysis using the equation proposed by Gitirana Jr. & Fredlund (2004). The hydraulic conductivity function, k(θ), was estimated using the equation proposed by Brooks & Corey (1964), based on the SWCC and the value of the saturated hydraulic conductivity, k s ( Figure 5 and Table 3). The estimation model was applied using a pore-size distribution index obtained for the first branch of the soil-water characteristic curve (Table 3). This first branch corresponds to the macropores and is assumed to account for most of the liquid water flow.

Stage 1
The specimen length has no effect on the moisture movement and equilibrium time for stage 1 because, with the specimen in the horizontal position, the problem is one-dimensional. Therefore, only one specimen length (L D = 1.0 m) was considered for the simulation and discussion of results related to stage 1. Figure 6 shows the changes in matric suction and degree of saturation over time at the top of the specimen. The presented curves are chosen as representative curves for all scenarios.
It can be seen that the hydraulic equilibrium process changes significantly depending on the soil type and initial remolding degree of saturation. The differences can be attributed mainly to the permeability and air-entry value of each material. Similar behavior has been observed by Bittelli & Flury (2009), Dexter et al. (2012 and Wang et al. (2015). The results show that materials with a lower degree of saturation, having a lower average permeability, require a longer period to reach equilibrium. This is particularly obvious for soil S1. The observed changes in degree of saturation for soil S1 mean that the equilibrium period in stage 1 is necessary as a preparation for the next testing stages.
Soil S2 presents no significant variation in matric suction and degree of saturation along the surface of the specimen,   where: ψ b1 is first air-entry value; ψ res1 is first residual suction; S res1 is first residual degree of saturation; ψ b2 is second air-entry value; S b is second degree of saturation; ψ res2 is second residual suction; S res2 is second residual degree of saturation; a is additional parameter; λ is controlling variable of the soil infiltration transition.
as shown in Figure 6b and 6d. However, the extremely small variations in matric suction and degree of saturation take a long time to develop. Regardless of the equilibrium time, the equilibrium process for soil S2 involves negligible changes in degree of saturation and should not affect the behavior of the specimen in the next testing stages. That means that stage 1 could be skipped for soil S1. Figure 7 summarizes the obtained equilibrium times for soils S1 and S2, respectively. For soil S1 and considering the initially saturated condition and H D = 0.5 m, the obtained equilibrium time is only 50 minutes. For S i = 75 and 50% the equilibrium times reach 26 and 48 hours, respectively. The obtained equilibrium times can be considered reasonable for the sample preparation in an RA. For soil S2 the equilibrium times may reach values that seem unfeasible, but stage 1 could be skipped as explained previously. Equilibrium times are not found in literature for similar experiments. However, for smaller scale experiments, such as pressure plate tests, the minimum equilibrium times recommended are 24, 48, and 72 hours for matric suctions below 100, 500, and 1500 kPa, respectively (American Society for Testing and Materials, 2016). Complete tests may take from a few days to several weeks (Bittelli & Flury, 2009;Dexter et al., 2012;Wang et al., 2015).
The effect of specimen thickness and initial degree of saturation is noticeable. The specimen thickness affects equilibrium times as a result of longer moisture movement paths. The initial degree of saturation, on the other hand, controls the initial and average unsaturated hydraulic conductivity. As a result, larger and drier specimens require significantly higher equilibrium times.
The main results from the analysis of stage 1 can be summarized as follows: a) the equilibrium process after specimen remolding must consider both the time required for moisture equilibrium and the magnitude of variations of degree of saturation; b) the higher the required equilibrium time, the smaller the predicted change in degree of saturation; and c) for soils S1 and S2 an equilibrium time of 48 hours would be adequate for all specimen thicknesses and initial degrees of saturation.

Stage 2
Due to the small influence of the specimen length perceived in stage 1, the next stage was also only analyzed for L D = 1.0 m. The perturbation times will be compared against the typical rainfall simulation time of one hour. A perturbation time longer than one hour implies that stage 3 (rainfall simulation) may follow stage 1 without any interference caused by specimen sloping. Low perturbation times are obviously undesirable and would produce heterogeneous surface conditions during stage 3 testing. It is interesting to note that previous studies (Egeli & Pulat, 2011;Cecconi et al., 2015;Montrasio et al., 2016) have completely neglected the possibility of a perturbation in the equilibrium state due to the sloping of the specimen. Figures 8 and 9 show that specimens of both soil types and with S i = 50% or lower would not suffer any surface disturbance after specimen sloping for a period of at least one hour. Specimens with S i = 75% would also not present significant moisture movement during a period of one hour, except for soil S1 with slopes of 30º (for thickness 0.15 m) and 45º.
The obtained results may be interpreted based on the SWCC of the soils analyzed (Figure 4). Relatively small variations of suction in soil S1 can cause significant changes in water content, depending on their initial condition. For example, a slope of 45º for L D = 1.0 m produces a maximum difference in elevation at the specimen surface of 71 cm, resulting in a maximum matric suction of 7.0 kPa. This value of matric suction, higher than the air-entry value of 2.0 kPa, is enough to remove water from the macropores.
All the simulations performed with soil S2 present extremely long perturbations times, due to the relatively low permeability, and the results have not been plotted. The perturbation due to sloping for soil S2 results in negligible changes in degree of saturation, similar to what was observed for stage 1. The behavior of soil S2 after specimen sloping can be justified by its air-entry value, equal to 90 kPa. The maximum suction values produced when the RA is inclined are smaller than the air-entry value, not allowing the redistribution of water in the specimen. Lower initial degrees of saturation of 75 and 50% obviously correspond to higher matric suctions. For these scenarios with lower degree of saturation, the negligible disturbance observed is due to the fact that large variations in matric suction would be required to produce significant changes in water content. This behavior is consistent with what is presented in most hydraulic conductivity models, such as Mualem (1976), van Genuchten (1980 and . The specimen condition with S i = 100% is not usual in infiltration/runoff RA tests with closed bottom boundary condition. The absence of basal flow would not allow any additional infiltration, resulting in C = 1 (full tank condition). However, the S i = 100% scenario was considered herein in order to establish a reference limiting condition. Figures 8 and 9 show for soil S1 and S i = 100% very early redistribution of pore-water pressure and degree of saturation, which would interfere with stage 3 testing results. For soil S2, only the most inclined  condition seems to present significant surface disturbance, as shown in Figure 8b.
In summary, for most of the typical testing conditions, considering usual field characteristics (S i ≤ 75%), there would be no interference of water redistribution processes due to specimen sloping for a period of one hour. For soil S2 there is no possibility of significant changes in degree of saturation due to specimen sloping. However, testing conditions for soil S1 must be considered with greater care and may require that stage 3 be limited to a shorter period of time, such as 30 minutes.

Stage 3
Considering that S i = 100% is not representative of usual RA testing conditions and that the behaviors for S i equal to 50 and 75% are similar, stage 3 was analyzed only for S i = 50%. Different specimen lengths, thicknesses, and slope angles were considered (Table 1). Figure 10 presents the behavior of soil S1, considering H D = 0.15 m and P P = 200 mm h -1 = 3.33 mm min -1 . The plots show a fast increase of the degree of saturation at the top of the specimen, reaching full saturation after approximately three minutes. A delayed and slightly slower increase of the degree of saturation is observed at the specimen base.
The generation of runoff can occur in two ways: a) when the effective precipitation (P P ) is greater than the surface infiltration capacity (I P ), regardless of the influence of the soil profile or stratification; or b) when the soil storage capacity has already been reached (full-reservoir). The first condition is related to the surface hydraulic conductivity and has been considered by Green & Ampt (1911) and Horton (1933), among others, using empirical or semi-empirical approaches. The second situation depends on the soil profile characteristics. Complex site conditions with highly conductive or impervious layers may produce different storage and internal redistribution conditions. An RA may artificially impose an undesired full-reservoir condition if the specimen thickness is small and the testing time is relatively long. It is important to note that the specimen thickness and the depth of the RA impervious bottom are often a matter of convenience. The RA dimensions need to consider the difficulties associated with preparing, transporting, and manipulating heavy specimens (Egeli & Pulat, 2011;Mendes, 2019). As a result, it is important to identify when and how the imposed bottom boundary condition affects the runoff evaluation.
The runoff curves shown in Figure 10 present three distinct behavior features over time. From the start of the test up to a time between 3 and 7 minutes there is no runoff generation as the unsaturated soil is able to infiltrate all the precipitation. Once runoff starts, it gradually increases and tends to an equilibrium value controlled by the hydraulic gradients produced that tend to constant values. However, a third factor controlling runoff is observed between times 25 and 30 minutes. These times correspond to the moment when the infiltration front reaches the impervious specimen bottom and can no longer advance. It is interesting to note that standard conceptual infiltration-runoff models, such as Green & Ampt (1911), Horton (1933), Philip (1957) and SCS (Soil Conservation Service, 1973) do not combine all of the mechanisms observed in the presented simulations. Figure 10. Stage 3 -pore-water pressure, degree of saturation, and point runoff for S i = 50% and H D = 0.15 m.

11/16
Full saturation at the specimen base coincides with a constant rate of maximum point runoff (R P ) at the three observed specimen points. The average value of maximum point runoff can be calculated based on the precipitation and surface slope, using Equation 7 and assuming no infiltration. For a value of α D of 15º and 45º the average values of runoff are 3.22 and 2.35 mm min -1 , respectively. Figure 10 shows that, under the imposed conditions, specimens with a thickness of only 0.15 m reach full saturation under 60 minutes, influencing the runoff (R P ) generation and consequently the C P value. The testing conditions become influenced by the specimen thickness and are only representative of specific field conditions where basal flow would be impeded by the presence of a low permeability layer at a depth of 0.15 m. Such dimensions for RA specimens may not produce representative runoff values for other field conditions. Similar studies (Greco et al., 2010;Moussouni et al., 2012;Madi et al., 2013;Montrasio et al., 2016) have adopted specimens with relatively small thicknesses (< 0.15 m). Therefore, the results presented in these studies should be interpreted with caution, because their boundary condition may not represent field conditions.  Table 1. Only the extreme values of L D = 0.5 and 2.0 m are shown. The influence of the bottom impervious boundary condition is no longer observed for the thicker specimens. It may also be observed that the length of the specimen (L D ) did not influence the test results and the values of point runoff (R P ). The bottom boundaries presented some increase in their degrees of saturation for some scenarios. However, the bottom boundary conditions did not seem to influence the moisture flow dynamics at the top surface and the runoff values. In summary, scale effects don't seem to arise from the selected values of L D , but may be an issue for small values of H D . It is important to emphasize that these observations are only valid for a testing period of 60 minutes. Longer or shorter testing periods would allow or required different testing specifications. Figures 10 to 12 show that point runoff is not uniform along the specimen surface. The non-uniformity is more pronounced for higher sloping angles. However, the central point value is approximately equal to the average value obtained from the two edge points. As a result, the total runoff obtained by collecting runoff produced over the entire surface is representative of average specimen conditions. Figure 13 presents the total surface runoff coefficient (C T ) for multiple specimen and RA conditions. The first obvious observation is that the total runoff coefficient is not a single value as typical assumed (Sen, 2008 13/16 but a dynamic parameter that evolves during a given precipitation event. For H D = 0.15 m, the C T values increases rapidly after around 25 minutes of rainfall because of the influence of the RA bottom boundary conditions. For H D ≥ 0.30 m, those boundary conditions became less important because, for the same α D value, C T values are practically the same.
The shape of a C T curve seems to be a good indicator of the quality of the testing specifications. Curves that present two humps, as those obtained for H D = 0.15, indicate that bottom boundary effects have influenced runoff measurements. The remaining curves presented in Figure 13 have shapes that correspond to adequate testing conditions.
The values of C T depend on the specimen slope angle Lower slope angles produce higher coefficients of runoff. This observation may seem counterintuitive at first. However, this behavior is a logical consequence of the larger fraction of the total precipitation values that reaches the ground surface when α D is lower. As more precipitation reaches the surface, the increase in the degree of saturation becomes faster, producing higher values of runoff.

CONCLUSIONS
Computer-aided design has been used to assess runoff testing conditions and evaluate the capabilities and limitations of a proposed experimental protocol. It was shown that the initial conditions of an RA specimen require careful examination. Hydraulic equilibrium prior to testing depends on the material properties, initial degree of saturation, and specimen dimensions. Attention must also be given to the possibility that moisture redistribution due to specimen sloping may overlap with water flow processes during runoff testing, altering the measured runoff.
Runoff curves may be affected by the impervious bottom of the RA. This undesirable boundary effect imposes an important limitation to the apparatus capability to represent field conditions. This study has shown that the minimum specimen thickness must be investigated prior and during experimental investigations.
Simulated runoff testing results were shown to be insensitive of specimen length, at least for the ranges considered herein. As a consequence, it's recommended that the specimen length in an RA be specified as a function of the capabilities and limitations of the rainfall simulation component of the system. Smaller surface areas are preferable to produce uniform artificial rainfalls.
The runoff coefficient was shown to depend on the surface slope. Flat slopes produce higher runoff due to the increase in the precipitation effectively applied to the horizontal projection of the ground surface. The lower slopes were also found to produce more uniform point runoff distributions, which is desirable in terms of testing data interpretation. However, the non-uniform point runoff distributions obtained when the specimen is tested using higher slopes can still be evaluated in an unbiased manner using the total runoff as it corresponds to the average value of a linearly distributed point runoff.
In conclusion, prior numerical simulations should always be performed when the soil hydraulic properties are available, in order to specify ideal testing duration and specimen thickness. Additionally, experimental runoff curves should be interpreted following the general rules presented herein and the testing specifications should be updated when boundary effects are observed in the runoff curves.