Shallow water table effects on water, sediment, and pesticide transport in vegetative filter strips – Part 1: nonuniform infiltration and soil water redistribution

Vegetation buffers like vegetative filter strips (VFS) are often used to protect water bodies from surface runoff pollution from disturbed areas. Their typical placement in bottomland often results in the presence of a seasonal shallow water table (WT) that can decrease soil infiltration and increase surface pollutant transport during a rainfall/runoff event. Simple and robust components of hydrological models are needed to analyse the impacts of WT in the landscape. To simulate VFS infiltration under realistic rainfall conditions with WT, we propose a generic infiltration solution (Shallow Water table INfiltration algorithm: SWINGO) based on a combination of approaches by Salvucci and Entekhabi (1995) and Chu (1997) with new integral formulae to calculate singular times (time of ponding, shift time, and time to soil profile saturation). The algorithm was tested successfully on 5 distinct soils both against Richards’s numerical solution and experimental data in terms of infiltration and soil moisture redistribution predictions, and applied to study the combined effects of varying WT depth, soil type, and rainfall intensity and duration. The results show the robustness of the algorithm and its ability to handle various soil hydraulic functions, and initial non-ponding conditions under unsteady rainfall. The effect of a WT on infiltration under ponded conditions was found effectively decoupled from surface infiltration/excess runoff processes for depths larger than 1.2 to 2 m, shallower for fine soils and shorter events. For non-ponded initial conditions, the influence of WT depth also varies with rainfall intensity. Also, we observed that soils with a marked air entry (bubbling pressure) exhibit a distinct behaviour with WT near the surface. The features and good performance of SWINGO support its coupling with an existing VFS model in the companion paper, where the potential effects of seasonal or permanent WTs on VFS pollutant transport and control are studied.


Introduction
The use of vegetative filter strips (VFSs) can reduce sediment and surface runoff (RO) pollutant (i.e., sediment, colloids, nutrients, pesticides, pathogens) movement into receiving water bodies.The dense vegetation-soil system reduces runoff pollutants in three ways by increasing (a) soil infiltration that reduces total runoff volume (and dissolved runoff pollutants), (b) surface roughness that reduces surface velocity and produces settling of sediment and sediment-bonded pollutants, and (c) contact between dissolved and particulate pollutants with the soil and vegetation surfaces that enhances their removal from runoff (Muscutt et al., 1993;Muñoz-Carpena et al., 1999, 2010;Dosskey, 2001;Fox et al., 2010;Yu et al., 2013;Lambretchs et al., 2014;Wu et al., 2014).The efficiency of VFS in trapping pollutants is heavily influenced by the highly variable spatial and temporal dynamics introduced by site-specific combinations of soil, climate, vegetation, and human land use.For the case of runoff pesticides, these influences have been recognised in multiple field stud-Published by Copernicus Publications on behalf of the European Geosciences Union.
ies (Lacas et al., 2005;Reichenberger et al., 2007;Poletika et al., 2009;Sabbagh et al., 2009).Other effects like hydraulic loading under concentrated flow conditions (Fox et al., 2010) or timing of the pesticide application (Sabbagh et al., 2013) can also result in reduced filter trapping efficiencies.As these systems are complex, the practice of using generic, simple regression equations relating the reduction efficiency of pollutants with VFS physical characteristics (i.e., length, slope) is often inadequate (Fox and Sabbagh, 2009).
Mechanistic understanding of VFS behavior has advanced significantly in the last 20 years, and numerical simulation tools are available to analyze this important best management practice (BMP) under upland field conditions where runoff is governed by excess rainfall and field inflow processes (Muñoz-Carpena et al., 1993, 1999;Abu-Zreig, 2001;Muñoz-Carpena and Parsons, 2004;Poletika et al., 2009;Sabbagh et al., 2009;Carluer et al., 2017).A recent linked mechanistic model has investigated multiple input factors and their relative importance and uncertainties on the predicted reduction of runoff, sediments, and pesticides (Fox et al., 2010;Lambrechts et al., 2014;Muñoz-Carpena et al., 2010, 2015).
However, because of their location near or at the riparian zone, VFS can at times be bounded by a seasonal shallow water table (WT) (Borin et al., 2004;Ohliger and Schulz, 2010).Examples of areas where these conditions exist either seasonally or on a more permanent basis are humid coastal flatland zones, floodplains near water bodies, and soils with limiting horizons resulting in perched WTs.Generally, capillary effects from a WT can reduce infiltration and increase subsequent runoff processes as well as have a major effect on contaminant transport to surface waters (Gillham, 1984).In spite of the potentially important environmental impacts of the presence of shallow water under VFS, there is a dearth of studies addressing this problem either experimentally or mechanistically.Several authors suggest the importance of this factor in VFS experimental studies (Lacas et al., 2005;Arora et al., 2010) or when designing or implementing this field BMP (Simpkins et al., 2002;Dosskey et al., 2006Dosskey et al., , 2011)), but they do not provide a mechanistic interpretation.Some authors suggest that the reduction of infiltration and VFS efficiency can be problematic for seasonal WT depths above 2 m, typical of hydromorphic soils (Dosskey et al., 2006(Dosskey et al., , 2011;;Lacas et al., 2012).As cited by Salvucci and Entekhabi (1995), the importance of accounting for areas of WT effects in water balance and runoff studies has been recognized for a long time, and specialized analysis and simulation approaches have been proposed by numerous authors (for example, Vachaud et al., 1974;Srivastava and Yeh, 1991;Salvucci and Entekhabi, 1995;Chu, 1997;Basha, 2000).
In spite of the ubiquity and importance of these areas and previous specialized analysis and modeling efforts, commonly used field and watershed hydrological models are limited when describing infiltration and soil water redistribution with WT (Beven, 1997, Liu et al., 2011).Among ex-isting simulation approaches, solutions to the fundamental Richards (1931) partial differential equation (RE) can describe the infiltration and redistribution of water in soil, including the specific case of when a system contains a WT.However, RE does not have a general analytical solution and its application in real-world systems requires computationally intensive numerical approximations that can result in mass-balance and instability errors in some cases (e.g., for coarse soils and highly dynamic boundary conditions) (Celia et al., 1990;Paniconi and Putti, 1994;Miller et al., 1998;Vogel et al., 2001;Ross, 2003;Seibert et al., 2003).As a result, soil infiltration is often modeled in field and watershed models using simpler physically-based approaches (Jury et al., 1991;Smith et al., 1993;Haan et al., 1994;Singh and Woolhiser, 2002;Talbot and Ogden, 2008;Ogden et al., 2015).One of the most often used approaches in hydrologic modeling is the Green-Ampt (1911) model adjusted to account for variable rainfall (Mein and Larson, 1973;Chu, 1978;Skaggs and Khaleel, 1982).The model has the advantages of being computationally efficient and that its parameters can be directly estimated from physical measurements or derived indirectly from soil texture (Rawls et al., 1982(Rawls et al., , 1983)).However, the limitation of the original Green-Ampt model is that it assumes isotropic soil with uniform initial moisture content and saturated "piston" infiltration.Even with these nonrealistic assumptions, if effectively parameterized, this method still generates useful and reliable results compared with other numerical and approximated approaches (Skaggs et al., 1969;Mein and Larson, 1973).Considering its advantages, Bouwer (1969) highlighted the utility of this method when taking into account the computational trade-offs with RE solutions.
Extensions of the Green-Ampt model beyond its initial assumptions have enabled its application to other natural infiltration cases, such as nonuniform soil profiles (Bouwer, 1969;Beven, 1984) and multistorm infiltration and redistribution (Ogden and Saghafian, 1997;Smith et al. 2002;Gowdish and Muñoz-Carpena, 2009).A particularly important case where an extension of the original assumption of the Green-Ampt model is necessary is when there is a WT.In general, depth-averaged soil moisture values in traditional infiltration equations like Green-Ampt (i.e., semi-infinite, uniform initial soil moisture) overpredict infiltration estimations when the soil is bounded by a WT.This is due to the difficulty in obtaining an equivalent initial uniform soil water content that effectively represents the real nonuniform water content condition with WT (Salvucci and Entekhabi, 1995;Chu, 1997).Recently, Liu et al. (2011) presented a modification to Craig et al.'s (2010) nondimensional form of the Green-Ampt model to account for the presence of a WT.Although this modification is shown to provide acceptable results compared with a RE solution for a range of WT depths, the method assumes an initial uniform soil water content profile, and its performance relies on an empirical correction between RE and standard Green-Ampt results.Alternatively, previous works (Childs, 1960;Holmes and Colville, 1970;Duke, 1972) have suggested describing the soil water redistribution over a WT as an equilibrium hydrostatic condition (Fig. 1).This approach assumes a linear relationship of soil matric potential (h, [L]) and soil depth (z, [L], positive downwards from the surface), whereby the nonuniform water content of the soil (θ , L 3 L −3 ) is described by the soil water characteristic curve, θ = θ (h) (Jury et al., 1991): where L [L] is the depth of the fixed shallow water table below the soil surface (i.e., the distance from the surface).
Based on these initial and boundary hydrostatic equilibrium conditions, Chu (1997) proposed an incremental calculation technique to evaluate infiltration into ponded soils with a WT.This calculation relies on Bouwer's (1969) expression of the Green-Ampt equation that accounts for infiltration of water into a nonuniform soil as follows: where t [T] is time since the beginning of the event, θ s [L 3 L −3 ] is the saturation water content, f [LT −1 ] is the rate of surface infiltration, and z F [L] is the wetting front depth.Following Neuman (1976) and Chu (1997), after substitution of the equilibrium condition (Eq.1), the cumulative (F , [L]) and instantaneous infiltration (f , [LT −1 ]) can be calculated as follows: where the subscript "p" denotes under-ponding or "capacity", i.e., when the flux at the surface is not limited by available water and is therefore maximum for each time; K s and K(h) [LT −1 ] represents the soil saturated and unsaturated hydraulic conductivity, respectively.Chu (1997) proposed the solution to Eqs. ( 2)-( 4) using sufficiently small increments of z, z = z i − z i−1 .If an initial value of F 1 and f 1 for the first z (from the surface to a small depth) is known, then successive values of time (t i = t i−1 + t) for each z can be approximated by substituting Eq. (3) into Eq.(2) as follows: .
(5) Chu (1997) further proposed that a valid initial step could be obtained by assuming standard Green-Ampt conditions (i.e., piston flow) from the surface, hydrostatic equilibrium of the surface water content with the WT (θ o ), and calculating the suction at the wetting front (S av ) as follows (Bouwer, 1964): Vachaud et al. (1974) was able to use experimental data to test the solution of this equation successfully.However, their experimental data did not allow enough time to determine how the model would respond when the wetting front reaches L.
An elegant and useful approximate solution to ponded infiltration with the WT was proposed by Salvucci and Entekhabi (1995).Their solution is based on the assumptions of initial hydrostatic equilibrium and uses Philip (1957) integral approximation of RE (Fig. 1).advantageous, as it describes not only the infiltration but also soil water redistribution during infiltration, and the characteristics of the wetting front as it moves towards the WT during long events.In addition, the method assumes a more realistic piecewise linear wetting front with a variable slope during infiltration (α in Fig. 1).This algorithm was successful when compared with RE solution for three different soil types and when tested with the soil moisture profile data from Vachaud and Thony's (1971) experiments.However, the applicability of the algorithm for coupling with commonly used hydrological models is limited as it requires ponded conditions, Brooks and Corey's soil water function (Appendix Eq.A1), and similarly to the original Green-Ampt it requires an implicit solution.
The overall objective of this work and its companion paper (Lauvernet and Muñoz-Carpena, 2018) is to analyze the impact of the presence of a WT on VFS efficiency.In this first paper, we will expand the Green-Ampt-based infiltration solution to soils bounded by a WT under variable rainfall with no initial ponding.We accomplish this by combining Salvucci and Entekhabi (1995) and Chu (1997) approaches with a generic solution technique and developing novel integral formulae to calculate the singular times (time to ponding, t p , shift time, t 0 , and time to column saturation, t w ) for soils with no initial ponding.We assess the ability of the simplified method to accurately predict surface infiltration and water content predictions for a variety of soils compared with RE numerical solutions and previously pub-lished experimental data.An illustrative example of calculation during an unsteady rainfall event is also presented along with examples of applications of the proposed algorithm to analyze the effects of WT depth.In the companion paper, we couple the new shallow water infiltration algorithm with an existing VFS numerical model (VFSMOD, Muñoz-Carpena et al., 1999, 2010, 2015) and analyze the effects on runoff, sediment, and pesticide removal efficiency.
2 Proposed algorithm 2.1 Infiltration rate in soils bounded by a WT with a nonponded initial state and subject to constant rainfall In general, the infiltration rate (f , LT −1 ) of a WT-bounded soil with uniform rainfall rate (i, LT −1 ) and no initial surface ponding will have a similar profile to the example shown in Fig. 2a, described by the following: The identification of three singular times during the infiltration calculations is necessary for a solution to Eq. ( 7).These singular times are time to reach ponding (t p ) (which depends on the shift time, t 0 , described later) and time to column saturation (t w ), when the wetting front approaches the WT capillary fringe at depth z w (see Fig. 1).The effective saturation depth z w relies on L and soil air entry pressure (h b ), z w = L − h b (z w ≥ 0; z w = 0 when L < h b , i.e., the soil is effectively saturated by the capillary fringe).Often, h b is set at 0 (i.e., z w = L), even if some of the soil characteristic functions take the air entry pressure into account (Brooks and Corey, 1964;Clapp and Hornberger, 1978).At t w , the soil column is saturated and the rate of infiltration sharply drops to f w , or i if i < f w (Fig. 2a).The value of t w depends on L and the slope of K(h) (Salvucci and Entekhabi, 1995).If the WT is very shallow, the time to saturation t w can occur before the time to ponding.Salvucci and Enthekabi (1995) and Liu et al. (2011) initially proposed that the infiltration rate is equal to f w = K s , when t ≥ t w , meaning that the vertical hydraulic gradient at the initial WT is 1.However, in most field situations when the wetting front has reached the WT, the profile's hydraulic gradient is less than 1 and the proposed solution might overestimate the final infiltration rate.Instead, another solution is to consider that for t ≥ t w the infiltration flow at the surface (Q f ) is controlled by lateral drainage flow (Q L ) at the downslope boundary of the simulated soil elementary volume (Fig. 1b), applicable to floodplain conditions typical of VFS.If we consider that the soil profile is saturated at t ≥ t w , with an effective saturation depth z w = L − h b , following Dupuit-Forchheimer assumptions (Van Hoorn and Van Der Molen, 1973) the discharge (Q f = Q L ) can be estimated as follows: where K sh is the lateral (horizontal) soil saturated hydraulic conductivity, w and b are the width (VFS dimension perpendicular to the flow) and length (VFS dimension in the flow direction) of the VFS surface area, and S is the slope of the initial WT.In hillslope hydrological modeling S is typically assumed to equal soil surface slope (S o ) (Beven and Kirkby, 1979;Vertessy et al., 1993).If the position of the infiltration elementary volume is close to a draining stream where S> S o , Eq. ( 8) may underestimate the infiltration rate and a 2-dimensional drainage approach like Hooghoudt (1940) equation should be used instead (Kao et al., 2001;Ritzema, 1994;van Schilfgaarde, 1957).In the algorithm developed here, the two options for the boundary condition are implemented: with "lateral drainage" (Eq.8) and Vachaud's "vertical drainage" (f w = K s ).In some field conditions, a mixture of both end-time boundary conditions might be expected, requiring empirical weighing of the two conditions.

Calculation of singular time points
Following Mein and Larson (1973), time to ponding t p is the time for f p = i (intersection of the curves in Fig. 2a), typically when the surface water content is equal to saturation (Fig. 2c).At t = t p the equivalent wetting front depth (z p ) can be calculated by equating Eqs. ( 4) and (7): Since Eq. ( 9) is implicit in z p ≥ 0, it can be solved for each time step by defining the function G p : R → R and its derivative G p so that the root z p ∈[0, z w ] (i.e., G p (z p ) = 0) is the wetting front depth at t p , where z p can be obtained applying a bracketed Newton-Raphson algorithm (Press et al., 1992) obtaining the following: where k is the Newton-Raphson iteration level, and ε the error tolerance (here ε = 10 −8 ).From Eq. (3) at t = t p and z = z p , and F p = i • t p we obtain the following: Next, to ensure that F p (Eq. 3) and F = i • t p match at the intersection of the two curves on t = t p (Fig. 2b), an abscissa translation (shift time, t 0 ) is applied to F p (Mein and Larson, 1973).Setting z = z p on Eqs. ( 2) and (3) yields t 0 as follows: Lastly, t w is determined by calculating the integral Eq. ( 2) at z F = z w = L − h b (Fig. 1) and adjusting for t p and t o , and using Eq. ( 3), the cumulative infiltration at t w is determined by the following: The value of t w is equivalent to the nondimensional time X c proposed by Liu et al. (2011) that relies on the empirical error correction between RE solution and the Green-Ampt model.However, here t w (Eq.14) is calculated analytically for the more general case of nonuniform soil water content.The solution of Salvucci and Entekhabi (1995) can be simplified by setting the wetting front slope to zero (i.e., a horizontal front (α = 0) at the depth z F , Fig. 1).This approach reduces the solution, making it analogous to Eq. ( 2), which was employed by Bouwer (1969) in his explanation of the Green-Ampt model's applicability.For initial nonponding conditions, the equation becomes As the wetting front travels deeper into the soil, α could increase, contingent on the type of soil (e.g., α is larger for fine soils).However, as the wetting front approaches WT, the pore space available for infiltration is small, which limits the error of the calculations (Salvucci and Entekhabi, 1995).This assumption is tested in Sect.2.4.For a given time t, to solve for z F we specify the implicit function of Eq. ( 16) G: R → R and its derivative G , so that the root z F ∈[z i−1 , z w ] of the function G is equal to the new depth of the wetting front, In summary, for each time increment the proposed algorithm computes the depth of the wetting front (Eq.17), F (Eqs. 3, 15) and f (Eqs.4, 7 and 8) using the singular times auxiliary Eqs. ( 12)-( 14).A bracketing step in the Newton-Raphson algorithm is necessary, as the function G is undefined outside its physical range (z p < z F < z w ) (Press et al., 1992).The proposed algorithm is generic in that it can be used with any soil hydraulic functions like those of Gardner (1958), van Genuchten (1980) or Brooks and Corey (1964) (Appendix A) if numerical integration is used.Here, we used a Gaussian-quadrature integration scheme (Abramowitz and Stegun, 1972;Press et al., 1992).

Infiltration of soils with a WT and variable rainfall without initial ponding
For real VFS field situations, unsteady rainfall without initial soil ponding must be considered.Nonuniform rainfall is described by a hyetograph as a series of constant rainfall periods j (i.e., i = i j for t j < t < t j +1 ).The runoff produced by excess infiltration (i.e., Hortonian) and WT saturation (i.e., Dunne) are then determined at each time by water balance at the surface without accounting for evaporation during the rain event (Chu, 1997): where is the increment for that rainfall period, P and RO [L] are cumulative precipitation and runoff (excess rainfall), respectively, and when present s is the surface storage (0 < s < s max ) that acts as a reservoir that must be filled (s = s max ) before runoff is generated (Chu, 1978;Skaggs and Khaleel, 1982).For each period, if there is excess at the surface ( P − F > 0), the excess is first distributed to fill up the surface storage ( s ≤ s max −s) and the remainder (if any) to runoff ( RO = P − F − s ≥ 0).
For nonponding conditions at the beginning of the event, t p and t 0 must be calculated (Eqs.12-14), otherwise if initial ponding is present t p = t 0 = 0.If during a rainfall period the surface storage becomes zero, and if the new i j +1 of the following period is larger than the infiltration rate at the end of the last period, t p (and t 0 ) must be recomputed (Eqs.12-14) for the subsequent rainfall event (Chu, 1978;Skaggs and Khaleel, 1982).Also, each time t p , and t 0 are calculated, t w has to be recalculated.
To allow for predictions of soil water content redistribution during the event (Fig. 1) and to maintain mass balance during infiltration for alternating periods of ponding and nonponding conditions, it is necessary to track the "effective" position of the wetting front z F for periods with no ponding.To do this, the value of z F must satisfy the total cumulative infiltration amount at every time step (Fig. 1) and F (Eq. 3) becomes implicit in z F .As before, the root z F ∈[z j −i , z w ] (z j −1 is the wetting front depth at the previous time step) of the function G F : R → R and its derivative is as follows: The wetting front depth estimates provided by the algorithm are key in many hydrological applications where the aim is to simulate the potential for direct contamination of the WT by pollutants.
The next section provides an illustrative application of the full algorithm (herein referred to as SWINGO: Shallow Water table INfiltration alGOrithm) under unsteady rainfall conditions, typical in VFS settings (see the Supplement for coding details, source code, inputs, and outputs).Comparison of normalized infiltration rates (f/K s ) for soils without initial ponding described in Table 1, with vertical drainage (Vachaud) bottom boundary (f w ) conditions.Lines represent the SWINGO simplified model results.Symbols represents Richards equation numerical solution.The rainfall rate i was selected based on the saturated hydraulic conductivity (K s ) for each soil to ensure ponding at the surface.RE (Celia et al., 1990) using Nofziger and Wu's (2003) CHEMFLO-2000 model.We used four soils that represented a variety of attributes.The Brooks and Corey soil water attributes and hydraulic conductivity curves (Table 1) were used for the initial soil description, and this description was later compared with van Genuchten parameters yielding similar results (results not shown).The first three soils represent typical clay, silty loam, and sandy loam soils with a 1.50 m deep WT (Salvucci and Entekhabi, 1995).The fourth soil corresponds to a fine sandy soil experimentally studied by Vachaud and Thony (1971) with a WT at 1.01 m.
The soil water initial condition in CHEMFLO-2000 was set to hydrostatic equilibrium with a WT (Eq.1).The bottom boundary condition was set to a fixed matric potential h(z = L) = 0, to be representative of a WT at depth L. To simulate rainfall, the top boundary condition is set to a mixed type boundary with the flux density equal to the specified rainfall rate and the critical matric potential equalling zero (Nofziger and Wu, 2003).To allow for the development of distinct t p and t w values during the simulation, the constant rates of rainfall were chosen based on the soil texture.This selection was done utilizing a ratio of i/K s = 6 for the fine soils (clay and silty loam) and i/K s = 2 for the coarse soils, corresponding to the sandy loam and fine sandy soils studied by Vachaud and Thony (1971).
The comparison of the relative infiltration rates (f/K s ) calculated by RE (symbols) and the proposed SWINGO (lines) for the case of vertical drainage end boundary conditions (f w = K s ) is shown in Fig. 3.The performance of the algorithm is similar to RE for all soils studied.The median efficiency coefficients C eff (Nash and Sutcliffe, 1970) 1 with vertical drainage (Vachaud) bottom boundary (f w ) conditions.ranged from 0.927-0.9997,with the highest values being for clay and yielding statistically acceptable models at a 0.01 level of significance (Ritter and Muñoz-Carpena, 2013) (Table 1).For the same clay soil with ponded conditions and a WT, Salvucci and Entekhabi (1995) reported errors of approximately 5 % at time t w , at the point when the wetting front reaches the WT saturation (z w ), and the infiltration rate switches to the saturated hydraulic conductivity f w = K s (f/K s = 1).Smaller differences (1 % for clay and sandy loam and 3 % for the rest) were found between both solutions in our tests.These observations indicate that the simplification (horizontal wetting front, α = 0) did not affect the predictive ability of the infiltration rate.A crucial pattern to notice is that the estimates of time to ponding acquired across our tested soil types and normalized rates of rainfall closely matched the outputs of the RE solution.Our results also indicate that the use of the nonuniform integral equations (Eqs.9-12) effectively limit errors in the t p estimation that sometimes occur when utilizing the Green-Ampt model (Barry et al., 1996).
Figure 4 displays the cumulative infiltration and the depth of the wetting front determined using Eqs.( 20)-( 21) for the vertical drainage boundary condition for the cases from Table 1.Similar to the infiltration curves, z F values exhibited a plateau at t w as they reach column saturation (Fig. 4b), corresponding to the capillary fringe at a depth of z F = z w = L−h b (Fig. 1), and therefore are not equal to the depth of the WT (fine sand: L = 1.01 m; other soil types: L = 1.50 m).
As the simplified approach is able to produce reliable z F predictions, it also allows for the depiction of the redistribution of the soil water content during infiltration.We display the predictions of soil water (Fig. 5) calculated by the proposed algorithm (dashed lines) compared with the outputs of the RE solution (solid lines) for the nonponding numerical test examples used previously.The simplified model is able to identify the midpoint of the wetting front depth at all time points.Additionally, our simplification of including the horizontal wetting front (α = 0) generates an accurate prediction of soil water at earlier time points for all soil types, but this prediction decays somewhat at later time points when approaching column saturation for fine soils.The model does not degrade at later time points for the sandy soil type when it matches a horizontal wetting front redistribution.As mentioned previously, because of the smaller pore space near column saturation, the mass errors generated by nonzero slopes stay negligible.The infiltration mass balance error at the end of the simulation (Fig. 4a) ranges from 3 to 8 %.This range of error values is deemed satisfactory, as these errors are the summation of approximation errors of both the infiltration and redistribution of soil moisture generated during the simulation.
A quick comparison of execution times between CHEM-FLO and SWINGO for the cases in Fig. 5 yielded small reduction of 1-5 s with SWINGO (CPU: 1.6 GHz Intel Core 2 Duo).These results are machine, computer, and compiler dependent, where a CHEMFLO finite-difference solution is implemented in Java computer language (run in Oracle ® jre-8u144) and contains a graphical user interface not intended for optimized simulations.SWINGO was implemented as a command line application in Fortran (Intel ® Fortran Compiler v17.0.4).Admittedly, the differences will likely be smaller with optimized code and new developments of Richards equation implementations (e.g., de Maet et al., 2014).However, these small differences will likely be compounded in the context of throughput simulations where the algorithm will be used in some applications.For example, the model coupled in the companion paper (VFSMOD) is used in current long-term pesticide regulatory assessments (30-year daily time steps in the USA or 10-year daily time steps in the EU) (Muñoz-Carpena et al., 2010, 2015).Considering ∼ 1/3 to ∼ 2/3 of days with rainfall runoff, the model would be run between ∼ 3000 and ∼ 7000 times for a typical 30 years.assessment.Under this type of throughput applica-  1.
tions condition, even a marginal time improvement can prove advantageous.In addition to marginal speed benefits, the proposed algorithm is robust (from the direct integral solution) and has physical consistency with the original VFSMOD that uses the extension of Green-Ampt for unsteady rainfall conditions (Chu, 1978;Skaggs and Khaleel, 1982) without the presence of a shallow water table.

Experimental testing
The physics of the model were tested in a second step using experimental data from Vachaud et al. (1974) and Chu (1997).The data collected in the laboratory represents infiltration under ponded conditions in a vertical column of fine sand soil with a WT at 0.925 m depth.To demonstrate the generality of the proposed algorithm, the Vachaud et al. (1974) measured soil hydraulic characteristics were fitted to van Genuchten soil water characteristic and related unsaturated hydraulic conductivity function based on Mualem (1976) simplification (vG : vG), and the latter was  also fitted to Gardner function (vG : Grd) (see Appendix A and soil parameters in Table 1).The goodness of fit of these hydraulic functions (inset of Fig. 6) shows a small improvement of the K(h) function for Gardner over that of van Genuchten-Mualem against the experimental data.
The simulated relative infiltration rates obtained with the proposed algorithm matched the observed data well (C eff = 0.913-0.942,RMSE for f = 5.07 × 10 −6 -6.20 × 10 −6 m s −1 ), yielding statistically acceptable models at α = 1 % for vG : Grd and α = 5 % for vG : vG combinations (Ritter and Muñoz-Carpena, 2013) (Table 1).The main differences observed between SWINGO solutions with vG : Grd or vG : Grd soil water functions are near the time when the wetting front depth approaches the WT, with a small advance (∼ 0.02 h) introduced by the vG : Grd option.These small differences are related to the slope of the wetting front being different than 0, especially close to the intersection with the WT at the end of the event (Fig. 5).Note also that in this experimental case no observed data were available for comparison at the time when the wetting front reached the WT.
In all, these results provide not only a test of the simplified model against experimental data, but also illustrate its robustness and flexibility to handle other soil hydraulic functions.

Illustration for unsteady rainfall conditions
The use of SWINGO to simulate realistic unsteady rainfall conditions is presented for a storm composed of four rainfall periods: , and i 4 = 0.25 cm (5 < t ≤ 6.9 h) (Table 2 and Fig. 7).The soil is clay (Table 1) with bottom vertical drainage boundary  beginning of the event the soil is not ponded and is in equilibrium with the WT at L = 150 cm below the surface, so max(z F ) = z w = L − h b = 0.6 m (Fig. 7).For the initial period, we calculate first the time to ponding with Eqs. ( 9)-( 12) and ( 19) (t p = 4657.2s = 1.29 h), the corresponding t 0 (2319 s = 0.64 h) with Eq. ( 13), and the time to reach the WT t w (16 100 s = 4.47 h) with Eq. ( 14).Since the t w is higher than the rainfall period and t p lower than the rainfall period, infiltration is equal to the rainfall rate (f = i 1 ; 0 < t < t p1 ) before ponding.After ponding it follows the infiltration capacity curve described by the solution of Eqs. ( 16)-( 17).
At the beginning of the second rain period, since the new rainfall rate is less than the infiltration rate at the end of the previous period (i 2 = 0.25 cm h −1 < f p = 0.52 cm h −1 ) and t w is still beyond the period, the infiltration rate equals the new rainfall rate (f = i 2 ).At the beginning of the third period, the new rainfall rate is larger than the corresponding potential infiltration rate at that time (i 3 =1 cm h −1 > f p = 0.44 cm h −1 ) and ponding starts again immediately such that the new t p = t 3 (15 000 s = 4.2 h, beginning of the new rainfall period), and t 0 (13 764 s = 3.82 h) and t w (18 500 s = 5.14 h) are recalculated.Since t w is beyond the period, the infiltration is maintained at capacity for the duration of this rainfall period.For the last period, the rainfall rate is lower than the ending infiltration capacity for the last period (i 4 = 0.25 cm h −1 < f p = 0.34 cm h −1 ), and infiltration is initially set to the rainfall rate.However, since t w is within this period, the soil saturates when the water front reaches the WT depth (t ≥ t w ), and this results in saturated vertical drainage flow with unit hydraulic gradient f = f w = K s (Eq.7) until the end of the storm.The values of the wetting front position (z F ) in Table 2 are calculated from the solution of Eq. ( 17) during infiltration capacity (ponding) periods, and the equivalent depths described by Eq. ( 19) during nonponding periods.Similarly, cumulative totals are calculated with Eqs. ( 3) or ( 15), and excess rainfall amounts are calculated with the surface mass balance Eqs. ( 18) for every time step.

Evaluation of WT effects on infiltration under conditions of ponding and nonponding
Figure 8 presents the effect of the WT depth variation (L = 0-200 cm) and event duration (0.5 < D < 6.0 h) on cumulative infiltration under ponding conditions for the soils in Table 1.The two end-time boundary conditions are compared: f w vertical (a-d) and f w lateral (e-h) with S o = 0.02, b = 1 m, and K sh = K s (Eq.8 and Table 1).For the conditions tested it is possible to identify three clearly defined regions (denoted I, II, and III in Fig. 8) based on the influence of the WT depth on the cumulative infiltration.Region I (left, shaded in Fig. 8) represents the WT near the surface, i.e., when it is within the capillary fringe area L < h b (Fig. 1).The position of the WT in this region does not affect infiltration since the soil column is already saturated regardless of L with F = D • f w .Next, Region II (clear background in Fig. 8a-d) is the most sensitive to variations of WT depth, located between L = h b and a limit depth (L = 125-180 cm) where the variation of F is small (slope less than 0.2 %).This limit depends on the shape of the soil water characteristic curve for each soil.Finally, Region III represents a region where surface infiltration can be considered effectively decoupled from the presence of the WT.Next, the robustness and physical behavior of the algorithm under nonponded initial conditions was tested with  1).
different rainfall rates (i = 0.1-20 cm h −1 ), event durations (D = 1-12 h), and WT depths (L = 0-400 cm). Figure 9a-d summarizes the results for D = 6 h and the vertical drainage boundary condition (f w = K s ).Two main effects are identified.Firstly, as expected F is insensitive to changes in L for rainfall intensities lower than K s , when f = i (no ponding) and F = D • K s .Notice that this effect, although present, is not visible in the clay soil (Fig. 9b) since the K s is below the first isoline.Secondly, for rainfall rates above K s , the sensitivity to L varies by soil, depending on h b and the time to ponding values for each rainfall rate (Eq.12).As in the ponding case, the soil column is saturated when L ≤ h b , and there is no sensitivity below this depth.In finer, less permeable soils (Fig. 9a-b) ponding happens earlier for the same rainfall rate i, resulting in an increased sensitivity to L with lower rainfall rates.For the lateral drainage boundary condition, results are similar for the finer soils (Fig. 9e-f), but much more sensitive to WT depth and rainfall rate values for more permeable soils (Fig. 9g-h).Importantly, since excess rainfall runoff is complementary to F (Eq. 18), these results also quantify the important influence that the combined effects of WT, soil type, and rainfall intensity can have on surface runoff flow and transport processes in the VFS.
Figure 9. Change in cumulative infiltration (F ) as a function of rainfall rate (i) and water table depth (L) under nonponded initial conditions after a rainfall duration D = 6 h for the four types of soils in Table 1 and two end drainage bottom boundary conditions (f w ): (a-d) vertical and (e-h) lateral (S o = 0.02, b = 1 m, and here K sh = K s from Table 1).The isolines describe the change of F with water table depth for the same rainfall rate (i); h b is the bubbling pressure (capillary fringe) for each soil type (Table 1).

Summary and conclusions
Limitations in current modeling approaches hamper the evaluation of the effects of WTs on soil infiltration and runoff.A promising way to overcome these issues is by utilizing simplified yet realistic specialized algorithms in conjunction with available hydrological models to evaluate the impact of WTs in the environment.Previously, Salvucci and Entekhabi (1995) and Chu (1997) recommended the use of Green-Ampt implicit integral equations to examine infiltration into ponded soils with WT.We developed and assessed a simplified generic algorithm that is appropriate for coupling with available hydrological models, in particular the study of WT effects on VFS runoff pollution control performance.The proposed SWINGO algorithm is generic -it can utilize any configuration of soil hydraulic functions -and can be operated under nonponded, ponded, and realistic variable rainfall conditions to determine runoff (excess rainfall), infiltration, and soil water redistribution during the event.
SWINGO performed well (C eff from 0.91 to 0.99) in comparison with the RE solution and using experimental data on five representative soils.The algorithm also was able to successfully describe the soil water redistribution during the simulated event.These useful and reliable predictions indicate that the proposed approach incorporating a horizontal slope of the wetting front is suitable for most real-world applications.Through an application of our proposed SWINGO algorithm, we showed the sensitivity of the infiltration and excess runoff to the depth of the WT, the length and intensity of the rainfall event, the soil type, and drainage bottom conditions.
Some of the limitations of the proposed algorithm are the assumptions of a homogeneous soil profile and horizontal wetting front for fine soils.Future research is recommended to determine the general validity of the assumption of a hydrostatic equilibrium and the proposed computation of singular points during the infiltration episode.Additional experimental testing of the model should be conducted using data collected under various controlled and natural conditions (especially during events long enough for the wetting front to reach the WT).As in real soils a mixture of both end-time lateral and vertical boundary conditions might be expected, these field studies could also investigate site and event characteristics for which these boundary conditions might be dominant or have relative influence.
As SWINGO was accurate, fast and robust when analyzing a variety of conditions, it is appropriate to couple with currently available hydrological models to gauge the influence of the presence of WTs on other landscape processes beyond the simulation of filter strips.The proposed integral equation has broader relevance as a step forward in improving the science of hydrologic modeling in general in many other settings, e.g., to study shallow water table effects on surface runoff, infiltration, flooding, transport, ecological, and land use processes.
The dynamic coupling with VFS overland flow and sediment and pesticide transport processes is developed in the companion paper (Lauvernet and Muñoz-Carpena, 2018).Global sensitivity and uncertainty analysis of the coupled model is conducted to identify important input factors and their interactions that will provide better understanding of the fundamental processes controlling VFS efficiency under WT conditions and guide users to select effective parameters for practical applications.
Data availability.Model, code, data and additional materials are publicly accessible through the supplement provided in this paper.

3
Figure3.Comparison of normalized infiltration rates (f/K s ) for soils without initial ponding described in Table1, with vertical drainage(Vachaud)  bottom boundary (f w ) conditions.Lines represent the SWINGO simplified model results.Symbols represents Richards equation numerical solution.The rainfall rate i was selected based on the saturated hydraulic conductivity (K s ) for each soil to ensure ponding at the surface.

Figure 4 .
Figure 4. Comparison of (a) cumulative infiltration (F ) and (b) wetting front depth (z F ) movement results.Lines represent the SWINGO simplified model and symbols represent Richards equation numerical solution for soils without initial ponding in Table1with vertical drainage(Vachaud)  bottom boundary (f w ) conditions.

Figure 5 .
Figure 5.Comparison of soil water (θ) redistribution between Richards equation numerical solution (solid lines) and the SWINGO simplified model (dashed lines) during infiltration without initial ponding and with vertical drainage (Vachaud) bottom boundary condition (f w ) for soils in Table1.

Figure 6 .
Figure 6.Comparison of the simplified and RE results against Vachaud et al. (1974) experimental data set (figure body), and fitting of soil water characteristics to different equations (inset); vG and Grd represent, respectively, the van Genuchten and Gardner's soil characteristic curves used to parametrize the simplified and RE models (seeTable 1 for details).

Figure 7 .
Figure 7. Simulation of an unsteady rainfall event on the clay soil in initial equilibrium with a shallow water table at 150 cm depth, nonponded conditions, and vertical drainage (Vachaud) bottom boundary conditions (f w ): (a) infiltration (f ) and rainfall rates (i, subindices in i 1 to i 4 represent the rainfall periods within the hyetograph); (b) cumulative rainfall (P ), infiltration (F ), excess runoff (RO), and wetting front depth (z F ) during the event.

Figure 8 .
Figure8.Effect of water table depth (L) on cumulative infiltration (F , represented by isolines) for distinct soils under initial ponding and different durations of infiltration events (D) for four types of soils and two end drainage bottom boundary conditions (f w ): (a-d) vertical and (e-h) lateral (S o = 0.02, b = 1 m, and here K sh = K s from Table1).
Figure 1.Conceptual depiction of infiltration and soil water redistribution for soils with shallow water table for (a) time before wetting front reaches the water table and (b) time after the wetting front reaches the water table (t ≥ t w ), where surface infiltration flow (Q f ) is limited by lateral Boussinesq subsurface flow (Q L ).See explanation of symbols in the text.
Conceptual curves of (a) infiltration rate, f ; (b) cumulative infiltration, F ; and (c) soil water redistribution, θ , under shallow water table, for soil without initial ponding, and constant rainfall rate (i) conditions.The singular times for ponding (t p ), shifting (t 0 ) and to reach column saturation (t w ), and final infiltration rate (f w ) after the wetting front reaches the water table (t ≥ t w ) are represented.
R. Muñoz-Carpena et al.: Shallow water table effects on water, sediment, and pesticide transport 2.3 Infiltration capacity algorithm after surface ponding

Table 2 .
Infiltration and excess runoff calculations for an illustrative unsteady rainfall event on a clay soil with no initial ponding at equilibrium with a shallow water table at 150 cm depth, s max = 0, and end vertical boundary conditions.The "+" sign in the first column represents any time right after the beginning of that time step.