A Transformation from a reference element to an element in physical space and vice versa and streamline calculation with MHFEM A.1 Transformation from a reference element to an element in physical space . . A.2 Transformation from an element in the physical space to a reference element A.3 Streamline calculation with MHFEM . . . . . . . . . . . . . . .

List of Figures

. Homogeneous transmissivity case without barrier. 22 2.7 Piezometric head and streamlines (magenta lines) computed with a leakage coefficient of 6 m -1 d -1 . Homogeneous transmissivity case with the barrier. 22 2.8 Piezometric head and streamlines (magenta lines) computed with a leakage coefficient of 6 m -1 d -1 . Homogeneous transmissivity case without barrier. 24 2.9 Piezometric head and streamlines (magenta lines) computed with a leakage coefficient of 6 m -1 d -1 . Heterogeneous transmissivity case with the barrier. . ).Homogeneous transmissivity case. . . 3.16 Piezometric head and streamlines (blue lines) computed with a leakage coefficient of 9 m -1 d -1 and approximated streamlines (green lines) computed with (C λ , δ C λ )=(6 m -1 d -1 ,3 m -1 d -1 ). Homogeneous transmissivity, case with the barrier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.17 Piezometric head and streamlines (blue lines) computed with a leakage coefficient of 9 m -1 d -1 and approximated streamlines (green lines) computed with (C λ , δ C λ )=(6 m -1 d -1 ,3 m -1 d -1 ). Heterogeneous transmissivity, case without barrier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Piezometric head and streamlines (blue lines) computed with a leakage coefficient of 9 m -1 d -1 and approximated streamlines (green lines) computed with (C λ , δ C λ )=(6 m -1 d -1 ,3 m -1 d -1 ). Heterogeneous transmissivity, case with the barrier. .

Abbreviations and definitions

General introduction

Throughout the world, groundwater supplies are hidden vital resources that are facing rising pressure owing to pollution and over-consumption from anthropogenic activities. Pollutions putting groundwater at risk notably include: discharge of waste and wastewater onto or into the ground, use of chemicals such as fertilizers and pesticides, spreading of slurry, manure and abattoir wastes, poor storage of solvents, petroleum products by using above or underground storage tanks. These last facilities represent a significant contamination risk to groundwater. Mitigation measures can be used to prevent soil and water pollution from waste disposal, land fill sites, septic or chemical storage tanks. Among them, drains and impervious barriers may be set up. A drainage system is assumed to be designed to convey all potentially contaminated water and spills of fuel to suitable collection points for disposal or treatment. Besides, specific measures can be installed to prevent groundwater pollution. The construction of facilities such as constructed wetlands (CWs) for the treatment of urban rainy weather can be a promising wastewater treatment system where chemical and biological processes are similar to those found in natural wetlands. CW systems are increasingly considered as a natural and rustic treatment process by local authorities in developing and developed countries. Their significant efficiency may be evaluated regarding their designs and dimensioning by using computational simulation tools. This thesis focuses mainly on two cases of study: 1) the modeling of the drain-aquifer exchanges with an evaluation of the sensitivity of the underground flows with respect to the drainage exchange coefficient 2) the modeling of the flow in porous media variably saturated within a CW using inverse approach parameter identification. This dissertation comprises two parts each of which includes two chapters organized as: Introduction, State of the art, Methodology, Results and discussions, and Conclusions. In addition, seven useful appendices provided at end of this thesis. These contain materials relevant to mathematical functions, discretization of equations, detailed statistical data and results whose direct inclusions in the main body would obscure the manuscript's readability.

In part I, in a case of underground groundwater storage site in which drains and an impervious barrier are implemented, computational simulations are used to assess the effectiveness of these installations. Under the Dupuit-Forchheimer assumption, in the second chapter, the groundwater flow and the leakage drain/aquifer interactions are implemented in a conforming finite element method (FEM) and a mixed hybrid finite element method (MHFEM) in a horizontal two-Chapter 1 -General introduction dimensional domain modeling a regional aquifer below chemical storage tanks. As a complementary part for Chapter 2, in Chapter 3, a sensitivity analysis of the piezometric head, velocity field and streamlines is carried out with respect to the leakage coefficient to provide a new insight into groundwater waterbody exchanges. It moreover allows for a careful analysis of the differences between FEM and MHFEM implementations for both an ideal homogeneous hydraulic conductivity over the domain, and a heterogeneous one.

In the second part, in Chapter 4, a flow modeling in vertical variably saturated CWs is carried out by implementing the Richards model by means of a MHFEM especially well adapted to the simulation of heterogeneous media. A particular attention is brought to boundary condition modeling (surface ponding or evaporation) to be able to tackle different sequences of rainfall-runoff events as well as the soil oven-dry conditions.

Finally, in Chapter 5, for a proper parameter identification, large field data-sets would be needed. As van Genuchten-Mualem (vGM) soil hydrodynamic parameters depend on water content and are notoriously difficult to be measured, their estimation is subject to considerable experimental and numerical studies. Meticulous attention is brought to the sensitivity analysis performed with respect to the vGM parameters that reveals a predominant influence of the shape parameters and the saturated conductivity of the filter on the piezometric heads, during saturation and desaturation. As these are usually not available, notably due to the randomness of the storm events, we thus propose a simple, robust and low-cost numerical method for the inverse modeling of the soil hydrodynamic properties. The calibration of hydrodynamic parameters for subsurface CWs is a sensitive process and remains to be a challenging task since unsaturated flow modeling involves highly non-linear equations. To that end, a data assimilation technique is implemented by applying automatic differentiation (AD) to augment computer codes with derivative computations. Identification experiments are conducted by comparing measured and computed piezometric head by means of the least square objective function. The temporal variability of hydrodynamic parameter is then assessed and analyzed.

Part I

Modeling of drain-aquifer exchanges in saturated media

French Abstract

Les outils de modélisation des eaux souterraines sont intéressants pour étudier le risque de contamination souterraine par les fuites des réservoirs de stockage et pour évaluer la performance des systèmes de drainage conçus pour prévenir la pollution des eaux souterraines. Le chapitre 2 s'applique à modéliser l'impact d'un drain mis en place dans un aquifère non confiné pour contrôler les flux d'eau afin de minimiser le risque de migration des polluants. Dans le chapitre 3, deux modèles numériques sont mis en oeuvre afin d'être étudiés et comparés. Le premier modèle est basé sur une méthode d'éléments finis (FEM), le deuxième est basé sur une méthode d'éléments finis mixtes hybrides (MH-FEM). La méthode MHFEM conduit à un système d'équations dont les inconnues sont les niveaux piézométriques moyens et les flux à travers les facettes, alors que les inconnues sont les niveaux piézométriques aux noeuds du maillage dans le modèle FEM. On rappelle que le modèle MHFEM peut conduire à un champ de vitesse de meilleure qualité grâce à l'hypothèse de continuité de la composante normale de la vitesse utilisée dans l'approche mixte. Deux cas considerés comme homogène et hétérogène sont étudiés pour un site hébergeant sur une véritable plate-forme industrielle chimique située entre un canal de navigation à l'ouest et une rivière à l'est.

Les flux calculés au niveau du canal, de la rivière et du drain ont un comportement similaire mais des valeurs distinctes en raison de la différence dans la modélisation des échanges de masse d'eau. Dans tous les cas, le drain capte plus d'eau en utilisant la FEM. Les écarts entre les solutions FEM et les solutions MHFEM sont difficiles à évaluer à partir de calculs de modèles. De la même manière, les simulations de l'écoulement avec conductivité hydraulique hétérogène naturelle sans barrière, ont été réalisées et comparées avec celui avec la barrière. Cela confirme que la barrière imperméable augmente l'efficacité du drainage. Les flux calculés au niveau du canal, de la rivière et du drain, ainsi que l'erreur de bilan de masse estimée comme la somme du débit d'entrée du canal et du débit de sortie de la rivière et des drains montre que les deux codes sont conservateurs. En outre, il est à remarquer que les drains captent beaucoup plus d'eau dans les cas FEM. Ceci est moins considérable dans le cas hétérogène.

En fournissant beaucoup plus d'informations sur les hauteurs piézométriques et le champ de vitesse, une analyse de sensibilité révèle l'influence d'une perturbation du coefficient d'échange de manière claire. Le chapitre 3 est consacré à une question classique en sciences de l'environnement : l'évaluation de la sensibilité des « sorties » d'un modèle numérique aux variations des paramètres d'entrée utilisés dans ce même modèle. Ces paramètres peuvent être des paramètres descripteurs du milieu physique ou peuvent être des conditions initiales ou des conditions aux limites. La différentiation automatique (AD) constitue un ensemble de techniques conçues pour « augmenter » les codes numériques avec des calculs de dérivées. Les codes peuvent ainsi être utilisés pour étudier la sensibilité de différentes variables par rapport à différents paramètres. Parmi les outils logiciels disponibles, nous avons choisi Tapenade qui offre toute la généralité de la AD et fournit une interface utilisateur facilitant l'implémentation de cette analyse de sensibilité sur les deux codes FEM et MHFEM. Il a été nécessaire d'effectuer une restructuration des codes à différencier pour qu'ils soient utilisables par Tapenade en vue d'implémenter la AD. Une méthodologie générale est proposée pour différencier les codes de l'équipe Mecaflu dans le paragraphe "Minimal effort AD strategy of existing codes" de la thèse. Les codes différenciés en mode linéaire tangent sont validés en utilisant du test de Taylor effectué sur les hauteurs piézométriques et le champ de vitesse. Nous avons choisi d'étudier la sensibilité des résultats des simulations (hauteur piézométrique, vitesse, ligne de courant) par rapport au coefficient d'échange drain-nappe. Cette analyse de sensibilité est conduite pour les deux modèles numériques utilisés pour résoudre des problèmes posés (FEM et MHFEM) sur des domaines de conductivité hydraulique homogène d'une part, et hétérogène d'autre part. Les résultats montrent clairement le rôle de la barrière imperméable car les sensibilités par rapport au coefficient d'échange sont très faibles en dehors de la zone protégée par la barrière. On observe aussi que le niveau piézométrique calculé par MHFEM est plus sensible que celui de FEM, notamment à proximité du drain. Bien que les drains captent moins d'eau en utilisant la MHFEM, on peut voir que les hauteurs piézométriques calculée par MHFEM sont les plus sensibles près du drain. On note également que les lignes de courant reconstruites en utilisant la AD sont en bon accord avec celles calculées avec les codes non différenciés. A notre connaissance, ces résultats sont nouveaux dans la modélisation de l'écoulement des eaux souterraines. Ils ont été publiés dans [START_REF] Moezzibadi | Sensitivity of groundwater flow with respect to the drain-aquifer leakage coefficient[END_REF].

CHAPTER 2

Drain-aquifer exchanges in saturated porous media

Introduction

Mitigation measures may be used to prevent soil and water pollution from waste disposal, land fill sites, septic or chemical storage tanks. Among them, drains and impervious barriers may be set up. The efficiency of this technique can be evaluated by means of groundwater modeling tools. The groundwater flow and the leakage drain-aquifer interactions are implemented in a conforming finite element method (FEM) and a mixed hybrid FEM (MHFEM) in a horizontal two-dimensional domain modeling regional aquifer below chemical storage tanks. Computations are performed with both an ideal homogeneous hydraulic conductivity and a realistic heterogeneous one.

State of the art

Groundwater contamination commonly occurs due to human activities whereby chemicals or wastes may be released to the environment, either intentionally or accidentally. Groundwater can be exposed to these chemicals by means of infiltration from surface runoff in the landfill. Pollutions putting groundwater at risk notably include: discharge of waste and wastewater onto or into aquifers, use of chemicals such as fertilizers and pesticides, spreading of slurry, manure, and abattoir wastes, poor storage of solvents and petroleum products by using above-ground or underground storage tanks. These last facilities represent a significant contamination risk for groundwater if appropriate mitigation measures, including surface and subsurface drainage, are not designed in their conception and construction [Environment-Agency, 2017]. The drainage system should be designed to convey all potentially contaminated water and spills of fuel to suitable collection points for disposal or treatment. The Groundwater Daughter Directive (2006/118/EC) clarifies the requirements for measures to prevent or limit inputs of pollutants into groundwater. Factors to be considered when carrying out risk assessment for tank leaking into the ground are identified in Environment-Agency [2017].

In order to confront subsurface pollution or to address the treatment, it is essential to understand the mutual inter-influence of surface waters and groundwaters. Furthermore, with growing use of water resources and increasing uncertainties in water supply, investigating groundwater and surface water as an integrated system is crucially needed.

Early, [START_REF] Rushton | Possible mechanisms for leakage between aquifers and rivers[END_REF] studied the leakage between aquifers and rivers by using an idealized one-dimensional problem. More generally, groundwater modeling tools allow for the computation of the exchanges of water between streams and aquifers and their influences on the quality and quantity of water within both domains (see [START_REF] Fleckenstein | River-aquifer interactions, geologic heterogeneity, and low-flow management[END_REF]; [START_REF] Ruehl | Differential gauging and tracer tests resolve seepage fluxes in a strongly-losing stream[END_REF]; [START_REF] Doppler | Field evidence of a dynamic leakage coefficient for modelling river-aquifer interactions[END_REF]; [START_REF] Sanz | Modeling aquifer-river interactions under the influence of groundwater abstraction in the mancha oriental system (se spain)[END_REF] and the references therein, for instance). Throughout the world, groundwater supplies are hidden vital resources that are facing rising pressure owing to pollution and overconsumption from anthropogenic activities.

The exchanges between an aquifer and a river, which are generally complex and dependent on many physical factors, are of great concerns in hydrological studies regarding the quantification of groundwater contributions to rivers [START_REF] Ellis | Quantifying urban river-aquifer fluid exchange processes: A multiscale problem[END_REF] and the investigation of its ecological significance. Numerous approaches and techniques are described in the literature by [START_REF] Kalbus | Measuring methods for groundwater-surface water interactions: a review[END_REF] to estimate fluxes at the groundwater-surface water interface. Low river flows are commonly controlled by exchanges between groundwater and surface water, the magnitude of which is controlled by hydraulic characteristics of the aquifer and the aquitard materials overlying the aquifer [START_REF] Fleckenstein | River-aquifer interactions, geologic heterogeneity, and low-flow management[END_REF]. [START_REF] Brunner | Disconnected surface water and groundwater: From theory to practice[END_REF] discussed the influence of conceptual assumptions on simulation results of the interaction between disconnected streams and groundwater. [START_REF] Wang | A quantitative analysis of hydraulic interaction processes in stream-aquifer systems[END_REF] showed the importance of the transition from connection to disconnection for stream-aquifer interactions and demonstrated that the hydraulic connectedness of the stream-aquifer system can reach a critical disconnection state depending on the horizontal and vertical hydraulic gradient at the free water surface. [START_REF] Ellis | Quantifying urban river-aquifer fluid exchange processes: A multiscale problem[END_REF] highlighted the significance of a range of fluid exchange processes and pathways at the urban river-aquifer interface in spatial and temporal scales. The effect of water level fluctuations on stream-aquifer exchanges and the influence on the fluxes of contaminants were investigated by a few authors [START_REF] Zachara | Persistence of uranium groundwater plumes: Contrasting mechanisms at two doe sites in the groundwater-river interaction zone[END_REF][START_REF] Cloutier | Controls of groundwater floodwave propagation in a gravelly floodplain[END_REF]. Recently, [START_REF] Baratelli | Estimation of stream-aquifer exchanges at regional scale using a distributed model: Sensitivity to in-stream water level fluctuations, riverbed elevation and roughness[END_REF] studied the effects of in-stream water level fluctuations on the stream-aquifer exchanges for a regional hydro-system.

Depending on the hydrological properties of the area under study, a river gains water from or loses water to the aquifer. Such exchanges are frequently modeled considering a linear relationship between the exchange rate and the difference between the river head and the groundwater head [START_REF] Rushton | Possible mechanisms for leakage between aquifers and rivers[END_REF]. More precisely, the exchange flux per unit area of riverbed is modeled using the leakage concept. According to [START_REF] Bear | Hydraulics of groundwater[END_REF], mainly horizontal flows are observed at a distance from the river of about 1.5 times the aquifer thickness. In a steady-state unconfined aquifer, groundwater flows between the drain and the aquifer may be thus modeled under the Dupuit-Forchheimer assumption by considering a horizontal two-dimensional spatial domain and by neglecting vertical flow components presented by [START_REF] Kinzelbach | Groundwater modelling. an introduction with sample programs in basic. developments in water sciences[END_REF].

Besides field investigations which are essential in order to understand the physical environments, numerical models provide significant contributions to bring new insights into various environmental processes. Groundwater modeling tools are of interest to investigate the risk of underground contamination by leakage from storage tanks and to evaluate the performance of drainage systems often designed to prevent groundwater pollution and to provide a modeling framework managing the interaction between river and groundwater. [START_REF] Ruf | Modeling the interaction between groundwater and river flow in an active alpine floodplain ecosystem[END_REF] developed a model using MODFLOW-2000, which is a tool based on finite difference method, to couple the hydrodynamic and the groundwater model and to understand the interaction between hydrological and vegetation dynamics in a floodplain. The importance of the spatial pattern of infiltration and exfiltration rates is highlighted for the spatial distribution of the riverine habitat. [START_REF] Du | Integrated hydraulic modeling of groundwater flow and river-aquifer exchanges in the lower valley of Var River[END_REF] set up a hydraulic model with FE-FLOW software by considering precipitation, evapotranspiration, groundwater extraction and river-aquifer exchanges and eventually could calibrate the transfer rates in the riverbed along the river.

Context of the study

This chapter discusses the impact of a drain set up in an unconfined aquifer from a modeling point of view. The objectives are to evaluate the efficiency of the drainage in the control of water fluxes as a mitigation measure to lessen the risk of pollutant migration into the aquifer. The modeling of the drain-aquifer interaction is carried out here by considering the drain as a river, that is by applying methods developed for the modeling of the classical river-aquifer interaction.

The groundwater flow and leakage interactions between the drain and the aquifer are implemented by two numerical methods: a conforming finite element method (FEM) and a mixed hybrid finite element method (MHFEM). The FEM modeling of such exchanges is classical [START_REF] Kinzelbach | Groundwater modelling. an introduction with sample programs in basic. developments in water sciences[END_REF]. To the best of our knowledge, the implementation of waterbody-aquifer exchanges in a mixed formulation has not yet been published. A comparison between these two discrete approaches is then performed in terms of computed piezometric head, velocity field, and streamlines considering an engineering case study. The considered case (see Figure 2.1a) deals with a regional aquifer located below an industrial platform with chemical storage tanks under rain leaching and water infiltration. A real chemical industrial platform located between a channel and a river is considered as a case study. Two drains are used between the channel and the river in order to catch the polluted water due to the discharges from the storage tank into the aquifer. An impervious barrier may be set up around the area influenced by industrial discharges in order to improve the drainage. Computations are performed with both an ideal homogeneous hydraulic conductivity over the domain and a realistic heterogeneous one.

River-aquifer exchanges occur principally in three ways [START_REF] Winter | The design and use of a hydraulic potentiomanometer for direct measurement of differences in hydraulic head between groundwater and surface water[END_REF]. In the general case, the aquifer-waterbody exchange Q e is modeled by

Q e =        λ (h d -h db) , for h ≤ h db , λ (h d -h), for h db < h ≤ h d , λ (h -h d ), for h d < h, (2.1)
where h d and h db are the water level in the water-body and the level of water-body bottom (L), respectively, and λ is the leakage coefficient (T -1 ) (see Figure 2.1b).

The first case happens only if the drain is disconnected from the aquifer, when an unsaturated zone occurs underneath the drain. When then drain is entirely connected with the underlying aquifer, the flow between the drain and the aquifer is fully saturated and the infiltration rate varies linearly with the changes in the groundwater table. In the drain case (see Figure 2.1) the exchange term Q e represents the exfiltration from the aquifer to the artificial drain only, that is

Q e = λ (h -h d ), for h > h d , (2.
2) The direction of the exchange depends on the difference of the two head levels. The leakage coefficient λ depends on the permeability, the thickness and the leakage area of the drain bed. A major issue in groundwater modeling is thus its parameterization. Classically [START_REF] Kinzelbach | Groundwater modelling. an introduction with sample programs in basic. developments in water sciences[END_REF], the leakage coefficient is chosen as being constant in space and time.

The mathematical equations modeling the behavior of a slightly compressible fluid in an undeformable porous media are based on the Darcy law and the mass continuity equation.

To solve the underground flow problems, the FE and MHFE methods are considered, among the best known numerical methods. The system equations can be discretized to elements, using the Galerkin method and by means of interpolation functions (FEM). The approximate solution h is computed by means of a linear combination:

h(x, y,t) = N ∑ i=1 h i (t)w i (x, y), (2.3)
where N is the number of elements, h i is the unknown head at the node i of the mesh, w i is the basis scalar function related to the node i,

w i (x, y) = M ∑ G=1 Φ G i (x, y), (2.4)
and Φ G i is the unit function related to the node i of the triangle G (Figure 2.2a).

On the other hand, using the MHFEM makes it possible to calculate, contrarily to the velocity field obtained by the conformed finite elements, a new velocity field whose normal component is continuous at the interface of two elements. This treatment greatly improves the quality of the velocity field over the whole domain. Indeed, in this method, Darcy's flux on each triangle G is approximated as:

q G = 3 ∑ j=1 Q G,E j w j , (2.5)
where the vector q G belongs to the lowest order Raviart-Thomas space [START_REF] Raviart | A mixed finite method for the second order elliptic problems: Mathematical aspects of the finite element methods[END_REF], Q G,E j is the water flux over the edge E j of G and w j is the basis vector, for

j = 1, 2, 3 (see Figure 2.2b).
Finally, the streamlines are computed by means of either a classical interpolation method [START_REF] Cordes | Continuous groundwater velocity fields and path lines in linear, bilinear, and trilinear finite elements[END_REF] for FEM and the property for MHFEM explained in Appendix A, Section A.3.

Two dimensional hydrodynamic modeling

Buckingham-Darcy law

In 1856, Darcy formulated an empirical law that relates in a proportional manner the flow

Q (L 3 T -1
) through a porous medium is proportional to the flow section A (L 2 ) and the difference of piezometric heads ∆H = H 1 -H 2 (L) and in an inversely proportional manner to x,y)

( i K K i k j x y h x k j i y w 3 w 1 w 2 E 1 E E 3 2 K B B G K B B B B G G (a) (b) Figure 2.2: (a): FEM unit function Φ G i (x,y) (b): MHFEM basis vectors ω i for i=1,2,3 .
the distance (L) between two points H 1 and H 2 . For a one dimensional flow with an incompressible homogeneous fluid in a saturated porous medium with a hydraulic conductivity

K (L T -1
), Darcy's law can be written as

Q A = K H 1 -H 2 L ,
(2.6) Darcy velocity can then be expressed as a function of the pressure gradient and the gravity [START_REF] Bear | Hydraulics of groundwater[END_REF]:

q = - k µ (∇p + ρg∇z), (2.7)
where • q is the velocity vector of Darcy (L T -1 ).

• k is the intrinsic permeability tensor of the porous medium (L 2 ),

• p is the pressure

(M L -1 T 2 ),
• z is the elevation head * (L),

• ρ is the volumetric mass density (L T -3 ),

• µ is the dynamic viscosity (M L -1 T -1 ).

When the volumetric density is constant, Equation (2.7) in horizontal flows is reduced to: q = -K∇h, (2.8) where h is the pressure head (L) and K is the hydraulic conductivity tensor (L/T ).

Mass continuity equation

The mass continuity equation expresses the principle of conservation of the mass. For an unsaturated flow regime, it is written in the form:

S ∂ h ∂t + ∇.q + Q = 0 in Ω, (2.9)
* The vertical coordinate is defined positive upward.

where S is the storage coefficient which is equal to the effective porosity when the compressibility can be neglected compared to the storage related to the movement of the water table . The flow region is denoted by Ω and the internal source/sink term Q describes injection/pumping wells and groundwater recharge (L/T ). In a steady-state flow, the charge and the water content are constant and Equation (2.9) can be reduced to the following form:

∇.q + Q = 0 in Ω, (2.10)

Hydrodynamic equation

For a saturated flow, under the Dupuit-Forchheimer assumption, combining and integrating Darcy's law (2.8) and the mass continuity Equation (2.9) yield the general formulation

S ∂ h ∂t + ∇.( T∇h) + Q = 0 in Ω,
(2.11) where h is the piezometric head (L), Ω is the flow region considered as a 2D unconfined aquifer, t is the time variable (T ), T is the transmissivity tensor (L 2 /T ). Equation (2.11) is usually subject to Neumann, Dirichlet or mixed boundary conditions (see Subsection 2.4.4). Introducing the exchange expression Q e (L/T )(infiltration/exfiltration) in Equation (2.11), it can be expressed as:

S ∂ h ∂t + ∇.( T∇h) + Q e + Q = 0 in Ω, (2.12)
In an unconfined aquifer, the transmissivity tensor T (L 2 /T ) is equal to M(h) × K where M(h) is the aquifer thickness which varies with the water table elevation h (see Figure 2.1):

M(h) = h -h b (2.13
) where h b is the elevation of the aquifer impervious bottom (L) and the phreatic surface is always above the arbitrary impervious bed. The storage coefficient S is set equal to the effective porosity n e for an unconfined aquifer.

Boundary conditions

Let ∂ Ω D and ∂ Ω N be boundaries of Ω such that ∂ Ω D ∩ ∂ Ω N = / 0 and ∂ Ω D ∪ ∂ Ω N = ∂ Ω. Classical initial and boundary conditions are        h(x, y, 0) = h 0 (x, y) in Ω, h(x, y,t) = h 1 (x, y,t) on ∂ Ω D , K(h(x, y,t) -h b ) ∂ h ∂ n (x, y,t) = Q n on ∂ Ω N .
(2.14)

In other words, a Dirichlet condition and a Neumann condition are imposed on the boundaries ∂ Ω D and ∂ Ω N , respectively.

Discretisation methods

Flow domain Ω is defined as a two-dimensional (2D), and is space-discretized into triangular elements G. Thus, the elements are composed of three edges denominated E i (∀i = 1, 2, 3). Dirichlet (Ω D ) or Neumann (Ω N ) are imposed as boundary conditions.

The two FE and MHFE numerical methods implemented in this study, manage the fluxes in a different manner, (see Figure 2.3), since the exchanges with the drain are taken into account at the nodes (triangle vertices) in the FEM modeling, while they are assigned to the edges in the MHFEM formulation. For the sake of generality, the drains are implemented following Equation (2.1). Note that in the discrete form the leakage coefficient λ (L -1 ) depends on the edge length |E| following:

λ = C λ |E|, (2.15) where C λ (L -1 T -1
) is the leakage coefficient per drain meter.

Finite Element method

Using the Galerkin method, Equation (2.12) can be solved by the FEM [START_REF] Kinzelbach | Groundwater modelling. an introduction with sample programs in basic. developments in water sciences[END_REF] based on linear triangular elements and an implicit finite difference scheme for the time integration.

Substituting Equation (2.3) in Equation (2.12) and with the integration extending over the domain Ω, we obtain a system of N linear differential equations with respect to time in the N unknown piezometric heads h i (t):

Ω ∂ ∂ x (T ∂ ∂ x N ∑ j=1 h j w j ) + ∂ ∂ y (T ∂ ∂ y N ∑ j=1 h j w j ) + S ∂ ∂t ( N ∑ j=1 h j w j ) + Q + w i dxdy = 0, ∀i = 1, ..., N.
(2.16) Using Green's theorem, the first two terms can be integrated by parts. This yields:

Ω -T ∂ ∂ x ( N ∑ j=1 h j w j ) ∂ w i ∂ x -T ∂ ∂ y ( N ∑ j=1 h j w j ) ∂ w i ∂ y + S ∂ ∂t ( N ∑ j=1 ∂ h j ∂t w j )w i + Qw i dxdy + ∂ Ω q n w i ds = 0, ∀i = 1, ..., N,
(2.17) where ds is the differential of path-length along the boundary and q n is the flow across the boundary per unit length of it, expressed as:

q n = N ∑ j=1 T ( ∂ w j ∂ x + ∂ w j ∂ y ). n ∂ Ω h j ,
with n ∂ Ω being the unit vector normal to the boundary ∂ Ω. Equation 2.17 can be written in the form of (2.18) where

N ∑ j=1 P i j h j + N ∑ j=1 R i j ∂ h j ∂t -F i = 0 or Ph + R dh dt -F = 0,
P i j = Ω T ∂ w i ∂ x ∂ w j ∂ x + T ∂ w i ∂ y ∂ w j ∂ y dxdy, R i j = Ω Sw i w j dxdy, F i = Ω (Q e + Q)w i dxdy + ∂ Ω Q n w i ds.
Equation (2.18) constitutes a time-dependent system of ordinary differential equations. A difference method is used to solve it in (0, T ):

(P)(h(t ′ )) + (R) (h(t + ∆t) -h(t)) ∆t = (F), (2.19)
After rearranging, the fully implicit scheme (t ′ = t + ∆t) yields the equation system

(R) ∆t + (P) (h(t + ∆t)) = (R) ∆t (h(t)) + (F).
(2.20)

Mixed Hybrid Finite Element method

The flow-transport coupling implemented in numerical simulations in porous media requires a good approximation of the velocity. The mixed finite element method provides an answer to this type of problem by approaching the piezometric head and the velocity simultaneously. The mixed approach was presented for the first time in the field of potential flows by [START_REF] Meissner | A mixed finite element model for use in potential flow problems[END_REF]. Equation (2.11) being discretized by means of a mixed approach, computes the state variable and its gradient in a simultaneous manner. This method leads to a system of equations whose unknowns are the mean piezometric heads per mesh and the fluxes through the edges (for mixed elements of order one). The disadvantage of this approach is that the matrix, symmetric undefined positive, associated with the system is difficult to solve. The mixed hybrid formulation, as it is presented by [START_REF] Arnold | Mixed and non conforming finite element methods: Implementation,postprocessing and error estimates[END_REF] and [START_REF] Chavent | Mathematical models and finite elements for reservoir simulation[END_REF], makes it possible to eliminate this disadvantage by choosing the unknowns of the system as the mean piezometric head for each edge of elements. The matrix associated with the mixed hybrid system is positively symmetric as it is defined in the classical methods. The resolution of this type of system is very efficient with the preconditioned conjugate gradient method.

In this subsection, the mixed hybrid approximation applied to the resolution of the flow in saturated porous medium is presented.

Continuity of pressures and fluxes and boundary conditions

Considering T h G,E i as the mean piezometric head of each edge i which belongs to the element G, the continuity of the piezometric head at edges, as well as Dirichlet boundary conditions of is expressed by

T h G,E i = T h G ′ ,E i , (2.
21) The continuity of fluxes between two adjacent elements is given by

Q G,E i + Q G ′ ,E i = 0 , (2.22
) where E i is the common edge between two elements, G and G ′ . This equation is valid for all the interior edges E i (∀i = 1, 2, 3) of the domain Ω.

Dirichlet boundary conditions: Dirichlet boundary conditions, provided by a piezometric head equality on the considered edges, are represented as:

T h G,E = T h D,E , ∀E ⊂ ∂ Ω D .
(2.23)

Neumann
Q G,E = Q N,E , ∀E ⊂ ∂ Ω N .
(2.24)

Variational formulation of Darcy's law

Darcy's law can be written as: K-1 . q = -∇h.

(2.25)

By multiplying each member of Equation (2.25) by a test function s, integrating on the element G and using the formula of Green, it yields to the following equation:

G (K -1 • q) • s = - G ∇ h• s = G h ∇ • s - ∂ G h s • n G .
(2.26)

By replacing K -1 , h and q by their approximation on the element G, and by taking successively the basic functions w i as test function, the mixed formulation is obtained as:

G (T G -1 . q G ). w i = h G G ∇. w i - n f ∑ j=1 T h G,E j E j w i . n G, j , ∀i = 1, ..., n f (2.27)
where n f is the number of the edges of the mesh.

On each element the vector function q G has the following properties [START_REF] Chavent | A unified physical presentation of mixed, mixed-hybrid finite elements and usual finite differences for the determination of velocities in waterflow problems[END_REF]]:

• ∇ q G is constant over the element G.

• q G n G,E i is constant over the edge E i of the triangle,∀i = 1, 2, 3, where n G,E i is the normal unit vector exterior to the edge E i . • q G is perfectly determined by knowing the flux through the edges.

Darcy's flux integrated over the aquifer thickness, namely q = -T.∇h, is approximated on each triangle G by Equation (2.5). The basis vector w j , see Figure 2.4, verifies

E j w i n G,E j = δ i, j , ∀i = 1, .., 3 (2.28)
where δ i, j is the Kronecker symbol and n G,E j is the exterior normal unit vector to E j . Functions w i (i = 1, 2, 3) correspond to a vector q G having a unitary flux through the edge E i and null flux through the other edges. Thus, with the MHFEM, the normal component q G is continuous from G to the adjacent element G ′ . On each element, the approximation is such that ∇ q G is constant over the element G and q. n G,E i is constant over the edges E i of the triangle.

In the media with an isotropic hydraulic conductivity distribution, using the fact that on the one hand ∇. w i = 1 and considering Equation (2.5) on the other hand, Equation (2.27) is rewritten as

n f ∑ j=1 Q G,E j G (T -1 • w j ) • w i = h G -T h G,E i , ∀ i = 1, .., n f (2.29)
The n f × n f symmetric elementary matrix B G associated with the element G is defined as:

B G = [B G,i j ] , where B G,i j = (T -1 • w j ) • w i
(2.30) Thus, the variational formulation of Darcy's law is:

n f ∑ j=1 Q G,E j B G,i j = h G -T h G,E i , ∀ i = 1, .., n f (2.31)
The elementary matrix equation is deduced as:

B G Q G = h G DIV T G -T h G,E i , ∀ i = 1, .., n f (2.32)
where

DIV T G =        1 1 . . 1        , Q G =        Q G,E 1 Q G,E 2 . . Q G,E n f        , T h G =        T h G,E 1 T h G,E 2 . . T h G,E n f       
The matrix B G being invertible and using the continuity of the fluxes (2.21), Equation (2.32) can be expressed as

Q G,E = h G α G,E -∑ E ′ ⊂∂ G B -1 G,EE ′ • T h E ′ , where α G,E = ∑ E ′ ⊂∂ G B -1 G,EE ′ (2.33) and ∑ E ′ ⊂∂ G
is the sum of the edges E ′ belonging to the element G.

Matrix form of the continuity of flux

Substituting Equation (2.33) in Equation (2.22), the expression for each edge inside the domain Ω is obtained as following:

a G α G,E h G + a G ′ α G ′ ,E h G ′ -a G ∑ E ′ ⊂ ∂ G B -1 G,EE ′ T h E ′ -a G ′ ∑ E ′ ⊂ ∂ G ′ B -1 G ′ ,EE ′ T h E ′ = 0.
(2.34) Using the relation (2.33) in Equation (2.24) which expresses the Neumann boundary conditions, the following equation is obtained as:

a G α G,E h G -a G ∑ E ′ ⊂ ∂ G B -1 G,EE ′ T h E ′ -Q N,E = 0, ∀ E ⊂ ∂ Ω N (2.35)
Equations (2.34) and (2.35) can be expressed by a matrix relation expressing the continuity of the normal component of the velocity vector as Dh -RTh -V -I = 0 (2.36) where Th is the mean piezometric head of each edge which is not known, h is the mean piezometric head of each element and the matrices are defined as:

D = [D EG ] n f r,nm such that D EG = α G,E if E ⊂ ∂ G, 0 if E ⊂ ∂ G, R = [R EE ′ ] n f r,n f r such that R EE ′ = ∑ E, E ′ ⊂∂ G B -1 G,EE ′ V = [V E ] n f r such that V E = Q N,E ∀ E ⊂ ∂ Ω N 0 ∀ E ⊂ ∂ Ω N I = [I E ] n f r such that I E = ∑ G⊃E ∑ E ′ ⊂∂ G B -1 G,EE ′ T h E ′ ∀ E ′ ⊂ ∂ Ω D where ∑ E, E ′ ⊂∂ G
is the sum on the elements G which contains both the edges E and E ′ , n f r is the number of edges over the domain Ω where the pressure head has not been imposed and nm is the number of elements in Ω.

Discretization of the continuity equation

Equation of the continuity (2.26) can be written as the following general form:

S(x, y) ∂ h ∂t + ∇ • q = f (x, y,t) on Ω, for t ∈ ] 0,T [ (2.37)
The domain Ω is discretized into the elements G. Thus, Equation (2.37) is valid for each element G:

G S(x, y) ∂ h ∂t + G ∇ • q = G f (x, y,t) (2.38)
The temporal space is discretized in N time intervals. The notation T h n E means that T h E is considered at time n∆t.

Replacing h and q by their approximations on the element G and using temporal discretization, the Equation (2.38) is obtained as:

G S G h n+1 G -h n G ∆t + G ∇ • q n+1 G = G f n+1 (2.39)
where h G is constant on the element G and

• ∇. q G = 1 |G| is constant on G, |G| is the surface of the element.
• S G is an approximation of S(x, y) constant on G.

Considering F n+1 G as an approximation of G f n+1 being constant on the element G and by integrating the equation (2.39) on the element G, the following equation is obtained:

|G| S G h n+1 G -h n G ∆t + ∑ E ⊂ ∂ G Q n+1 G,E -F n+1 G = 0 (2.40)
The elementary equation 2.33 gives directly:

∑ E⊂ ∂ G Q n+1 G,E = (α G h n+1 G -∑ E ′ ⊂ ∂ G α G,E ′ T h n+1 G,E ′ ) where α G = ∑ E⊂ ∂ G α G,E (2.41)
Substituting the expression (2.41) in Equation (2.40) yields:

|G| S G h n+1 G -h n G ∆t + a G α G h n+1 G -∑ E ′ ⊂ ∂ G α G,E ′ T h n+1 E ′ -F n+1 G =0 ∀ G ∈ Ω (2.

42) Considering the intermediate variables as following:

γ G = α G ∆ t S G |G| and β G = γ G 1 + γ G ∀ G ∈ Ω
the equation can be expressed as:

h n+1 G -(1-β G ) h n G -β G ∑ E ′ ⊂∂ G α G,E ′ α G T h n+1 E ′ - β G F n+1 G α G = 0. (2.43)
Writing the equation (2.43) in matrix form, a second equation of state (mass balance equation) is obtained:

h n+1 -M h n -NTh n+1 -K F n+1 -H n+1 = 0 (2.44)
where F is a sink-source vector and the matrices defined as below:

M = [M GG ′ ] nm,nm where M GG ′ = 1 -β G if G = G ′ 0 if G = G ′ K = [K GG ′ ] nm, nm where K GG ′ = β G α G if G = G 0 if G = G ′ N = [N GE ] nm,n f r where N GE = β G α G,E α G if E ⊂ ∂ G 0 if E ⊂ ∂ G H = [H G ] nm where H G = ∑ E ⊂ (∂ G ∩ ∂ Ω D ) β G α G,E α G • Th E Combining the Equation (2.36) and (2.44) yields the MHFEM system (R -DN -R e )Th n+1 = DSh n + DKF n+1 + DH n+1 -V n+1 -I n+1 -Q e (2.45)
where Th is the vector of the piezometric head for non-Dirichlet edges, h is the vector of piezometric head for the triangles.

Drain distribution: Matrices related to the drain contributions are computed as follows

R e = [R e EE ′ ] ∈ IR N e\∂ Ω D ×N e\∂ Ω D such that R e EE ′ = -λ , 0, if E = E ′ and Th n E > h f ,E ∀E ⊂ drain, else where,
and

Q e = [Q e E ] ∈ IR N e\∂ Ω D such that Q e E =    λ h r,E , if Th n E > h f ,n ∀E ⊂ drain, λ (h r,E -h f ,E ), if Th n E < h f ,n ∀E ⊂ drain, 0, ∀E ⊂ drain,
where λ is the leakage coefficient, h r,E is the level of the drain water free surface of edge E and h f ,E is the drain bed elevation of edge E.

Numerical results and discussions

Case study

The considered case corresponds to a real site in France (Figure 2.5(a)). The exact location is not given for confidentiality reasons. It is based on a real chemical industrial platform with an area of 1 km 2 located between a navigation channel to the west and a river to the east. The water flows from the channel (height: 221.5 m) to the river (height: 217 m).

The underlying aquifer is subject to drainage for environmental protection. The domain (Figure 2.5(a)), comprises two drains (red lines) laid out to a elevation of 217.1 m to 217.4 m in a north-south alignment to catch and to evacuate the groundwater flowing below the possible contamination sources. Four cases are investigated by considering the domain as homogeneous or heterogeneous, regarding the hydraulic conductivity distribution, for each an impervious barrier either is accounted in or not. In the first homogeneous case, in which

Approximated piezometric heads and streamlines

Computed piezometric head and streamlines are displayed in Figures 2.6, 2.7, 2.8 and 2.9. Piezometric head isolines are distributed between the channel (upstream head of 221.5 m) and the river (downstream head of 217 m).

Homogeneous transmissivity cases: In the homogeneous cases, the groundwater smoothly flows from the channel to the river both in the southern (Y<4000 m) and northern zones (Y>5500 m). Considering the case without barrier, in the middle zone (4000 m<Y<5500 m), some (see Figure 2.6, streamlines 5, 6, 7 and 8) but not all of the fluid particles are caught by the drains. This effect is corroborated by the drawdown of the piezometric head to the drains. It can be mentioned that in the FEM the streamline number 9 is more deviated in comparison with the one in the MHFEM. Contrarily to the case without barrier, all the particles (see Figure 2.7, streamlines 5, 6, 7, 8 and 9) delimited by the impervious barrier is captured by the drains in the case with the barrier. Although the two discretization methods yield the same flow pattern, the drawdown is greater in the FEM case. This is an expected result since the FEM and MHFEM formulations manage exchange fluxes between the aquifer and the drain in a different manner (see Subsection 2.5). In the homogeneous case with the barrier (see Figure 2.7) contrarily to the case without barrier (see Figure 2.6), inwhich the streamline number 9 is not captured by the drains, all the particles in the middle zone are caught.

Table 2.1 reports fluxes computed at the channel, the river and the drain boundaries as well as the Mass Balance Error (MBE) estimated as the sum of the inlet flow rate at the channel and the outlet flow rate at the river and the drains. Null values for the MBE prove that the two codes are conservative. Table 2.1 moreover shows that the drains catch much more water in the FEM cases with homogeneous hydraulic conductivity distribution. This is less noticeable in the other cases due to the heterogeneity.

Heterogeneous transmissivity cases:

A similar behavior, in Figures 2.8 and 2.9, may be observed in the northern and southern zones in the heterogeneous transmissivity cases, whatever the discretization method. Streamlines are slightly different in the middle zone. Near to the drain, in the case where the streamlines are confined by the almost impervious barrier, all the particles are caught by the drains (see Figure 2.9, streamlines 5, 6, 7, 8 and 9) which is not the case for the one with no barrier (see Figure 2.8, streamline number 9 is not captured). This confirms that the impervious barrier increases the efficiency of drainage. Flux values at the channel, river and drain boundaries are reported in Table 2

Approximated velocity fields

Velocity fields computed in the homogeneous transmissivity case without barrier are plotted in Figure 2.10(a). As expected from the difference in the exchange modelings (Subsection 2.5), velocity vectors are larger close to the northern drain in the FEM because more water is caught. One also observes a difference in the orientation of the velocity vectors at the east of the drain indicating that the FEM is able to catch water in this zone, contrarily to the MHFEM.

In the case with the barrier carried out by the MHFEM (see Figure 2.10(d)), the barrier caused a significant change in the flow velocity field after the drains in which, contrarily to the case without barrier (see Figure 2.10(b)), the orientations of the velocities are toward the drains. This confirms the efficiency of the impervious barrier to enforce the drainage. The velocities in the homogeneous case with the barrier are greater in comparison to the one without barrier due to the difference in the hydraulic conductivity explained in Subsection 2.6.1.

Comparing Figures 2.10 and 2.11, velocities are smaller in the heterogeneous cases. The main reason to explain this behavior is that the hydraulic conductivity is subject to high variations in the domain (including small values). This heterogeneity has an important impact on the flow as it reduces its velocity in a significant manner. Setting the barrier reinforces very slightly the prohibition of the flow from going down to the river in the heterogeneous cases, contrarily to the homogeneous ones, in which, it improves significantly the drainage and the velocity field in the vicinity of the drains.

Conclusions

The groundwater flow and drain-aquifer interactions are implemented in a conforming FEM and a MHFEM considering the drain as a river. The drain-aquifer exchange flux is modeled by means of the leakage concept and Darcy's law where the leakage coefficient depends on the permeability, the thickness, and the leakage area of the drain bed. In this chapter, the leakage coefficient per meter of drain is assumed to be constant. Numerical experiments have been carried out using ideal homogeneous hydraulic conductivity over the domain, and actual heterogeneous case study. Fluxes computed at the channel, the river, and the drain boundaries have similar behavior but distinct values due to the difference in the modeling of water-body exchanges. In all the cases, the drain catches more water using FEM. Finally, the importance of the role of the impervious barrier is highlighted.

Introduction

The potential variations of any model related to the value of parameters or assumptions may be investigated in order to have minimal errors as well as optimal solutions. In the presence of uncertainty, the sensitivity analysis techniques help address these issues by understanding relationships between input and output variables and thus, the robustness of the results of a model can be verified. Sensitivity analysis may be used for various reasons, such as decision-making or development of recommendations, communication, increasing understanding or quantification of system, and model development. In model development, it can be used for the purposes of model validation or accuracy, simplification, calibration, and even to identify important parameter for further studies [START_REF] Pannell | Sensitivity analysis of normative economic models: theoretical framework and practical strategies[END_REF]. The importance of sensitivity analysis can be highlighted for calibrating models with large number of parameters which can facilitate the calibration by providing the priority of parameters identifications.

In addition, for a deterministic model, a possible way to tackle the problems to evaluate the influence of the uncertainty of input data on the uncertainty of output data is to impose a perturbation of the input values and to observe the resulting perturbations in the output variables. Contrarily to the deterministic sensitivity analysis, in which model input is specified as multiple point estimates and are varied manually, in the probabilistic one, model inputs are specified as a distribution and thus, multiple parameters are varied simultaneously and the model may be run many times [START_REF] Scheidegger | Deterministic and statistical characterization of porous media and computational methods of analysis[END_REF].

The hydraulic conductivity tensor, the effective porosity and the leakage coefficient are modeling parameters, the values of which are subject to uncertainties which may be considered when choosing approaches to study groundwater-surface water interactions. General sources of these uncertainties may be errors and approximations in input data measurement, parameter values, model structure and model solution algorithms. A classical issue in environmental sciences is the evaluation of the sensitivity of the model outputs to the input parameters [START_REF] Pianosi | Sensitivity analysis of environmental models: A systematic review with practical workflow[END_REF]. A sensitivity analysis is supposed to determine the change in model output values that results from modest variations in model input values such as parameters which may have a physical or a geometrical meaning, or some initial or boundary conditions, for instance.

In the literature, various studies have investigated the sensitivity of waterbody-aquifer exchange to modeling parameters as a key to understand and to evaluate the ecological structure of groundwater-surface water interactions. As a measure of the model response to perturbations, sensitivity analysis is often acknowledged as an important element in the modeling process [START_REF] Helton | Sensitivity analysis of model output: Samo 2004[END_REF]. [START_REF] Du | Integrated hydraulic modeling of groundwater flow and river-aquifer exchanges in the lower valley of Var River[END_REF] performed the sensitivity of the groundwater level with respect to the values of infiltration and exfiltration. This study showed a significant influence of the infiltration when the groundwater level is low and the predominant impact of exfiltration after the head level peaks. Furthermore, flux estimates based on the Darcy equation are known to be inaccurately approximated and measured [START_REF] Kalbus | Measuring methods for groundwater-surface water interactions: a review[END_REF]. [START_REF] Sheets | Ground-water modeling of pumping effects near regional ground-water divides and river/aquifer systems-results and implications of numerical experiments[END_REF] examined the sensitivity of hydraulic heads to changes in hydraulic conductivity, riverbed hydraulic conductance, areal recharge and the general head boundary conductance as well as pumping rate of wells in a steady-state model. Among the statistical ones, a Monte Carlo method was used to study the variability of the river-aquifer seepage flow to the spatial variability in the aquifer-saturated hydraulic conductivity [START_REF] Bruen | Sensitivity of stream-aquifer seepage to spatial variability of the saturated hydraulic conductivity of the aquifer[END_REF]. Furthermore, for model calibration purposes, sensitivities were computed from the differentiation of groundwater flow equations with respect to the hydraulic conductivity, the aquifer thickness, or the aquifer recharge rate [START_REF] Mazzilli | Sensitivity analysis of two-dimensional steady-state aquifer flow equations. implications for groundwater flow model calibration and validation[END_REF].

In the context of drain-aquifer exchanges of Chapter 2, a pending question is: 'how much are the simulation results sensitive to the value of the leakage coefficient?' Indeed, neither an empirical formula exists to calculate this coefficient, nor a device to measure it directly. The uncertainty on the leakage coefficient in this specific context has thus to be evaluated. In this chapter and in the paper authored by [START_REF] Moezzibadi | Sensitivity of groundwater flow with respect to the drain-aquifer leakage coefficient[END_REF], sensitivities of head, velocity field, and streamlines are computed with respect to the leakage coefficient to provide a new insight into the groundwater modeling, notably to grasp the differences between FEM and MHFEM models in an accurate manner. From a technical implementation point of view, the automatic differentiation (AD) [START_REF] Griewank | Evaluating Derivatives[END_REF][START_REF] Hascoet | The tapenade automatic differentiation tool: Principles, model, and specification[END_REF] of the numerical code is performed in the so-called tangent linear (TL) model [START_REF] Elizondo | Automatic versus manual model differentiation to compute sensitivities and solve non-linear inverse problems[END_REF][START_REF] Charpentier | Local sensitivity analysis of a numerical model of volcanic plinian columns through automatic differentiation[END_REF]. This chapter is devoted to sensitivity methods and presents AD basics provided in Subsections 3.2 sensitivity analysis, 3.3 the AD basics, and finally, 3.3.1 AD application to FEM and MHFEM numerical codes (Subsection 3.4).

Sensitivity analysis

To emphasize the generality of the proposed method, let us consider the general model M (P) = x, the inputs and output of which are the modeling parameters P and the state variable x, respectively.

Given a model M (P), a sensitivity analysis measures the impact of any small perturbation δ P of its input parameter set P on the model response x = M (P). Some qualitative information may be deduced from repeated evaluations of the model M (P + δ P) with different values of δ P . However, when possible, sensitivity analysis based on derivative computations should be preferred as it generally provides more valuable qualitative and quantitative information on the model [START_REF] Elizondo | Automatic versus manual model differentiation to compute sensitivities and solve non-linear inverse problems[END_REF][START_REF] Charpentier | Local sensitivity analysis of a numerical model of volcanic plinian columns through automatic differentiation[END_REF].

Assuming the model is differentiable, its sensitivity s M (P 0 , δ P) with respect to P may be evaluated with the parameter set P = P 0 in the direction of perturbation δ P following

s M (P 0 , δ P) = ∂ M ∂ P (P 0 ).δ P = |M (P + ωδ P ) -M (P)| ω|δ P | , (3.1)
where ω is a small parameter. Within simulation codes, sensitivity computations may be carried out by means of a finite difference method to obtain approximate derivatives, or a differentiation of the discrete equations of the model, or a differentiation of the numerical code implementing the discrete equations. The interest is three-fold. First, such an analysis provides a qualitative and quantitative insight into the physical behavior of the model and how it is impacted by a change on one of its inputs. Second, it allows the assessment of the uncertainties in the parameters of the model and how the outputs may be altered by some parameter misfit. Third, accurate sensitivity computations are a necessary prerequisite for parameter identification methods involving a gradient computation.

The correctness of the sensitivity process is checked with a classical Taylor test. This compares the sensitivity result obtained with AD to first order finite difference approximations following

r ω (M (P)) = s M (P 0 , δ P) M D (P, δ P ) . (3.2)
Theoretically, this ratio should tend linearly towards 1 as ω tends to zero. In practice, the subtraction of too close floating-point numbers yields a cancellation error that dominates finite difference truncation error for smaller ω.

Automatic differentiation

AD [START_REF] Griewank | Evaluating Derivatives[END_REF][START_REF] Naumann | The Art of Differentiating Computer Programs[END_REF] is a set of techniques designed to augment computer codes with derivative computations.

Within AD, a computer code may be viewed as a sequence of statements run in a prescribed order. Given the input data, the execution flow includes information about how the code starts, the actual order of execution of the statements and how it terminates. This execution flow is no more than a large composition of arithmetic operations and intrinsic functions.

The differentiation in Tangent Linear mode (TL) consists in the differentiation of this compound function by applying the chain rule and classical rules such as "the derivative of a sum is the sum of the derivatives". AD may be applied to very large codes [START_REF] Charpentier | The mesodif package for gradient computations with the atmospheric model meso-nh[END_REF]. Note that the differentiation in the so-called adjoint mode, frequently used for identification purposes, is beyond the scope of this study.

AD relies on two kind of software. On the one hand, source transformation tools such as Tapenade [START_REF] Hascoet | The tapenade automatic differentiation tool: Principles, model, and specification[END_REF] or Adifor [START_REF] Bischof | Adifor-generating derivative codes from fortran programs[END_REF] are able to generate source codes containing derivative statements. On the other hand, operator overloading libraries such as Adol-C [START_REF] Griewank | Algorithm 755: Adol-c: A package for the automatic differentiation of algorithms written in c/c++[END_REF], Rapsodia [START_REF] Charpentier | Fast higher-order derivative tensors with rapsodia[END_REF] and Arbogast [START_REF] Charpentier | Higher order automatic differentiation for special functions with modular c[END_REF]] may be used at compile time to propagate derivatives at runtime. In this study, we use the Tapenade software [START_REF] Hascoet | The tapenade automatic differentiation tool: Principles, model, and specification[END_REF]. For small codes (less than 3000 lines), the web interface gives access to the "Tapenade On-line AD Engine". Larger codes require to download the Tapenade software for a local installation. Differentiation may be then performed using the graphical user interface. The user provides the source of its code, the name of the top routine to be differentiated, the "independent" input variables (modeling parameters of interest) and, optionally, the "dependent" output variables (state vector, for instance). When the provided information is coherent, Tapenade generates a tangent linear code differentiated with respect to the independent variables. At runtime, this linear code propagates one (or several) direction(s) of perturbation to compute the dependent variables and their sensitivity(ies).

The key-point of a successful AD is to be confident in the tool. From a practical point of view, this signifies that "differentiation errors" arising when using an AD tool are generally user's ones.

Given a code and differentiation instructions, Tapenade first proceeds to a syntactic analysis of the code structure and statements (parser phase). The parser checks the correctness of the statements. It notably finds out type incoherences, potential source of errors (equality test on real numbers), and I/O of active variables (independent and dependent), for instance. If these exist, parser errors and warnings have to be cured (or at least understood) before the actual differentiation phase. Very little effort is needed to obtain a differentiated code with the Tapenade AD tool (a few seconds to a few hours, depending on the quantity of "errors" in the source code).

Minimal effort AD strategy for existing codes

In this subsection a minimal effort AD strategy to obtain a tangent linear code for sensitivity calculation is proposed. The present study takes advantage of two existing groundwater flow codes, namely GW_FEM and GW_MHFEM, written in Fortran. Their general implementation scheme is similar (see Figure 3.1) that is: -Read-input: reading the case study data (modeling parameters P, boundary conditions, mesh, simulation duration, ...),-(MH)FEM solver: solving the FEM or MHFEM discrete problem, and -Write-output: writing output data x (piezometric head, velocity field, streamlines).

The differentiation may be performed with respect to the parameter set P, here to the two leakage coefficients (one per drain, see Figure 2.1a). As described in Figure 3.1, the original user code GW_FEM may be automatically differentiated with a very few effort. In the preprocessing stage (blue statements), the program GW_FEM becomes the top subroutine to be differentiated, the arguments of which should relate the independent variable P to the dependent variable x. One can notice that the AD tool cannot consider P as an actual independent variable since its value is read in the original program. The solution is to consider a null input vectorial variable deltaP in the argument of the top routine and to add it to the parameter P after reading. This has no effect on the simulation results. AD is carried out with respect to the independent variable deltaP.

In the tangent linear mode, the source to source AD tool, Tapenade, generates the differentiated routine GW_FEM_D (green statements), the arguments of which are deltaP and x as well as their derivatives deltaP_D and x_D. One can observe that AD differentiates computing routines such as FEM_solver, while it does not differentiate I/O routines, namely read_input(P) and write_output(x). A write statement for the tangent linear variable xd should be thus added to the main routine afterwards (magenta statements to save the sensitivities in a file).

The tangent linear code is evaluated at point (P,deltaP) where P is read in a file and deltaP=0 in the direction of perturbation deltaP_D (usually a canonical vector of the This tool generates the differentiated codes (see Figure 3.2). Some details about the codes (number of routines, number of lines) together with run-times are provided in Table 3.1 for the original user codes and the differentiated ones. It can be observed that the FEM code is faster than the MHFEM code because the number of unknowns is greater in the MHFEM case. The ratios between the TL codes and their respective original code are lower than the theoretical bound of 4 related to the differentiation of the division operator [START_REF] Morgenstern | How to compute fast a function and all its derivatives: A variation on the theorem of baurstrassen[END_REF]. This is an excellent result since the TL codes comprise the original code statements for trajectory computations.

Recycling the linear solver

One can notice that the MHFEM solver comprises an iterative bi-conjugate gradient method for the solution of the general linear system Ax = b, the differentiation of which is more efficient and accurate by calculating the tangent linear variable ẋ following

A ẋ = ḃ -Ȧx.
At this stage, we already know A and x by running direct routine. The variables Ȧ and ḃ are the output variables of the Tapenade tangent linear routines. Next, we call the linear solver in order to obtain ẋ.

Numerical results and sensitivity computations

The sensitivity of the piezometric head, velocity fields and streamlines with respect to the leakage coefficient applied to a real case-study provides a new insight into groundwater waterbody exchanges. It moreover allows for a careful analysis of the differences between the FEM and MHFEM modelings.

Computations are performed, considering the same case study mentioned in Section 2.6.1, with the leakage coefficient C λ = 6 m -1 d -1 for the two drains. Sensitivities are computed with respect to a perturbation of 50% of the leakage coefficient of the northern drain, that is δ C λ = 3 m -1 d -1 , for a deeper understanding of the underground flow.

Once the Taylor test is validated, piezometric head, velocity field, streamlines and their sensitivities to the leakage coefficient value are computed in the homogeneous and the heterogeneous transmissivity cases, without and with the impervious barrier, using the GW_FEM and GW_MHFEM codes.

Validation of the differentiated codes

In the present case, the GW_FEM and GW_MHFEM codes are differentiated with respect to the leakage coefficient of the northern drain. Taylor tests are carried out by considering the position of the particle at the end of the second streamline from the top of Figure 2.6. One observes the expected behaviors for the two ratios. The tangent linear codes are correct. Taylor tests were performed on the head and the velocity field too.

Piezometric head sensitivity

Sensitivity isolines accurately describe how and where the drain catches water. Figure 3.3 displays head sensitivity isolines computed in the homogeneous transmissivity case without barrier. As expected, a positive perturbation of 50% in the exchange capacity between the aquifer and the north drain results in a negative piezometric head sensitivity indicating a stronger drawdown.

Although the drain catches less water using the MHFEM code, it can be seen that the piezometric head computed by MHFEM is more sensitive near to the drain since the closed domain delineated by the isoline -0.04 m is larger in the MHFEM. This is confirmed in Table 3.3 by ratios between sensitivities and computed fluxes, FEM_D/FEM and MHFEM_D/MHFEM, of 19.57% and 25.20%, respectively. In the homogeneous case, in presence of the impervious barrier, the difference of the and3.6). Again, the role of the impervious barrier in limiting the sensitivities inside the protected area is obviously noticeable. Even trends are identical, it can be noticed in Table 3.4 that the ratio between sensitivities and computed fluxes computed by MHFEM at the drain (6.1%) is twice larger than the ones computed by FEM (2.7%) in the case without barrier (9.0% for the MHFEM and 4.4% for the FEM in the case with the barrier).

Velocity field and sensitivity

Velocity fields and their sensitivity computed in the homogeneous and heterogeneous transmissivity cases are shown in Figures 3.7 to 3.14.

Considering the homogeneous case and comparing the numerical methods (see Figures 3.7 and 3.8, 3.9 and 3.10), it can be seen that the sensitivities of the velocity fields are almost similar along the drain, but differ at the drain ends. This is in agreement with the two drain modelings. Particularly, considering Figure 3.8, contrarily to the velocity vectors after the drains are toward the river, the sensitivities in the same region that are mostly greater than the velocities of direct code, are toward the drains. Thus, it is inferred that perturbation of 50% of the leakage coefficient and increasing the drainage changes totally the direction of the velocity field from the river to the drains.

In the heterogeneous case, velocities are smaller (see Figures 3.11,3.12,3.13 and 3.14) and sensitivities are even smaller. Important differences in the vector lengths of the two velocity fields are noticeable at the northern end of the drain because the MH-FEM_D/MHFEM ratio is here twice larger than FEM_D/FEM ratio in Table 3.3.

Streamlines sensitivity

The streamlines computed using C λ = 9 m -1 d -1 (blue lines) and approximated streamlines (green lines) are compared. They are built from heads and velocity fields and their sensitivity computed with the TL codes with parameter and perturbation values

(C λ , δ C λ ) = (6 m -1 d -1 ,3 m -1 d -1 ).
In all cases, the approximated streamlines are close to the streamlines computed with C λ = 9 m -1 d -1 . In the homogeneous transmissitivity cases, one can notice some erratic behavior on the streamline number 5 of Figures 3.15 and 3.16. Although originated from the south-west of the drains, this particle is impacted by the north drain while running to the river. The erratic behavior occurs at a piezometric head maximum. This behavior is probably also favored by the small size of the element in this part of the computational domain. The particle hesitates between the drain and the river, finally reaching the river.

(C λ , δ C λ )=(6 m -1 d -1 ,3 m -1 d -1
).Homogeneous transmissivity case.

C λ , δ C λ )=(6 m -1 d -1 ,3 m -1 d -1
). Homogeneous transmissivity, case with the barrier. .17: Piezometric head and streamlines (blue lines) computed with a leakage coefficient of 9 m -1 d -1 and approximated streamlines (green lines) computed with

(C λ , δ C λ )=(6 m -1 d -1 ,3 m -1 d -1
). Heterogeneous transmissivity, case without barrier.

Conclusions

A new method has been proposed here to elaborate a sensitivity analysis which can provide a new insight regarding the behavior of two numerical FE methods implemented for river-aquifer interaction.

Numerical experiments have been carried out using an ideal homogeneous hydraulic conductivity over the domain, and an actual heterogeneous case study to measure the sensitivity of the computed piezometric head, velocity field, and streamlines to a perturbation of the leakage coefficient value. The proposed implementation of sensitivity analysis through AD, which needs a minimal effort, is described with details as it may serve for other applications.

Taking in account the direct results in chapter 2, discrepancies between the FEM solutions and the MHFEM solutions are difficult to evaluate from model computations. Hence, the accurate sensitivity analysis we propose, provides much more information on piezometric head and velocity field, revealing the influence of a perturbation of the leakage coefficient in a clear manner, whatever the case study. The role of the impervious barrier is clearly stated since sensitivities are 0 outside the protected area. To the best of our knowledge, the streamline reconstruction results are new in groundwater flow modeling. Although some erratic behavior is observed, the first-order approximated streamlines are very close to the computed streamlines.

It is also shown that the dependency of the state variables on leakage coefficient in MHFEM is more significant than FEM whereas the drains in FEM catches more water

(C λ , δ C λ )=(6 m -1 d -1 ,3 m -1 d -1
). Heterogeneous transmissivity, case with the barrier. than ones in MHFEM.

Consequently, model results indicate that surface-groundwater interactions may need sophisticated approachs to be able estimate the leakage coefficient which has been proven to be an important component of simulations for such exchanges.

In comparing with other sensitivity analysis methods, AD may allow us to approximate the sensitivities even though if the direct solution does not converge and may avoid simulations which are computationally expensive.

Part II

Modeling and parameter identification of flow in a variably saturated porous medium

French Abstract

Les rejets urbains par de temps de pluie sont épisodiques, et peuvent avoir des effets de nature chronique liés à la répétition des phénomènes. Ils peuvent modifier les différentes composantes des milieux récepteurs : altérations physiques des écoulements et de la morphodynamique, altérations chimiques par apports de matières en suspension, de fertilisants ou de micropolluants, altérations biologiques par sélection des espèces, toxicité et bio-accumulation.

Les Zones Humides Artificielles (ZHA) sont de plus en plus utilisées comme installations rustiques de traitement d'eau. Elles sont connues pour leur capacité de traitement des eaux usées et leur faible coût d'installation et d'exploitation, leur facilité d'utilisation ainsi que leur bonne intégration dans le paysage. Néanmoins, leur efficacité varie dans le temps en raison de problèmes de fonctionnement tels que le colmatage dû à des charges excessives en matières en suspension ou à un dimensionnement inadéquat. La modélisation numérique des flux peut être utilisée pour évaluer les performances de filtration et améliorer la conception des ZHA. Les rejets urbains en temps de pluie ont été diagnostiqués comme étant la source principale provoquant la dégradation de la qualité du cours d'eau urbain appelé Ostwaldergraben et situé dans la ville de Strasbourg. Notre cas d'étude est la zone humide artificielle de l'Ostwaldergraben (constituée de trois couches différentes : une couche de sable, une couche de transition et une couche de drainage) engendrant de par sa conception un écoulement vertical en régime transitoire lié aux événements de précipitation en surface. Des périodes sèches prolongées peuvent potentiellement avoir des impacts sur les macrophytes ou sur les micro-organismes, ce qui peut induire une diminution significative de l'efficacité du traitement. Ainsi, pour éviter le stress hydrique des plantes, une zone saturée est mise en oeuvre dans la partie inférieure du filtre.

Le chapitre 4 s'intéresse à la modélisation par la méthode MHFEM des écoulements au sein de ces filtres plantés. Lorsqu'un filtre à écoulement vertical n'est pas saturé, l'infiltration a lieu au niveau du point d'alimentation, ce qui n'utilise qu'une faible partie du filtre. En conséquence, des concentrations plus importantes en sortie sont observées au début des évènements. Pour éviter l'amplification de ce phénomène sur le terrain, un soin particulier a été porté à la conception des systèmes de distribution en surface. Pour en rendre compte dans le modèle, un contrôle des conditions aux limites au sommet du filtre a été mis en place pour traiter les faibles débits d'alimentation ainsi que les débits d'alimentation plus importants. Le modèle prédictif de van Genuchten Mualem (vGM), modifié par l'introduction d'une pression d'entrée d'air pour améliorer son efficacité à proximité de la saturation, a été choisi pour décrire les propriétés hydrodynamiques du sol. Une telle modélisation implique Part II -General introduction de connaître les paramètres hydrodynamiques dits de van Genuchten-Mualem qui sont notoirement difficiles à mesurer. Une comparaison des résultats obtenus avec le code MHFEM, d'une part, et le logiciel commercial Hydrus, d'autre part, lors de la simulation de différents événements pluvieux montre que ceux-ci sont semblables.

Le chapitre 5 présente le processus d'optimisation que nous avons mis au point pour caler les paramètres du modèle vGM. Nous discutons en particulier une approche de modélisation inverse simple, robuste et peu coûteuse, déterministe et stochastique, pour l'identification des propriétés hydrodynamiques du sol à partir des mesures piézométriques obtenues lors des précipitations successives. L'optimisation inverse sur la minimisation d'une fonction "objective" est mise en place. Les écarts entre les données observées et les hauteurs piézométriques calculées sont évalués à l'aide de la fonction des moindres carrés. L'erreur quadratique moyenne et le critère de Nash-Sutcliffe sont largement utilisés pour le calage et l'évaluation de modèles hydrologiques avec des données observées. Compte tenu du fait que, dans un calage avec une approche déterministe, l'influence de l'estimation des valeurs initiales du paramètre sur les paramètres calibrés ne peut être négligeable, le choix d'un jeu de paramètres initial approprié permet de faire converger rapidement les estimations et d'obtenir le jeu paramètres recherché.

L'attention apportée dans cette thèse à la sensibilité du modèle aux paramètres vGM révèle une influence prédominante des paramètres de forme et de la conductivité à saturation du filtre. Les hauteurs piézométriques évaluées avec les paramètres optimisés pour les périodes d'événements pluviométriques montrent clairement que le fait d'ignorer l'effet d'hystérésis peut être la cause principale des discordances de ces résultats avec les hauteurs piézométriques observées. En divisant les périodes en sous-périodes, le processus du calage pendant les sous-périodes d'alimentation comporte plus de difficultés que celui de drainage. Cela peut être dû à la procédure de traitement des conditions aux limites pendant l'infiltration qui dépend de la saturation de la colonne de sol. La mise en oeuvre de l'optimisation du gradient pour ces sous-périodes fournit des courbes mieux ajustées. D'autre part, les courbes ajustées pendant les sous-périodes de drainage sont excellentes, quelle que soit la modélisation numérique. La variabilité temporelle des paramètres hydrodynamiques est évaluée et analysée pour montrer l'efficacité de la modélisation et l'effet de l'hystérésis (see [START_REF] Feddes | Modelling soil water dynamics in the unsaturated zone -state of the art[END_REF]) sur l'évolution de ces paramètres pendant les évènements de drainage et de précipitation en utilisant le critère de Nash-Sutcliffe (see [START_REF] Moriasi | Model evaluation guidelines for systematic quantification of accuracy in watershed simulations[END_REF]). Les paramètres calés lors des sous-périodes de drainage et d'alimentation qui sont fournies par les optimisations de gradient et stochastique, montrent principalement les mêmes variations temporelles. Cependant, quelques discordances entre les paramètres identifiés par les deux méthodes mises en oeuvre peuvent être dues au risque de piégeage dans les minima locaux en utilisant la méthode de gradient ou un nombre insuffisant de paramètres aléatoires traités par le calage stochastique. La variabilité temporelle des paramètres de vGM au long de la période devrait être associée à la dynamique du processus d'écoulement de l'eau à travers de la ZHA. En fait, la saisonnalité devrait modifier la dynamique de la végétation qui pourrait avoir des impacts sur le débit de l'eau. La période de printemps peut fournir l'activité racinaire la plus élevée qui doit augmenter la teneur en eau résiduelle. Par conséquent, les résultats montrent la plus grande valeur de la teneur en eau résiduelle pendant la période de printemps. Après plusieurs périodes de fonctionnement de la ZHA, les matières en suspension accumulés sur le dessus du filtre peuvent réduire la conductivité hydraulique de la couche supérieure. Néanmoins, il reste difficile d'analyser la tendance dynamique des paramètres de forme n et α.

CHAPTER 4

Variably saturated vertical flow constructed wetland modeling

Introduction

The present chapter provides a brief review on the effects of urban discharges and pesticides on environment and human health, potential threats and follows by presenting several guidelines. In order to minimize runoff pollution, constructed wetlands, whose efficiency and significant role on water treatment has been recently proven, are introduced as an efficient remediation. Thus, a short history about the utilization of constructed wetlands in water treatment is explored. The next section concerns an overview of the application of constructed wetlands for the treatment of discharges. As a modeling approach, a MH-FEM two-dimensional simulation with different numerical modeling techniques applied to tackle the Richards equation is presented. Consequently, this chapter covers the temporal increment and switching boundary condition procedures with a particular attention to the latter. Furthermore, several test cases are simulated as a validation of the model and the results are compared with those existing in the literature and HYDRUS modeling software.

Finally, the model is implemented for a typical CW during multiple-rainfall events employing different methods in order to select the appropriate model for the inverse modeling which is carried out in Chapter 5.

State of the art

Urban runoff pollution and threats to ecosystems

In urban areas, storm-water quality depends significantly on many factors such as land use features, number of dry days before a rain event, intensity and duration of rainfall events, topography, soil pollutant deposits and other factors which make the characterization of urban runoff much more complex than that of wastewater treatment [START_REF] Lee | Characterization of urban stormwater runoff[END_REF].

Increasing impervious surfaces can rise serious environmental issues by reducing infiltration and favoring surface run-off. This is commonly referred to as degradation of urban aquatic ecosystems caused in part by increased volumes and peaks of storm-water runoff and higher water pollutant concentrations contained in it [START_REF] Bernhardt | Restoring streams in an urbanizing world[END_REF][START_REF] Peters | Effects of urbanization on stream water quality in the city of atlanta, georgia, usa[END_REF][START_REF] Daly | Modelling of stormwater biofilters under random hydrologic variability: a case study of a car park at monash university, victoria (australia)[END_REF].

Urban wet weather discharges are episodic, but may have long term effects including physical and chemical alterations on the receiving media, depending on the number of events. They can impact different components of media: physical alterations (modification of flows, morphodynamics), chemical alterations (contributions of suspended matter, fertilizers, micro-pollutants), biological alterations (species selection, toxicity, bioaccumulation) and the hydro-system as a whole (eutrophication) [START_REF] Beck | Transient pollution events: Acute risks to the aquatic environment[END_REF][START_REF] Brent | A method for the toxicity assessment of wet weather events[END_REF]. These different types of effects have consequences on the different components of the receiving environments, particularly in the case of small rivers. These effects may be associated with changes of the hydraulic characteristics, an increase in suspended solids, fertilizers and micropollutants. The urban wet weather discharges can also generate selection mechanisms, toxicity and bioaccumulation in living organisms, or even influence on the ecosystem as a whole by promoting its eutrophication. The rainwater toxicity also depends on where it comes from. [START_REF] Bay | Water quality impacts of stormwater discharges to santa monica bay[END_REF] showed that waters from highly urbanized watersheds are much more toxic to marine environment than largely undeveloped watersheds.

Pollution sources

The concentrations of pollutants in rain water are extremely low, and generally the rainwater quality is good before it arrives at the ground level. The sources of the pollutants found in urban runoff are generally categorized with respect to the location and the use of the site. These may be due to the atmospheric pollution, the leaching of dry weather deposits, accumulation of fallouts on watersheds, the erosion of urban materials or the restoration of suspended pollutants in sewerage networks (see Figure 4.1). During periods of dry weather, a large number of pollutants are accumulated on roofs, roads, car parks, surfaces. Contaminants can either be dissolved or attached on particles. The concentrations of urban discharge pollutants depend on factors such as the intensity of the rain, the amount of runoff, the surface material and the nature of activities which are carried out on or near surface. Consequently, the most important factor is the distance traveled by the flow before the runoff can be collected in sewerage systems.

Threats to biotopes

Many studies have investigated the impact of urban discharges during rainy episodes as these pose serious risks to the environment in terms of the physico-chemical quality of the biotopes. It remains one of the major threats to surface water quality [START_REF] Chocat | La pollution due aux rejets urbains de temps de pluie : impacts sur les milieux rècepteurs[END_REF][START_REF] Parent-Raoult | Etude en systèmes artificiels de laboratoire des effets de Rejets Urbains de Temps de Pluie sur les communautès pèriphytiques : influence de facteurs d'èxposition[END_REF] (see Figure 4.2).

Pesticides may enter water bodies through diffuse and point sources. Among diffuse pollution, [START_REF] Reichenberger | Mitigation strategies to reduce pesticide inputs into ground-and surface water and their effectiveness; a review[END_REF] indicated surface runoff and erosion, leaching and drainage represent the major pathways.

The relationship between the pesticide exposure and probability of cancer development has been the subject of numbers of studies [START_REF] Brito | Pesticide exposure and cancer among rural workers from bariri, sao paulo state, brazil[END_REF][START_REF] O'leary | Pesticide exposure and risk of breast cancer: a nested case-control study of residentially stable women living on long island[END_REF]. The vulnerability in the reproductive system is surveyed by [START_REF] Petrelli | Environmental risk factors and male fertility and reproduction[END_REF] and the potential risk due to pesticide contamination in aquatic ecosystems has been assessed notably by [START_REF] Belfroid | Relative risks of transformation products of pesticides for aquatic ecosystems[END_REF], [START_REF] Yamaguchi | Concentrations and hazard assessment of pcbs, organochlorine pesticides and mercury in fish species from the upper thames: River pollution and its potential effects on top predators[END_REF], [START_REF] Capkin | Water quality and fish size affect toxicity of endosulfan, an organochlorine pesticide, to rainbow trout[END_REF], [START_REF] Houdart | Assessing pesticide pollution risk: from field to watershed[END_REF] and [START_REF] Van Den Brink | Effects of a herbicide-insecticide mixture in freshwater microcosms: Risk assessment and ecological effect chain[END_REF]. Furthermore, different procedures to evaluate and to analyze the environmental impacts of pesticides have been developed [START_REF] Levitan | how to" and "why":: assessing the enviro-social impacts of pesticides[END_REF][START_REF] Falconer | Pesticide environmental indicators and environmental policy[END_REF][START_REF] Finizio | Environmental risk assessment for pesticides[END_REF][START_REF] Bues | Assessing the environmental impacts of pesticides used on processing tomato crops[END_REF] [START_REF] Ritter | Sources, pathways, and relative risks of contaminants in surface water and groundwater: A perspective prepared for the walkerton inquiry[END_REF] 2004; De [START_REF] De Schampheleire | Risk assessment of pesticide spray drift damage in belgium[END_REF][START_REF] Hansen | Risk assessment of emerging contaminants in aquatic systems[END_REF][START_REF] Guérit | Environmental risk assessment: A critical approach of the european tgd in an in situ application[END_REF].

Guidelines

Several guidelines or standards values for maximum residue levels in food, soil, atmosphere and water by different regulatory agencies or organizations been have established in order to protect human health, environment and biodervisity [WHO., 1997; US-EPA-600/R, 2000; [START_REF] Hamilton | Regulatory limits for pesticide residues in water (iupac technical report)[END_REF][START_REF] Menard | Relevance of integrating agricultural practices in pesticide dietary intake indicator[END_REF]. Among them, the European Water Framework Directive [Commision Directive, 2008] introduced quality objectives to improve the global quality of water bodies across Europe to achieve good status in surface water and groundwater. Among the pollution sources which should be reduced, contaminants from storm-water runoffs are of great concern, since they could negatively alter the habitat for aquatic fauna if they are not treated properly. Depending on the type of catchment -urban, rural or industrial -the impact of storm-water runoffs on the receiving water body might change due to the large variety of pollutants involved [START_REF] Zgheib | Priority pollutants in urban stormwater: Part 1 -case of separate storm sewers[END_REF][START_REF] Ladislas | Floating treatment wetlands for heavy metal removal in highway stormwater ponds[END_REF][START_REF] Schmitt | Constructed wetlands treating stormwater from separate sewer networks in a residential strasbourg urban catchment area: Micropollutant removal and fate[END_REF].

As a result, storm-water managers are invited to apply the best practice environmental storm-water management regarding operational and environmental conditions of the treatment systems on removal efficiency [START_REF] Daly | Modelling of stormwater biofilters under random hydrologic variability: a case study of a car park at monash university, victoria (australia)[END_REF].

Constructed wetlands

4.2.2.1 Stormwater CW as an alternative solution to water treatment system Stormwater Constructed Wetlands (SCWs) are considered as a sustainable and promising option, whose performance, cost and resources utilization can complement or replace conventional water treatment [START_REF] Tack | Contaminants in natural and constructed wetlands: Pollutant dynamics and control[END_REF][START_REF] Arias | Feasibility of using constructed treatment wetlands for municipal wastewater treatment in the bogota savannah, colombia[END_REF][START_REF] Zhang | de cas Résumé Des analyses de sensibilité et des estimations de paramètres sont étudiées sur deux études de cas de transfert de masse en milieu poreux[END_REF].

Storm-water CWs are engineered systems to manage flood peaks and they are designed to favour natural treatment processes by involving wetland vegetation and soil, and their microbial flora. They are designed to take advantage of many of the processes that occur in natural wetlands, while being a more controlled environment [START_REF] Vymazal | Wastewater Treatment in Constructed Wetlands with Horizontal Sub-Surface Flow[END_REF]. Storm-water is diverted to the system and, by flowing vertically through the soil media, is subjected to several physical and chemical treatment processes, such as sedimentation, fine filtration, adsorption and biological uptake [START_REF] Davis | Bioretention technology: Overview of current practice and future needs[END_REF]. Among SCW systems, planted filters with horizontal or vertical flow have long demonstrated their effectiveness in the treatment of stormwater [START_REF] Albalawneh | Efficiency of a horizontal sub-surface flow constructed wetland treatment system in an arid area[END_REF]. The abilities of a SCW to improve water quality are widely recognized, and their efficiency at reducing the ecological impact of urban runoff by reducing suspended solids and micropollutants has been reported in several studies [START_REF] Ellis | of Hydrological Sciences, International Association; of Geodesy, International Union, and Assembly, Geophysics. General. Impacts of Urban Growth on Surface Water and Groundwater Quality[END_REF][START_REF] Huett | Nitrogen and phosphorus removal from plant nursery runoff in vegetated and unvegetated subsurface flow wetlands[END_REF][START_REF] Khan | Use of constructed wetland for the removal of heavy metals from industrial wastewater[END_REF][START_REF] Tang | Processes impacting on benzene removal in vertical-flow constructed wetlands[END_REF].

At present, constructed wetlands are used for a wide variety of pollution, including agricultural and industrial wastewater, various runoff waters and landfill leachate. Their design includes decisions about the amount of parallel flow paths and compartmentation. The removal efficiency of the wetland is a function of the surface area (length multiplied by width), while the cross-sectional area (width multiplied by depth) determines the maximum possible flow. Generally, a surface area of about 5 to 10 m 2 per person equivalent is required [US-EPA-600/R, 2000].

The capacity of constructed wetlands to retain pesticides is achieved through the process of sorption either by plant or by sediment material improving the quality of water. Generally, CWs have been considered as black boxes in which the measurements of the concentrations of influent and effluent pollutants are challenging or impossible. To understand and take advantages of CWs efficiently, this "black-box" which includes the operations of CWs has to be investigated [START_REF] Dordio | Wetlands: Water Living Filters?[END_REF]. Several studies has challenged the efficiency of wetlands, regarding insecticide retention, by monitoring input and output measurements [START_REF] Schulz | Effectiveness of a constructed wetland for retention of nonpoint-source pesticide pollution in the lourens river catchment, south africa[END_REF][START_REF] Moore | Mitigation of chlorpyrifos runoff using constructed wetlands[END_REF][START_REF] Avila | Capacity of a horizontal subsurface flow constructed wetland system for the removal of emerging pollutants: An injection experiment[END_REF]. By investigating the capacity of a CW in removing two particular pesticides, parathion and omethoate, [START_REF] Cheng | Xenobiotics removal from polluted water by a multifunctional constructed wetland[END_REF] observed that the decontaminating of synthetic waters from pesticides of related chemical structures was performed efficiently in spite of a low removal of herbicides. The biofilters are assumed to be designed with a combination of vegetation and soil filter media effective in removing nitrogen, as suggested in the literature [START_REF] Read | Plant traits that enhance pollutant removal from stormwater in biofiltration systems[END_REF][START_REF] Fawb | Guidelines for Filter Media in Biofiltration Systems. Facility for Advancing Water Biofiltration[END_REF]. [START_REF] Daly | Modelling of stormwater biofilters under random hydrologic variability: a case study of a car park at monash university, victoria (australia)[END_REF] assessed the hydrologic and nitrogen treat-ment performance of alternative applications of biofiltration for storm-water management across a range of climates.

Types of CWs

Water can move through a CW wetland system along several parallel paths and directions, considering both surface and subsurface flow processes. Factors such as the nature of the wetland vegetation, the surface topography and the properties of the wetland sediments influence the movement of water through a wetland.

Basically, CWs are categorized in three types. Their designs depend on the treatment target and the amount and the quality of the influent (see Figure 4.3). The direction of the flow within the porous medium (vertical or horizontal) influences very significantly the oxygen renewal of the porous medium, and therefore the type of microorganisms which are present in.

Surface Flow Constructed Wetlands (SFCWs), similar to natural wetlands, consist of a submerged layer of sand or gravel in which plants are rooted in. They may have one or more shallow basins with an impervious barrier to prevent seepage to groundwaters. SFCWs are less expensive to be constructed, they require large land area, especially if nitrogen or phosphorus removal are required.

Horizontal Sub-surface Flow Constructed Wetlands (HFCWs) have been used primarily for the treatment of municipal or domestic wastewater. There is a variety of industrial wastewaters which have been treated in HF constructed wetlands [START_REF] Wallace | Treatment of cheese-processing waste using subsurface-flow wetlands[END_REF][START_REF] Yang | Treatments of oil-refinery and steel-mill wastewaters by mesocosm constructed wetland systems[END_REF][START_REF] Davison | Aspects of design, structure, performance and operation of reed bedseight years' experience in northeastern new south wales, australia[END_REF][START_REF] Mbuligwe | Comparative treatment of dye-rich wastewater in engineered wetland systems (ewss) vegetated with different plants[END_REF][START_REF] Dotro | Treatment of chromium-bearing wastewaters with constructed wetlands[END_REF]. The quality of industrial wastewaters varies widely, with many wastewaters having very high concentrations of pollutants.

A Vertical Flow Constructed Wetland (VFCW) is a planted filter bed that is drained with a drainage pipe at its bottom. Contaminated water is distributed across the top surface by means of a mechanical dosing system. The main difference between a vertical and a horizontal constructed wetland, beyond the direction of the flow path, is the aerobic conditions. A VFCW should be more aerobic and thus intermittently loadings improve the oxygen diffusion on the filter bed [START_REF] Vymazal | Constructed wetlands for wastewater treatment: Five decades of experience[END_REF].

Vertical flow reed beds are basins filled with layers of gravels of different grain sizes superimposed, covered with a layer of sand in which the macrophytes are planted. The VFCWs are filled with layers of gravel or sand of different particle sizes depending on the quality of the wastewater to be treated. These are then distributed uniformly on the filter. The input effluent then percolates through the porous medium and then goes under physical, chemical and biological treatments. The treated effluent is collected on a layer of pebbles and recovered by drains [START_REF] Molle | How to treat raw sewage with constructed wetlands: an overview of the french systems[END_REF].

VFCW systems, which have commonly high removal efficiency, are recently used more to reduce the required land area. Many studies of VFCW were conducted in order to find the maximum removal efficiencies regarding different kinds of media [START_REF] Korkusuz | Comparison of the treatment performances of blast furnace slag-based and gravel-based vertical flow wetlands operated identically for domestic wastewater treatment in turkey[END_REF], different bed depth [START_REF] Taniguchi | Evaluation of extremely shallow vertical subsurface flow constructed wetland for nutrient removal[END_REF] and different loading regimes [START_REF] Kantawanichkul | Effect of dosing regime on nitrification in a subsurface vertical flow treatment wetland system[END_REF].

The plant roots in the filters also have an important role. The most used plants in Europe are Phragmites Australis. In addition to the obvious aesthetic aspect, the plants provide a slight assimilation of nutrients and metals and provide oxygen to the porous medium, developing specific microorganisms around rhizomes or roots.

The vegetation has a positive impact on the slowing down of the clogging phenomena of the porous medium. Regarding the substrate chosen for the porous medium, it has an absolutely fundamental role, since it is responsible for sedimentation, filtration, sorption and chemical precipitation. The composition of sand is also important and influences for instance the retention of phosphorus. The thickness of the filter layer is important to allow the purification mechanisms to develop, but it can generally be limited to a few tens of centimeters.

The choice of the substrate is fundamental and must be adapted according to the objectives and the mechanisms of elimination of different pollutants, such as filtration, sorption and precipitation [START_REF] Garcia | Contaminant removal processes in subsurface-flow constructed wetlands: A review[END_REF].

The type granulates used has influences both in terms of hydrodynamics and in the treatment of pollutants. Indeed, [START_REF] Stottmeister | Effects of plants and microorganisms in constructed wetlands for wastewater treatment[END_REF] indicated that for filters planted with reeds, the main parameter that influences hydrodynamics is the size of particles of the substrate. The finer the sand, the lower its hydraulic conductivity. The characterization of soil is necessary as the knowledge of the properties of modeling parameters, such as hydraulic conductivity, and land-use effects on these properties are of great importance for an efficient soil and water management.

Precipitations and performance of CWs

VFCWs are accumulative systems and the hydraulic overloads might penalize the filter longevity. Hence, the intervals between feeding and resting periods might be planned accordingly.

Rainfall intermittency and dry periods have an important role in the improvement of such systems. The periods of dry weather allows for the mineralization of the upper layer of the porous medium. However, extended periods of dry weather may affect the quality of treatment of the rainy event. In spite of this, the purification efficiency of the filter can be recovered in a short time when they are fed again. Thus, the feeding-rest periods are of great importance to guarantee the growth of the biomass on the filter media to maintain aerobic conditions. In addition, urban discharges caused by runoff can result chemical, physical and biological characteristics changes of the receiving environment. It can considerably modify the flow as morphodynamics, generating strong suspended solids, fertilizers or micropollutant and can cause the emergence of a phenomenon of species selection, toxicity, bioaccumulation and eutrophication. The smaller the surface of the receiving area, the greater the effect of rainwater pollution. The water flow is more polluted at the beginning of the rainy event because it carries the pollutants accumulated in dry weather. This is the principle of first flow which depends on several factors, including the basin area, its shape, the intensity of the rain, the amount of rain, the period of precipitation, the area of the impervious surfaces and drainage systems. The accumulation of these contaminants contributes to the amplification of the siltation phenomenon because of the abundance of pollutants being carried by precipitations. Eventually, the level of the water will tend to increase gradually as the total surface of the watercourse decreases due to erosion of the watercourse.

During rainy events, ponding may appear at top surface of the filter and hence, infiltration velocity may increase. Oxygen supply can be restricted for a longer time due to more intense storm events producing longer and heavier ponding. Long periods of oxygen restriction can considerably affect filter performances regarding nitrification and organic matter degradation [START_REF] Forquet | Diphasic modelling of vertical flow filter[END_REF]. Indeed, hydraulic overloads may be physically acceptable on filters, however, constant ponding effects the biological activities in a negative manner due to the deficiency of oxygen renewal. Although organic and hydraulic loads in VFCWs may vary significantly, hydraulic load shall not exceed certain defined limits. To avoid the negative effects of prolonged dry periods, and thus their impact on plants and microorganisms, it is recommended to set a permanently saturated layer at the bottom of the filter. However, additional ventilation is required to prevent odor nuisance and maintain aerobic conditions. This saturated zone has little or no impact on the quality of treatment. [START_REF] Molle | Effect of reeds and feeding operations on hydraulic behaviour of vertical flow constructed wetlands under hydraulic overloads[END_REF] found that after each feeding period on a plant running for one year, three or four days of rest (without any water infiltration) were necessary to restore the infiltration rates observed at the beginning of the feeding period. Thus, under hydraulic overloads, we can assume that the rest periods are essential to occur more in order to minimize the ponding times and thus to avoid the risk regarding the lack of oxygen. In a VFCW system, feeding is programmed to be greater than infiltration capacity in order to provide a sufficient water distribution onto the filter surface. This ponding means that the hydraulic flow able to pass through the filter is limited. However, this is not necessarily the case on a newly-installed treatment plant where the infiltration material is still clean [START_REF] Molle | Filtres plantés de roseaux: Limites hydrauliques et retention du phosphore[END_REF]. During spring, summer, and fall, surface ponding may be effectively eliminated due to evapotranspiration demands [Army Corps of Engineers, 1976].

Context of the study

Several issues concerning the design, construction and operation of CWs must be considered in order to increase the potential treatment efficiency and to avoid malfunctioning problems such as clogging due to excessive loads of suspended solids, water stress of macrophytes and microorganisms due to extended dry periods and inadequate CWs sizing. Because of the complexity of the processes in CWs, all the available design guidelines and rules are based on empirical rules of thumb [START_REF] Brix | The use of vertical flow constructed wetlands for on-site treatment of domestic wastewater: New danish guidelines[END_REF]. Parameters of existing design guidelines and rules have been derived from experiments under specific conditions. Therefore, designing CWs using these parameters is limited to these conditions. The predictive ability of simulation tools is of great interest to evaluate filtering characteristics of a CW Langergraber [2011]. Furthermore, they can be used to conduct a reliable sensitivity analysis and eventually to provide an appropriate prediction ability outside the calibration range. Once the filter has reached saturation, infiltration is dictated by the permeability of the media and the pressure forces exerted by ponding depth. For modeling VFCWs with intermittent feeding, a transient variably-saturated flow model has to be considered. These systems are generally highly dynamic, due to the intermittent loading, and thereby, more numerical difficulties are generated. Models capable to assess the performance of VFCWs can be the Richards equation which is known to perform well to describe variably-saturated flow. The spatial discretization of the Richards equation is commonly performed using the finite difference method [START_REF] Van Dam | Numerical simulation of infiltration, evaporation and shallow groundwater levels with the richards equation[END_REF][START_REF] Brunone | Numerical simulations of one-dimensional infiltration into layered soils with the richards equation using different estimates of the interlayer conductivity[END_REF] or finite elements [START_REF] Pan | Finite element methods for modeling water flow in variably saturated porous media: Numerical oscillation and mass-distributed schemes[END_REF][START_REF] Lehmann | Comparison of iterative methods for improved solutions of the fluid flow equation in partially saturated porous media[END_REF][START_REF] Diersch | On the primary variable switching technique for simulating unsaturatedsaturated flows[END_REF]] method.

Considering a numerical method suitably applied to Richards model, here we have focused on MHFEM the essential idea of which is to approximate simultaneously the pressure and its gradient (MFE) and to enforce the continuity of the normal component (hybridization) [START_REF] Mosé | Application of the mixed hybrid finite element approximation in a groundwater flow model: Luxury or necessity?[END_REF]. In addition to the advantages of the MHFEM comparing to the FEM which are explained in Chapter 2, it allows to conserve mass locally [START_REF] Hoteit | The maximum principle violations of the mixedhybrid finite-element method applied to diffusion equations[END_REF]. Moreover, MHFEM can be efficiently applied to flow with general irregular grids and highly heterogeneous permeability. This heterogeneity is due to both the heterogeneous sediment distribution and the non-uniform water content in the constructed wetland.

For a better convergence behavior in a heterogeneous medium, the applied technique switches between the mixed-form and the pressure-head based form of the Richards equation [START_REF] Wanko | Contribution to numerical modeling of water flow in variably saturated, heterogeneous porous media[END_REF]. Modeling issues arise when the soil reaches oven-dry conditions. A particular attention should also be brought to boundary condition modeling (surface ponding or evaporation) to be able to tackle different sequences of rainfall-runoff events.

A mass lumping scheme is used in order to avoid oscillation problems related to the discrete expression of the term representing mass variation in the volume by separating the stationary and accumulation parts of the flux (terms with a time derivative). Nonphysical oscillations are of great importance when modeling flows in variably saturated porous media, particularly for water infiltration problems in initially dry soils. Furthermore, the lumped formulation considerably reduces the CPU time compared with the standard one. In addition, to improve the convergence, a transformation of the variables is adapted to the MHFEM approximations [START_REF] Younes | A new mass lumping scheme for the mixed hybrid finite element method[END_REF].

A short explanation of the techniques used to solve the resulting systems of equations and the convergence criteria are presented in Appendix C followed by the procedure used for temporal discretization and procedure for time stepping adjustment at each time level. This discretization procedure is inspired by the thesis authored by [START_REF] Padilla | Modeling and Optimization of Biological Mitigation Processes in Porous Matrices Contaminated by Pesticides: Towards a New Functionality of Stormwater Basins[END_REF].

For validation purposes, several test cases including infiltration under intensive rain at a dry soil, high evaporation at a wet soil and [START_REF] Pan | A transformed pressure head-based approach to solve richards' equation for variably saturated soils[END_REF] are considered. Finally, the flow in a three-layered CW is simulated using the vGM parameters by literature during multi-rainfall events. The results are compared with those approximated by HYDRUS (a software based on a finite element discretization).

Hydrodynamic modeling 4.4.1 Governing equations

The soil can be assimilated to a porous medium formed by solid aggregates (mineral and organic elements) between which there are empty spaces that can be filled with gas or liquid. A detailed presentation of the modeling of mass transfer in porous medium is given as below.

The Darcy's law (Equation (2.7)) in the vertical flows, when the volumetric density is essential to be constant, can be written as:

q = -K(h)∇(h + z). (4.1)
Therefore, for an unsaturated flow with an incompressible fluid, through porous media, the hydrodynamic equation is obtained by associating the continuity equation with the Buckingham-Darcy's law in which the hydraulic conductivity depends on the pressure h (L) head or the water content θ (L 3 L -3 ):

               ∂ θ ∂t + S s S w (θ ) ∂ h ∂t + ∇. q = f in Ω, q = -K(h)∇(h + z) in Ω, h(x, y, 0) = h 0 (x, y) in Ω, h(x, y,t) = h 1 (x, y,t) on Ω D , K(h(x, y,t)) ∂ h ∂ n (x, y,t) = Q n on Ω N , (4.2)
where t is the time variable (T ), Ω is the flow region considered as a 2D domain, S s the specific storage (L -1 ) and S w the degree of saturation (-). K is the unsaturated hydraulic conductivity (LT -1 ). The hydraulic conductivity varies under the influence of capillary and gravity forces depending on saturation of soil. As the water content decreases, the gravity has less control on water infiltration.

Variable transformation

Some researches have attributed convergence difficulty to the presence of sharp wetting fronts [START_REF] Diersch | On the primary variable switching technique for simulating unsaturatedsaturated flows[END_REF][START_REF] Williams | An evaluation of temporally adaptive transformation approaches for solving richards' equation[END_REF]. To smear the wetting front and to improve the convergence, the pressure head variable h is transformed into a new dependent variable ĥ [START_REF] Pan | A transformed pressure head-based approach to solve richards' equation for variably saturated soils[END_REF].

ĥ = h 1+κh , for h < he, h , for h ≥ he, (4.3)
where

• κ a is a constant independent of both the K(h) and C(h) relationship. [START_REF] Pan | A transformed pressure head-based approach to solve richards' equation for variably saturated soils[END_REF] suggest a value ∼ = -0.04 cm -1 or -4 m -1 .

• ĥ is transformed pressure head (L).

Then the Darcy's law can be expressed as

q = -K ∂ h ∂ ĥ .
∂ ĥ

∂ h ∇(h + z) = -K ∂ ĥ ∂ h .∇ ĥ -K∇z. (4.4)
Using the inverse equation of the transformed variable: h = ĥ 1-κ ĥ , for ĥ < he, ĥ , for ĥ ≥ he.

(4.5)

A transformed hydraulic conductivity is defined as K:

K = K ∂ h ∂ ĥ    K 1-κ ĥ 2 , for h < he, K , for h ≥ he, (4.6)
or in terms of the variable h, it can be rewritten as:

K = K 1 + κh 2 , for h < he, K , for h ≥ he. (4.7)
Thus, equation (4.4) can be expressed as: q = -K∇ ĥ -K∇z . (4.8)

Total porosity

The total porosity of a porous medium represents the density of the pores that can be occupied by a liquid or gaseous fluid. It can be expressed by the ratio between the volume of the void spaces and the total volume of the porous medium. The soil system is composed of three phases of solid, liquid and gas:

φ = V l +V g V s +V l +V g , (4.9)
where φ is the total porosity (-), V s , V l and V g are the the volumes of the solid phase, the liquid phase and the gaseous phase respectively (L 3 ).

Effective porosity

Effective porosity n e is the ratio of the volume of the pores which are interconnected in the rock to the total volume of this rock. Generally it can be defined as the portion of the soil or rock through which chemicals move, or that portion of the media that contributes to flow [START_REF] Fetter | Contaminant hydrogeology[END_REF]. The effective porosity is estimated by laboratory or field tracer test, but often they are relatively expensive and time consuming [START_REF] Stephens | A comparison of estimated and calculated effective porosity[END_REF].

n e = q v , (4.10) where q is the specific discharge and v is the mean velocity of a conservative tracer (seepage velocity).

Hydraulic conductivity

Hydraulic permeability or conductivity is the parameter relating the flow velocity to the hydraulic gradient in Darcy's law. It characterizes the ability of the medium to allow a fluid to pass through it under the effect of a gradient of pressure head. The conductivity can be estimated by the relationship:

K = K s K r K A ,
(4.11) where K s is the saturated hydraulic conductivity (LT -1 ), K A is a dimensionless anisotropy tensor and K r is the relative hydraulic conductivity function. In the case of a twodimensional and anisotropic medium, K A is defined by a symmetric matrix:

K A = K xx K xz K zx K zz .
Supposing that hydraulic conductivity tensor coincide with the coordinate axis x and z then

K A = K X 0 0 K Z .
An isotropy ratio a = K Z K X can be defined as a relationship between the components of the hydraulic conductivity matrix. Depending on the pore structure, this ratio can be equal to one (K A = I = 1 0 0 1 ). In pore structures with a preferred orientation in the plane, the values of K x and K z are not identical. This indicates that the medium is anisotrope. The relation between θ and h can be described by the van Genuchten-Mualem (vGM) model.

The Brooks and Corey model shows a discontinuity of the derivative of the content at the air-entry pressure. This mathematical modification leads the van Genuchten model to be continuously differentiable. However, fine textured soils may be less well represented by this model, especially near to saturation. Therefore, modifications have been proposed to introduce the equivalent of an air inlet value (h e ), but with restricted values compared to those of Brooks and Corey model [START_REF] Vogel | Effect of the shape of the soil hydraulic functions near saturation on variably-saturated flow predictions[END_REF]. The dimensionless relative hydraulic conductivity function K r is given by the modified vGM model [START_REF] Ippisch | Validity limits for the van Genuchten-Mualem model and implications for parameter estimation and numerical simulation[END_REF]:

K r =      Se τ 1-1-(SeSc) 1/mv mv 1-(1-Sc 1/mv ) mv 2 , for Se < 1, 1 , for Se ≥ 1, (4.12)
where Se is the effective saturation, which is approximated by the expression:

Se = 1 S c [1 + |α v h| n v ] -m v , for h < h e , 1
, for h ≥ h e , and S c is the saturation at the cut-off point h e in the classical van Genuchten model which is given by: and: τ is is an empirical parameter for tortuosity (-), α v is a free parameter related to the mean pore size of the soil (L -1 ), n v is a free parameter related to the uniformity of the soil pore-size (-), h e is a free parameter, referred as the air entry value (L), m v is a parameter in the modified vGM expression given by m

S c = [1 + |α v h e | n v ] -m v ,
v = 1 -1 n v (-).
The maximum value obtained from the three hydraulic conductivity values at the edges for each element is considered as the average hydraulic conductivity in that element.

Effective Saturation

The dimensionless value Se, referred as effective saturation, was defined by van Genuchten [1980] as a normalized water content:

Se = θ -θ r θ s -θ r , (4.13)
where θ r is the residual water content (L 3 L -3 ), defined as the water content for which the gradient dθ dh becomes zero and θ s is the saturated water content (L 3 L -3 ).

Specific water capacity

The specific water capacity is defined as the slope of the soil moisture characteristic curve at given pressure head: C(h) = dθ dh . Therefore, specific water capacity is given by:

C(h) = d (Se [θ s -θ r ] + θ r ) dh = α v m v n v (θ s -θ r )|α v h| nv-1 S c [1+(α v h) nv ] mv+1
, for h < h e , 0

, for h ≥ h e .

(4.14)

Top boundary conditions

On a VFCW, the monitored infiltration on the filter depends on several factors such as infiltration rate at surface, the ponding depth and ponding time which may occur at the top of CW, and thus the oxygen content in porous media.

A particular attention is given to the switching boundary condition procedure, which is important for situations with ponded water layers or fluctuating groundwater close to the soil surface. This procedure which is adapted from Van Dam and Feddes [2000] (see Figure 4.4), switches from head (h top ) to flux (q top ) controlled boundary condition and vice versa depending on the saturation of the soil column at the top. This algorithm was adapted to be used with the mixed hybrid formulation and 2D flow problem. The inflow Q in into the soil column is calculated considering the fluxes positive if they are directed upward. Q in includes the flux (LT -1 ) at soil profile bottom q bot , the potential flux at the soil surface q top , root water extraction q root , and the total lateral flux to drains or ditches q drain during the time step. If Q in is positive, more water enters than leaves the Saturated ? soil column. Hence a head condition is used. If Q in is negative, the soil profile becomes unsaturated and a flux condition is imposed. If the soil column is unsaturated, a comparison between Q in and V air (the total air volume in the soil profile at the start of the time step) determines whether the soil column will remain unsaturated or becomes saturated during the time step. If the soil becomes saturated a head condition is used. In the case the soil remains unsaturated, the procedure distinguishes between evaporation and infiltration. The evaporation is limited to E max , a maximum soil evaporation (LT -1 ), which is in relationship with h atm , the soil water pressure head in equilibrium with the air humidity (L): h atm = 13.3 × 10 5 In(e act /e sat ) where e act and e sat are the actual and saturated vapor pressure, respectively ( M L.T 2 ). (h atm ≈ -2.75 × 10 5 cm) [Kroes and Van Dam, 2003]. In the case of infiltration, if the potential flux q top exceeds I max (the maximum infiltration flux at the soil surface (LT -1)) and K sat (the top soil saturated hydraulic conductivity (LT -1 )), a head boundary condition is prescribed as h pond , the height of water ponding (L) on the soil surface at t j which is approximated by:

Q in > 0 Q in > V air Q in < h pond q top < I max q top = q top h top = Q in h top = h pond h top = Q in -V air q top > 0 q top < I max and q top < -K s h top = h pond q top = q top q top > E max h top = h atm q top = q top
h n+1 pond = h n pond + (Q load -Q in f )∆t n , ( 4
.15) where Q load is the quantity of water loaded during a time lapse (Q load = Q top ) and Q in f is the quantity of actual flux infiltrating at the surface during the same time step.

Mass balance error and convergence criteria

Mass balance error: A global mass balance error ε MB can be computed based on the ability of the model to conserve mass [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF]. The mass balance measure MB is defined as: MB = total additional mass in the domain total net flux into the domain

= nt ∑ n=1 nm ∑ G=1 Mass G nt ∑ n=1 nm ∑ G=1 Flux G (4.16)
and ε MB = 1 -MB (4.17) where

Mass G =      |G| ∆t n Ĉn+1,m G + θ n+1,m G φ G S s n+1,m G ĥn+1,m+1 G -ĥn G , i f θ n+1,m G,E i φ G > tol f , |G| ∆t n θ n+1,m+1 G -θ n G + |G| ∆t n θ n+1,m G φ G S s n+1,m G ĥn+1,m+1 G -ĥn G , i f θ n+1,m G,E i φ G ≤ tol f , Flux G = α n+1,m G ĥn+1,m+1 G + K-1,n+1,m G K n+1,m G z G - 3 ∑ i=1 α n+1,m G,E i T ĥn+1,m+1 G,E i + K-1,n+1,m G K n+1,m G z G,E i -|G| f n+1 G .
A mass balance ratio can be computed for each time step MB n , with its corresponding error ε MB n = 1 -MB n . Locally, mass balance error is approximated as:

E 1G = ∆t n |G| 3 ∑ i=1 Q n+1,m+1 G,E i -f G ∆t n + 3 ∑ i=1 E 1G i , (4.18)
where

E 1Gi =        1 3 T Ĉn+1,m G,E i + S s n+1,m G T θ n+1,m G,E i φ G T ĥn+1,m+1 G,E i -T ĥn G,E i , i f T θ n+1 G,E i φ G > tol f , 1 3 T θ n+1,m+1 G,E i -T θ n G,E i + S s n+1,m G T θ n+1,m G,E i φ G T ĥn+1,m+1 G,E i -T ĥn G,E i , i f T θ n+1 G,E i φ G ≤ tol f , E 2G = ∆t n |G| (Mass G + Flux G ) . (4.19)
Kosugi [2008] concluded that is important to check the mass balance and solution convergence at each time step to determine the accuracy of the numerical scheme when using the discretization to simulate unsaturated water flow. It has to be noticed that, a low mass balance error is a necessary but not a sufficient condition to ensure accuracy to the solution.

Mass balance error can decrease even if the solutions do not converge [START_REF] Tocci | Accurate and economical solution of the pressure-head form of richards' equation by the method of lines[END_REF][START_REF] Kosugi | Comparison of three methods for discretizing the storage term of the richards equation[END_REF].

Numerical solution and convergence criteria: In this system of linear equations, the number of unknowns which are the water pressure traces (T h), is equal to the number of edges to which the pressure has not been imposed. The matrix associated to the hydrodynamics equations can be solved by the conjugate gradient method, preconditioned with an incomplete Cholesky decomposition using the Eisenstat procedure [START_REF] Eisenstat | Efficient implementation of a class of preconditioned conjugate gradient methods[END_REF].

The iteration process for the hydrodydamics using the mass condensation scheme is stopped when the following convergence criteria are met:

• The discrepancy between the approximated pressure head of the edge for two successive iteration levels should be smaller than an absolute iteration convergence tolerance. The smaller the values of tolerance, the more accurate the solution. However it makes the computational time increase. Taheri Shahraiyni and Ataie Ashtiani [2008] defined the range of this value between 0.001 cm and 1 cm.

T h n+1,m+1 G,E i -T h n+1,m G,E i ≤ tol a .
(4.20)

• The difference between the calculated values of water content between two succes-sive iteration levels should be smaller than a tolerance. A value of 0.0001 for this tolerance criterion is determined by [START_REF] Huang | A new convergence criterion for the modified picard iteration method to solve the variably saturated flow equation[END_REF]. This convergence criterion is specially for cases where water content changes dramatically with small changes in the pressure head.

T θ n+1,m+1 G,E i -T θ n+1,m G,E i ≤ tol c . (4.21)
• The iteration convergence test, which involves both absolute and relative error, reduces the number of iterations, particularly when the pressure head changes significantly but not the water content. Values for the relative tolerance generally in the range of 10 -2 to 10 -5 are often used depending upon the desired accuracy [START_REF] Kavetski | Adaptive time stepping and error control in a mass conservative numerical solution of the mixed form of richards equation[END_REF]. This mixed criterion will serve to

T h n+1,m+1 G,E i -T h n+1,m G,E i -tol r T h n+1,m+1 G,E i -tol a < 0 . (4.22)
For the standard MHFEM formulation, the iteration process is stopped when the relative residual norm is smaller than a relative tolerance predetermined by the user.

T h n,m+1 -T h n,m T h n,m+1 = n f ∑ i=1 T h i n,m+1 -T h i n,m 2 n f ∑ i=1 T h i n,m+1 2 ≤ tol r . (4.23) max( θ G φ G ) < tol indsat ? Pe hyd = 1 |G| ( 1 C G dK dh | G ) 2 1 |G| K G C G V = 1 |G| ( 1 C G dK dh | G ) 2 V = ( q G,x θ G ∆x ) 2 + ( q G,z θ G ∆z ) 2 ∆t max G,hyd = 1 V ∆t max G,hyd = 0.5|G| C G K G ∆t New = min(∆t Ghyd ) ∆t heuristic > ∆t new ? ∆t new > T S ? ∆t = ∆t heuristic ∆t = ∆t new ∆t = ∆t init ∆t + t n ≥ T S ?
Time step= 1 ?

N iter > Maximum N iter ? Heuristic method ∆t = T S -t n ∆t min ≤ ∆t ≤ ∆t max ? Co = ∆t V Fn hyd = Co

Numerical results

Ponding conditions

Infiltration under intensive rain at a dry soil is simulated for a sandy soil with initial groundwater level at -8320 mm depth. At the reference case [START_REF] Van Dam | Numerical simulation of infiltration, evaporation and shallow groundwater levels with the richards equation[END_REF], the hydraulic head gradient at the soil surface is large enough to absorb the infiltration rate of 1000 mm d -1 but at the time of simulation t= 0.008 d, the flux boundary condition is changed to a head one (h top = h pondmax = 0.0 mm, see Figure 4.4), and the infiltration rate starts to decrease gradually. The calculated infiltration rate approximated by the simulations carried out either without or with switching technique and variable transformation are shown in Figure 4.6 which are close to that of the reference [START_REF] Van Dam | Numerical simulation of infiltration, evaporation and shallow groundwater levels with the richards equation[END_REF].

The cumulative amount of infiltration are given in Table 4.1 which is nearly 4 cm, however, the cumulative potential infiltration equals 10 cm . Evaporation flux-controlled condition to head-controlled one (h top = h atm = -1377000 mm, see Figure 4.4) and the evaporation rate gradually decreases. After 5 days, the cumulative actual evaporation amounts of the simulations carried out by non-transformed variables is 11 mm, while it is approximated to be 1.02 cm for those obtained by variable transformation (see Table 4.2) and the cumulative potential evaporation equals 25 mm.

Switching technique and variable transformation

Several one-dimensional cases differing by characterization of the porous media (soil properties and simulation parameters are listed in Appendix D), the initial and boundary conditions (see Table 4.3). Indeed, the test cases 1. It has to be remarked that the results for these 12 cases are compared with those of [START_REF] Pan | A transformed pressure head-based approach to solve richards' equation for variably saturated soils[END_REF]. Two kinds of boundary conditions were used: -50,000 3.4722 × 10 -4 cm s -1 0 cm s -1 21600 Layered 1.2 -1,000 3.4722 × 10 -4 cm s -1 0 cm s -1 18000 Layered 1.3 -200 3.4722 × 10 -4 cm s -1 0 cm s -1 13680 Layered 2.1 -50,000 8.3333 × 10 -5 cm s -1 0 cm s -1 43200 Layered 2.2 -1,000 8.3333 × 10 -5 cm s -1 0 cm s • A constant head condition at the upper and lower ends of the second (uniform) soil profile was used.

The initial conditions used for each problem are given in Table 4.3.

Simulations were performed using the Richards equation on their h-based form, mixed-form or using a switching method between these two forms. A transformed pressure was introduced as the dependent variable and results were compared to those without using transformation of variable. The test cases 3.1 and 4.1 are related to the simulations in an arid soil imposing positive and negative pressures as Dirichlet boundary conditions, respectively. The test cases 3.2 and 4.2 simulate the infiltration in a semi-arid medium by imposing positive and negative Dirichlet boundary conditions, respectively. Similarly, the test cases 3.3 and 4.3 are related to the the infiltration in a wet soil by imposing Dirichlet boundary conditions.

Approximated depth distributions of pressure head and volumetric water content for

Pressure head and water content distributions are very similar to those given by [START_REF] Pan | A transformed pressure head-based approach to solve richards' equation for variably saturated soils[END_REF].

In order to analyze quantitatively the efficiency of the considered methods, three parameters indicating time steps, CPU time and mass balance error are distinguished and compared. A glance at Table 4.4 suffices to see that in all the test cases, except for test cases 4.1 and 4.2, the transformation the primary variables can reduce time steps. Based on these parameters and according to the established criteria of selection, the models in which the non-transformation of variable coupled to the h-based form of Richards equation with the exception of the test case 4.1 and 1.3.

Multiple rainfall-runoff events and drainage in a loamy soil

In order to investigate the spatial variation of groundwater responses to a transient top boundary condition in a 80 cm column of loamy soil, a single scenario of double precipitation events is simulated and the verification was also performed using HYDRUS 1D. The initial water piezometric head for all depths were specified as -100 cm. At the bottom of the soil, a "seepage" face boundary condition is chosen such that if the simulated piezometric head is greater than 22 cm the drainage will be activated. Soil properties and simulation parameters are listed in Appendix D.

In Figure 4.11, the piezometric head at different elevations of the column are compared to those of simulations obtained by HYDRUS. The results provided by MHFEM and HYDRUS are in a good agreement.

The water pressure head time-series at the top surface with run off are compared with the results performed with surface ponding. The maximum height of water ponding on the soil surface is considered as 5 cm. 4.5 the indicators time steps, CPU time and mass balance error can be compared for all methods. In the scenario within surface run off, and similarly within surface ponding, the transformation of primary variables can slightly reduce the CPU time implementing the h-based and mixed methods. It has to be remarked that the CPU times for simulations within surface runoff are about three times more than those within surface ponding.

Conclusions

Various test cases in a variably saturated porous media were simulated by employing different techniques using a MHFEM. These techniques, which are aimed to improve the hydrodynamic modeling, are compared together in terms of CPU time and mass balance error and for validation, they are compared to those of the relevant recent references. A promoted auto-selective boundary conditions procedure is implemented to be able to tackle transient flow regimes switching from Dirichlet to Neumenn boundary condition and vice versa. These techniques are also implemented for a double rainfall events scenario, taking in account both runoff and ponding, in a typical homogeneous VFCW. The results are very similar to those reported in the references and those provided by HYDRUS. The transformation of primary variables can decrease slightly the CPU time whilst the mixed form Richards equation, specifically in runoff surface, increases it considerably. Consequently, the choice of proper method depends on the aim of the study (inverse method, accuracy of the results, initial and boundary conditions).

Parameter identification-Inverse method Contents

Introduction

The calibration of hydrodynamic parameters for subsurface constructed wetlands (CWs) is a sensitive process since highly non-linear equations are involved in unsaturated flow modeling [START_REF] Ritter | Using inverse methods for estimating soil hydraulic properties from field data as an alternative to direct methods[END_REF].

As vGM soil hydrodynamic parameters depend on water content, their estimation is subject to considerable experimental and numerical studies. Moreover, many of the parameters can generally not be measured directly in the field. These can be estimated through calibration against a historical record of variations in specific measured data.

The purpose of this chapter is to quantify the vGM soil parameters by monitoring the discrepancies between field-based observations and numerical simulation and by minimizing the mean square error so that the model reproduces at best the observations. In particular, the sensitivity analysis performed with respect to the vGM parameters reveals a predominant influence of the shape parameters α, n and the saturated conductivity of the filter on the piezometric heads, during saturation and desaturation.

This chapter is briefly organized as follows. Fundamentals are reviewed in Section 5.2. A detailed description of the area (Ostwaldergraben, Strasbourg) and its instrumentations is provided is Section 5.3. The proposed calibration methodology for the soil hydrodynamic parameters are discussed in Section 5.4. The calibration process is performed for multirainfall event periods. The temporal variability of vGM parameters during these periods are evaluated and analyzed in Section 5.6.

State of the art

According to the explanation of [START_REF] Asme | Guide for verification and validation in computational solid mechanics[END_REF] guide for the representation of a physical system in the real world, a CW can be modeled in three different manners, from the general to the specific: (i) a conceptual model, (ii) a mathematical model, and (iii) a numerical one. Eventually, the computed numerical model results can be compared to available experimental data for model validation. In view of that, a calibration process may be implemented to configure the numerical model within an acceptable range of parameters (see Figure 5.1).

Although various unsaturated subsurface flow models in porous media are available, the versatility and the efficiency of these modeling tools depend considerably on the quality of the identification process for unsaturated flow parameters.

Nevertheless, the parameters of soil-water retention curves which are the key functions required in hydrological, environmental and ecological modeling, are difficult to measure. Hence, several studies investigated how to predict these parameters from basic soil properties using Pedotransfer functions (PTFs). While different PTFs show good accuracy and reliability, many of them have limited validity [START_REF] Teixeira | Application of Soil Physics in Environmental Analyses: Measuring, Modelling and Data Integration[END_REF].

Identification of vGM model

As an alternative, the vGM model (Equation (4.12)) is still one of the most frequently used models in hydrological applications for the prediction of hydraulic conductivities applied for various types of soils. The characterization of vGM parameters is either based on their geometry or obtained by empirical estimation. Lambot et al. [2004] tested the inverse modeling method on laboratory-scale flow experiments conducted on artificial sand and three undisturbed soil columns. This provided relatively good estimates for the water retention curves and also for the saturated conductivity when the moisture range explored was not too small. [START_REF] Thoma | Estimating unsaturated hydraulic functions for coarse sediment from a field-scale infiltration experiment[END_REF] indicated that vGM relationships can be applied to these coarse, conglomeratic soils to characterize unsaturated flow behavior even under high sustained infiltration rates, and by using four parameters of vGM and implementing Metropolis-Hastings search with multiple independent-chains, the results were well fitted with observed data. Furthermore, [START_REF] Dittmer | Simulation of a subsurface vertical flow constructed wetland for cso treatment[END_REF] fitted manually the parameters by matching the curves of simulated and measured effluent rates as functions of time. [START_REF] Younes | Sensitivity and identifiability of hydraulic and geophysical parameters from streaming potential signals in unsaturated porous media[END_REF] conducted global sensitivity analysis and parameter estimation procedures to evaluate the influence of vGM and geophysical parameters on the Streaming Potential (SP) signals and studied the identifiability of these parameters from SP measurements. [START_REF] Langergraber | Modeling of processes in subsurface flow constructed wetlands: A reviewall rights reserved. no part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system[END_REF] provided a survey of existing simulation tools for SCWs ensuring that measured data from SCWs can be matched. Using HYDRUS simulation tool when the hydraulic behavior of the system is well described, it was concluded that the influence of the hydraulic parameters of the filter material is much higher compared to the influence of the biokinetic model parameters. [START_REF] Fournel | Modeling constructed wetlands with variably saturated vertical subsurface-flow for urban stormwater treatment[END_REF] implemented HYDRUS to explore specific features of SCWs and added an additional conceptual layer at the outlet of the wetland to mimic the local head loss resulting from seepage boundary condition. Some characterization of the hydraulic properties of the system is also described in [START_REF] Morvannou | Characterizing hydraulic properties of filter material of a vertical flow constructed wetland[END_REF], where the direct laboratory methods and the inverse modeling from in situ measurements were used. The results showed that the direct laboratory measurements failed to reproduce the water content obtained from the laboratory experiments, while the inverse modeling showed a better reliability in the simulation.

Carrying out a hydrodynamic characterization of the infiltration material, a set of parameters (initial values) can be produced to be used in the inverse optimization module included in HYDRUS-1D. This calibration methodology implemented by [START_REF] Morvannou | Characterizing hydraulic properties of filter material of a vertical flow constructed wetland[END_REF] reduced the risk of non-convergence of the model, since the values obtained are likely to be close to provided optimal values. As a deterministic inverse approach, a gradient method may be used due to its efficiency in comparison with gradientless methods. The main disadvantage of this method is that it may obtain local minimums. Also, it demands unimodality, continuity and differentiability of the objective function.

For many real problems, generally the data cannot be known with certainty due to measurement errors. Hence, contrarily to deterministic optimization in which the data for a given problem are assumed to be known accurately, in stochastic optimization, uncertainty may be considered in the model. [START_REF] Maier | Calibration of hydraulic parameters for large-scale vertical flow constructed wetlands[END_REF] used an inverse method considering a global optimization based on a stochastic search strategy in which hydraulic calibration was carried out on the outflow rate measurement. Thus, soil hydraulic functions based on vGM coefficients and preferential flow characteristics in large pores at high saturity were obtained. They also used the information from multiple restarts of the optimization algorithm to determine suitable parameter ranges and reveal potential correlations. Furthermore, "intelligent" algorithms including Genetic Algorithm (used for complex problems where traditional techniques can not be applied) and Particle Swarm Optimization (a set of particles search for good solutions to a given optimization problem) have a reliable ability of global optimization and provide accurate estimations. Contrarily to gradient-type methods, these methods are generally not dependent on initial values. For instance, [START_REF] Abbaspour | Estimating unsaturated soil hydraulic parameters using ant colony optimization[END_REF] applied a global search algorithm for estimating model parameters in a lysimeter experiments based on optimization by a colony of ants.

As uncertainties in unsaturated hydraulic conductivity function is significantly large, stochastic approaches such as Latin hypercube sampling may be taken to efficiently minimize the input parameter combinations. The feasibility of using Latin hypercube sampling method to estimate vGM parameters from time-lapse vertical electrical soundings during a constant inflow infiltration was investigated by [START_REF] Farzamian | Constraining unsaturated hydraulic parameters using the latin hypercube sampling method and coupled hydrogeophysical approach[END_REF]. [START_REF] Lambot | A global multilevel coordinate search procedure for estimating the unsaturated soil hydraulic properties[END_REF] carried out hydraulic parameters identification during natural infiltration events using both heuristic and stochastic global optimization methods by means of Global Multilevel Coordinate Search optimization algorithm in sequential combination with the local Nelder-Mead Simplex algorithm.

In order to escape local minimums in deterministic methods and to avoid costly computations in stochastic methods, Simulated Annealing approaches are provided. It can be considered as a stochastic version of the descent optimization method. Indeed, this algorithm, instead of taking steps along the gradient at a given point, takes steps stochastically. Simulated Annealing Algorithm is usually used to approximate near optimal solution related to hard optimization problems. Besides, [START_REF] Friedlander | Hybrid deterministic-stochastic methods for data fitting[END_REF], by analyzing rate of convergence, investigated hybrid methods that can provide the benefits of both approaches. It was shown that by controlling the sample size in an incremental-gradient algorithm, the steady convergence rates of full-gradient methods can be improved.

Parameter variability and the hysteresis effect

The VFCW hydraulic characteristics depend on several factors including age of the system, season and climate conditions such as periodicity of rainfall events. Although the filter has to be designed to accept hydraulic overload, but minimizing the surface ponding is a key point to ensure sufficient oxygen renewal for treatment. Dynamic models thus emerge as valuable tools for describing system hydraulics in response to various number of influencing factors and in order to understand better the filter hydraulics which may lead to improve the filter design. Simulation tools may be useful for sensitivity analysis or CW dimensioning and design, for instance. Modeling combined with field data may also allow for an evaluation of the filter aging.

Modeling issues arise when the soil reaches oven-dry conditions and runoff or ponding occurs. Hence, a particular attention should also be brought to boundary condition modeling (surface ponding or evaporation) to be able to tackle different sequences of rainfallrunoff events. The implementation of a virtual layer at the bottom of the column with a low saturated hydraulic conductivity K s has been proven to effectively simulate the calibrated output rate [START_REF] Fournel | Modeling constructed wetlands with variably saturated vertical subsurface-flow for urban stormwater treatment[END_REF].

The present study focuses on the long-term hydrodynamic characterization by implementing a numerical model of a vertical flow stormwater constructed wetland (VFSCW). The basic scheme of the studied VFCW is presented in Figure 5.2 whose design is inspired by wastewater treatment guidelines [Groupe [START_REF] Macrophytes | Traitement des Eaux. Epuration des eaux usées par filtres plantés de macrophytes -recommandations techniques pour la conception et la réalisation[END_REF].

The hysteresis phenomenon is generally due to several reasons cited by [START_REF] Feddes | Modelling soil water dynamics in the unsaturated zone -state of the art[END_REF] such as swelling and shrinkage for fine grained clays which may be resulted by wetting and drying or thermal effects [START_REF] Nimmo | Unsaturated Zone Flow Processes[END_REF]. Figure 5.3 shows a typical example of hysteretic water retention in a soil. The main outer drying and wetting curves correspond The variations of vGM parameters and thus, the behavior of unsaturated soils, are greatly be influenced by hydraulic hysteresis. This phenomenon describes the dependency of capillary pressure saturation curves on the history of flow. Indeed, the relation between permeability and water content can be a function of water and soil matrix properties, as well as relation between water content and water head which is considerably affected by hysteresis. Assuming singled-values functions in the approximation of the hydrodynamic conductivity in order to characterize the hydrodynamic properties, unsaturated flow models often are commonly applied without taking in account the hysteresis effect due to a lack of practical experimental methods for measuring water retention behavior along wetting paths. Seemingly such simplifications would be acceptable for many flow applications, but in various cases it may be necessary to consider hysteresis in the soil hydraulic properties to generate more realistic results.

The influence of hydraulic hysteresis on the behavior of unsaturated soils due to rainfall-induced landsliding, shrinkage and swelling of expansive soils, may have a detrimental effect on the variations of the parameters. If hysteresis in such analyses are not accounted for the analyses of the behavior of unsaturated soils, it may lead to significant errors in prediction of solute movement and contaminant concentrations and thus to discrepancies between predicted and observed results [START_REF] Kool | Development and evaluation of closed-form expressions for hysteretic soil hydraulic properties[END_REF]. This phenomenon generally is due to several reasons cited by [START_REF] Feddes | Modelling soil water dynamics in the unsaturated zone -state of the art[END_REF] such as swelling and shrinkage for fine grained clays which may be resulted by wetting and drying or thermal effects [START_REF] Nimmo | Unsaturated Zone Flow Processes[END_REF]. Figure 5.3 shows a typical example of hysteretic water retention in a soil. [START_REF] Likos | Hysteresis and uncertainty in soil water-retention curve parameters[END_REF] performed laboratory experiments in order to measure hysteretic soil-water characteristic curves for a wide range of soil types and then assessed uncertainty in three simplifications adopted to estimate wetting path parameters from easily measured drying curves.

Sensitivity analysis and calibration

The sensitivity analysis performed with respect to the vGM parameters (see Subsection 5.5.4) reveals a predominant influence of the shape parameters α, n and the saturated conductivity of the filter as well as of the virtual layer on the piezometric heads, during saturation and desaturation [START_REF] Fournel | Modeling constructed wetlands with variably saturated vertical subsurface-flow for urban stormwater treatment[END_REF]. Thus, the local context parameters used for the sensitivity analysis were reduced to those parameters that most heavily influenced the flow produced.

If hysteresis, which is physically caused by the presence of entrapped air, is not accounted in the analyses of the behavior of unsaturated soils, it may lead to significant errors in prediction of solute movement and contaminant concentrations and thus to discrepancies between predicted and observed results [START_REF] Kool | Development and evaluation of closed-form expressions for hysteretic soil hydraulic properties[END_REF]. Consequently, primary wetting and draining branches of soil-water characteristics may be considered in calibration process for both empirical [START_REF] Gillham | Hydraulic properties of a porous medium: Measurement and empirical representation 1[END_REF] and theoretical models [START_REF] Mualem | A modified dependent-domain theory of hysteresis[END_REF].

In the procedure of calibration, the quality of the final model parameter values could be affected by many factors, such as the quality of the input data, the simplifications and errors inherent to the model structure, the "power" of the optimization algorithm, the estimation criteria, the objective functions and so on [START_REF] Madsen | Parameter estimation in distributed hydrological catchment modelling using automatic calibration with multiple objectives[END_REF].

Furthermore, non-uniqueness of a solution in an inverse problem depends highly on the considered soil response, and consequently on the experimental design. Assuming that the used model describes the soil system efficiently and the parameters are identifiable, an inverse problem without enough field experiments will be ill-posed. Thus, for proper parameter identification, large field datasets would be needed. As these are usually not available, notably due to the randomness of the storm events, thus a simple, robust and low-cost numerical method for the inverse modeling of the soil hydrodynamic properties is proposed. Such robust optimization techniques can be used when the parameters are known within certain bounds. Thus, the goal is to find a solution that is rather feasible for all data and optimal.

The hydrodynamic modeling is carried out by implementing the Richards model by means of a MHFEM adapted to the simulation of heterogeneous media (see [START_REF] Younes | A new formulation of the mixed finite element method for solving elliptic and parabolic pde with triangular elements[END_REF] and Chapter 4), and the van Genuchten-Mualem (vGM) parametrization. A particular attention is brought to the top boundary conditions -surface ponding or evaporation -to be able to model the sequences of rainfall-runoff events. MHFEM results after parameter estimation are compared to those of HYDRUS [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 monitoring of the VFSCW provides us information in terms of water exchanges, filtering abilities during feeding and drainage sub-periods, and ageing along time. Large datasets are needed to account for significant spatial and temporal variability in the evaluation of the variability of the vGM soil parameters hidden behind the hysteresis effect. On the one hand, this study relies on the automatic differentiation (AD) of the MHFEM code for sensitivity analysis and gradient-based data assimilation. On the other hand, a stochastic method is implemented to, eventually, determine a global optimum. The modeling efficiency, as defined by [START_REF] Moriasi | Model evaluation guidelines for systematic quantification of accuracy in watershed simulations[END_REF], and a statistical analysis are then evaluated for the different parameter sets to demonstrate the hysteresis effect.

Description of the study area

The field data were collected for a three-layered CW located in Strasbourg (Alsace, France) at the water edge of the urban water stream Ostwaldergraben, during several months. The total surface of the catchment area is 18000 m 2 in which the areas of total roofs and roads represent 13.8% and 15.6.% of the total area. The catchment is exclusively residential and drained by a sewerage network (see Figure 5.4). Before 2012, the water from the storm-water system was discharged directly into the water stream without undergoing any type of treatment process. This was one of the main causes of the poor quality of water of the Ostwaldergraben.

In order to improve the water quality of the water stream, a constructed wetland sys- tem designed by the Urban Community of Strasbourg * (CUS) and in partnership with the engineering science, computer science and imaging laboratory (ICUBE ) has been planned by installing a storm-water treatment system consisting of a sedimentation pond and a filtration media (Figure 5.5). The reader who wants more details concerning the SCW instrumentation could refer to [START_REF] Walaszek | Urban stormwater treatment by a constructed wetland: Seasonality impacts on hydraulic efficiency, physico-chemical behavior and heavy metal occurrence[END_REF].

The main characteristics of the Ostwaldergraben catchment are provided in Table 5.1.

Settling pond

The water in the pond comes only from rainwater, runoff, leaching and snow melt. The pond is in the form of a pear and is preferred over the rectangular shape to avoid areas of water stagnation. The pond contains two overflow sections in order to prevent the overflow of filtration media in the case of rising water level resulted from storm water. A beam is installed at the entrance of the pond functioning as a flow dispersion device to dampen the water force and thus to prevent erosion (see Figure 5.6).

Filtration bed

Schematically, the main components of a biofilter are the ponding zone, the filter media and the drainage layer, (see Figure 5.6). Water enters the system from the top, infiltrates through the filter media and starts percolating into the drainage layer when soil moisture reaches field capacity. Water in the drainage layer can infiltrate to the surrounding soil (exfiltration) or can be collected by a perforated under drain pipe [Melbourne [START_REF] Water | WSUD engineering procedures : stormwater[END_REF] and directed to a receiving water body or the urban drainage system. Because the catchment area draining to the biofilter is much larger than that of the biofilter itself (commonly approximately 50 times larger) [START_REF] Wong | Australian Runoff Quality: A Guide to Water Sensitive Urban Design[END_REF], the level of saturation of the filter media often becomes greater than field capacity. In these cases, water flows through the filter, generating outflow that can be seen during the same day as the inflow events. When inflow rates are larger than the ability of the filter to accommodate or percolate water, a pond will form on top of the system. The depth of this pond is limited (typically 0.2-0.4 m), thus, part of the inflow might overflow and bypass the system without being treated. Between inflow events, the water stored in the filter media is mainly lost via evapotranspiration.

Precipitation during the calibration periods

Datasets of precipitations and piezometric heads were measured during several months (from 16/04/2013 to 18/09/2013) which represent 10.5% of monitored days. For a better understanding of the rainwater treatment systems, the data collected by the sensors are analyzed statistically to deduce information on the cumulative water level, the dry period and the average intensity (see Table 5.2). In particular, the time series show that the groundwater level is significantly impacted by the season. Figure 5.7 shows the rain events characteristics. Overall, 48 rain events have been simulated during this study. In particular, one can notice that both the average intensity and the dry period have the highest variabilities (209% of variation coefficient) and the lowest variability is related to the water depth (150%). Besides, regarding the monthly rain characteristics analysis, the dry period during September has the highest variability (275% of variation coefficient). Detailed characteristics of sampled rainfalls are provided in Table F.1 in Appendix F.

Feeding device

Generally, in an unsaturated zone with a downward vertical flow, infiltration occurs near the feed point using only a small part of the infiltration bed. Consequently, higher concentrations are observed at the beginning of the events. In order to avoid the amplification of this phenomenon, special attention has been taken in the design of the distribution systems at the surface of the filtration bed for low flow rates and high flow rates.The system is directly fed with incoming raining events, creating a random pulse.

Considering the incoming flow, it can be noticed that the sedimentation pond and the filter act as moderators at the arrival of a rainy episode and ensure a smooth flow. For example, for particularly intense episodes, for a flow rate of 200 m 3 /h at the inlet, a flow rate of 40 m 3 /h at the floating weir and a flow rate of less than 10 m 3 /h at the outlet of the filtered was observed. Thus, the treatment system makes it possible to avoid shock effects in the receiving environment. In addition, this means that the floating weir is immediately reactive to the arrival of a rainy episode while the filter is slower to respond.

Transpiration rate

In general, moisture loss must be accounted in all hydrologic models including in a wetland in which the rate of removal of water stored within the wetland sediments is influenced by evapotranspiration . Thus, the experiment to estimate the transpiration rate was performed on 23rd July 2013 by monitoring the water level decrease in the filter after a dry period of 20 days and 25th July, 2013 in the field, following a day with a short and intense rainy episode. Considering weather conditions in July 2013 (mean temperature 22.1 °C, varying between 15.9 °C and 36.4 °C, total sunshine period during the month 294.1 h and cumulative rainfalls 9.7 mm). The measured transpiration value (approximated as 6.6 mm d -1 ) is probably overestimated the real mean rate value. As results obtained are not significantly different, it can be assumed that the extended dry weather period did not influence the transpiration rate of the filter. Yet evapotranspiration contributed to less than 1% of water loss in the system for all events [START_REF] Schmitt | Constructed wetlands treating stormwater from separate sewer networks in a residential strasbourg urban catchment area: Micropollutant removal and fate[END_REF].

The amount of biomass present on the filter planted with reeds calculated as 389 gDM/m 2 , is low compared to the data found in the literature for reed-planted filters on a real scale. However, at the time of the experiment, the system was only one year old, compared to the ones with the age of 10 years or more sectors found in the literature.

Calibration methodology

The conceptual layer always remains saturated, hence only the calibration of its saturated conductivity is required.

HYDRUS inverse method implementation

The HYDRUS 1D model may be used as a mechanistic model to characterize the soil hydrodynamic parameters in SCW [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]. HYDRUS 1D implements a Marquardt-Levenberg algorithm for the estimation of soil hydrodynamic parameters from measured transient flow from the data measured in the laboratory and the field data.

Model evaluation techniques

The agreement between a numerical model on the one hand, and a dataset on the other hand may be evaluated through error computations. Two model evaluation techniques, which are widely used for calibration and validation of hydrological models with observed data, are considered.

I) The mean squared error: The normalized mean squared error (MSE) is a criteria most widely used for calibration and evaluation of hydrological models with observed data.

MSE = 1 n . n ∑ t=1 (h t Sim -h t Obs ) 2 ∈ (0, +∞), (5.1)
where h t Sim and h t Obs are respectively the simulated and the observed values at n observation times. In a calibration, MSE is subject to minimization. It can be expressed in the units of the quantities that it compares when a unique variable is considered. However, this absolute dimensional criterion has major drawbacks. First, it cannot fairly compare MSEs computed for different variables (heads and velocities, for instance) as units are different. Second, it does not fairly compare MSEs for storm event datasets of different magnitude.

Note that a relative dimensionless error such as can be used.

rMSE = ∑ n i=1 (h i -h o i ) 2 ∑ n i=1 (h o i ) 2
∈ (0, +∞).

(5.2)

II) The Nash-Sutcliffe efficiency: The so-called Nash-Sutcliffe efficiency (NSE) is a dimensionless value obtained by dividing MSE by the variance of the observations and subtracting that ratio from 1.0. Obviously, it can be seen that NSE and MSE are inversely related:

NSE = 1 - ∑ n t=1 (h t Sim -h t Obs ) 2 ∑ n t=1 (h t Obs -µ Obs ) 2 = 1 - MSE σ 2 Obs ∈ (-∞, 1), (5.3)
where µ Obs is the mean and σ Obs is the standard deviation of the dataset.

A value of 1 for NSE indicates a perfect agreement between the model and the dataset. Guidelines were established for model evaluation based on the review results and projectspecific considerations. The interested reader is referred to review papers such as [START_REF] Harmel | Cumulative uncertainty in measured streamflow and water quality data for small watersheds[END_REF] and [START_REF] Moriasi | Model evaluation guidelines for systematic quantification of accuracy in watershed simulations[END_REF]. Reported values for model calibration and model validation are synthesized in [START_REF] Moriasi | Hydrologic and water quality models: Performance measures and criteria[END_REF].

In this chapter, MSE is used for calibration while model efficiency is evaluated from estimated vGM parameters through a NSE computation. Model performance is rated as "very good" if 0.75 < NSE < 1.00, as "good" if 0.65 < NSE ≤ 0.75, as "satisfactory" if 0.50 < NSE ≤ 0.65, and "unsatisfactory" otherwise.

Inverse modeling

The general inverse problem defined as Find p * ∈ P such that min

p∈P (J • M (p)) = (J • M )(p * ) (5.4)
is adapted to the paper's objective. The unknown vGM parameters p are searched in the set of admissible values P. The set that minimizes Equation (5.4) is denoted by p * . The objective function J • M is a compound function that applies MSE to compare the piezometric data to the simulation results of MHFEM computed by using the parameter set p. Such a minimization problem can be solved using either a gradient method or a stochastic one, depending on the abilities of the software or the availability of its source code. A basic stochastic method is also implemented in an attempt to find out the global minimum.

Sensitivity analysis

Among the set of modeling parameters, critical ones may be determined from a sensitivity analysis in order to downsize the number of unknowns for the inverse problem. From a computer point of view, a sensitivity analysis may be carried out through automatic differentiation (AD) explained in 3.3.

Algorithmic parameter estimation chart

The proposed parameter estimation methodology is summarized in Figure 5.8. The top level (red) tree is concerned with key modeling hypothesis set on the initial and boundary conditions, the mesh discretization, the soil hydrodynamic parameterization (the vGM model, for instance) and the computer hydraulic code (MHFEM or HYDRUS, for instance).

The green arrows of the chart are concerned with code differentiation and usages. Automatic differentiation with respect to vGM parameters is applied to the MHFEM code of Chapter 4 as described in Chapter 3. The resulting code is used for sensitivity analyses in order to decide on the most sensitive parameters, then in the identification by a gradient method of the vGM parameter values with respect to a given piezometric head dataset.

The blue loop organizes the temporal estimation of vGM parameters along the storm event sequence. It implements a deterministic gradient method as well as a stochastic one to avoid being trapped in local minima. Attention is paid to the quality of the parameter estimation by evaluating the model efficiency criterion (NSE).

Initially, the calibration process is performed for multi-rainfall event periods ignoring the stochastic method in order to avoid costly computations. Considering two specified groups of feeding and drainage sub-periods calibration, both calibration methods are applied for each sub-period. Finally, the temporal variability of vGM parameters during these periods are evaluated and discussed.

Initial and boundary conditions

Initial conditions : The initial pressure distribution in the field is set as a hydrostatic distribution which is obtained from the pressure head measured initially by a piezometer in the filter. Thus, the initial piezometric head at the lower and upper parts are approximated.

Top boundary conditions : Drainage occurring after every feeding event can improve the filter bed aeration and the aerobic degradation during the dry period. Modeling issues arise when the soil reaches oven-dry conditions and a special attention should be given to the switching boundary condition procedure, which is important for situations involving ponded water layers or fluctuating groundwater close to the soil surface. This procedure switches from head (dirichlet) to flux (Neumann) controlled boundary condition and vice versa depending on the saturation of the soil column at the top. The algorithm from [START_REF] Van Dam | Numerical simulation of infiltration, evaporation and shallow groundwater levels with the richards equation[END_REF] was adapted to be used with the mixed hybrid formulation and 2D flow problem [START_REF] Padilla | Modeling and Optimization of Biological Mitigation Processes in Porous Matrices Contaminated by Pesticides: Towards a New Functionality of Stormwater Basins[END_REF] and was improved as illustrated in Figure 4.4 in Chapter 4.

Lower boundary condition : At the bottom of the drainage layer, a "seepage" face boundary condition is chosen. This corresponds to the effectiveness of a flux at the lower boundary when the simulated piezometric head is greater than a threshold limit value.

Figure 5.9 displays the conceptual numerical domain, the different layers that have been discretized by triangle elements and boundary conditions. From top to bottom: 30 cm of sand layer; 34 cm of intermediate layer (gravel: 4-8 mm) and drainage layer (gravel and stone: 16-22.4 mm). Two aeration pipes installed between the upper part and the intermediate layer and also at the bottom of the sedimentation pond contribute to the aerobic conditions. As said previously, extended dry periods can potentially have an impact on macrophytes or on microorganisms, leading to potential treatment efficiency decrease. To avoid this, 16 cm of saturated zone has been lied at the bottom of the filter [START_REF] Schmitt | Constructed wetlands treating stormwater from separate sewer networks in a residential strasbourg urban catchment area: Micropollutant removal and fate[END_REF]. Observations of matrix pressure head at the elevation -64 cm were the experimental data used for inverse analysis of water flow experiments.

The computational domain Ω is discretized using a mesh containing 160 triangular elements with the same size (see Figure 5.10).

Model set-up and initial parameters

Practically, in a deterministic calibration, the influence of the guess for the initial values of the parameter on the final calibrated parameters cannot be neglected. Thus, choosing a suitable initial parameter set makes the calibration process converge quickly and leads to accurate parameters.

The hydraulic parameters of the different layers are determined from the model of vGM (Equation (4.12)). In total, seven parameters allow to describe the three hydraulic functions of the massif filter. The residual water content (θ r ), the saturated water content (θ s ) and the saturated hydraulic conductivity (K s ) have a simplified physical representation [START_REF] Ippisch | Validity limits for the van Genuchten-Mualem model and implications for parameter estimation and numerical simulation[END_REF]] called also as shape parameters (α v and n v ) are determined by adjustment.

A priori assumptions on vGM modeling parameters

The hydraulic conductivity through isotropic porous media may be approximated using the vGM parameters and the tortuosity matrix as described in Subsection 4.4.1.4. A set of seven parameters is used to model each layer of the filter.

Orders of magnitude for θ r and θ s are available in the literature for different kind of materials that can be implemented in a VFCW. Table 5.3 reports parameter values for the sand layer [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF] and the fine gravel layer [DanHan- [START_REF] Danhan-Cheng | Capillary effect on flow in the drainage layer of highway pavement[END_REF]. It also reports values for the third, virtual, layer [START_REF] Fournel | Modeling constructed wetlands with variably saturated vertical subsurface-flow for urban stormwater treatment[END_REF]. On the other hand, the shape parameters α v and n v are generally identified from data [START_REF] Ippisch | Validity limits for the van Genuchten-Mualem model and implications for parameter estimation and numerical simulation[END_REF]. These are used as initial guesses in the identification process.

Assumptions are usually made to reduce the number of modeling parameters to be identified. Following [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF], the matrix tortuosity parameter τ is set to 0.5 as an average value, while the air entry value h e is set to 0. The sensitivity analysis results presented in Subsection 5.5.4 allow for downsizing the set of parameters subject to identification.

Sensitivity analysis and temporal parameter estimation

The sensitivities of the MHFEM model to the changes in the various vGM parameters of the first layer and the saturated hydraulic conductivity of the last layer were examined. To that end, a perturbation of 50% was applied to the vGM parameters related to the first layer and to the hydraulic conductivity of the last layer.

Figure 5.11-a represents the sensitivity of water pressure to vGM parameters timeseries for the 1st period of May. The parameters identified as having the greatest influence are : the saturated water content θ sat and the shape parameterα for the first layer, the hydraulic conductivity K sat for the last layer. The parameters n, K sat and θ res for the filter layer can modify the water pressure profile, particularly during each rainfall event, thus, their influence cannot be ignored. For the sake of simplicity, Figure 5.11-b is presented as a zoom and simplified plot of a rainfall event. The temporal intervals in which each parameter can be estimated, are highlighted with different patterns. It has to be remarked that the major influence of the hydraulic conductivity of the virtual layer is revealed due to its high sensitivity comparing to the one for the first layer. This corresponds, physically, to the the role of drainage pipe. The sensitivity results with respect to those parameters that are not considerable (the sensitivity of water pressure with respect to the hydraulic conductivity of the intermediate layer is zero, for instance) are not presented here. Regarding the sensitivity analysis which highlighted the parameters with major influence, specifically the hydraulic conductivity for the virtual layer and considering the fact that the variations of properties due to the plant activities (decay, growth and death) and suspended solid accumulation are expected only for the first layer, the values for the parameters of the second layer are fixed through the calibration process.

Parameter estimation and model efficiency

Calibration of the sensitive parameters is obtained through inverse modeling implementing a deterministic (gradient) and a stochastic (simple random processing) method by using the experimental measures during 48 rainfall events. Consequently, the temporal estimation of the considered parameters and the efficiency of the two methods are discussed.

The calibration results are provided in Figures 5.12 to 5.17 which are categorized in three groups depending on the precipitation occurrence: multi-rainfall event periods (Subsection 5.6.1), feeding sub-periods (Subsection 5.6.2) and drainage sub-periods (Subsection 5.6.3). These figures can highlight the discrepancies between the calibrated model and the observed data qualitatively. The values of the corresponding vGM optimized parameters are provided in Tables 5.4 to 5.9.

Parameter estimation for multi rainfall events

Tables 5.4 and 5.5 show clearly that several parameter sets obtained for multi rainfall events provide unsatisfactory NSEs. Hence, these periods are divided into 20 feeding and drainage calibration sub-periods as illustrated in Figures 5.12 and 5.13. Table 5.4: Multi rainfall events. Optimized parameters by MHFEM (gradient method) and NSEs: very good (green), good (blue), satisfactory (orange), unsatisfactory (red).

α 1 n 1 θ r1 θ s1 K 1 K 3 MSE NSE Month (cm -1 ) (-) (-) (-) (cm.

Parameter estimation during feeding sub-periods

Feeding periods are, here, defined as temporal ranges corresponding to piezometric head data series with increasing value. Figures 5.14 and 5.15 plot observed and piezometric heads computed from estimated parameters reported in Tables 5.6 and 5.7. Among the eleven simulated feeding sub-periods, six optimized set of parameters provide good and very good NSEs. One observes in Figure 5.14-c that the model predicts drainage subperiods, mostly when runoff stops. This disagreement is monitored by NSE computations, see Table 5.6. This indicates that the calibration process -modeling, inverse modeling and gradient method -often fails at identifying the parameters. A well known issue is that a gradient method can be trapped in a local minimum. However, the stochastic inverse method involving 1000 parameter sets did not succeed in a better manner, for any of these events. In the present case, the failure explanations are mainly within the modeling, and 0 9 / 0 0 9 / 0 2 0 9 / 0 3 0 9 / 0 4 0 9 / 0 5 0 9 / 0 6 0 9 / 0 7 0 9 / 0 8 0 9 / 0 9 0 9 / 1 0 0 9 / 1 1 0 9 / 1 2 0 9 / 1 3 0 9 / 1 4 0 9 / 1 5 0 9 / 1 6 0 9 / 1 7 0 9 / 1 8 notably the account for ponding and infiltration that depend on the saturation of the soil on the top of the filter. Other assumptions can also interfere. For instance, the threedimensional behavior of the water flow may be predominant during the feeding process since surface flow may occur near to the top layer.

Obviously, no trend on the evolution of the parameter value can be deduced from Table 5.6 since some of the NSE values are less than 0.5. Crossing the NSE values from Table 5.6 with the rain event characteristics, the lowest NSE values are obtained only during May. These events are remarkably identified by the highest average water depth and the highest average intensity (see Figure 5.7). Comparing these results (see Table 5.6) with those obtained during the multi-rain events (see Table 5.4), an improvement of the calibration is obtained by a global increase in the NSE values for some events. Drainage sub-periods are, here, defined as temporal ranges corresponding to piezometric head data series with decreasing value. Modeling parameters identified for a total of ten different drainage sub-periods are reported in Tables 5.8 and 5.9. Their variations along time can be better observed in Figure 5.19-I for the gradient method (blue ▽) and the stochastic method (red ).

Hydraulic conductivities

The virtual layer hydraulic conductivities (K 3 ) at the bottom of the filter provides similar trends for both the gradient and stochastic methods. From April to September, K 3 was decreasing, suggesting that the flowing capacity of water through the virtual layer is reducing along the months (Figure 5.19-I-f). This phenomenon might be due to the migration of fine particles, suspended solid or root debris from the top layer to the bottom that house the drainage pipe. The trend of the top layer hydraulic conductivity (K 1 ) is completely different to the third layer one (see Figure 5.19-I-e). There is no evidence of a reduction in infiltration capacity of the top layer which is likely to accumulate particles as well as the root system. In fact some phenomena with opposite impact on the top layer hydraulic conductivity might occurred along the months. The filtration process on the top layer might reduce the infiltration capacity but stalling due to shear stresses during water flow through porous media for high rain event intensity as well as the seasonal dynamic of vegetation could increase the permeability.

Water content Residual and saturated water content displayed opposite trend (see . Process dynamic previously highlighted inside the planted porous media could explain this result. While residual water content is closely associated with the capillary forces due to accumulated fine particles and roots dynamic, saturated water content should be impacted by porous media reorganization due to plant growth, settlement or swelling of sand or gravel layers.

vGM shape parameters Even if gradient and stochastic methods provide similar trends for top layer α and n parameters, there is no clear explanation for their variabilities. Physically, the drainage process occurs through vertical flow in porous media when gravitational forces are higher than capillary ones. Hence the porous media changes from near saturated condition characterized by the lack of shape parameters toward unsaturated one that is sensitive to the shape parameters. Furthermore, it can be noticed that the parameter α provides the highest variability of estimated parameters regardless the numerical method (see Tables 5.8 and 5.9). Hence it will be clearly difficult to catch the meaning of this parameter along multi-rainfall event.

Several conclusions arise. First, model efficiencies (NSEs) are evaluated as very good for all these drainage events, whatever the inverse method is. Second, parameters identified by one or the other inverse method are generally close, but some differences may exist even for NSEs above 0.95. This assesses the temporal variability of top layer vGM parameters from one drainage period to another.

Hysteresis effect

The variability of the return periods of the events accentuates the phenomenon of hysteresis (see Figure 5.7) since the saturation condition of the porous media could be different before each rainfall event (see Figure 5.3). Hence, the ability of the model to simulate a large range of event will be assessed. Results discussed in the previous paragraphs highlight the fact that MHFEM model partly failed to estimate properly the vGM parameters for hydrodynamic characterization through vertical flow stormwater constructed wetland when considering the multi-rain events. After dividing the water flow process with feeding and drainage sub-periods, better fits are obtained. The temporal variability of vGM parameters during all sub-periods is presented in Figure 5.18. Note that including those estimations which are rated as satisfactory and unsatisfactory may make the interpretation difficult and thus we had to ignore them in the investigation.

In this paragraph, differences with estimated parameters are deeper discussed in order to understand the hysteresis effect regarding the evolutions of parameters which are calibrated in a very good or good manner.

For very good parameter identifications (NSE > 0.75), a Principal Component Analysis (PCA) for the optimized vGM parameters is carried out and presented in Figure 5.20. The individuals are the sub-periods (feeding and drainage) and the variables are the estimated parameters. The first principal component axis (33,3%) is negatively well-correlated with the saturated hydraulic conductivity of the first and third layers. The second principal component axis (26.1%) is in a negative correlation with the shape parameters α and n (see Figure 5.20). Indeed, the bigger the residual water content value, the smaller the α and n values. Sand material has clearly small n and α values, knowing that the α value is related to the inverse of the air entry pressure value and the n value depends on the width of the pore size distribution. These evidences are in agreement with the results provided by the PCA.

The feeding sub-periods and the drainage ones are significantly separated (see Figure 5.20). The PCA highlights that the estimated parameters of the feeding sub-periods are clearly different from those of the drainage sub-periods due to the hysteresis effect. The NSE values without considering the feeding and drainage sub-periods would be reduced enormously.

In addition, the Shapiro-Wilk test proves that among all the estimated parameters, the saturated hydraulic conductivity of the first layer (w = 0.8837, p-value= 0.05384), the saturated hydraulic conductivity of the third layer (w = 0.9077, p-value= 0.1249) and the saturated water content (w = 0.9588, p-value= 0.6720) are normally distributed. The homogeneity of variances are verified by the Levene's test. Thus, the analysis of variance explains that the influence of the feeding and drainage sub-periods on:

• the estimated value of the saturated water content is significant (Pr(> F) = 0.0095); • the first layer saturated hydraulic conductivity is considerable (Pr(> F) = 0.0010); • the last layer saturated conductivity is not significant because the outlet boundary condition governs the hydrodynamic process in the third layer.

Moreover, a temporal estimation of parameters in a very good manner during feeding and drainage sub-periods can bring out useful information about long term SCW hydrodynamic variations. Figure 5.19-II demonstrates the temporal variations of vGM parameters during these feeding (blue ) and drainage (blue ▽) sub-periods.

Conclusions

This study discusses numerical approaches to investigate the temporal variability of vGM parameters in a VFSCW during multi-rainfall periods considering the hysteresis effect. The results show that the characteristics of the filter layer alter along time. Applying a model, several issues such as water accumulation due to biomass and plants growth, presence of organic matter, and the relying of total suspended solids in this layer have to be considered as the main causes of these changes. The pressure head curves provided by simulations using the calibrated parameters sets are not perfectly in accordance with observed data and the NSE evaluations for some periods, particularly those with the high average rain event intensity, were not ranked as satisfactory. Thus, the study had to be focused on the impact of the hysteresis effect on the calibration process.

The estimated parameters during both drainage and feeding sub-periods which are provided by gradient and stochastic optimizations, demonstrate dominantly the same temporal variations. However, a few discordances between the parameters identified by the two implemented methods may be strongly due to either the risk of trapping in local minima using a gradient method or an insufficient number of random parameters processed by the stochastic calibration. Furthermore, several optimizations during feeding sub-periods are unsatisfactory using either MHFEM or HYDRUS. Both are not able to properly reproduce some of the observed wetting periods. In addition, errors in observed data may also produce these mismatches. Carrying out a statistical analysis for the different parameter sets showed that the estimated parameters of the feeding sub-periods are clearly different from those of the drainage sub-periods which evidences the hysteresis effect. CHAPTER 6

Conclusions and perspectives

Sensitivity analyses and parameter identification for water flow through variably saturated porous media are applied to two remediation facilities. Both studies are carried out using automatic differentiation, as an efficient tool for computing the derivatives and Jacobians of functions, performed in tangent linear mode.

In Chapter 2, the influence of drainage on the primary variables and the streamlines of particles was studied for an aquifer with both homogeneous and heterogeneous hydraulic conductivity in a conforming FEM and a MHFEM implementation for the drain-aquifer interactions. It was clearly demonstrated that for both numerical methods the fluxes approximated at all the boundaries, including the drains, have similar behavior but distinct values due to the difference in the modeling of water-body exchanges.

A sensitivity analysis was carried out in Chapter 3. For further investigations in the two discretization methods dealing with the drainage process in the aquifer, an important conclusion from this sensitivity analysis is that the dependency of the state variables on the leakage coefficient in MHFEM is more significant than in FEM. Besides, the drains in FEM catches more water than ones in MHFEM. To the best of our knowledge, streamline reconstruction results presented in the thesis are new in groundwater flow modeling.

The Chapter 4 describes different techniques implemented into the MHFEM in a variably saturated medium to improve the hydrodynamic modeling. In particular, a promoted auto-selective boundary condition procedure was needed to be able to model transient flow regimes switching from Dirichlet to Neumann boundary condition and inverse.

In chapter 5, once the proper technique for MHFEM was chosen, the differentiated numerical code was adapted to be implemented in the optimization process. Several sets of vGM parameters were identified using pressure head data. Those parameter sets which are calibrated in a very good manner provide information to investigate their temporal variability in a VFSCW along time. The important role of the hysteresis effect is highlighted by distinguishing the feeding and drainage sub-periods. A principal component analysis for the different parameter sets showed that the estimated parameters of the feeding subperiods are clearly different from those of the drainage sub-periods, which evidences the hysteresis effect. The results show that the characteristics of the filter layer alter along time which can be explained by the fact that applying a model, several issues such as water accumulation due to biomass and plants growth, presence of organic matter, and the relying of total suspended Chapter 6 -Conclusions and perspectives solids in this layer have to be considered.

As a further work, efforts on the validation of the procedure of switching boundary condition have to be continued particularly for the extremely high rainfall events. Although the model is calibrated during several periods, the model may be validated through the prediction of the evolution vGM parameters during several years by referring to some extra collected data from experimental site. As automatic differentiation has proven its efficiency for sensitivity analysis and parameter estimation for hydrodynamic modeling, in addition, the same processes for modeling of the flow coupled with reactive transport will be carried out to obtain information in terms of ageing of the filter along time.

l 1 (0, 0) F 2D --→ L(x L1 , z L1 ) l 2 (1, 0) F 2D --→ L(x L2 , z L2 ) l 3 (0, 1) F 2D --→ L(x L3 , z L3 )
The analytic epressions of x(u, v) and z(u, v) are:

x (u, v) = (N l 1 , N l 2 , N l 3 )   x L 1 x L 2 x L 3   , z (u, v) = (N l 1 , N l 2 , N l 3 )   z L 1 z L 2 z L 3  
where the base function of interpolations are equal to the geometric transformation functions:

N l 1 = m Ĝ,l 1 = 1 -u -v , N l 2 = m Ĝ,l 2 = u , N l 3 = m Ĝ,l 3 = v
So then the nodal coordinates x(u, v) and z(u, v) and can be expressed as:

x (u, v) = (1 -u -v) x L 1 + (u) x L 2 + (v) x L 3 = x L 1 + ua + vb z (u, v) = (1 -u -v) z L 1 + (u) z L 2 + (v) z L 3 = z L 1 + uc + vd where: a = x L2 -x L1 , b = x L3 -x L1 c = z L2 -z L1 , d = z L3 -z L1
A vector transformation can be written as:

→ w i (x, z) = J → ∧ w i (u,v) det J
The Jacobian matrix associated to the transformation F2D is defined by:

J 2D = ∂ x ∂ u ∂ x ∂ v ∂ z ∂ u ∂ z ∂ v = x L 2 -x L 1 x L 3 -x L 1 z L 2 -z L 1 z L 3 -z L 1 = a b c d
and its determinant, denoted as det J = adbc.

The base functions associated to the edges E i are defined over ∧ G:

→ ∧ w 1 = u v -1 , → ∧ w 1 = u -1 v , → ∧ w 1 = u v .
Therefore, the vector transformation can be written as:

→ w i (x, z) =   J 11 ∧ w iu +J 12 ∧ w iv det J J 21 ∧ w iu +J 22 ∧ w iv det J  
where J i j is the i j component of the Jacobian matrix J.

w 1 =   ∧ w 1u J 11 + ∧ w 1v J 1,2 det J ∧ w 1u J 21 + ∧ w 1v J 2,2 det J   = ua+bv-b ad-bc uc+vd-d ad-bc → w 2 =   ∧ w 2u J 11 + ∧ w 2v J 1,2 det J ∧ w 2u J 21 + ∧ w 2v J 2,2 det J   = (u-1)a+bv ad-bc (u-1)c+vd ad-bc → w 3 =   ∧ w 3u J 11 + ∧ w 3v J 1,2 det J ∧ w 3u J 21 + ∧ w 3v J 2,2 det J   = ua+bv ad-bc uc+vd ad-bc
A.2 Transformation from an element in the physical space to a reference element

The transformation from an element in the physical space to a reference element by an integral of a scalar function can be expressed as

G f (x, z)dxdz = Ĝ f (x (u, v) , z (u, v)) det Jdvdu
Thus, the components of the matrix B G can be calculated as

B i j = G w ix L xx w jx + w iz L zz w jz dxdz = Ĝ J 11 ∧ w iu +J 12 ∧ w iv L xx det J J 11 ∧ w ju +J 12 ∧ w jv det J + J 21 ∧ w iu +J 22 ∧ w iv L zz det J J 21 ∧ w ju +J 22 ∧ w jv det J det Jdvdu = 1 u=0 1-u v=0 w ix kx w jx + w iz kz w jz det Jdvdu B 11 = 1 u=0 1-u v=0 w 1x kx w 1x + w 1z kz w 1z det Jdvdu B 11 = 1 12 3b 2 k z +3d 2 k x -3abk z -3cdk x +a 2 k z +c 2 k x k z k x (ad-bc) B 12 = 1 u=0 1-u v=0 w 1x kx w 2x + w 1z kz w 2z det Jdvdu B 12 = -1 12 b 2 k z +d 2 k x -3abk z -3cdk x +a 2 k z +c 2 k x k z k x (ad-cb) B 13 = 1 u=0 1-u v=0 w 1x kx w 3x + w 1z kz w 3z det Jdvdu B 13 = -1 12 b 2 k z +d 2 k x +abk z +cdk x -a 2 k z -c 2 k x k z k x (ad-cb) B 21 = 1 u=0 1-u v=0 w 2x kx w 1x + w 2z kz w 1z det Jdvdu B 21 = -1 12 b 2 k z +d 2 k x -3abk z -3cdk x +a 2 k z +c 2 k x k z k x (ad-cb) B 22 = 1 u=0 1-u v=0 w 2x kx w 2x + w 2z kz w 2z det Jdvdu B 22 = 1 12 b 2 k z +d 2 k x -3abk z -3cdk x +3a 2 k z +3c 2 k x k z k x (ad-cb) B 23 = 1 u=0 1-u v=0 w 2x kx w 3x + w 2z kz w 3z det Jdvdu B 23 = 1 12 b 2 k z +d 2 k x -abk z -cdk x -a 2 k z -c 2 k x k z k x (ad-cb) B 31 = 1 u=0 1-u v=0 w 3x kx w 1x + w 3z kz w 1z det Jdvdu B 31 = -1 12 b 2 k z +d 2 k x +abk z +cdk x -a 2 k z -c 2 k x k z k x (ad-cb) B 32 = 1 u=0 1-u v=0 w 3x kx w 2x + w 3z kz w 2z det Jdvdu B 32 = 1 12 b 2 k z +d 2 k x -abk z -cdk x -a 2 k z -c 2 k x k z k x (ad-cb) B 33 = 1 u=0 1-u v=0 w 3x kx w 3x + w 3z kz w 3z det Jdvdu B 33 = 1 12 b 2 k z +d 2 k x +abk z +cdk x +a 2 k z +c 2 k x k z k x (ad-cb)
and the matrix inverse B -1

G deno = a 2 k z -abk z + d 2 k x -cdk x + b 2 k z + c 2 k x (ad -bc) B -1 K11 = 2 -a 3 bk z 2 + c 2 d 2 k x 2 -c 3 dk x 2 + a 2 b 2 k z 2 + 3b 2 c 2 k x k z + 3a 2 d 2 k x k z +2a 2 c 2 k x k z + a 4 k z 2 + c 4 k x 2 -4abcdk x k z -abc 2 k x k z -a 2 cdk x k z /deno B -1 K12 = -2 2a 3 bk z 2 -c 2 d 2 k x 2 + c 3 dk x 2 -a 2 b 2 k z 2 + ab 3 k z 2 + cd 3 k x 2 -2b 2 c 2 k x k z -2a 2 d 2 k x k z + b 2 cdk x k z + abd 2 k x k z + 2abcdk x k z + abc 2 k x k z + a 2 cdk x k z /deno B -1 K13 = -2   -2a 3 bk z 2 + 2c 2 d 2 k x 2 -2c 3 dk x 2 + 2a 2 b 2 k z 2 -ab 3 k z 2 -cd 3 k x 2 -b 2 c 2 k x k z -a 2 d 2 k x k z + 2a 2 c 2 k x k z + a 4 k z 2 + c 4 k x 2 -b 2 cdk x k z -abd 2 k x k z + 6abcdk x k z -2abc 2 k x k z -2a 2 cdk x k z   /deno B -1 21 = -2   a 3 bk z 2 -c 2 d 2 k x 2 + c 3 dk x 2 -a 2 b 2 k z 2 + ab 3 k z 2 + cd 3 k x 2 -2b 2 c 2 k x k z -2a 2 d 2 k x k z + b 2 cdk x k z + abd 2 k x k z +2abcdk x k z + abc 2 k x k z + a 2 cdk x   /deno B -1 K22 = 2 c 2 d 2 k x 2 + a 2 b 2 k z 2 -ab 3 k z 2 -cd 3 k x 2 + 2b 2 d 2 k x k z + 3b 2 c 2 k x k z +3a 2 d 2 k x k z + b 4 k z 2 + d 4 k x 2 -b 2 cdk x k z -abd 2 k x k z -4abcdk x k z /deno B -1 K23 = 2   a 3 bk z 2 -2c 2 d 2 k x 2 + c 3 dk x 2 -2a 2 b 2 k z 2 + 2ab 3 k z 2 + 2cd 3 k x 2 -2b 2 d 2 kxk z + c 2 b 2 k x k z + a 2 d 2 k x k z -b 4 k z 2 -d 4 k x 2 + 2b 2 cdk x k z +2abd 2 k x k z -6abcdk x k z + abc 2 k x k z + a 2 cdk x k z   /deno B -1 K31 = -2   -2a 3 bk z 2 + 2c 2 d 2 k x 2 -2c 3 dk x 2 + 2a 2 b 2 k z 2 -ab 3 k z 2 -cd 3 k x 2 -b 2 c 2 k x k z -a 2 d 2 k x k z + 2a 2 c 2 k x k z + a 4 k z 2 + c 4 k x 2 -b 2 cdk x k z -abd 2 k x k z + 6abcdk x k z -2abc 2 k x k z -2a 2 cdk x k z   /den B -1 K32 = 2   a 3 bk z 2 -2c 2 d 2 k x 2 + c 3 dk x 2 -2a 2 b 2 k z 2 + 2ab 3 k z 2 + 2cd 3 k x 2 -2b 2 d 2 kxk z + c 2 b 2 k x k z + a 2 d 2 k x k z -b 4 k z 2 -d 4 k x 2 + 2b 2 cdk x k z +2abd 2 k x k z -6abcdk x k z + abc 2 k x k z + a 2 cdk x k z   /deno B -1 K33 = 2   -3a 3 bk z 2 + 4c 2 d 2 k x 2 -3c 3 dk x 2 + 4a 2 b 2 k z 2 -3&b 3 k z 2 -3cd 3 k x 2 +2b 2 d 2 k x k z + 4b 2 c 2 k x k z + 4a 2 d 2 k x k z + 2a 2 c 2 k z k x + b 4 k z 2 + d 4 k x 2 +c 4 k x 2 -3b 2 cdk x k z -3abd 2 k x k z -3abc 2 k x k z -3a 2 cdk x k z   /deno

A.3 Streamline calculation with MHFEM

Providing the initial location of the paricles, in each element, the trajectory of the particle is calculated and the output edge and the coordinates of the exit point. Knowing the output edge of the element, the edge-element adjacency matrix is used to determine the new element, into which the particle enters.

In the two-dimensional case with triangular elements, for each element G, the fluxes through the edges are known, as well as the coordinates of the entry point of the particle. The velocity field at any point of the element is given by the relation (2.5). Since the functions w i are complex on the real element G, the problem is transferred to the element of reference Ĝ [START_REF] Kaasschieter | Mixed finite elements for accurate particle tracking in saturated groundwater flow[END_REF]. The basic functions ˆ w i are perfectly determined on the element Ĝ (see Section A.1):

ˆ w 1 = u v -1 ˆ w 2 = u -1 v ˆ w 3 = u v (A.1)
The coordinates of the entry point in the element of reference (u in , v in ) are determined using the function F -1 .F where F is the bijective transformation that allows to pass from the element of reference to the real element (see Section A.1). The components of the velocity in the element of reference are written as:

q Ĝ,u = du dt = Q G,T u -Q G,2 q Ĝ,v = dv dt = Q G,T v -Q G,1 (A.2) where Q G,T = 3 ∑ j=1 Q G, j
Knowing the initial position of the particle in the element of reference (coordinates of the input point of the particle in the element of reference):

u(0) = u in v(0) = v in
By solving the system of differential equations (A.2), the trajectory of the particle in the element can be obtained [START_REF] Pollock | Semianalytical computation of path lines for finite-difference models[END_REF]:

u(t) = e Q G,T t u in + (1-e Q K,T t )Q K,2 Q K,T v(t) = e Q G,T t v in + (1-e Q K,T t )Q K,1 Q K,T if Q G,T = 0 u(t) = u in -t Q G,2 v(t) = v in -t Q G,1 if Q G,T = 0 (A.3)
If Q T is negative, a drainage point is present and the particle converges to the point with coordinates of (

Q G,2 Q G,T , Q G,1 Q G,T
) in the element of reference (it is needed to to apply F to have these coordinates in the real coordinate space). Note that if the coordinates of this point are outside of the element of reference, the particle only pass through the element without converging. In a steady state, Q G,T is zero in the absence of a well-source term in an element. Thus, the velocity q in each element is then constant and the trajectory of the particle is rectilinear.

Then, the times that the trajectory of the particle meets the other two edges of the triangle, is determined. It has to be considered that the particle can not exit from the edge from which it enters (continuity of the normal component of the velocity vector q at the interface of two elements) and it exits through the edge with the smallest positive time. In addition to the point of exit, a series of points belonging to the trajectory in the element of reference is calculated. The function F is applied to all these points, in order to trace the trajectory of the particle in the real element.

APPENDIX C

Discretization of Darcy's law over an element

The hydraulic conductivity tensor K is assumed invertible. Then, integrating over the element G , using the base function u G , Equation (4.8) can be represented as

G u G K-1 q = - G u G ∇ ĥ - G u G K-1 K∇z (C.1)
An approximation can be obtained under the application of the product rule of divergence:

G u G K-1 q = - G ∇ ĥ u G -ĥ∇ u G - G K-1 K ∇ z u G -z∇ u G (C.2)
Applying Green's theorem to the right hand side of Equation (C.2) and considering K constant over the element G, it can be obtained as:

G K-1 q G u G = G ĥG ∇ u G - i⊂∂ G T h G,E i u G n G,E i + K-1 G K G G z G ∇ u G -K-1 G K G i⊂∂ G z G,E i u G n G,E i (C.3)
where ∂ G is the edge of the element G, which is composed of three edges E j (∀ j = 1, 2, 3) and n G,E i is the exterior normal unit vector to E j . T h G,E i and KG,E i are the approximations of the mean transformed pressure head (L) and the mean transformed hydraulic conductivity (LT -1 ) at the edge E i of the triangular element respectively. Z G,E i is the elevation head at the center of the edge E i .

Using the basis functions, Equation (C.3) can be rewritten as:

G K-1 G → q G → w i = G ĥG ∇ → w i - i⊂∂ G T ĥG,E i → w i → n G,E i + K-1 G K G G z G ∇ → w i -K-1 G K G i⊂∂ G z G,E i → w i → n G,E i (C.4)
So that, using expression (2.5) in Equation (C.4), we obtain

3 ∑ j=1 Q G,E j G K-1 G → w j → w i = ĥG G ∇ → w i -T ĥG,E i i⊂∂ G → w i → n G,E i + K-1 G K G z G G ∇ → w i -K-1 G K G z G,E i i⊂∂ G → w i → n G,E i (C.
Q G,E BG = ĥG DIV T G + K-1 G K G z G DIV T G -T ĥG,E -K-1 G K G z G,E (
C.8) by defining matrices as below: Hence equation (C.9) can be expressed as:

Q G,E =   Q G,E 1 Q G,E 2 Q G,E 3   , DIV T G =   1 1 1   , T ĥG,E 3 =   T ĥG,E 1 T ĥG,E 2 T ĥG,E 3   z G,E =   z G,E 1 z G,E 2 z G,E 3   , BG =   BG,
Q G,E i = αG,E i ĥG + 3 ∑ j=1 B-1 G,i, j K-1 G K G z G - 3 ∑ j=1 B-1 G,i, j T ĥG,E j - 3 ∑ j=1 B-1 G,i, j K-1 G K G z G,E j ∀i = 1,
Q G,E i = αG,E i ĥG + 3 ∑ j=1 B-1 G,i, j K-1 G K G z G -∑ E j ⊂∂ Ω D B-1 G,i, j T ĥG,E j -∑ E j ⊂∂ Ω D B-1 G,i, j K-1 G K G z G,E j -ÎG,E i ∀i = 1, 2, 3 (C.11)
Neumann boundary conditions: It can be represented by the equality

Q G,E i = Q N G,E i ∀E i ⊂ ∂ Ω N (C.12) where Q N G,E i = αG,E i ĥG + 3 ∑ j=1 B-1 G,i, j K-1 G K G z G -∑ E j ⊂∂ Ω D B-1 G,i, j T ĥG,E j -∑ E j ⊂∂ Ω D B-1 G,i, j K-1 G K G z G,E j -ÎG,E i ∀E i ⊂ ∂ Ω N (C.13)
Hydraulic gradient boundary condition: A unit hydraulic gradient boundary condition also known as free drainage boundary condition, which is equal to the hydraulic conductivity for any particular pressure head at given time, can be used to allow water to leave a flow domain under non-forced drainage conditions:

Q N G,E i = K G,E i (C.
ĥ = ĥG nm L = LG nm R = RE,E ′ n f ,n f T ĥ = T ĥE n f Ĵ = ĴE n f Î = ÎE n f V = VE n f DE,G = αG,E i f E ⊂ ∂ G 0 i f E ⊂ ∂ G ĥG = ĥG ∀G ⊂ Ω LG = K-1 G K G z G ∀G ⊂ Ω RE,E ′ = ∑ G⊃(E and E ′ ) B-1 G,E,E ′ ∀E ⊂ ∂ Ω D , ∀E ′ ⊂ ∂ Ω D T ĥE = T ĥG,E ∀E ⊂ ∂ Ω D ĴE = K-1 G K G z G,E ∀E ⊂ ∂ Ω D ÎE = ∑ G⊃E I G,E ∀E ⊂ ∂ Ω D VE = Q N G,E ∀E ⊂ ∂ Ω N 0 ∀E ⊂ ∂ Ω N

C.3 Mass conservation

Using the standard pressure based form of the Richards equation, the equation of mass conservation is given by 16) which is valid in each element G and by applying the chain rule it can be rewritten as (C.20) where F G = G f (x, y,t) = |G| f G , and |G| is the area of the element G.

G C ∂ h ∂t + G S w S s ∂ h ∂t + G ∇ → q = G f (x, z,t) ∀ G over Ω for t ⊂]0, T [ (C.
G C ∂ h ∂ ĥ ∂ ĥ ∂t + G S s θ φ ∂ h ∂ ĥ ∂ ĥ ∂t + G ∇ → q = G f (x,
Base on Raviart-Thomas properties, ĥG and

∇ → q G = 1 |G| ∑ E i ⊂∂ G Q G,E

C.3.1 Time discretization and linearization

The mass balance equation can be time discretizated using different time discretization approximations. Comparing the implicit finite difference schemes, Taheri Shahraiyni and Ataie Ashtiani [2008] concluded that among the different numerical schemes (Fully implicit, Crank-Nicolson and Runge-Kutta schemes), the fully implicit scheme is better than other schemes for numerical solution of pressure based Richards equation.

Solutions with poor mass balance and associated poor accuracy can be improved by the numerical approach proposed by [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF] using a method based on a fully implicit (backward Euler) time approximation applied to the mixed form of the Richards equation.

Pressure based formulation and a backward Euler time discretization has shown to produce mass balance errors in several cases because of the manner in which the time derivative ∂ θ ∂t is approximated as C ∂ h ∂t . Even if these terms are mathematically equivalent in the continuous partial differential equation, their discrete analogues are not equivalent. This inequality is aggravated by the highly nonlinear nature of C. [START_REF] Kirkland | Algorithms for solving richards' equation for variably saturated soils[END_REF] also recommended the use of the fully implicit formulation as the use of a Crank-Nicholson scheme on the closely related mixed form of Richards equation fails to reduce truncation error, producing potential instabilities.

Thus, a fully implicit numerical scheme is used to express the time discretization for the Equations (C.21) and (C.23), resulting the following equations for the pressure head based and the mixed form of the Richards equation, respectively

Ĉn+1 G + θ n+1 G φ G S s n+1 G ĥn+1 G -ĥn G = ∆t n |G| |G| f n+1 G - 3 ∑ i=1 Q n+1 G,E i (C.24) θ n+1 G -θ n G + θ n+1 G φ G S s n+1 G ĥn+1 G -ĥn G = ∆t n |G| |G| f n+1 G - 3 ∑ i=1 Q n+1 G,E i (C.25)
The linearization of the Richards equation is based on Picard iterative process which consists of constructing a sequence of functions, which will get closer to the desired solution and results symmetry of the final system of matrix equations.

The linearized pressure head based equation is expressed as

Ĉn+1,m G + θ n+1,m G φ G S s n+1,m G ĥn+1,m+1 G -ĥn G = ∆t n |G| |G| f n+1,m+1 G - 3 ∑ i=1 Q n+1,m+1
G,E i G,E i (C.27) where m is defined as the iteration index. Using Celia ′ s approach [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF] the mass balance problem can be eliminated by directly approximating the temporal term ∂ θ ∂t with its algebraic analog [START_REF] Clement | A physically based, two-dimensional, finite-difference algorithm for modeling variably saturated flow[END_REF] . θ n+1,m+1 G in a truncated Taylor serie is expanded with respect to h . Neglecting all terms higher than linear, Celia ′ s function can be representated as:

θ n+1,m+1 G = θ n+1,m G + ∂ θ G ∂ h G n+1,m h n+1,m+1

C.3.2 Switching technique

Using switching technique, the primary variable is changed at each iteration, thus a better convergence behavior is achieved compared to both the mixed-form and pressure-head based form of the Richards equation [START_REF] Diersch | On the primary variable switching technique for simulating unsaturatedsaturated flows[END_REF][START_REF] Hao | A mass-conservative switching method for simulating saturatedunsaturated flow[END_REF].

If θ n+1,m φ G > tol f (0 ≤ tol f ≤ 1) then pressure head is used as primary variable, if not a mixed-form of the Richards equation is implemented. For the problems dealing with unsaturated water flow using MHFE formulations, significant unphysical oscillations are exhibited due to the time-dependent terms in the mass matrix (M) in the non-diagonal coefficients of the matrix. To avoid this phenomenon in the numerical solution, a mass condensation scheme was used [START_REF] Belfort | Modélisation des écoulements en milieux poreux non saturés par la méthode des elements finis mixtes hybrids[END_REF]. The expression of flux at each edge is defined in terms of stationary and transitory flow regimes [START_REF] Younes | A new mass lumping scheme for the mixed hybrid finite element method[END_REF].

F = FG nm NG,E =    β G αn+1,m G,E αn+1,m G i f E ⊂ ∂ G 0 i f E ⊂ ∂ G PG,E = β G αn+1,m G i f G = G ′ 0 i f G = G ′ ĤG = ∑ E⊂(∂ G∩∂ Ω D ) β G αn+1,m G,E T h G,E + K-1 G K G z G,E αn+1,m G βG,E = β G i f G = G ′ 0 i f G = G ′ FG = |G| f G G ⊂ Ω Û = ÛG nm ÛG =        (1 -β G ) ĥn G i f θ n+1,m G φ G > tol f
Using pressure-head based form of the Richards equation (

T θ n+1,m G,E i φ G > tol f ): Q G,E i = Q G,E i + |G| 3 f G - |G| 3 TC G,E i + S sG T θ G,E i φ G ∂ T h G,E i ∂t ∀i = 1, 2, 3 (C.35)
Using the mixed form of the Richards equation ( 2, 3 (C.36) where Q G,E i is the flux corresponding to the stationary problem without the sink/source term over the element G with an area of |G| and T h G,E i , TC G,E i , T θ G,E i represents traces of pressure, specific water capacity and water content, respectively, over the edge E i .

T θ n+1,m G,E i φ G ≤ tol f ): Q G,E i = Q G,E i + |G| 3 f G - |G| 3 ∂ T θ G,E i ∂t + S sG T θ G,E i φ G ∂ T h G,E i ∂t ∀i = 1,
In the equation (C.21), considering a stationary problem ( ∂ ĥ ∂t = 0), without the sink / source term over the element G and by calculating the average pressure over the element as a function of traces of pressure, as below ĥG = ∑ E j ⊂∂ Ω D αG,E j T ĥG,E j + ∑ E j ⊂∂ Ω D αG,E j T ĥG,E j αG ∀i = 1, 2, 3

+ 3 ∑ j=1 αG,E j K-1 G K G z G,E j αG -K-1 G K G z G (C.37)
Then the term Q G,E i is approximated as below

QG,E i = αG,E i ∑ E j ⊂∂ Ω D αG,E j T ĥG,E j + ∑ E j ⊂∂ Ω D αG,E j T ĥG,E j αG + αG,E i 3 ∑ j=1 αG,E j K-1 G K G z G,E j αG -∑ E j ⊂∂ Ω D B-1 G,i, j T ĥG,E j -∑ E j ⊂∂ Ω D B-1 G,i, j K-1 G K G z G,E j -I G,E i ∀i = 1, 2, 3 (C.38)
The expressions considering the transient flow regime obtained as Hence The flux at each edge is calculated by:

S sG T θ G,E i φ G ∂ T h G,E i ∂t = S sG T θ G,E i φ G ∂ T h G,E i ∂ T ĥG,E i ∂ T ĥG,E i ∂t = S sG T θG,E i φ G ∂ T ĥG,E i ∂t (C.39) TC G,E i ∂ T h G,E i ∂t = TC G,E i ∂ T h G,E i ∂ T ĥG,E i ∂ T ĥG,E i ∂t = T ĈG,E i ∂ T ĥG,E i ∂t (C.
Q n+1,m+1 G,E i =     ∑ E j ⊂∂ Ω D αn+1,m G,E i αn+1,m G,E j αn+1,m G -∑ E j ⊂∂ Ω D B-1,n+1,m G,i, j     T ĥn+1,m+1 G,E j + |G| 3 f G - |G| 3∆t T Ĉn+1,m G,E i + S s n+1,m G T θ n+1,m G,E i φ G T ĥn+1,m+1 G,E i +Y 0 G,E i +Y 1 G,E i +Y 2 G,E i (C.43)
where:

Y 0 G,E i =     ∑ E j ⊂∂ Ω D αn+1,m G,E i αn+1,m G,E j αn+1,m G -∑ E j ⊂∂ Ω D B-1,n+1,m G,i, j     T ĥG,E j D Y 1 G,E i =      3 ∑ j=1 αn+1,m G,E i αn+1,m G,E j αn+1,m G - 3 ∑ j=1 B-1,n+1,m G,i, j      K-1,n+1,m G K n+1,m G z G,E j Y 2 G,E i =        |G| 3∆t T Ĉn+1,m G,E i + S s n+1,m G T θ n+1,m G,E i φ G T ĥn G,E i i f T θ n+1 G,E i φ G > tol f -|G| 3∆t T θ n+1,m G,E i -T θ n G,E i -S s n+1,m G T θ n+1,m G,E i φ G T ĥn G,E i -T Ĉn+1,m G,E i T ĥn+1,m G,E i i f T θ n+1 G,E i φ G ≤ tol f
System of equations using mass lumping scheme: The system of linear equations governing the hydrodynamics using the Lumped MHFE formulation is given as: αG,E α′

G,E αG -B-1 G,E,E ′ ∀E ′ ⊂ ∂ Ω D Ŷ1 = Ŷ1 E n f Ŷ1 E = ∑ G⊃E Y 0 G,E +Y 1 G,E + |G| 3 f G ∀E ⊂ ∂ Ω D Ŷ2 = Ŷ2 E n f Ŷ2 E = ∑ G⊃E Y 2 G,E ∀E ⊂ ∂ Ω D
In this system of equations, as the time-dependent terms appears only in the diagonal coefficients of the matrix M, it would be no more oscillation problems.

As [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF] proposed the modified Picard iteration technique, this system can be re-expressed in terms of the increment in iteration ∆T ĥn+1,m

G,E = T h n+1,m+1

G,E

X2 E = ∑ G⊃E    |G| 3∆t T Ĉn+1,m G,E T ĥn+1,m G,E -T ĥn G,E i f T θ n+1 G,E n G > tol f |G| 3∆t T θ n+1,m G,E -T θ n G,E i f T θ n+1 G,E n G ≤ tol f ∀E ⊂ ∂ Ω D X3 = X3E,E n f ,n f X3E,E ′ = ∑ G⊃(EandE ′ )
αG,E α′

G,E αG -B-1 G,E,E ′ ∀E ⊂ ∂ Ω D ∀E ′ ⊂ ∂ Ω D APPENDIX D
Soil parameters used in the modified Mualem-van Genuchten model APPENDIX F

Characteristics of precipitations

  2.1 (a): Domain of the case-study and top view of the aquifer (b): Drainaquifer exchanges and cross sectional view of the aquifer . . . . . . . . . . 9 2.2 (a): FEM unit function Φ G i (x,y) (b): MHFEM basis vectors ω i for i=1,2,3 . . 11 2.3 FEM and MHFEM unknowns at the drain. . . . . . . . . . . . . . . . . . . 13 2.4 Triangular meshes and MHFEM basis vectors w 1 , w 2 , w 3 on the triangle G . 16 2.5 (a): Computational domain. (b): Distribution of the hydraulic conductivity in the heterogeneous case (m/d). . . . . . . . . . . . . . . . . . . . . . . . 20 2.6 Piezometric head and streamlines (magenta lines) computed with a leakage coefficient of 6 m -1 d -1

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.10 Velocity field near to the drains with the leakage coefficient C λ = 6 m -1 d -1 . Homogeneous transmissivity. (a): FEM, without barrier case. (b): MHFEM, without barrier case. (c): FEM, with the barrier case (d): MHFEM, with the barrier case. . . . . . . . . . . . . . . . . . . . . . . . 26 2.11 Velocity field near to the drains with the leakage coefficient C λ = 6 m -1 d -1 . Heterogeneous transmissivity. (a): FEM, without barrier case. (b): MHFEM, without barrier case. (c): FEM, with the barrier case (d): MHFEM, with the barrier case. . . . . . . . . . . . . . . . . . . . . . . . 27 3.1 Differentiation process of GW_FEM. Black: Original user code. Blue: User pre-process. Green: AD generation. Magenta: User post-process. . . . . . . 33 3.2 Differentiation of the MHFEM code by Tapenade. . . . . . . . . . . . . . . 34 3.3 Piezometric head sensitivity isolines to a perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . Homogeneous transmissivity, case without barrier. Left: FEM. Right: MHFEM. . . . . . . . . . . . . . . . . . . . . . 37 3.4 Piezometric head sensitivity isolines to a perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . Homogeneous transmissivity case with the barrier. Left: FEM. Right: MHFEM. . . . . . . . . . . . . . . . . . . . . . 37 3.5 Piezometric head sensitivity isolines to a perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . Heterogeneous transmissivity case without barrier. Left: FEM. Right: MHFEM. . . . . . . . . . . . . . . . . . . . . . 38 3.6 Piezometric head sensitivity isolines to a perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . Heterogeneous transmissivity case with the barrier. Left: FEM. Right: MHFEM. . . . . . . . . . . . . . . . . . . . 38 vii 3.7 Velocity field and its sensitivity to perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . FEM, Homogeneous transmissivity, case without barrier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Velocity field and its sensitivity to perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . MHFEM, Homogeneous transmissivity, case without barrier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Velocity field and its sensitivity to perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . FEM, Homogeneous transmissivity, case with the barrier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Velocity field and its sensitivity to perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . MHFEM, Homogeneous transmissivity, case with the barrier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Velocity field and its sensitivity to perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . FEM, Heterogeneous transmissivity, case without barrier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Velocity field and its sensitivity to perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . MHFEM, Heterogeneous transmissivity, case without barrier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Velocity field and its sensitivity to perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . FEM, Heterogeneous transmissivity, case with the barrier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Velocity field and its sensitivity to perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . MHFEM, Heterogeneous transmissivity, case with the barrier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 Piezometric head and streamlines (blue lines) computed with a leakage coefficient of 9 m -1 d -1 and approximated streamlines (green lines) computed with (C λ , δ C λ )=(6 m -1 d -1 ,3 m -1 d -1

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Sources of pollutants in urban runoff (Reference: Urban Community of Lyon -Water Department) . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Potential pathways of surface (blue) and ground water (orange) contamination [Ritter et al., 2002] . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Schematic figures of different types of constructed wetland: (a) Surface flow CW (b) Horizontal sub-surface flow CW (c) Vertical flow CW (Modified from Tilley et al. [2008]) . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Procedure to select head (h top ) or flux (q top ) top boundary condition (the flux leaving an element is considered positive) . . . . . . . . . . . . . . . . 4.5 Procedure for time stepping adjustment . . . . . . . . . . . . . . . . . . . 4.6 Infiltration rate -Van Dam and Feddes [2000] . . . . . . . . . . . . . . . . 4.7 Evaporation rate -Van Dam and Feddes [2000] . . . . . . . . . . . . . . . 4.8 (a) Layered soil profile (b) uniform soil profile (from Pan and Wierenga [1995]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Comparison between the approximated water pressure head / water content and Pan and Wierenga [1995] test cases with Neumann boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Comparison between the approximated water pressure head / water content and Pan and Wierenga [1995] test cases with Dirichlet boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Water pressure head time-series with surface run off-HYDRUS and MH-FEM approximations at different elevations . . . . . . . . . . . . . . . . . 4.12 Water pressure head time-series with surface run off and surface ponding-HYDRUS and MHFEM approximations at the surface . . . . . . . . . . . 5.1 Evolution of numerical simulation process and model validation using calibration (adapted from Xiong et al. [2009]) . . . . . . . . . . . . . . . . . 5.2 Conceptual diagram of the filteration bed in a CW . . . . . . . . . . . . . . 5.3 Hysteresis in the moisture characteristic . . . . . . . . . . . . . . . . . . . 5.4 Geographic map of the studied area and the presented sections . . . . . . . 5.5 Ostwaldergraben site, the sedimentation pond (front) and the constructed wetland (behind) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Schematic diagram of a vertical flow filter (cross-section view) . . . . . . . 5.7 Characteristics of simulated rain events. . . . . . . . . . . . . . . . . . . . 5.8 Modeling and temporal identification chart. Top red dash-dotted tree: Modeling hypothesis. Green solid lined subchart: Automatic differentiation and related features. Blue dashed loop: Temporal estimation of hydrodynamic parameter variability. . . . . . . . . . . . . . . . . . . . . . . 5.9 Numerical diagram of the CW : rainfall (R), inflow (Qin), evapotranspiration rate (ET) and the outflow from drainage pipe (Qout) . . . . . . . . . . 5.10 Numerical diagram of the SCW : (a) General scheme (b) 1D domain (c) 2D MHFEM computational domain . . . . . . . . . . . . . . . . . . . . . 5.11 Sensitivity of water pressure to vGM parameters time-series. (a) Sensitivities during the 1st period of May (b) Calibration strategy and temporal estimation intervals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.12 Approximated water pressure head time-series at the surface during periods: (a) April (b) May 1 (c) May 2 . . . . . . . . . . . . . . . . . . . . . . List of Tables 2.1 Comparison between the computed flux values (m 2 d -1 ) for the two methods in the homogeneous cases. . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Comparison between the computed flux values (m 2 d -1 ) for the two methods in the heterogeneous cases. . . . . . . . . . . . . . . . . . . . . . . . 3.1 Information about the original and tangent linear (TL) codes. . . . . . . . . 3.2 Taylor tests for GW_FEM and GW_MHFEM. . . . . . . . . . . . . . . . . . . . . 3.3 Comparison between the computed flux values in the homogeneous cases (m 2 d -1 ) and the sensitivity of the fluxes for the two methods. . . . . . . . 3.4 Comparison between the computed flux values in the heterogeneous cases (m 2 d -1 ) and the sensitivity of the fluxes for the two methods. . . . . . . . 4.1 Cumulative actual infiltration under intensive rain at a dry soil . . . . . . . 4.2 Cumulative actual evaporation under high evaporation demand . . . . . . . 4.3 Initial and Boundary Conditions, Simulation Times, and Profile Types. . . 4.4 Time steps, CPU times and mass balance errors for 12 cases with or without variable transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Time steps, CPU times and mass balance errors for a multiple rainfallrunoff event with or without variable transformation . . . . . . . . . . . . . 5.1 Characteristics of the Ostwaldergraben catchment area . . . . . . . . . . . 5.2 Statistics of precipitation during the calibration periods . . . . . . . . . . . 5.3 Hydraulic properties reported in the literature. . . . . . . . . . . . . . . . . 5.4 Multi rainfall events. Optimized parameters by MHFEM (gradient method) and NSEs: very good (green), good (blue), satisfactory (orange), unsatisfactory (red). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Multi rainfall events. Optimized parameters by HYDRUS (stochastic method) and NSEs: very good (green), good (blue), satisfactory (orange), unsatisfactory (red). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Feeding periods. Optimized parameters (gradient method) and NSEs. . . . 5.7 Feeding periods. Optimized parameters (stochastic method) and NSEs. . . 5.8 Drainage periods. Optimized parameters (gradient method) and NSEs. . . 5.9 Drainage periods. Optimized parameters (stochastic method) and NSEs. . . D.1 Soil parameters used in the modified Mualem-van Genuchten model . . . . F.1 Characteristics of sampled rainfalls from April 2013 to September 2013 . .

κ

  Abbreviations and definitions

Figure 2

 2 Figure 2.1: (a): Domain of the case-study and top view of the aquifer (b):Drain-aquifer exchanges and cross sectional view of the aquifer

Figure 2

 2 Figure 2.3: FEM and MHFEM unknowns at the drain.

Figure 2

 2 Figure 2.4: Triangular meshes and MHFEM basis vectors w 1 , w 2 , w 3 on the triangle G

Figure 2 . 5 :

 25 Figure 2.5: (a): Computational domain. (b): Distribution of the hydraulic conductivity in the heterogeneous case (m/d).

  Figure 2.6: Piezometric head and streamlines (magenta lines) computed with a leakage coefficient of 6 m -1 d -1 . Homogeneous transmissivity case without barrier.

Figure 2

 2 Figure 2.7: Piezometric head and streamlines (magenta lines) computed with a leakage coefficient of 6 m -1 d -1 . Homogeneous transmissivity case with the barrier.

Figure 2

 2 Figure 2.8: Piezometric head and streamlines (magenta lines) computed with a leakage coefficient of 6 m -1 d -1 . Homogeneous transmissivity case without barrier.

Figure 2 . 9 :

 29 Figure2.9: Piezometric head and streamlines (magenta lines) computed with a leakage coefficient of 6 m -1 d -1 . Heterogeneous transmissivity case with the barrier.

Figure 2

 2 Figure 2.10: Velocity field near to the drains with the leakage coefficient C λ = 6 m -1 d -1 . Homogeneous transmissivity. (a): FEM, without barrier case. (b): MHFEM, without barrier case. (c): FEM, with the barrier case (d): MHFEM, with the barrier case.

Figure 2

 2 Figure 2.11: Velocity field near to the drains with the leakage coefficient C λ = 6 m -1 d -1 . Heterogeneous transmissivity. (a): FEM, without barrier case. (b): MHFEM, without barrier case. (c): FEM, with the barrier case (d): MHFEM, with the barrier case.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2 Sensitivity analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.3 Automatic differentiation . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.3.1 Minimal effort AD strategy for existing codes . . . . . . . . . . . 33 3.4 Numerical results and sensitivity computations . . . . . . . . . . . . . . . 35 3.4.1 Validation of the differentiated codes . . . . . . . . . . . . . . . . 36 3.4.2 Piezometric head sensitivity . . . . . . . . . . . . . . . . . . . . 36 3.4.3 Velocity field and sensitivity . . . . . . . . . . . . . . . . . . . . 39 3.4.4 Streamlines sensitivity . . . . . . . . . . . . . . . . . . . . . . . 40 3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

Figure 3 . 1 :

 31 Figure 3.1: Differentiation process of GW_FEM. Black: Original user code. Blue: User pre-process. Green: AD generation. Magenta: User post-process.

  Figure 3.2: Differentiation of the MHFEM code by Tapenade.

Figure 3 . 3 :Figure 3 Figure 3 . 5 :

 33335 Figure 3.3: Piezometric head sensitivity isolines to a perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . Homogeneous transmissivity, case without barrier. Left: FEM. Right: MHFEM.

Figure 3

 3 Figure 3.6: Piezometric head sensitivity isolines to a perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . Heterogeneous transmissivity case with the barrier. Left: FEM. Right: MHFEM.

Figure 3

 3 Figure 3.7: Velocity field and its sensitivity to perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . FEM, Homogeneous transmissivity, case without barrier.

Figure 3 Figure 3

 33 Figure 3.8: Velocity field and its sensitivity to perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . MHFEM, Homogeneous transmissivity, case without barrier.

Figure 3 Figure 3

 33 Figure 3.10: Velocity field and its sensitivity to perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . MHFEM, Homogeneous transmissivity, case with the barrier.

Figure 3

 3 Figure3.12: Velocity field and its sensitivity to perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . MHFEM, Heterogeneous transmissivity, case without barrier.

Figure 3 .Figure 3

 33 Figure 3.13: Velocity field and its sensitivity to perturbation of 50% of the leakage coefficient C λ = 6 m -1 d -1 . FEM, Heterogeneous transmissivity, case with the barrier.

Figure 3

 3 Figure3.15: Piezometric head and streamlines (blue lines) computed with a leakage coefficient of 9 m -1 d -1 and approximated streamlines (green lines) computed with(C λ , δ C λ )=(6 m -1 d -1 ,3 m -1 d -1).Homogeneous transmissivity case.

Figure 3

 3 Figure3.16: Piezometric head and streamlines (blue lines) computed with a leakage coefficient of 9 m -1 d -1 and approximated streamlines (green lines) computed with (C λ , δ C λ )=(6 m -1 d -1 ,3 m -1 d -1 ). Homogeneous transmissivity, case with the barrier.

  Figure3.17: Piezometric head and streamlines (blue lines) computed with a leakage coefficient of 9 m -1 d -1 and approximated streamlines (green lines) computed with(C λ , δ C λ )=(6 m -1 d -1 ,3 m -1 d -1). Heterogeneous transmissivity, case without barrier.

Figure 3

 3 Figure3.18: Piezometric head and streamlines (blue lines) computed with a leakage coefficient of 9 m -1 d -1 and approximated streamlines (green lines) computed with(C λ , δ C λ )=(6 m -1 d -1 ,3 m -1 d -1). Heterogeneous transmissivity, case with the barrier.

Figure 4

 4 Figure 4.1: Sources of pollutants in urban runoff (Reference: Urban Community of Lyon -Water Department)

  Figure 4.3: Schematic figures of different types of constructed wetland: (a) Surface flow CW (b) Horizontal sub-surface flow CW (c) Vertical flow CW (Modified from Tilley et al. [2008])

Figure 4

 4 Figure 4.4: Procedure to select head (h top ) or flux (q top ) top boundary condition (the flux leaving an element is considered positive)

  heuristic = ∆t n ∆t heuristic = ∆t n *Factor smaller than 1 ∆t heuristic = ∆t n *

Figure

  Figure 4.5: Procedure for time stepping adjustment

Figure 4 .

 4 Figure 4.7 shows the simulated actual evaporation rate of sand. Initially the potential soil water flux is large enough to meet the potential soil evaporation rate. In cases simulated without variable transformation, at t=1.1 d the upper boundary condition changes from

Figure 4

 4 Figure 4.8: (a) Layered soil profile (b) uniform soil profile (from Pan and Wierenga [1995])

Figure 4 Figure 4

 44 Figure4.9: Comparison between the approximated water pressure head / water content and[START_REF] Pan | A transformed pressure head-based approach to solve richards' equation for variably saturated soils[END_REF] test cases with Neumann boundary conditions.

Figure 4

 4 Figure 4.11: Water pressure head time-series with surface run off-HYDRUS and MHFEM approximations at different elevations

Figure 4 .

 4 Figure 4.12: Water pressure head time-series with surface run off and surface ponding-HYDRUS and MHFEM approximations at the surface

Figure 5 . 1 :

 51 Figure 5.1: Evolution of numerical simulation process and model validation using calibration (adapted from Xiong et al. [2009])

Figure 5 .Figure 5 . 3 :

 553 Figure 5.2: Conceptual diagram of the filteration bed in a CW

Figure 5

 5 Figure 5.4: Geographic map of the studied area and the presented sections

Figure 5 . 5 :

 55 Figure 5.5: Ostwaldergraben site, the sedimentation pond (front) and the constructed wetland (behind)

Figure 5 . 6 :

 56 Figure 5.6: Schematic diagram of a vertical flow filter (cross-section view)

Figure 5

 5 Figure 5.7: Characteristics of simulated rain events.

Figure 5 . 8 :

 58 Figure 5.8: Modeling and temporal identification chart. Top red dash-dotted tree: Modeling hypothesis. Green solid lined subchart: Automatic differentiation and related features. Blue dashed loop: Temporal estimation of hydrodynamic parameter variability.

Figure 5 .

 5 Figure 5.9: Numerical diagram of the CW : rainfall (R), inflow (Qin), evapotranspiration rate (ET) and the outflow from drainage pipe (Qout)

Figure 5 .

 5 Figure 5.10: Numerical diagram of the SCW : (a) General scheme (b) 1D domain (c) 2D MHFEM computational domain

Figure 5 .

 5 Figure 5.11: Sensitivity of water pressure to vGM parameters time-series. (a) Sensitivities during the 1st period of May (b) Calibration strategy and temporal estimation intervals.

Figure 5 .

 5 Figure 5.12: Approximated water pressure head time-series at the surface during periods: (a) April (b) May 1 (c) May 2

Figure 5 .

 5 Figure 5.13: Approximated water pressure head time-series at the surface during periods: (a) July (b) August (c) September

Figure 5 .Figure 5 .

 55 Figure 5.14: Approximated water pressure head time-series during feeding, (a): Sub-period 1 (b): Sub-period 3 (c): Sub-period 6 (d): Sub-period 7 (e): Sub-period 8 (f): Sub-period 9 (g): Sub-period 11

Figure 5 .Figure 5 .

 55 Figure 5.16: Approximated water pressure head time-series during drainage, (a): Sub-period 2 (b): Sub-period 4 (c): Sub-period 5 (d): Sub-period 9 (e): Sub-period 10 (f): Sub-period 12

Figure 5 .

 5 Figure 5.18: Temporal variability of vGM parameters during all sub-periods

  Variations of 1 s1 along the sequence of drainage events

Figure 5 .FeedFigure 5 .

 55 Figure 5.19: (I)Temporal variability of vGM parameters during the drainage periods and (II) Hysteresis effect for some feeding-drainage events: (a) α 1 ; (b) n 1 ; (c) θ r1 ; (d) θ s1 ; (e) K 1 ; (f) K 3 .

→

  Figure B.1: Homogeneous-hydraulic conductivity equal to 333 m d -1 without barrier and drainage with a leakage coefficient of 6 m -1 d -1 . a: Piezometric head and streamlines (magenta lines) . b: Velocity field (blue vectors) near to the drains (dotted red lines).

  be represented in a matrix form as:

A

  Dirichlet condition Ω D and a Neumann condition Ω N are imposed on the boundaries ∂ Ω D and ∂ Ω N , respectively.Dirichlet boundary conditions:As it is explained in chapter 2, the boundary conditions of Dirichlet are represented as:ÎG,E i = ∑ E j ⊂∂ Ω D B-1 G,i, j T ĥG,E j + K-1 G K G z G,E j ∀i = 1,2, 3 (C.10)

  14) C.2 Matrix form of the continuity flux Dĥ + D L -RT ĥ -R Ĵ -Î -V = 0 (C.15) where: D = DE,G n f ,nm

  sum over the elements G containing the set of internal edges E and E ′ and the sum over the elements G containing the edge E.

  i are constant over the element G. Thus the balance equation when using the standard pressure based form of the Richards equation is expressed as:E i = |G| f G ∀ G over Ω Standard pressure based form of richards equation for t ⊂]0, T [ (C.21)The standard pressure based form of the Richards equation can lead to large mass balance errors; To improve the properties with respect to accurate mass conservative solutions, the mixed form (Equation (C.22)) is provided as below: have convergence difficulties for dry initial conditions. By using approximation for each element:|G| ∂ θ G ∂t + |G| θG φ E i = |G| f G ∀ G over ΩMixed form of richards equation for t ⊂]0, T [(C.23)

C. 4 N

 4 Approximation and the matrix form of the average pressure By multiplying the linearized equation by the term Ĉn+1,m of the average pressure is expressed by ĥn+1,m+1 -NT ĥn+1,m+1 -N Ĵn+1,m -P Fn+1,m+1 + β Ln+1,m -Ĥ -= NG,E nm,n f P = PG,G ′ nm,nm Ĥ = ĤG nm β = βG,G ′ nm,nm

  of equations using mixed hybrid finite element method Combining the two equations of state (C.15) and (C.33) and considering T h as a principal unknown yields the system of hybrid mixed equations. Rn+1,m -Dn+1,m N T ĥn+1,m+1 = Dn+1,m Û -Dn+1,m β Ln+1,m + Dn+1,m N Ĵn+1,m + Dn+1,m P Fn+1,m+1 + Dn+1,m Ĥ + Dn+1,m Ln+1,m -Rn+1,m Ĵn+1,m -În+1,m -V (C.34) Considering ( Rn+1,m -Dn+1,m N) as matrix M, which is positive definite symmetric,

  whereX = V -Ŷ1 + X1 + X2 -X3 T h n+1,m and X1 = X1E n f

  Sensitivity of the model M with respect to P 0 in the direction of pertur-

	xiii	Abbreviations and definitions
	e act E max e sat h d h atm h db h Obs h pond h Sim he K A K r K s m M(h) m v MB n n n e n f n f r nm p Q in Q load q top r bation δ P Actual vapor pressure Maximum soil evaporation (LT -1 ) Saturated vapor pressure Water level in water-body (L) Soil water pressure head in equilibrium with the air humidity (L) Level of water-body bottom (L) Observed piezometric head value Height of water ponding on the soil surface Simulated piezometric head value Free parameter in the modified Mualem-van Genuchten model referred as the air entry value (L) Dimensionless anisotropy tensor Relative hydraulic conductivity function Saturated hydraulic conductivity (LT -1 ) Picard iteration level Aquifer thickness (L) which varies with the water table elevation h Parameter (-) in the modified Mualem-van Genuchten model Mass balance ratio for each time step Time level iteration Effective porosity Number of the edges of the mesh Number of edges over the domain Ω where pressure head is not imposed Number of elements over the domain Ω Inflow into the soil column Quantity of water loaded during a time lapse Potential flux at the soil surface S s Specific storage (L -1 ) S w Degree of saturation (-) Sc Saturation at the cut-off point he Se Effective saturation tol a Tolerance for discrepancy between the approximated pressure head of the edge for two successive iteration levels tol c Tolerance for the calculated values of water content between two succes-sive iteration levels tol f Tolerance for the switching procedure tol r Relative tolerance for the iteration convergence test v Mean velocity of a conservative tracer V g Volume of the gaseous phase V l Volume of the liquid phase V s Volume of the solid phase V δ i, j Kronecker symbol AD Automatic Differentiation ET Evapotranspiration Ex Exflitration FEM Finite Element Method MHFEM Mixed Hybrid Finite Element Method MSE Mean Squared Error NSE Nash-Sutcliffe Efficiency PCA Principal Component Analysis TL Tangent Linear VFSCW Vertical Flow Stormwater Constructed Wetland vGM van Genuchten-Mualem T h G,E element G CHAPTER 1

T -1 ) per meter

E j

Edge of the element G for j = 1, 2, 3 * A set of parameters minimizing the objective function

Q e Exchange flow q bot Flux at soil profile bottom (LT -1 ) Q G,E j Water flux (L 2 T -1

) over the edge E j belonging to the element G Q in f Quantity of actual flux infiltrating at the surface ω Ratio that compares the sensitivity obtained with Automatic Differentiation to first order finite difference approximations s M (P 0 , δ P) i Mean piezometric head (L) of each edge ∀i = 1, 2, 3 which belongs to the air Total air volume in the soil profile at the start of the time step

Z G,E i

Elevation head (L) at the center of the edge E i of the element G

  L'objectif est d'évaluer l'influence du coefficient d'échange dans la modélisation d'une mesure d'atténuation du drainage. Les modèles mathématiques décrivant l'écoulement de l'eau dans un milieu poreux indéformable

sont basés sur la loi de Darcy et l'équation de continuité de masse. La modélisation des échanges drain-aquifère est étudiée en adaptant le modèle classique des échanges nappe-rivière. Plus précisément, le flux d'échange est modélisé à l'aide du concept d'échange présenté par

[START_REF] Kinzelbach | Groundwater modelling. an introduction with sample programs in basic. developments in water sciences[END_REF]

et basé sur une relation linéaire entre le niveau d'eau dans la rivière et le niveau d'eau dans l'aquifère.

  boundary conditions: Neumann boundary condition is used where a prescribed flow across the bounding edges is known. It can be represented by the equality

Table 2 .

 2 1: Comparison between the computed flux values (m 2 d -1 ) for the two methods in the homogeneous cases.

			Without barrier			With the barrier
	Methods Channel River	Drains MBE Channel River	Drains MBE
	FEM	144.9	-113.2 -31.7	0.0	269.5	-173.5 -96.0	0.0
	MHFEM 138.7	-118.5 -20.2	0.0	244.7	-173.1 -71.6	0.0

Table 2 .

 2 2: Comparison between the computed flux values (m 2 d -1 ) for the two methods in the heterogeneous cases.

			Without barrier			With the barrier
	Methods Channel River	Drains MBE Channel River	Drains MBE
	FEM	315.6	-293.4 -22.2	0.0	300.5	-280.0 -20.5	0.0
	MHFEM 306.5	-286.7 -19.8	0.0	289.0	-271.7 -17.3	0.0

Table 3 .

 3 1: Information about the original and tangent linear (TL) codes. This differentiation method has been successfully applied to both the FEM and the MHFEM codes.

		Original code		TL code	Runtime Ratio
	Code	Routines Lines Runtime Routines Lines Runtime TL/Original
	FEM	36	2489 0.51 s	24	3245 0.61 s	1.19
	MHFEM 26	1491 7.44 s	18	1677 14.6 s	1.96
	same dimension as P).

Table 3 .

 3 2: Taylor tests for GW_FEM and GW_MHFEM.

	ω	r ω (GW_FEM) r ω (GW_MHFEM)
	10 1	0.186	0.216
	10 0	0.186	0.741
	10 -1 0.959	0.967
	10 -2 0.996	0.997
	10 -3 0.999	1.000
	10 -4 0.962	1.004
	10 -5 1.040	1.095
	10 -6 1.183	0.816
	10 -7 15.740	3.504

Table 3 .

 3 3: Comparison between the computed flux values in the homogeneous cases (m 2 d -1 ) and the sensitivity of the fluxes for the two methods. head sensitivity approximated by the two methods is less remarkable. A glance to see Figure3.4 suffices to observe that the almost impervious barrier is efficient since a perturbation of the leakage coefficient has no impact outside the protected area. Likewise, the computed flux sensitivity at the river is almost 0 in Table3.3.

			Without barrier			With the barrier
	Methods	Channel River	Drains MBE Channel River	Drains MBE
	FEM	144.9	-113.2 -31.7	0.0	269.5	-173.5 -96.0	0.0
	FEM_D	2.9	3.3	-6.2	0.0	4.0	0.0	-4.0	0.0
	FEM_D/FEM	2.0%	-2.9% 19.6% -	1.4%	0.0%	4.1%	-
	MHFEM	138.7	-118.5 -20.2	0.0	244.7	-173.1 -71.6	0.0
	MHFEM_D	2.3	2.7	-5.0	0.0	10.7	0.0	-10.7	0.0
	MHFEM_D/MHFEM 1.7%	-2.3% 25.2% -	4.4%	0.0%	14.9% -
	Same calculations were performed in the heterogeneous transmissivity cases (see Fig-
	ures 3.5

piezometric

Table 3 .

 3 4: Comparison between the computed flux values in the heterogeneous cases (m 2 d -1 ) and the sensitivity of the fluxes for the two methods.

			Without barrier			With the barrier
	Methods	Channel River	Drains MBE Channel River	Drains MBE
	FEM	315.6	-293.4 -22.2	0.0	300.5	-280.0 -20.5	0.0
	FEM_D	0.3	0.3	-0.6	0.0	0.908	0.0	-0.9	0.0
	FEM_D/FEM	0.1%	0.1%	2.7%	-	0.03%	0.0%	4.4%	-
	MHFEM	306.5	-286.7 -19.8	0.0	289.0	-271.7 -17.3	0.0
	MHFEM_D	0.6	0.6	-1.2	0.0	1.5	0.0	-1.5	0.0
	MHFEM_D/MHFEM 0.2%	-0.2% 6.1%	-	0.1%	0.0%	9.0%	-

  al.,

	Ground water	Urban and
	infiltration	rural runoff
	Direct discharge to water bodies	Potential pathways of contamination	Naturally occurring substances
	Infiltration of pesticide, fertilizers	Injection of contaminants

Figure 4.2: Potential pathways of surface (blue) and ground water (orange) contamination

Table 4 .

 4 3: Initial and Boundary Conditions, Simulation Times, and Profile Types.

	Case Initial Pressure (cm) Upper Boundary	Lower Boundary Simulation Time (s) Profile Type
	1.1

Table 4 .

 4 4: Time steps, CPU times and mass balance errors for 12 cases with or without variable transformation Neumann boundary condition by imposing low and high flux, respectively. The test cases 1.2 and 2.2 consider low and high rate infiltration, respectively, in a semi-arid. The test cases 1.3 and 2.3 simulate low and high rate infiltration, respectively, in a wet soil.

	Model			h-based			Mixed			Switch
	Case	κ (cm -1 ) Time steps CPU (s) εMB	Time steps CPU (s) εMB	Time steps CPU (s) εMB
	1.1	0.0	5465	26.01	1.2 × 10 -12	4910	22.93	2.7 × 10 -9	5234	24.10	-5.5 × 10 -10
	1.1	-0.04	4580	20.57	-1.4 × 10 -10	4448	20.16	1.5 × 10 -10	4427	20.19	-1.9 × 10 -10
	1.2	0.0	2666	9.41	1.7 × 10 -14	2608	9.23	-3.2 × 10 -8	2571	9.07	-7.2 × 10 -9
	1.2	-0.04	2525	9.00	3.0 × 10 -14	2522	9.01	7.6 × 10 -9	2520	8.97	-3.8 × 10 -10
	1.3	0.0	1785	6.47	2.0 × 10 -14	1770	6.44	-1.8 × 10 -9	1772	6.45	-1.7 × 10 -9
	1.3	-0.04	1779	6.52	2.0 × 10 -15	1746	6.34	1.1 × 10 -8	1735	6.47	-7.6 × 10 -10
	2.1	0.00	4647	21.17	1.9 × 10 -12	4405	19.11	-4.6 × 10 -10 4433	19.41	-4.5 × 10 -10
	2.1	-0.04	4266	18.31	-1.1 × 10 -9	4244	18.38	-1.0 × 10 -9	4244	18.38	-1.0 × 10 -9
	2.2	0.0	2142	7.22	1.0 × 10 -14	2077	6.84	-3.8 × 10 -8	2077	6.90	-1.4 × 10 -8
	2.2	-0.04	2018	6.60	-7.9 × 10 -10	2017	6.63	3.3 × 10 -9	2018	6.66	-7.9 × 10 -10
	2.3	0.00	972	3.32	-2.8 × 10 -14	933	3.20	-2.0 × 10 -8	933	3.20	-1.6 × 10 -8
	2.3	-0.04	921	3.02	-1.39 × 10 -13 922	3.00	-6.2 × 10 -10 922	3.01	-2.6 × 10 -9
	3.1	0.00	7288	38.01	-1.3 × 10 -13	7590	40.37	1.1 × 10 -9	8064	44.86	-3.7 × 10 -10
	3.1	-0.04	6588	11.14	-2.76 × 10 -11 6283	36.15	5.5 × 10 -12	6040	35.93	-6.5 × 10 -13
	3.2	0.00	3401	18.28	-2.1 × 10 -14	3153	16.27	3.7 × 10 -7	3225	16.42	-1.9 × 10 -8
	3.2	-0.04	2158	11.06	2.9 × 10 -13	2149	11.06	-1.9 × 10 -8	2219	11.29	-9.7 × 10 -10
	3.3	0.0	2667	14.12	4.4 × 10 -12	2745	14.51	-1.0 × 10 -7	2680	14.13	-2.1 × 10 -13
	3.3	-0.04	1843	9.71	8.0 × 10 -14	1748	9.48	9.3 × 10 -8	1911	9.99	-1.4 × 10 -10
	4.1	0.00	1963	5.00	1 × 10 -15	1959	5.07	-1.1 × 10 -5	1963	5.10	-1.3 × 10 -5
	4.1	-0.04	2031	5.01	1.5 × 10 -12	2031	5.16	-2.8 × 10 -5	2031	5.11	-2.9 × 10 -5
	4.2	0.00	2052	6.50	-1.2 × 10 -14	1974	6.40	-8.5 × 10 -8	1979	6.37	-7.4 × 10 -8
	4.2	-0.04	2005	6.00	1.5 × 10 -13	2004	6.17	-1.07 × 10 -8 2004	5.92	-1.1 × 10 -8
	4.3	0.00	1857	4.85	2.2 × 10 -16	1852	4.86	-8.5 × 10 -8	1866	4.92	-1.4 × 10 -8
	4.3	-0.04	1857	4.75	9.3 × 10 -14	1843	4.73	-6.6 × 10 -12 1857	4.79	-4.8 × 10 -8

moist (h = -200 cm), very dry (h = -50000 cm) and intermediate (h = -1000 cm) initial conditions are presented in Figures 4.9 and 4.10.

The test cases 1.1 and 2.1 show the infiltration in arid soil with a

Table 4 .

 4

	Model			h-based			Mixed			Switch
	Case	κ (cm -1 ) Time steps CPU (s) εMB	Time steps CPU (s) εMB	Time steps CPU (s) εMB
	Surface runoff	0.0	41121	115.65	3.408 × 10 -13 68614	205.09	1.7 × 10 -11	43422	121.58	-3.0 × 10 -6
		-0.04	40941	112.73	4.2 × 10 -13	68205	201.92	2.1 × 10 -10	43806	124.95	8.1 × 10 -6
	Surface ponding 0.0	21443	46.17	-3.2 × 10 -12 21666	48.72	-2.2 × 10 -8 21411	45.52	1.47 × 10 -5
		-0.04	21491	45.68	-1.1 × 10 -12 21534	46.34	-3.6 × 10 -9 21439	45.72	-1.1 × 10 -5

5: Time steps, CPU times and mass balance errors for a multiple rainfall-runoff event with or without variable transformation

Table 5 .

 5 1: Characteristics of the Ostwaldergraben catchment area

		Total surface	18000
		Active surface	5200 m 2
		Percentage of runoff	0,29
	Catchment area	Total roof area (% of total area)	2 470 m 2 (13,8%)
		Area of the total road (% of total area) 2 780 m 2 (15,6%)
		Type	Residential
		Accommodation (% of surface area)	Suburban houses (100%)

Table 5 .

 5 2: Statistics of precipitation during the calibration periods

		Average dry	Number of	Percetage of wet	Maximum	Minimum	Maximum
	Month	period	rainfall	periods to the	Rainfall	rainfall	intensity
		duration (day)	events	total rainy days	duration (min)	duration (min)	(cm/min)
	April	1.41	4	30.82%	2119	3	0.645
	May	0.60	9	33.81%	761	17	0.312
	July	1.97	2	3.16%	270	31	0.216
	August	1.61	2	6.26%	360	120	0.311
	September	0.80	11	25.94%	276	34	0.654

Table 5 .

 5 3: Hydraulic properties reported in the literature. determined by experimental methods. Particularly, the order of magnitude θ r and θ s are given in the literature for different types of environment (clay, silt, sand, gravel,...). The hydraulic parameters of van Genuchten law

	Layer i	References	α (cm -1 )	n (-)	θ r (-)	θ s (-)	K s (cm.min -1 )
	1	Celia et al. [1990]	0.0335 2.0 0.102 0.368	0.5532
	2	DanHan-Cheng et al. [2012]	1.0	2.19 0.005 0.42	2.0
	3	Fournel et al. [2013]	0.036 1.56 0.078 0.43	0.00823
	and can be

Table 5 .

 5 5: Multi rainfall events. Optimized parameters by HYDRUS (stochastic method) and NSEs: very good (green), good (blue), satisfactory (orange), unsatisfactory (red).

	min -1 ) (cm 2 )	(-)

Table 5 .

 5 6: Feeding periods. Optimized parameters (gradient method) and NSEs.

	Period	Start α 1 date (cm -1 ) (-) n 1	θ r1 (-)	θ s1 K 1 (-) (cm.min -1 ) (cm 2 ) MSE NSE K 3 (-)
	01	16/04 .010	2.31 .013 .05 .009 .0008	.48	.97
	03	25/04 .924	2.77 .007 .15 .028 .0034 33.86	.84
	06	02/05 .151	2.63 .110 .48 .346 .0087 11.55	-.35
	07	03/05 .877	4.94 .013 .81 .420 .0116 15.49	.35
	08	04/05 .889	4.95 .028 .71 .422 .0071	8.72	-.03
	09	26/05 .137	2.02 .005 .35 .028 .0038 27.72	.58
	11	29/05 .026	1.73 .013 .90 .026 .0100 32.03	.34
	13	02/07 .387	3.04 .092 .44 .027 .0005	.77	.98
	16	24/08 .062	8.14 .012 .24 .019 .0011	1.50	.84
	17	26/08 .045	8.49 .016 .23 .005 .0008	.03	.99
	20	08/09 .035	3.81 .029 .29 .010 .0021 12.78	.72

Table 5 .

 5 7: Feeding periods. Optimized parameters (stochastic method) and NSEs.

	Period	Start α 1 date (cm -1 ) (-) n 1	θ r1 (-)	θ s1 K 1 (-) (cm.min -1 ) K 3	MSE NSE (cm 2 ) (-)
	01	16/04 0.044	3.94 .003 .07 0.016 .0220	7.92	.55
	03	25/04 1.118	3.14 .064 .50 0.570 .0166 137.53	.33
	06	02/05 0.206	2.25 .007 .47 0.481 .0192	5.14	.40
	07	03/05 0.582	4.82 .09	.90 0.522 .0087	17.27	.28
	08	4/05 0.119	5.89 .199 .89 1.401 .0129	4.52	.47
	09	26/05 0.159	7.30 .034 .34 0.026 .0042	23.63	.64
	11	29/05 0.010	5.82 .076 .18 0.014 .0054	22.01	.55
	13	02/07 0.305	2.78 .096 .41 0.022 .0004	10.72	.67
	16	26/08 0.105	6.45 .009 .40 0.028 .0007	5.34	.42
	17	17/08 0.112	5.84 .006 .72 1.347 .0008	2.25	.68
	20	08/09 0.042	4.49 .029 .20 0.010 .0021	13.79	.69

Table 5 .

 5 8: Drainage periods. Optimized parameters (gradient method) and NSEs.

	Period	Start α 1 date (cm -1 ) (-) n 1	θ r1 (-)	θ s1 K 1 (-) (cm.min -1 ) (cm 2 ) MSE NSE K 3 (-)
	17/04 .005	1.68 .250 .30 0.557 .0044 3.642 .836
	28/04 .500	2.80 .190 .48 0.148 .0025 0.119 998
	01/05 .011	2.00 .080 .39 0.550 .0027 0.198 .934
	05/05 .412	2.89 .024 .64 0.347 .0038 0.293 .990
	27/05 .590	3.85 .068 .65 0.250 .0028 0.497 .981
	30/05 .067	1.48 .010 .74 0.832 .0008 0.063 .993
	03/07 .498	5.10 .010 .95 1.106 .0024 0.458 .996
	21/08 .034	2.00 .102 .37 0.553 .0014 0.028 .989
	27/08 .053	2.23 .089 .41 0.550 .0011 0.010 .999
	01/09 .003	2.03 .020 .23 0.550 .0014 0.280 .945
	Mean	.217	2.61 .102 .52 0.544 .0023
	Dev.		.248	1.12 .084 .22 0.275 .0012
	Min		.003	1.48 .010 .23 0.148 .0008
	Max		.590	5.10 .250 .95 1.106 .0044

Table 5 .

 5 9: Drainage periods. Optimized parameters (stochastic method) and NSEs.

	Period	Start α 1 date (cm -1 ) (-) n 1	θ r1 (-)	θ s1 K 1 (-) (cm.min -1 ) (cm 2 ) MSE NSE K 3 (-)
	17/04 .005	1.01 .453 .76 0.648 .0038 3.897 .825
	28/04 .416	2.40 .160 .58 0.223 .0035 0.308 .996
	01/05 .052	1.61 .144 .27 1.080 .0033 0.235 .921
	05/05 .446	2.59 .182 .69 0.893 .0034 0.253 .992
	27/05 .110	5.21 .168 .93 0.814 .0029 1.177 .954
	30/05 .061	2.82 .016 .70 0.678 .0012 0.203 .978
	03/07 .402	4.09 .079 .92 0.206 .0021 1.054 .990
	21/08 .049	5.17 .080 .68 0.139 .0014 0.025 .990
	27/08 .060	3.07 .082 .39 0.550 .0011 0.016 .998
	01/09 .004	2.28 .173 .40 0.970 .0014 0.288 .943
	Mean	.161	3.03 .154 .63 0.620 .0024
	Dev.		.183	1.40 .118 .22 0.336 .0010
	Min		.004	1.01 .016 .27 0.139 .0011
	Max		.446	5.21 .453 .93 1.080 .0038

Table D

 D .1: Soil parameters used in the modified Mualem-van Genuchten model Material α (1/cm) n (-) K (cm/sec) θ r (-)

						θ s (-)	τ	h e (cm)	Reference
	Soil 1	0.0280	2.2390 0.006262 0.0286 0.3658	0.5	0	Pan and Wierenga [1995]
	Soil 2	0.0104	1.3954 0.0001516 0.1060 0.4686	0.5	0	Pan and Wierenga [1995]
	Soil 3	0.1020	2.0000	0.00922	0.1020 0.3680	0.5	0	Pan and Wierenga [1995]
	Sand	0.0249	1.507 0.0002025	0.01	0.43 -0.140	0	Van Dam and Feddes [2000]
	Loam	0.036	1.89	0.0002888 0.078	0.43	0.5	0	Simunek et al. [1998]

Table F

 F .1: Characteristics of sampled rainfalls from April 2013 to September 2013

	Start date event	Start time	Duration (min)	Depth (mm)	Average intensity (mm/h)	Dry period duration (d)	Return period
	16/04/2013	02:32:00 AM	7	0.4	3.4	3.46	<1 week
	26/04/2013	03:10:00 PM	24	0.4	1.0	6.14	<1 week
	26/04/2013	04:55:00 PM	279	15.4	3.3	0.06	3 to 6 months
	26/04/2013	11:25:00 PM	19	0.4	1.3	0.08	<1 week
	27/04/2013	12:53:00 AM	928	23.0	1.5	0.05	3 to 6 months
	27/04/2013	05:57:00 PM	217	5.6	1.5	0.07	1 to 2 weeks
	27/04/2013	10:09:00 PM	62	1.2	1.2	0.02	<1 week
	02/05/2013	06:56:00 PM	49	4.0	4.9	2.01	2 weeks to 1 month
	03/05/2013	05:26:00 PM	108	4.2	2.3	0.90	1 to 2 weeks
	04/05/2013	08:51:00 AM	175	2.6	0.9	0.57	<1 week
	04/05/2013	06:08:00 PM	71	1.8	9.8	0.27	2 weeks
	26/05/2013	10:57:00 PM	629	13.4	1.3	0.34	1 to 1,5 month
	29/05/2013	03:12:00 AM	928	23.0	1.5	0.05	3 to 6 months
	29/05/2013	10:10:00 PM	300	4.1	1.0	0.05	1 to 2 weeks
	02/07/2013	12:43:00 PM	274	8.9	1.9	3.69	1 month
	02/07/2013	11:04:00 PM	28	0.8	1.8	0.25	<1 week
	24/08/2013	12:29:00 PM	66	9.8	4.7	4.88	1,5 to 3 months
	25/08/2013	04:52:00 AM	16	0.4	1.5	0.59	<1 week
	26/08/2013	03:15:00 AM	260	10.0	2.3	0.92	1 to 1,5 months
	26/08/2013	09:13:00 AM	66	0.8	0.7	0.07	<1 week
	07/09/2013	06:06:00 AM	59	1.2	1.2	11.82	<1 week
	08/09/2013	01:19:00 AM	150	31.8	12.7	0.76	5 to 10 years
	08/09/2013	06:38:00 AM	30	0.4	0.8	0.12	<1 week
	08/09/2013	09:49:00 PM	60	3.8	1.9	0.61	1 to 2 weeks
	09/09/2013	01:21:00 AM	22	0.4	1.1	0.06	<1 week
	09/09/2013	07:37:00 PM	10	0.4	2.4	0.75	<1 week
	10/09/2013	02:43:00 AM	8	0.4	3.0	0.29	<1 week
	10/09/2013	03:24:00 AM	9	0.4	2.7	0.02	<1 week
	10/09/2013	04:28:00 AM	98	5.2	3.2	0.04	2 weeks to 1 month
	10/09/2013	10:20:00 AM	13	0.6	2.8	0.18	<1 week
	10/09/2013	04:24:00 PM	6	0.4	4.0	0.24	<1 week
	11/09/2013	11:07:00 AM	12	0.4	2.0	0.78	<1 week
	11/09/2013	12:04:00 PM	60	1.4	1.4	0.03	<1 week
	11/09/2013	04:50:00 PM	10	0.6	3.6	0.16	<1 week
	11/09/2013	06:16:00 PM	18	0.4	1.3	0.05	<1 week
	11/09/2013	08:27:00 PM	208	11.4	3.3	0.08	1,5 to 3 months
	12/09/2013	03:19:00 PM	3	0.4	8.0	0.64	<1 week
	12/09/2013	04:25:00 PM	24	0.6	1.5	0.04	<1 week
	12/09/2013	08:44:00 PM	28	0.4	0.9	0.16	<1 week
	14/09/2013	08:51:00 AM	17	2.2	4.9	1.49	1 to 2 weeks
	14/09/2013	05:17:00 PM	10	9.8	58.8	0.33	1 to 2 years
	14/09/2013	07:30:00 PM	41	0.6	0.9	0.09	<1 week
	15/09/2013	03:46:00 AM	15	0.4	1.6	0.32	<1 week
	16/09/2013	04:54:00 PM	27	3.4	7.6	1.54	2 weeks to 1 month
	17/09/2013	08:08:00 PM	301	7.6	1.5	1.12	2 weeks to 1 month
	18/09/2013	04:23:00 AM	86	1.4	1.0	0.13	<1 week
	18/09/2013	04:04:00 PM	5	0.4	4.8	0.43	<1 week
	18/09/2013	06:43:00 PM	25	3.2	7.7	0.11	2 weeks to 1 month

* Nowadays Eurométropole de Strasbourg

Chapter 5 -Parameter identification-Inverse method

Acknowledgements

I would like to take the opportunity to thank everyone who supported me throughout the

APPENDIX A

Transformation from a reference element to an element in physical space and vice versa and streamline calculation with MHFEM

A.1 Transformation from a reference element to an element in physical space The system represents a function F 2D , that maps a point (u,v) in an uv-coordinates system into a point (x, y) in a xy-coordinate system.

APPENDIX E

The differentiated code generated by Tapenade (MHFEM)

Publications and communications

Publications :