Khalil Alastal

Rachid Ababou

Moving Multi-Front (MMF): A generalized Green-Ampt approach for vertical unsaturated flows

A Moving Multi-Front (MMF) method is developed and tested for solving the Richards equation governing unsaturated flow in vertical homogeneous porous columns. The MMF model is a gridless method. It can be viewed as a generalization of the Green-Ampt piston flow approach, which models flow as a single abrupt moving front separating saturated and dry regions during infiltration. The MMF model further generalizes this concept, using a parametrization of unsaturated water retention and conductivity-pressure curves. It reduces the Richards PDE to an ODE system governing "M" moving front positions. Different tests are developed to validate this approach for 1D transient downwards/upwards flows submitted to constant and time-varying pressure boundary conditions. They include: (i) infiltration to deep water tables; (ii) infiltration to shallow water tables; (iii) capillary rise from fixed water tables; (iv) gradual water table rise (partially saturated column with evolving pressure condition at bottom). The MMF results are compared favorably to finely discretized fixed grid solutions of Richards PDE. Analyses of error and accuracy show satisfactory results in terms of water content profiles and boundary fluxes (e.g. infiltration rate).

Introduction

Vertical water flows through porous media like soils are of great interest for hydrologists and geotechnical scientists and engineers. In general, these water flows take place in the presence of air, and they are modeled based on a generalized Darcy approach (Darcy-Buckingham) where the water content, as well as the pore pressure, are taken into account in the unsaturated porous medium. Unsaturated flows include important phenomena such as downward infiltration, and capillary rise due for instance to a rise of the water table. Infiltration is a downward water flow process through which water penetrates the soil surface and flows towards drier soil due to capillary gradients and gravitational force. Capillary rise involves a rise of moisture driven upwards by capillary gradients and partially counteracted by the downwards gravity force.

Researchers have developed a number of analytical, semi-analytical, and numerical models or "solvers" to simulate infiltration and, more generally, unsaturated flow processes. The Darcy-based Richards equation [START_REF] Richards | Capillary conduction of liquids through porous mediums[END_REF] is still the most frequently used equation to describe unsaturated flow, and it can be generalized to deal with partially saturated as well as unsaturated porous media, and to include other effects like compressibility. The Richards equation is, mathematically, a quasi-linear PDE (Partial Differential Equation) with nonlinear coefficients. Numerical solvers have been developed for solving the partially saturated/unsaturated Richards PDE with implicit discretization in time, Finite Volume (Bigflow 3D) [START_REF] Ababou | BIGFLOW: A Numerical Code for Simulating Flow in Variably Saturated, Heterogeneous Geologic Media (Theory and User[END_REF] or Galerkin Finite Element (Hydrus 1D/2D/3D) [START_REF] Simunek | The HYDRUS-1D Software Package for Simulating the One-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media[END_REF] discretization in space, and iterative methods to linearize and solve the resulting nonlinear system of algebraic equations (modified Picard iterations in Bigflow 3D; Newton-Raphson iterations in Hydrus 1D/2D/3D). These two codes will be used here for benchmarking.

The highly nonlinear nature of the Richards PDE has attracted the attention of many researchers with the aim of improving or adapting nonlinear iteration methods, along with time-stepping methods. See for instance [START_REF] Islam | Numerical tests of the lookup table method in solving Richards' equation for infiltration and drainage in heterogeneous soils[END_REF] and references therein. Specifically, in [START_REF] Islam | Numerical tests of the lookup table method in solving Richards' equation for infiltration and drainage in heterogeneous soils[END_REF], the authors propose a "lookup table method" to enhance the performance of Picard and Newton iterations for highly nonlinear constitutive relationships by sampling adaptively the pressure (h) at points where the capacity curve C(h) varies most rapidly, in order to avoid divergence of the Picard or Newton solvers.

On the other hand, due to the complexity and non-linear nature of the Richards PDE, analytical and semi-analytical solutions to infiltration-type problems were developed under special constraints and limitations, usually for 1D vertical flow. In spite of their limitations, the analytical models are useful in several ways: for approximate analyses of real infiltration processes, but also, for evaluating the accuracy of

K s s F F
The GA model has been the focus of many interests because of its simplicity and adaptability for a variety of applications in soil and watershed hydrology. Other works in the literature extended or analyzed the GA approximation in various ways, e.g.: infiltration under time varying conditions [START_REF] Chu | Infiltration during an unsteady rain[END_REF][START_REF] Warrick | Infiltration under variable[END_REF]; layered and heterogeneous soils [START_REF] Liu | Green-Ampt model for layered soils with nonuniform initial water content under unsteady infiltration[END_REF][START_REF] Selker | Green and Ampt infiltration into soils of variable pore size with depth[END_REF][START_REF] Yanwei | Improvement and verification of the Green-Ampt model for sand-layered soil[END_REF] and [START_REF] Kacimov | Green-Ampt one-dimensional infiltration from a ponded surface into a heterogeneous soil[END_REF]; sloping soil surfaces [START_REF] Chen | Green-Ampt infiltration model for sloping surfaces[END_REF] and [START_REF] Gavin | A simple method to analyze infiltration into unsaturated soil slopes[END_REF]; air flow and air compression effects on infiltration [START_REF] Morel-Seytoux | Two-phase flows in porous media[END_REF]; non uniform initial moisture [START_REF] Liu | Green-Ampt model for layered soils with nonuniform initial water content under unsteady infiltration[END_REF]; moisture redistribution after infiltration [START_REF] Musy | Physique du sol, 6. Presses Polytechniques et Universitaires Romandes[END_REF]; horizontal imbibition instead of vertical infiltration [START_REF] Cao | Feasibility investigation of improving the modified green-ampt model for treatment of horizontal infiltration in soil[END_REF][START_REF] Prevedello | Generalization of the Green-Ampt theory for horizontal infiltration into homogeneous soils[END_REF][START_REF] Zhang | Modification and evaluation of Green-Ampt model: Dynamic capillary pressure and broken-line wetting profile[END_REF]. The latter work by [START_REF] Zhang | Modification and evaluation of Green-Ampt model: Dynamic capillary pressure and broken-line wetting profile[END_REF] also proposes a parametric extension of GA based on a piecewise linear approximation of moisture profiles with two linear zones (for horizontal flow). Similarly, [START_REF] Meng | Infiltration Simulation with Improved Green-Ampt Model Coupled with the Wet Zone Partition Function[END_REF] present an extended GA model for vertical flow comprising a partition between a saturated zone and an unsaturated zone with assumed linear pressure profile in the unsaturated zone (the wet/unsaturated partition function is fitted empirically using a neural network). In addition, the paper by [START_REF] Ali | Green-Ampt approximations: A comprehensive analysis[END_REF] presents an extensive review and performance analyses of various "explicit" approximations of the classical GA solution, in comparison with the exact "implicit" GA equation.

At this point, in comparison with the GA piston flow model, it is useful to consider a technical mathematical issue concerning the spacetime discretization methods used for solving numerically the nonlinear unsaturated flow equations (Richards PDE). For example, in the numerical based method like finite elements (Hydrus), or finite volumes (Bigflow), the nonlinear Richards PDE (Partial Differential Equation) is discretized in both space and time. In contrast, the GA equation avoids discretizing space (z) by assuming that the vertical moisture/pressure profile is a step function with only one sharp wetting front located at position Z(t), at any given time t. [START_REF] Warrick | Infiltration under variable[END_REF] used a slightly modified form of the GA infiltration equation under a boundary condition of time-variable ponded depth at soil surface. They proposed to discretize the GA model in time with relatively coarse time steps, and to reformulate it as a classical GA model under piecewise constant boundary conditions over time. [START_REF] Lee | A MATLAB method of lines template for transport equations[END_REF] applied the Method of Lines (MoL) for unsaturated flow problems governed by Richards' PDE. This method is based on a fixed discretization of the spatial domain while keeping the time domain continuous. This results in a system of nonlinear Ordinary Differential Equations (ODE) that can be easily solved by available numerical time integrator schemes. [START_REF] Matthews | Modeling flow through a one-dimensional multi-layered soil profile using the Method of Lines[END_REF] further considered the MoL for multi-layered soils and proposed methods to handle the discontinuity of water content and hydraulic conductivity at interfaces between soil layers. [START_REF] Talbot | A method for computing infiltration and redistribution in as follows: a discretized moisture content domain[END_REF] proposed a solution method for 1D vertical infiltration and redistribution based on the GA approximation, where they discretize the water content domain (θ) into equal bins Δθ instead of discretizing the vertical spatial dimension. This is equivalent to discretizing the water content profiles θ(z) into vertical bins of width Δθ. The vertical water movement is decomposed at every time step into a two-step process: the so-called "infiltration step" is followed by a "redistribution step", the latter being based on capillary forces. [START_REF] Ogden | A new general 1-D vadose zone flow solution method[END_REF] subsequently built on the finite water-content discretization of [START_REF] Talbot | A method for computing infiltration and redistribution in as follows: a discretized moisture content domain[END_REF], however this time they used the method of lines (MoL), and they approximate the partial derivative q ( / ) by neglecting a term that contributes to capillary suction gra- dient [ [START_REF] Ogden | A new general 1-D vadose zone flow solution method[END_REF], their Eq. ( 10)]. They also criticized the previous paper by [START_REF] Talbot | A method for computing infiltration and redistribution in as follows: a discretized moisture content domain[END_REF], stating that the redistribution sub-step used in infiltration simulation was incorrect, and they renamed and reformulated this sub-step as "capillary relaxation", which moves water from low to high capillary bin(s) or region(s).

It has come to our attention, while revising this paper, that two particle-based schemes for solving the Richards equation were tested in [START_REF] Beaudoin | Une méthode particulaire pour résoudre l'équation de Richards (A particle method for solving Richards equation)[END_REF]. Both schemes are Lagrangian (the MMF method can also be qualified as Lagrangian), but each of them is different from MMF: the first scheme is a "Particle Strength Exchange" or "Weighted Particle Method" (where each particle has variable mass m(t)); and the second scheme is based on the concept of "Diffusion Velocity", or rather "pseudo-velocity" (it is based on a transformation from nonlinear capillary diffusion to strongly nonlinear convection). [START_REF] Haq | Meshless method of lines for the numerical solution of generalized Kuramoto-Sivashinsky equation[END_REF] present a "meshless Method of Lines" for solving a 1D nonlinear fourth-order PDE known as the Kuramoto-Sivashinsky equation. Their method could potentially be adapted to solve the nonlinear second-order Richards PDE. In spite of its name, the method is based on a fixed grid of nodes x ( ) i , and then, on an approximate linear decomposition of the unknown u x t ( , ) on a basis of functions

x { ( )} j with time-dependent coefficients t { ( )} j = j N ( 1, , ). The unknowns u x t ( , )
i are nodal, although the unknown coefficients t { ( )} j themselves are indeed meshless, and these coefficients are finally used to solve a system of ODE's for the nodal unknowns u t ( )

i = u x t ( , )
i by the MoL approach evoked earlier.

The objective of this paper is to propose a new quasi-analytical efficient alternative model named Moving Multi-Fronts (MMF). We will focus here on a one-dimensional (1D) vertical system taking into account both saturated and unsaturated flow dynamics in the porous domain. We will propose to improve on the classical G-A approach by discretizing the state variables themselves (the pressure head and consequently the water contents and the hydraulic conductivity along the soil profile) such that the nonlinear behavior can easily be accounted for. The resulting nonlinear ODEs can then be solved either by numerical integration or, say, by explicit Runge-Kutta finite differences in time. 1 The MMF tracks movable multi-fronts instead of a single GA 1 Note: in all such cases, when the resulting equations boil down to single numerical solutions (validation, benchmarking).

The literature involves many analytical and quasi-analytical models that have been developed over the past years to tackle infiltration problem (Green and A mpt, 1911;[START_REF] Horton | The Role of infiltration in the hydrologic cycle[END_REF][START_REF] Philip | The theory of infiltration: 1. the infiltration equation and its solution[END_REF][START_REF] Morel-Seytoux | Derivation of an equation of infiltration[END_REF][START_REF] Smith | A parameter-efficient hydrologic infiltration model[END_REF][START_REF] Warrick | An analytical solution to Richards' equation ponding depths of water[END_REF][START_REF] Ghotbi | Infiltration in unsaturated soils -An analytical approach[END_REF][START_REF] Parlange | Analytical approximation to the solutions of Richards' equation with applications to infiltration, ponding, and time compression approximation[END_REF][START_REF] Srivastava | Analytical solutions for one-dimensional, transient infiltration toward the water table in homogeneous and layered soils[END_REF][START_REF] Su | Approximate solution of a one-dimensional soil water infiltration equation based on the Brooks-Corey model[END_REF][START_REF] Triadis | Analytical model of infiltration under constant-concentration boundary conditions[END_REF][START_REF] Wang | An Analytical Solution for One-Dimensional Water Infiltration and Redistribution in Unsaturated Soil[END_REF][START_REF] Wu | Analytical solution to 1D coupled water infiltration and deformation in unsaturated soils[END_REF]; among many others. [START_REF] Morbidelli | Rainfall infiltration modeling: a review[END_REF] present a review on infiltration modeling, accommodating constant or variable rainfall rate, ponding, and post-ponding infiltration, for homogeneous or 2-layered soils. They review several such models, some of them related to (or extensions of) the classical Green-Ampt approach (see references therein). The review paper by [START_REF] Morbidelli | Rainfall infiltration modeling: a review[END_REF] also points out other infiltration models, including the semi-analytical model of [START_REF] Corradini | Modeling local infiltration for a two-layered soil under complex rainfall patterns[END_REF] which models infiltration/redistribution in a 2-layer soil. Their model can be considered as a simplified conceptual model of the Richards equation. It yields a system of ODE's governing suctions vs. time at interface and/or boundary points. When the two layers are identical, their model reduces to the Green-Ampt infiltration model under zero ponding depth.

The Green-Ampt (GA) piston flow model [START_REF] Green | Studies on soil phyics[END_REF]) is one of the most simplified analytical and physically based models that was developed initially for water infiltration under ponded conditions into initially dry homogeneous soils, far above the water table. The piston flow assumption in the GA model approximates the soil moisture distribution with an abrupt moving front separating a wetted (saturated) region above the front, and a dry region below it. This model implies a constant hydraulic conductivity ( ) and water content ( ) within the wetted zone, and it assumes that the abrupt wetting front can be characterized by a constant suction driving the flow [START_REF] Bouwer | Unsaturated flow in ground-water hydraulics[END_REF][START_REF] Bouwer | Rapid field measurement of air entry value and hydraulic conductivity of soil as significant parameters in flow system analysis[END_REF][START_REF] Bouwer | Infiltration of water into nonuniform soil[END_REF]. The model assumes that the front suction is a parameter that characterizes the soil, and does not depend on other factors such as initial water content and top boundary condition.

• There exists a well-defined "wetting front" separating the fully sa- turated and the totally dry regions, as shown in Fig. 1, where the dashed line represents the GA approximation.

• The wetting front is assumed to be characterized by some effective, constant suction head F , or pressure head = h F F (considered to be a characteristic parameter of the soil).

The GA model approximation divides the soil column into two zones, each zone having a priori its own hydraulic conductivity and water content (K , ):

• Top zone: this zone is totally saturated according GA approximation; it extends from the ground surface ( = z 0) to the moving wetting front

[ = z Z t ( ) F ]; in this zone we have = K K S , = S .
• Bottom zone: this zone is the totally dry region extending from the moving wetting front [ = z Z t ( )

F
] up to infinite depth; it is a semiinfinite zone; the flux is null in this zone since the initial soil is totally dry: = K 0, = d , where d is also denoted r ("residual water content").

Let us now develop the governing equations under these assumptions, using Darcy's flux-gradient law and mass conservation principles.

Global Darcy equation and GA infiltration rate

First, one develops a global approximation of Darcy's law = q z t K h z g ( , )

[( / ) ] for the entire saturated region. This is ob- tained by integrating q z t ( , ) vertically, at any fixed time t, through the saturated region, i.e., from the ponded surface at = z 0 to the wetting front at depth = z Z t ( )

F

. This procedure necessarily assumes the existence of an "effective" front pressure h F or front suction F (para- meter). In terms of pressure head "h", and with "z" downwards, the global version of the pressure gradient

h z / is = h z h H Z t / ( ) / () . F p F
Also, since the "transmission zone" is saturated in the GA model, we have = K K S in that zone. This yields finally the GA relation for the global integrated flux q t ( ) or infiltration rate f t ( ), in terms of pressure heads:

= = f t q t K h H Z t g ( ) ( ) ( ) S F p F (1a)
or equivalently, in terms of suction heads ( = = h and h )

F F : = = + + f t q t K H Z t g ( ) ( ) ( ) S F p F (1b)
Note that H p represents the positive pressure condition due to ponded water depth at soil surface (H 0 p

). The other variables and parameters are defined as follows:

• f t m s ( )[ / ]: is the infiltration rate at a given time t;

• K m s [ / ]
S : is the saturated hydraulic conductivity of the saturated transition zone;

• Z t m ( )[ ] F
is the depth of the moving wetting front at time t (it is also the vertical extent of the saturated zone);

• m [ ] F is the capillary suction head at the wetting front;

• g [dimensionless] represents normalized gravity = g g g /

Z

, with = + g 1 if the Oz axis is directed downwards (as is the case here).

Note that the global Darcy flux density q t ( ) (obtained by integration through the saturated region) is assumed to represent also the infiltration rate at soil surface, denoted f t ( ) [m/s]. The above relations for the single front GA infiltration rate can also be generalized in terms of the gravity vector g as follows, simply by changing the Oz axis direction:

• Let = + g
1 in the present case, for water infiltration, with z directed downwards.

• Let = g 1, with z directed upwards, e.g. to treat capillary rise or water table rise.

• Let = g 0, e.g. for horizontal flow model and/or in the absence of gravity.

Global mass balance and closure of GA infiltration equations

The previous relation of Eq. ( 1) is not closed because it contains two unknowns, the infiltration rate f t ( ) and the wetting front depth Z t ( )

F
. The closure of this problem is obtained by applying a global mass balance principle at the scale of the transmission zone. Thus, neglecting the flux in the initially dry zone just below the wetting front

( = q K 0 d d
), the depth of the wetting front Z t [ ( )]

F can be related to the cumulative infiltration F t [ ( )] by the following mass conservation principle:

= F t Z t ( ) ( ) ( ) S d F (2)
Taking the derivative of the cumulative infiltration F t ( ) with re- spect to time yields:

= = f t dF dt dZ dt ( ) ( ) S d F
(3) (footnote continued) integrals or to integrable 1rst order ODE's, we refer to these as quasi-analytical or semi-analytical solutions.

moving front, in similar Lagrangian manner.

In comparison with the "meshless MoL" of [START_REF] Haq | Meshless method of lines for the numerical solution of generalized Kuramoto-Sivashinsky equation[END_REF],the MMF method proposed here is completely meshless. On the other hand, in comparison to the method of "discretized moisture content domain" by [START_REF] Talbot | A method for computing infiltration and redistribution in as follows: a discretized moisture content domain[END_REF] and [START_REF] Ogden | A new general 1-D vadose zone flow solution method[END_REF], the MMF method to be presented here does not require the two steps process used in their model: all processes in the present MMF model (capillary gradients and gravity forces) are treated at once in a fully coupled manner. Further-more, their model also neglects some contributions of the non-grav-itational part containing [D( )∂θ/∂z; D is the soil-water diffusivity], whereas the MMF model does not. Finally, the present MMF can handle as easily co-gravity and anti-gravity flows. Examples including down-wards infiltration, capillary rise, sudden and gradual rise or drop of water table, will be presented to demonstrate the capabilities of the proposed model. Uniform and a nonuniform initial soil water profiles, constant or time dependent pressure head boundary conditions, will also be considered.

Method: the Moving Multi-Front model

The classical Green-Ampt (GA) model: single moving front

The classical GA model is a simplified representation of vertical water infiltration in a dry soil, away from the water table (deep water table). The GA model assumes a saturated piston type flow into the dry soil, i.e., flow is modeled as the displacement of a single sharp wetting front into a uniformly dry homogeneous soil. The front abruptly sepa-rates two distinct regions, a fully saturated region upstream (above), and a very dry region downstream (below). The sharp wetting front is pulled downwards by gravity and capillary suction forces, as shown in Fig. 1. This model is also called "single front" model because, as in the classical GA approach, it is based on the movement of a single front F Zt() .

Accordingly, based on the usual Green-Ampt piston flow hy-potheses, we assume that:

where:

• F t m ( )[ ] is the cumulative infiltration at time t;
• S is the saturated water content (or porosity); • d is the "dry" initial water content;

If the soil is initially very dry, then we have ideally d r (the initially dry water content is the residual water content of the soil). In that case the initial conductivity of the soil vanishes (

= K K ( ) 0 d d
), and for this reason, the flux q d below the wetting front is negligible or null at all times. Note also that Eq. ( 3) can also be obtained by applying mass conservation around the moving front [Z t ( ) F ] as follows:

= q t q t dZ dt ( ) ( ) ( ) top bot top bot F (4a) where q t f t ( ) ( ) top , q t ( ) 0 bot , = top S , and = bot d r . The ap- proximate equalities become exact if = d r exactly, implying = q t ( ) 0 bot , = q t f t ( ) ( ) top , whence: = f t dZ dt ( ) ( ) S d F (4b)
Now, finally combining Eq. ( 1) and Eq. ( 3) leads to the governing equation of the GA model:

= + + dZ dt K H Z t ( ) () 1 F S S d F p F (5)
Eq. ( 5) is a nonlinear Ordinary Differential Equation (ODE). It can be solved for the unknown wetting front depth Z t ( ) F using semi-analytical procedures which are briefly described in Appendix 1 for the sake of completeness. Alternatively, Eq. ( 5) could also be solved directly numerically in discrete time t ( ) n using any available robust ODE solver based for instance on explicit Runge-Kutta methods, or more general combinations of explicit/implicit temporal discretization methods, which are available for example in Matlab.

Once Z t ( )

F
is determined, eq.1 or eq.4b can be used to obtain the infiltration rate f t ( ). Finally, this single front approach can be considered also for other unsaturated flow processes, such as upwards flow, and other flow conditions, such as redistribution. However, the single front approximation is not accurate in general, except if the soil is very coarsegrained, with a narrow range of pore sizes, and initially very dry outside the wetted zone. Otherwise, the single front approach is inaccurate.

The moving Multi-Fronts (MMF) model

In this article, the approximate Green-Ampt (GA) model is generalized to the new "Moving Multi-Front" (MMF) model in order to deal efficiently and accurately with 1D vertical flows in partially saturated/ unsaturated porous media, taking into account both gravity and capillary effects on the process as realistically as possible.

In the new MMF model, the 1D flow process governed by Richards' PDE is represented more accurately. For example, the water content and pressure profiles cannot be reduced to the movement of a single "sharp front" in cases such as water infiltration in a fine-grained soil, and/or an initially wet soil, and/or, in the presence of a shallow water table. The main idea underlying the MMF approach is to extend the classical GA model by introducing more moving fronts in order to allow a better representation of the vertical spatial distribution of state variables (vertical moisture and pressure profiles). Thus, at any given time "t", water content z ( ) will varied with depth in a more realistic step-by-step/or "gradual" fashion, and similarly for pressure head h z ( ) and hydraulic conductivity K z K h z ( ) ( ( )). Therefore, in the MMF method, the unsaturated profiles K z ( ) and z ( ) are approximated by M moving fronts: this is shown schematically ). At left: 1D porous column, with a sharp wetting front separating the fully saturated and the totally dry regions. At right: an instantaneous water content profile (solid line), and the corresponding Green-Ampt approximation (dashed line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) in Fig. 2 for the moisture profile z ( ). With this parametrization of the multi-front model, we expect that the accuracy of the results will be improved as the soil column is divided (approximated) by a larger number of moving "fronts" over depth. Accordingly, a set of iso-values or "fronts" of the state variables (such as suction F ) will be selected within a range that depends on initial/boundary conditions, and then tracked through time based on their positions Z t ( , )

F F
(position of a given iso-value or front F ). Given a time-discretization, this formulation implies that we are looking at each new time step for the new position

+ Z t t ( )
F n n of the state variable's iso-value, or "front". Clearly, no spatial discretization is required.

For simplicity and clarity, we describe now the MMF model specifically for the case of downwards infiltration in a very dry soil. Other cases of 1D flow are treated with the same MMF model, by changing initial and boundary conditions. (See test problems in Section 4).

Assume we have defined M fronts, each characterized by its own suction value (subscript "F " refers to "Front"):

= { , , , } F F F F M _ 1 2
The main objective of the MMF model is to find the positions of the fronts Z t ( )

i F at each time step = Z t Z t Z t Z t Z t ( ) { ( ), ( ), ( ), , ( )} F F F F F M _ 1 2 3
According to Fig. 2, the multi-front approximation divides the 1D column into M + 1 zones:2

• Zones 1 to M: these zones are comprised between the ground surface ( = z 0) and the last (Mth) moving front [ = z Z t ( )] FM . They are the M transition zones from the totally saturated zone to the totally dry one. The parameters of each zone (number i) are K K i S , i S .

• Zone M + 1: this is the totally dry region; it is a semi-infinite dry zone below the Mth moving front [ = z Z t ( )]

FM

; the flux in this zone is null, and the parameters are:

= + K 0 M 1 , = + M d 1 .
The flow in each zone (i.e. between each pair of successive fronts) is governed by Darcy's equation, while a mass conservation equation is satisfied around each moving front.

Let us now develop the governing equations, using Darcy's fluxgradient law and mass conservation principles. We will use the suction head z t ( , ) as the main state variable, recalling that pressure head h z t ( , ) is just the opposite of suction head ( = h ). The advantage of using suction is that it is positive in the unsaturated zone, which is our main concern here. Recall also that "z" indicates depth (downwards).

Darcy's law applied between pairs of moving fronts

Applying Darcy's equation for zone i, where ( i M 1):

= = + q t K h z g K Z t Z t ( ) ( ) ( ) 1 i i i F F F F 1/2 1/2 1/2 i i i i 1 1 (6)
where quantities q t ( )

i 1/2
are mid-frontal fluxes, defined at intermediate positions between fronts Z t ( )

Fi 1 and Z t ( ) Fi . Equivalently, each flux q t ( ) i 1/2
approximates the flux in the mobile zone Z t Z t [ (), ()]

F F i i 1
. Similarly, quantities like K i 1/2 are mid-frontal or zonal hydraulic conductivities. Note that K i 1/2 does not depend on time, because it is a function solely of the two iso-values ( ,

F F i i 1
) through the conductivity-suction curve K ( ).

In the case of infiltration in a dry soil, the mid-front flux is assumed to be null through the bottom dry zone, i.e., between the next-to-last front M and the last front M + 1 (this emulates the case of a semiinfinite, initially dry soil column):

= + q t ( ) 0 M 1/2 .

Mass conservation (around each moving front)

Moreover, applying mass conservation equation around each of the M moving fronts Z t ( )

Fi

, we obtain:

= + + q t q t dZ dt ( ) ( ) ( ) i i i i i 1/2 1/2 1/2 1/2 (7)
where again q t ( )

i 1/2
represents a mid-frontal flux, and i 1/2 is the midfrontal water content between fronts Z t ( )

i 1 and Z t ( ) i
. Note that i 1/2 does not depend on time because it is a function solely of the two isovalues ( ,

F F i i 1
) through the water retention curve ( ). Recall that the bottom mid-front flux is null:

= + q t ( ) 0 M 1/2
. In addition, for the same reasons, the bottom water content is equal to the initial dry water content, that is:

= + M d 1/2
, where we usually take d r (residual moisture) for infiltration in a very dry semi-infinite soil column, as explained earlier.

The mass conservation Eq. ( 7) combined with Darcy's Eq. ( 6), applied to the moving fronts, yields a system of 1rst order nonlinear Ordinary Differential Equations (ODE's), or more generally here, a system of Algebraic Differential Equations. This system can be solved numerically for the unknown front positions

Z t { ( )} Fi subject to initial condition Z { (0)} Fi .
The numerical solution of such a system can be accomplished in discrete time t ( ) n using robust ODE solvers based for instance on explicit Runge-Kutta methods, or more general combinations of explicit/implicit temporal discretization methods. A suite of robust and reliable

F1 , Z t ( ) F2 , Z t ( ) F3 , Z t ( ) F 4
] where the suction head values are F1 , F2 , F3 , F 4 respectively. Right: instantaneous water content profile with MMF for M = 4 fronts (dashed line).

ODE solvers is available in Matlab, with possible choices including stiff systems, and also, systems that may require computationally more intensive implicit methods. In this article, an explicit stiff Matlab ODE solver function was used (mainly "ODE15S") in order to solve the system and obtain the evolving front positions Z t ( )

Fi

.

Once the evolution of front positions Z t [ ( )]

Fi is determined, the Darcy Eq. ( 6) can be applied to calculate fluxes between the moving fronts, including the top flux [q t ( )

1 1/2
] which represents the infiltration rate [ f t ( )]. The "bottom" flux is null for the infiltration problem in dry soil, but it can be non null and may need to be calculated, e.g. in the case of initially wet soils and/or in the presence of a shallow water table. Indeed, various types of unsaturated flow problems will be treated with the MMF model (see next Section 4). For instance, in the case of infiltration to a shallow water table, the position of front = 0 is fixed at bottom (because the water table is fixed at some finite depth, and because it corresponds physically to iso-value = 0). In that case, the bottom mid-frontal flux + q t ( )

M 1/2
is not null, and it can be computed as a function of time using the mid-frontal Darcy Eq. (6).

Parametrization of the MMF model (generalized Green-Ampt)

As in the classical GA model, the performance of the proposed MMF models depends on a suitable estimation of the model parameters. These parameters are the constant suction head [ Fi ] at each front "F i ", the "mid-front" unsaturated hydraulic conductivity [K i 1/2 ], and the "mid-front" water content [ i 1/2 ] between successive fronts. For clarity, we start with the parametrization of the simpler single front model (GA) before developing the parametrization of the multi-front model (MMF).

Parametrization of the single front model (classical GA model)

In the classical GA or single front model, the main parameter is the front suction [ F ], which needs to be related to soil hydraulic characteristics. Several expressions for F have been suggested in the literature, e.g. [START_REF] Panikar | Suction head at wet front in unsaturated-flow problems-A new definition[END_REF] [START_REF] Musy | Physique du sol, 6. Presses Polytechniques et Universitaires Romandes[END_REF]. Here we select the following two expressions:

1) The first expression for F , proposed here, corresponds to the point of the inflection of the water retention curve h ( ) or ( ) at which the specific moisture capacity [

= C h h h ( ) ( )/
] is maximum, as shown in Fig. 3. The suction head at this inflection point can be treated as a global capillary length scale of the porous medium. For the Van Genuchten model F , it can be easily calculated analytically (e.g. [START_REF] Ababou | Approaches to Large Scale Unsaturated Flow in Heterogeneous, Stratified, and Fractured Geologic Media[END_REF] and [START_REF] Alastal | Partially saturated oscillatory flow in a sandy beach (numerical modeling[END_REF])):

= n 1 1 1 F n 1/ (8)
where and n are the Van Genuchten parameters of the water retention curve ( ).

1) The second proposed estimate of wetting front suction F is the following, modified from [START_REF] Neuman | Wetting front pressure head in the infiltration model of Green and Ampt[END_REF], and originally proposed by Bouwer [START_REF] Bouwer | Unsaturated flow in ground-water hydraulics[END_REF][START_REF] Bouwer | Rapid field measurement of air entry value and hydraulic conductivity of soil as significant parameters in flow system analysis[END_REF] under the name "critical pressure":

= K d ( ) F r 0 (9) where = K K K ( ) ( )/
r S is the relative conductivity curve vs. suction. This estimation of F is illustrated graphically in Fig. 4. In the case of the exponential "Gardner" conductivity curve

= × K exp ( ) ( ) r G , the resulting front suction is simply = 1/ F G .

Parametrization of the multi-front model

The parameters [K i 1/2 ] and [ i 1/2 ] of the multi-front model can be determined by different methods. We propose here two possible classes of parametrization methods (the first method is the one used in Section 4 on Test Problems):

1) The first parametrization method is a simple approach based on discretizing the water content and the hydraulic conductivity K into (M + 1) frontal values i and K i defining M equally spaced segments. Alternatively, a variant of this method consists in discretizing directly the suction head i and then obtain the corresponding i and K i . After this discretization, average "mid-frontal"

i i 1 and K K [ , ] i i 1
, which do not depend on time but only on the corresponding suction iso-values i . In the four test problems of Section 4, the selected averaging method is the arithmetic mean for and the harmonic mean for K . Accordingly, the mid-frontal values are calculated as follows:

= + = + K K K ( ) / 2 {( )/2} . i i i i i i 1/2 1 1/2 1 1 1 1
The second method uses an integral criterion to calculate the fronts suctions i based on pre-selected values of mid-frontal conductivities K i 1/2 . This approach generalizes the method described in Eq. ( 9), where the front suction F of the single front "GA" model is obtained from an integral of the relative hydraulic conductivity curve K ( ) r . For multiple fronts, a generalization of this method consists in maintaining equal areas under the curve K ( ) r . An illustration of this method for M = 4 fronts is shown in Fig. 5. This parametrization method keeps the area under the curve K ( ) r equal to the area under the discrete front model approximation defined by M discrete frontal suction values and the corresponding frontal conductivities K i . For example, to find (or adjust) the front suction , the area under the K ) r curve between and is set equal to the area under the discrete approximation between and . This step is repeated between each two pairs of intermediate suctions in order to obtain the corresponding suction pressure at each front . Note: the notation in Fig. 5 corresponds to inter-frontal or mid-front suction values (also denoted 1/2 in the MMF equations) at which inter-frontal water contents ( i 1/2 ) and hydraulic conductivities (K i 1/2 ) are also calculated.

Test problems, validation, and performance of MMF model

To evaluate the proposed MMF model, the solutions obtained from the gridless MMF system of ODE's are compared to other solutions obtained by classical numerical solvers of the Richards PDE. In order to provide benchmark tests for our gridless MMF model, our choice here is to use extremely refined space-time grids in the classical Richards PDE solvers in order to use their numerical solutions as surrogate "true" solutions. Typically, the number M of fronts in our tests will be one order of magnitude less in the MMF model than the number of grid nodes used in the Richards PDE solvers. Specifically, in this article we compare our model to both Hydrus1D and BigFlow3D codes.

Hydrus1D is a well-documented and tested numerical model for simulating the 1D movement of water, heat, and multiple solutes in variably saturated/unsaturated porous media. It numerically solves the Richards' PDE for saturated-unsaturated water flow, coupled to Fickianbased advection dispersion equations for heat and solute transport. Hydrus1D is based on an implicit discretization in time, and 1D Galerkin Finite Elements discretization in space [START_REF] Simunek | The HYDRUS-1D Software Package for Simulating the One-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media[END_REF]. The linearization of the space-time discretized system is handled with Newton-Raphson iterations in Hydrus1D. Hydrus 1D also solves advective-diffusive/dispersive solute and several other transport processes coupled to unsaturated water flow.

Bigflow3D is a 3D finite volume flow code that has been widely documented, tested and benchmarked [START_REF] Ababou | BIGFLOW: A Numerical Code for Simulating Flow in Variably Saturated, Heterogeneous Geologic Media (Theory and User[END_REF][START_REF] Ababou | Testing procedures for spatially distributed flow models[END_REF]. Bigflow3D solves the 3D Richards PDE based on a generalized Darcy-type flux law, with a mixed pressure/moisture formulation of mass conservation. It uses an implicit discretization in time, and Finite Volume discretization in space. The linearization of the discretized space-time system is handled with a modified Picard iteration loop in BigFlow3D. It is capable of simulating various types of partially saturated flows and solves other flow problems such as planar groundwater flow, seawater intrusion, and surface flow governed by the diffusive wave equation.

There are several reasons for using two different classical codes for comparisons with the new MMF approach: first, we do not wish to rely entirely on a single code; secondly, there is some value added to these comparisons given that the two selected codes are quite different, Finite Elements vs. Finite Volumes (in fact, our auxiliary benchmark tests show that Hydrus & BigFlow's responses are undistinguishable, although presenting these auxiliary results is not the purpose of this paper); thirdly, there also minor technical reasons for using one code or the other depending on post-processing needs.

To sum up, the validation of the MMF model proposed in this study is based on comparisons with the classical "fixed grid" numerical solvers Hydrus1D and BigFlow3D. These "classical" Richards PDE solvers are considered as surrogates for "exact" solutions when used with extreme refinement in both space and time. A fine regular mesh Δz will be used here, corresponding to about one thousand fixed grid points in both solvers, and the time step Δt will be kept extremely small as well (e.g. one second or less for the Fine Sand tests). The criteria to be used for comparisons and validation of the MMF method are the transient water content and/or pressure head profiles ( z t ( , ); h z t ( , )), and also, the temporal evolution of boundary fluxes (top or bottom boundary flux depending on the proposed test).

Richards' equation & constitutive relations

The Richards PDE [START_REF] Richards | Capillary conduction of liquids through porous mediums[END_REF] is obtained by the combination of Darcy's law with mass balance equation. The 1D mixed form of Richards' PDE, describing vertical water flow in variably saturated/unsaturated homogeneous soils, is expressed as follows: In the case of unsaturated flow, both the hydraulic conductivity K ( ) and water content ( ) are functions of pressure head h ( ). Various models have been suggested to describe the constitutive relationships [ h ( ), K ( In this article, the unsaturated hydrodynamic properties of the soils are described with the soil water retention model h ( ) of van Genuchten [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils 1[END_REF], in combination with the functional model of hydraulic conductivity K ( ) proposed by [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF].

= h t K h h z gK h ( ) ( ) ( ) (10)
The water retention function is given by:

= + + h h ( ) [1 ( ) ] r s r n m ( 11
)
where r is the residual water content; s the saturated water content; m [ ] 1 is a scaling factor that represents an inverse capillary length scale; n is a dimensionless real number exponent, or shape parameter, related theoretically to pore size distribution (usually, > n 1); and m is also an exponent, related to n by the relation = m n 1 (1/ ). Inserting the van Genuchten curve h ( ) in Mualem's functional model K ( ) yields the unsaturated hydraulic conductivity-pressure curve K h ( ):

= + + K h K h h ( ) (1 ( ) ) 1 1 1 (1 ( ) ) S n n m 2 m 2 (12) where K m s [ / ]
S is the saturated hydraulic conductivity of the soil. For more details on unsaturated constitutive relations, see Ababou (2018) [Vol.1: Sec. 4.2, Sec. 4.3, Sec. 5.4].

Soil properties

Two different soils, namely Fine Sand (FS) and Guelph Loam (GL), were used in this study to evaluate the performance of the MMF model. The fine sand and the loam are presented in [START_REF] Alastal | Partially saturated oscillatory flow in a sandy beach (numerical modeling[END_REF]. The soil water retention curves are presented in Fig. 6. The curves indicated that the fine sand has a sharp transition from the saturated to the dry water content (Green-Ampt soil behavior). In contrast with the Guelph loam which have gradual transition behavior. The soil hydrodynamic parameters are summarized in Table 1 based on van Genuchten model.

Benchmark tests: results and discussion

The following tests were used to validate the MMF model and to show its applicability for a broad range of cases (flow regimes and soil types).

The selected tests are as follows:

• Infiltration in a dry soil, i.e., towards a deep water table: the purpose is to test and validate the MMF model for a prototype case, for which the classical GA model was theoretically designed (downwards infiltration in an initially dry semi-infinite soil column). In all the tests, the performance of the MMF model was analyzed mainly in terms of the temporal evolutions of two variables: the water content profiles [ z t ( , )], and the boundary flux [q t ( )].

Infiltration to deep water table (Test 1)

The test considered here is analogous to semi-infinite infiltration of the classical GA and it is characterized by monotonic water content profiles. This test treated 1D water infiltration towards deep water table to ensure that the soil is initially dry. Under infiltration, the top boundary condition is subjected to saturation and the wetting front started to move downwards into the previously dry soil, where the simulation should be stop far before the front reaches the bottom water table . Here is a summary of the initial and boundary conditions:

Initial condition (t = 0): = h z h ( , 0) 0 (depending on the soil as explained below)

Boundary conditions:

• Top BC (z = 0): = h t (0, ) 0, which corresponds to saturation.

• Bottom BC (z = L):

= h L t h ( , )
0 , which is compatible with the initial condition. This validation test was applied for the two soils (Guelph Loam [GL], and Fine Sand [FS]) with the properties presented earlier in Table 1. The chosen value of initial pressure head is = h 1 0 m for the GL soil, corresponding to a moderately dry water content [ = 0.414]; and = h 0.5 0 m for the FS soil, corresponding to almost totally dry water content [ = 0.033]. See Table 1 for comparison of these initial water contents with the saturated and residual water contents of each soil.

The column length is 1 m deep (L = 1 m). The vertical domain of the numerical PDE solver is discretized into fine discretization of 1001 nodes (Δz = 0.001 m) to ensure a better reference for comparison results, while only 30 moving fronts are used in MMF model to prove the robustness of the MMF even with a moderate number of fronts. The results show the downwards evolution of water content profiles, as well as the infiltration rate vs. time, as shown in Fig. 7 for the Guelph Loam and Fig. 8 for the Fine Sand.

Comparison between the meshless MMF model with M = 30 fronts, and the numerical solver with N = 1001 fixed nodes, shows an excellent agreement between the two models for the evolution of water content profiles and for the infiltration rate as well, and this, for each of the two soils tested here (Fine Sand and Guelph Loam).

Infiltration to shallow water table (Test 2)

This test deals with water infiltration to a shallow water table. The flow regime in this case differs from the previous one by two distinct characteristics: (i) the initial condition is non-uniform along the column both for moisture z ( ) and pressure h z ( ); and (ii) the state variable profiles obtained for this flow test are non-monotonic in space (vertically).

Consider a homogeneous soil column of length L = 1 m of Guelph Loam (GL), with the hydraulic properties described earlier in Table 1. The initial pressure head distribution h z ( , 0) is hydrostatic, implying that h z ( , 0) is a linear function of "z" starting from h = 0 at the water table to h = -L at ground surface. The lower and upper boundary conditions are zero pressure head. These initial and boundary conditions are summarized below.

• Initial condition (t = 0):

= h z z L ( , 0)
where "z" is depth measured from the ground surface.

• Boundary conditions: Top BC (z = 0): = h t (0, ) 0, corresponding to a constant saturation of ground surface Bottom BC (z = L): = h L t ( , ) 0, corresponding to a fixed water table at depth z = L. Note: this bottom BC is compatible with the initial condition at z = L.

In the numerical simulation performed with the Hydrus1D code, the soil column is discretized into 1001 nodes in order to obtain a finely discretized reference solution, while we use only 30 moving fronts in the MMF model.

Fig. 9 shows a good fit between the results of the two methods with respect to the evolution of water content profiles as well as the infiltration rate vs. time.

The agreement between the two models can be enhanced by increasing the number of fronts in the MMF model. A quantitative analysis of global error in terms of the number of fronts M for this test is presented separately in Section 5.

It should be noted that, in this test, the soil column is partially saturated. It contains a saturated zone that increases with time, until the column reaches full saturation at steady state (large times). As a consequence, some pressure fronts should disappear with time when the corresponding pressure iso-values no longer exist. To solve this problem, an auxiliary numerical scheme was implemented in the ODE solver of the MMF model to monitor this process and to manage the emergence and disappearance of pressure iso-values ("fronts"). This auxiliary numerical scheme includes two steps: (1) detection and removal of any successive fronts

+ Z t Z t [ ( ), ( )] i i 1
whose positions may be approaching each other (a threshold value δZ was pre-defined in the script); and (2) redistribution of the remaining front positions through spatial interpolation techniques.

Capillary rise from a fixed water table (Test 3)

In this test, we analyze unsaturated capillary rise in a soil or porous medium. This test, in contrast with the previous tests, represents an anti-gravity flow problem. A soil column of length L = 1 m, made up of the Guelph Loam soil (GL), is assumed to be relatively dry initially (see initial condition below). A zero-pressure head is suddenly imposed at the bottom boundary, which corresponds to imposing a fixed water table at depth = z L for all times > t 0. The initial and boundary con- ditions of this validation test are as follows:

• Initial condition = t ( 0): = h z m ( , 0)
1 , corresponding to initial water content = 0.414.

• Boundary conditions:

Top BC = z ( 0): = h t m (0, )
1 , which is compatible with the initial condition.

Bottom BC = z L ( ): = h L t ( , ) 0, which corresponds to a fixed water table at = z L. Note: the initial pressure condition = h z m ( , 0) 1 corresponds to a water content θ = 0.414 and a degree of saturation S = 64.9%; this may seem fairly wet, but in fact the corresponding initial conductivity of the soil is 1.81 E-7 m/s, about 20 times smaller than the saturated conductivity Ks = 3.66E-6 m/s; therefore, the soil is initially relatively dry, in the sense that it is initially poorly conductive.

In order to verify the accuracy of our MMF model, a numerical benchmark simulation was implemented using the finite element code Hydrus1D with N = 1001 nodes, versus the present MMF model with • Infiltration to shallow water table: the purpose is to test and validate the MMF model for non-uniform initial water content along the column, and such that the vertical profiles of water content and the hydraulic conductivity are non-monotonic.

• Capillary rise problem (to test and validate the MMF model for antigravity flow problem).

• Gradual rise of water table (to test and validate the MMF model for variable pressure head boundary condition). M = 30 moving fronts. Fig. 10 shows an excellent agreement between the two models, both in terms of moisture and pressure profiles (top plots) and in terms of flux evolution q(t) (bottom plot). Note that q(t) represents the upward flux through the water table during the capillary rise process. The agreement is remarkable given that only 30 moving fronts were used in the MMF method.

Gradual rise of a water table (Test 4)

The objective of this test is to validate the MMF model against timevarying pressure head boundary condition. In addition, unlike the previous tests, the proposed test involves a partially saturated/unsaturated soil column.

A semi-infinite fine sand column (FS), with properties defined in Table 1, is initially hydrostatic, with a water table located at elevation 0.5 m measured from the bottom of the column (in this test, the vertical coordinate is elevation rather than depth). The water table is forced upward by gradually rising the bottom pressure h t ( ) 0 from = h 0 + 0.5 m to = h 0 + 1.0 m. For the purpose of this test, we select this time varying "entry pressure" function In this benchmark test, four output variables are selected to compare the two models: the water content profiles z t ( , ), pressure head profiles h z t ( , ), the evolution of water table elevation Z t ( ) WT , and the flux q t ( ) 0 calculated at the bottom boundary.

Performance and accuracy

In this section, the numerical performance of the proposed MMF model is studied in two different ways: (1) through global analyses of error norms (moisture content and surface flux) vs. number of fronts (M); and (2) through direct analyses of surface flux f(t) [m/s] for several values of the number of fronts (M = 1, 2, 4, 8) where the case M = 1 is taken to mean the classical Green-Ampt model in the case of Test 1 (infiltration in the case of deep water table).

Numerical analyses of error norms

In this section, the gridless MMF model is evaluated in terms of its accuracy. An appropriate global error norm (ε) is defined, calculated, and then analyzed as a function of the number M of moving fronts. Since analytical solutions are not generally available for this type of nonlinear transient flow problem, given the highly nonlinear soil properties involved, we use here the classical numerical solvers with a large number of nodes (e.g. N = 1001 nodes) as a surrogate for the exact solution.

With this in mind, the error is calculated for two types of output variables: the space-time distribution of water content z t ( , ), and the time-dependent boundary flux (e.g. infiltration rate f(t)). We choose to implement error norm criteria for Test 2 of the previous section, i.e., water infiltration to a shallow water table, Guelph Loam soil (Section 4.3.2: Test 2).

For the water content, the global error is calculated by applying a discrete form of the space-time Euclidian L2-norm, where integration is replaced by discrete summation over both space and time, as follows:

= z t z t z t Z t 1 [ ( , ) ( , )] i i n n n i NUM i n MMF i n i n 2 1/2 (13)
where:

• NUM : is the volumetric water content resulted from the fine re- solution numerical solver with large number of nodes (N = 1001).

• MMF : is the volumetric water content resulted from the MMF model with a moderate number of moving fronts (here M ≪ 1000).

It is worth mentioning that, in order to calculate the above error, we interpolate linearly the moving front solution MMF between the fixed nodes of the classical numerical solver (Hydrus1D), since the positions of the MMF moving fronts do not coincide with the spatial discretization (fixed nodes) of the classical numerical solver.

For the boundary flux f t ( ), which corresponds to the infiltration rate for Test 2, the global error is calculated by applying a discrete form of the Euclidian L2-norm over time, as follows:

= t f t f t t 1 [ ( ) ( )] F n n n NUM n MMF n n 2 1/2 (14)
where:

• f NUM : is the boundary flux or infiltration rate resulting from the fine resolution numerical solver with a large number of fixed nodes (N = 1001).

• f MMF : is the boundary flux or infiltration rate resulting from the MMF model with a moderate number of moving fronts (here M ≪ 1000).

For water contents z t ( , ), the behavior of the normalized global error as a function of the number M of moving fronts is plotted in Fig. 12 on a log-log scale.

It is clear from the graph of Fig. 12 that the error decreases as a power function with the number of fronts M of the MMF model. Based on a log-log regression of the numerical error norm vs. M, it is found that the error norm behaves as follows:

= ×M m m 0.6387 [ / ] 1.92 3 3 (15)
The accuracy of the MMF model is characterized by the number of moving fronts (M), or by its inverse (1/M). Note that the 1/M parameter of the MMF model, is similar to the mesh size parameter Δz of classical models (e.g. linear Finite Elements for Hydrus1D, structured Finite Volumes for Bigflow3D). The conclusion from Fig. 12 and Eq. 15 is that, in terms of (1/M), the order of accuracy of the MMF method is close to 2. We can conclude in particular that the method is consistent and provides a good approximation of water content profiles.

One should also investigate further whether the estimated error norm (a global measure of error) is small relative to the calculated water content profiles z t ( , ). One way to answer this question is to compare the global error m m

[ / ] 3 3 to the range of moisture contents of infiltration flow Test 2 (Fig. 9). In this test, infiltration to shallow water table, the range of moisture contents is defined by the initial moisture at For the error on the infiltration rate, or boundary flux f t ( ), a similar power law behavior is obtained. The calculated global error F decreases as a power law with the number of moving fronts (M) of the MMF model, as can be seen in Fig. 13. Based again on a linear regression log-log scale, error F appears to follow the power law:

× M m 8.696 10 •[ / ] F 6 0 . 9 9 (16)
The conclusion from Fig. 13 and Eq. ( 16) is that, in terms of (1/M), the order of accuracy of the MMF method is about 1 for the flux variable. We can infer in particular that the method is consistent for flux calculations as well.

Technically, one should also investigate whether the flux error norm

Performance of MMF front in terms of infiltration rate f(t) vs. number of fronts

The temporal evolution of surface flux f(t) is a particularly important variable in many applications of soil hydrology. For this reason, we now focus on infiltration Test 1 (deep water table) in order to directly compare the surface flux f(t) [m/s] obtained for several values of the number of fronts (M = 1, 2, 4, 8), including the classical Green-Ampt model represented by the case M = 1. In that case (M = 1 front), the classical Green-Ampt solution is implemented here with the widely used front suction model F of Bouwer [START_REF] Bouwer | Unsaturated flow in ground-water hydraulics[END_REF][START_REF] Bouwer | Rapid field measurement of air entry value and hydraulic conductivity of soil as significant parameters in flow system analysis[END_REF] given by Eq. ( 9), and selecting K S and S for the values of mid-frontal conductivity and water content. Fig. 14 shows that the infiltration rate agreement between the MMF model and a Richards PDE solver (Hydrus 1D) improves rapidly with M, and becomes almost perfect for M = 8, which is remarkable. The figure also shows that the single front model (classical Green-Ampt) does not perform very well, as expected (especially for this type of loamy soil). As a reminder, a comparison for M = 30 fronts was also shown earlier in Fig. 7. In conclusion, the present MMF model (generalized Green-Ampt) performs quite well for just a relatively small number of fronts (M = 8) in terms of infiltration rate.

Conclusions and outlook

The Moving Multi-Front (MMF) model presented in this paper is essentially a meshless method for calculating vertical water flows in unsaturated or partially saturated homogeneous soil columns, under the combined influence of gravity and capillary gradients. The MMF model requires the solution of a nonlinear ODE system governing the positions Z t ( ) In addition, two error norms were studied: the error norm of the space-time distribution of moisture z t ( , ), and the error norm F of the time-variable infiltration rate or flux f t ( ). It was found that these error norms behave like a power of the number of fronts M: M ~with power ~1.92 for , and ~0.99 for F . It can be concluded, first, that the MMF method is numerically consistent, in the sense that the error norm decreases as M increases. Secondly, it is interesting to note that, in terms of M (1/ ), the order of accuracy of MMF is close to 2 for water content and 1 for flux. The same order of accuracy is obtained for instance with classical centered finite difference schemes applied to the Richards PDE with regular mesh size Δz: the order of accuracy is O z ( ) 2 for pressure or moisture, and O z ( ) 1 for the calculated mid-nodal flux (see truncation error analyses in [START_REF] Vauclin | Résolution numérique d'une équation de diffusion non linéaire: application à l'infiltration de l'eau dans les sols non saturés[END_REF], and [START_REF] Ababou | Testing procedures for spatially distributed flow models[END_REF] among others). Further analyses and rescaling of our calculated errors indicate that the relative errors on both moisture and flux are quite small; typically 1%, using only M ≈ 30 moving fronts in the MMF model.

In conclusion, the behavior of the MMF method is quite satisfactory. It appears from global errors that only a few tens of fronts are needed to obtain accurate solutions, as can be seen from direct comparisons with the classical solvers in Section 4.3 (Tests 1,2,3,4), and from the error norm plots of Section 5.1. In addition, the infiltration rates f(t) of Test 1 were systematically compared with the reference solution for M = 1,2,4,8, and these tests confirmed the previous results: infiltration rate f(t) is in very good agreement for all times, even for moderate number of fronts (e.g. M = 8).

Overall, the results obtained with the meshless MMF model, compared with finely discretized space-time solutions of the nonlinear Richards PDE, demonstrate the efficiency and reliability of the MMF. The MMF model has distinct advantages, including its ability to track directly a given pressure or water content in time, and its robustness in terms of pressure and moisture profiles.

Note: CPU times benchmark tests were not implemented because the classical FE or FV solvers (Hydrus, BigFlow) are coded in compiled languages (Fortran), while our MMF method at this stage is programmed in an interpreted language (Matlab). Interpreted languages are notoriously slower than compiled languages in terms of execution time, and the CPU times comparisons would be meaningless. Nevertheless, the performance results of MMF already suggest that CPU times could be significantly reduced with MMF once re-programmed in a compiled language (such CPU time tests should be run for equal precision of results).

Given these encouraging results, the MMF model is currently undergoing several extensions. Current work is devoted to generalizing the MMF method to include flux and pressure gradient boundary conditions, such as a prescribed Darcy flux q m s [ / ] z at ground surface (e.g. to simulate internal drainage under zero flux), and/or, a "gravitational" boundary condition ( = h z / 0) at some depth. The generalization to time-variable boundary conditions, including the switching from pressure to flux condition (and vice-versa), will open the way for multiple events such as storm/inter-storm scenarios. On the other hand, the case of multi-layered soils, involving one or more interfaces of material discontinuity between layers, is also being considered currently. The main issues with multi-front modeling of flow in a multi-layered soil have to do with a combination of several factors: thus, non-monotonic pressure profiles require detection and updates of the mobile pressure fronts; furthermore, when pressure fronts move through a material interface, the discontinuous pressure gradient must be dealt with, as well as the discontinuous moisture content (while pressure itself remains continuous). Some of these issues may be resolved by using more robust and adaptive ODE solvers (ongoing work).

Fig. 1 .

 1 Fig. 1. Green-Ampt piston flow approximation: single wetting front model, illustrated here for water infiltration in a dry soil (deep water table). At left: 1D porous column, with a sharp wetting front separating the fully saturated and the totally dry regions. At right: an instantaneous water content profile (solid line), and the corresponding Green-Ampt approximation (dashed line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .

 2 Fig. 2. Schematic diagram of the Moving Multi-Front (MMF) model approximation for water infiltration towards a deep water table, with M = 4 fronts (note: the "z" axis is downwards, and we show only four fronts for clarity). Left: the 1D porous column showing the ground surface and four successive moving fronts at depths [Z t ( )

Fig. 3 .

 3 Fig. 3. (a) Water retention curve ( ) showing its inflection point; (b) Capillary moisture capacity curve C ( ). The characteristic suction corresponding to the inflection point of ( ) and, equivalently, to the maximum soil moisture capacity C ( ), is one of the proposed estimates of front suction F for the single front GA model.

Fig. 4 .

 4 Fig. 4. Relative hydraulic conductivity curve and its "single-front model" approximation, showing the front suction value F according to Bouwer's criterion.

Fig. 5 .

 5 Fig. 5. Relative hydraulic conductivity curve with the MMF model approximation. In this example, only 4 fronts are shown. The front suctions are denoted F i , and the intermediate inter-frontal suctions are denoted

Fig. 6 .

 6 Fig. 6. retention curves h ( ) for the Fine Sand [FS] (lower red curve) and for the Guelph Loam [GL] (upper blue curve). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

  , the bottom boundary pressure gradually increases from = h 0 + 0.5 m to = h 0 + 1.0 m within 150 s.

  Fig. 11 shows an excellent fit between the MMF model with M = 30 moving fronts and the numerical solver of Richards' PDE (Bigflow3D code) with N = 1001 nodes, and a fine time step Δt ≈ 0.2 s. The agreement is excellent between the two methods (with only 30 moving fronts in the MMF).Recall that the flow is forced by the positive bottom pressure h t be seen from the bottom left plot of Fig.11. The bottom flux q t ( ) 0 was calculated in the MMF by using Darcy's law between the first two bottom fronts.

Fig. 7 .

 7 Fig. 7. Water infiltration to deep water table (Test 1) for the Guelph Loam [GL]. Left: water content profiles at different times [hours]. Right: temporal evolution of infiltration rate [times in seconds]. Solid curves represent the simulation with Hydrus1D using 1001 fixed nodes. Symbols represent the MMF model with M = 30 moving fronts. Dashed line represents initial water content [ = 0.414].

Fig. 8 .

 8 Fig. 8. Water infiltration to deep water table (Test 1) for the Fine Sand [FS]. Left: water content profiles at different times [minutes]. Right: temporal evolution of infiltration rate [times in seconds]. Solid curves represent the simulation with Hydrus1D using 1001 fixed nodes. Symbols represent the MMF model with M = 30 moving fronts. Dashed line represents the initial water content [ = 0.033].

Fig. 9 .

 9 Fig. 9. Water infiltration to shallow water table (Test 2) for the Guelph Loam [GL]. Top left: water content profiles at different times [minutes]. Top right: pressure head profiles at different times [minutes]. Bottom plot: temporal evolution of infiltration rate. Solid curves represent the classical numerical solution of Richards' PDE with Hydrus1D using N = 1001 nodes. Symbols represent the MMF model with M = 30 fronts. Dashed lines represent initial water content or pressure profiles.

Fig. 10 .

 10 Fig. 10. Capillary rise (Test 3) in a column of Guelph Loam [GL] until steady state. Top left plot: water content profiles at different times [hours]. Top right plot: pressure head profiles at different times [hours]; a hydrostatic steady state is reached at t = 24 h. Bottom plot: evolution of bottom flux q(t) at z = L [times in seconds]. The solid curves represent the numerical simulation with Hydrus1D using N = 1001 nodes, and the symbols represent the moving front positions of the MMF model with M = 30 fronts. Dashed lines represent the initial condition in the moisture and pressure plots.

F

  (a global measure of error) is small relative to the calculated flux f(t) at various times. The flux error norm obtained for M = 30 fronts is F 2.9769 E-07 m/s. It is quite small compared to the maximum flux rate f = 3.5 E-5 m/s observed at early times in Fig.9(Test2), and it remains reasonably small compared to the final steady flux, which is close to the saturated conductivity of the Guelph Loam,

Fig. 11 . 0 .

 110 Fig. 11. Partially saturated simulation of gradual water table rise (Test 4) in Fine Sand [FS]. Upper left plot: water content profiles z t ( , ) [times in seconds]. Upper right plot: pressure head profiles h z t ( , ) [times in seconds]; the blue line represents = h 0 (recall h 0 in the saturated zone and < h 0 in the unsaturated zone). Lower left plot: evolution of water table elevation vs. time [seconds]. Lower right plot: evolution of bottom flux q t ( ) 0 . Solid curves represent the numerical simulation (Bigflow3D) with N = 1001 nodes. Symbols represent the MMF model with M = 30 moving fronts. On the lower left plot, the dashed curve represents the pressure head h t ( ) 0 imposed at the bottom boundary. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fi

Fi

  Fig. 12. Global error norm versus the number of moving fronts M, for the water content variable z t ( , ). The dashed blue line is obtained by linear regression on a log-log scale. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Zones are inter-frontal regions; conductivity K and flux q are zonal or midfrontal quantities; accordingly the zonal quantities like the conductivity K of "Zone (i)" will be later subscripted by (i-1/2).

Acknowledgements

The first author would like to extend his sincerest thanks and appreciation to Welfare association, Palestine, for their support to his research visit to the IMFT, France.

Declaration of Competing Interest

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

Appendix 1

This appendix briefly describes the semi-analytical procedure for solving the governing equation of the classical Green-Ampt (GA) infiltration model, which corresponds also to the single-front model (with M = 1 front) in our generalized Moving Multi-Front method. Let us start with the GA governing equation developed in the text (Eq. 5, renamed here Eq. A1):

Eq. ( A1) is a nonlinear Ordinary Differential Equation (ODE). It can be solved for the unknown wetting front depth Z t ( ) F using either semianalytical or fully numerical procedures (ODE solvers). We focus here briefly on semi-analytical procedures.

1. Eq. (A1) can be directly integrated quasi-analytically in time, to obtain a non-differential equation F Z t t ( ( ); )

2. The short time and large time behavior of Z t ( ) F can be analyzed directly from this Eq. (A2) without recourse to numerical methods. Short and large times can distinguished based on the following characteristic time:

T CHAR plays the same role for Green-Ampt infiltration as the gravitational time T GRAV in Philip's infiltration theory. The short time solution (t T ) CHAR is dominated by capillary suction:

The