A fully-explicit discontinuous Galerkin hydrodynamic model for variably-saturated porous media*

De MAET T., HANERT E.

Earth and Life Institute, Environmental Sciences, Université catholique de Louvain, B-1348 Louvain-la-Neuve,

Belgium

Georges Lema?tre Centre for Earth and Climate Research, Université catholique de Louvain, B-1348 Louvainla-Neuve, Belgium, E-mail:thomas.demaet@uclouvain.be

VANCLOOSTER M.

Earth and Life Institute, Environmental Sciences, Université catholique de Louvain, B-1348 Louvain-la-Neuve,

Belgium

A fully-explicit discontinuous Galerkin hydrodynamic model for variably-saturated porous media*

De MAET T., HANERT E.

Earth and Life Institute, Environmental Sciences, Université catholique de Louvain, B-1348 Louvain-la-Neuve,

Belgium

Georges Lema?tre Centre for Earth and Climate Research, Université catholique de Louvain, B-1348 Louvainla-Neuve, Belgium, E-mail:thomas.demaet@uclouvain.be

VANCLOOSTER M.

Earth and Life Institute, Environmental Sciences, Université catholique de Louvain, B-1348 Louvain-la-Neuve,

Belgium

(Received April 23, 2014, Revised June 10, 2014)

Groundwater flows play a key role in the recharge of aquifers, the transport of solutes through subsurface systems or the control of surface runoff. Predicting these processes requires the use of groundwater models with their applicability directly linked to their accuracy and computational efficiency. In this paper, we present a new method to model water dynamics in variably- saturated porous media. Our model is based on a fully-explicit discontinuous-Galerkin formulation of the 3D Richards equation, which shows a perfect scaling on parallel architectures. We make use of an adapted jump penalty term for the discontinuous-Galerkin scheme and of a slope limiter algorithm to produce oscillation-free exactly conservative solutions. We show that such an approach is particularly well suited to infiltration fronts. The model results are in good agreement with the reference model Hydrus-1D and seem promising for large scale applications involving a coarse representation of saturated soil.

3D subsurface flow model, discontinuous Galerkin method, slope limiters, explicit time integration

Introduction

A good understanding of subsurface water dynamics is essential in many hydrological, environmental and engineering applications. However, predicting such dynamics is still a difficult task. The difficulty mainly comes from the heterogeneity of the soil properties, the nonlinearity of the flow process and the absence of fast techniques to measure the hydraulic properties everywhere, at the appropriate scale. Given these issues, modeling the water flow in heterogeneous and often partially saturated porous media quickly and robustly is still a challenge.

By using a continuum approach, the water flow in variably saturated heterogeneous porous media can be modeled with the Richards equation (RE). This equation concerns both the water saturated zone (SZ) and the unsaturated zone (UZ). In the former, the equation models an incompressible fluid, leading to a constant water content θ. The pressure head hthen reacts instantaneously to the boundary conditions. In the latter, the equation represents capillary physics and is complemented with the so-called retentioncurve equation, which links the two variables. It is usual to refer to theθ-form and the h-form of the RE for the equations obtained by incorporating the retention curve in the conservation law to keep eitherθ orh as the sole model variable, respectively. On the one hand, the θ-form is not valid in the SZ, but is known to be more efficient in the UZ. On the other hand, theh-form is valid everywhere, although it is not mass conservative once discretized in time. The issues associated with theθ-and h-forms are usually overcome by combining both variables into a mixed formulation. When only the SZ is considered, the RE reduces to the groundwater equation that has its own numerical issues such as those associated with thetreatment of the free surface[1,2].

Despite the abundance of subsurface-flow models, the development of efficient and accurate discretizations of the RE is still an active field of research. The presence of two very different dynamics (in the SZ and UZ) leads to strong nonlinearity in the constitutive relations and then in the coupling between parabolic and elliptic partial differential equations. Implicit schemes have difficulties to handle such nonlinearities resulting in a lack of convergence[3]. Moreover, as the system dynamics can be highly transient, adaptive time steps are often mandatory[4]. In infiltration cases, this results in very small time steps, well below what is needed to reach a sufficient level of accuracy. As a result, the efficiency of implicit models for the RE is generally sub-optimal when infiltration fronts occur.

The efficiency of implicit models is further impaired by their poor scaling on parallel architectures. The current trend in computer designs is indeed to increase the performances relying on parallel architectures instead of enhancing the computational power of each processor individually. Current subsurface flow models therefore have to take advantage of parallelism, and some steps in that direction have been made with ParFlow[5,6], ParSWMS[7,8]and DuMux[9]. These models however only achieve sub-optimal scaling for the RE, except for some fully unsaturated test cases[8]. This is due to the use of implicit solvers as they require a large amount of communication between computational nodes to solve linear systems. Unlike implicit solvers, explicit solvers require only one exchange of information between nodes per time step. As implementing an explicit solver is simple, one can easily achieve a super-linear scaling, i.e., a scaling better than 100%, thanks to the additional computer caches coming from additional resources. Moreover, explicit solvers do not require any global matrix linear solver or the computation of a Jacobian, which are complex and costly.

In this paper, we present a 3D model of RE based on the discontinuous Galerkin (DG) finite element method (FEM) and an explicit time integration scheme. The model relies on slope limiters to locally ensure the monotonicity of the solution. It scales optimally at least up to 64 processors. A special treatment on the interface term between elements allows the existence of physical discontinuities in the water content between different types of soils. The model is mass-conservative at the machine precision. In the next two sections we present the mathematical formulation of the model and the explicit DGFEM discretization. In the fourth section, we present 1D and 3D numerical results, which focus on different physical and numerical aspects of groundwater flows.

1. Mathematical formulation

whereK is the largest eigenvalue of K andτis a free parameter that controls the relaxation towards the steady state defined by Eq.(10). When the general diffusivity tensorKis simply a scaled identity matrix (i.e., for isotropic soil), it obviously reduces toK. In Eq.(11), the relaxation parameter τcan then be interpreted as the diffusivity. As such, it will constrain the stability of any explicit time discretization of that equation, i.e. as τincreases, the time step should decrease. Reciprocally, for a given time step, it is possible to determine the maximum value ofτfor ensuring stability. It should be noted that the approximation of Eq.(10) by Eq.(11) is only made for numerical purposes in order to deal with a system that includes a parabolic component and an elliptic component. In the SZ, it physically amounts to increase the value of specific storage Ss.

If the h-form of Eq.(1) is used in zones where an explicit time integration scheme is stable and Eq.(11) is used otherwise, we can easily combine them as follows:

Equation (12) is exact when Cm=C , which is verified in most of the UZ but not in the SZ (except if Ks/τ≤Ss, which is unlikely). In this case, Eq.(12) is mathematically identical to the originalh-form. For practical purpose, the subset of the UZ whereCm=C will be called the dry zone (DZ) and the subset where Cm=K/τwill be called the wet zone (WZ). Equation (12) works in both the UZ and SZ with an upper threshold for the diffusion coefficient equal toτ, as shown in Fig.1. Although the h-form of the equation is not mass-conservative once discretized in time, the whole algorithm is mass-conservative at the machine error precision (see Subsection 2.1 related to the mass conservation).

In the time-stepping algorithm described in the next section, the model equations are complemented with the following equation

whereλis a free parameter andh～is the value ofh obtained after the time integration of Eq.(12), refreshed at each time step before the resolution of Eq.(14). Equation (14) allows us to spread out, in a mass-conservative way, the possible overshoots that can be caused by the approximation made in Eq.(13).

A weak form of the model equations is obtained by multiplying Eqs.(1), (12) and (14) respectively by the test functionsu,vand w∈H1(?). Taking the volume integral over the domain ?and using the divergence theorem, we obtain the following weak formulations of RE:

2. Space and time discretizations

The model equations are now discretized in space by means of a discontinuous Galerkin FE scheme, and then in time with the Euler explicit time integration scheme.

2.1 Discontinuous Galerkin space discretization

Through partitioning the domain?intoN nonoverlapping elements ?ewith interface Γe, the model variablesθandh can be approximated as

where Ndis the total number of degrees of freedom (DOF) and φjare piecewise polynomials defined on each element?esuch that

where xiis the vector of coordinates for the nodei and ne=Nd/Nis the number of DOF per element. Here we use piecewise linear (P1) basis functions for both variables. Since the discrete solution can exhibit discontinuities between mesh elements, the following jump[·]and averaging {·}operators on the interface Γehave to be introduced:where ??eis the contour of the element domain ?e. The material properties are assumed constant by element.

The DGFEM uses a piecewise linear approximation that allows discontinuities between mesh elements. The approximation however becomes continuous when the solution is smooth. Discontinuities between mesh elements usually appear when the spatial resolution is insufficient to represent sharp gradients such as the one appearing in infiltration fronts. In this case, the size of the jumps between the elements gives direct information of the local spatial error. To stabilize the diffusion term and enforce a weak continuity constraint between elements, penalty termsPand Piθhave to be added to Eqs.(21), (22) and (23)[12]. The continuity constraint onh can be expressed as

with n0the order of the FE approximation (in our case n0=1),Kscthe normal flux-oriented scalar conductivity,J the water flux, and lea characteristic length of the two adjacent elements. It should be noted that the exact solution for the pressure headhis continuous, while the water contentθcan be discontinuous between different soil horizons. One of the advantages of DGFEM is to have the possibility to represent these discontinuities. However, imposing a continuity constraint on θin a classical way would tend to smooth them out. The continuity constraint is thus only imposed onh with Eq.(29) to indirectly stabilize the Eq.(21). For the Eq.(23), the continuity constraint Piθis expressed as

where σdis similar toσbut with Kscreplaced by λ. The presence of fθ()weakly ensures similar jumps between the couplesθh+/θh-and f()θf(wàn)θ(). A special treatment, described in Subsection 2.5 related to the slope limiters, is needed at interfaces between media with different hydrological properties. To increase the stability, the mass matrices Mijin Eqs.(21) and (22) have been lumped. HenceCm,jis only present on the diagonal, which allows us to avoid inverting the productCm,jMijat each time step. Following other model designs, we use nodal element values forK.

An element is considered affiliated to the SZ, the WZ or the DZ in the following priority order: (1) if one node of the element is saturated (i.e.,h≥0), the element belongs to the SZ, (2) if one node of the element has the value Cm＜C, the element belongs to the WZ, (3) otherwise the element belongs to the DZ.

2.2 Mass conservation

As was mentioned above, the pureh-form of the RE is the easiest to solve but it is not mass-conservative. This is due to the time discretization of the Cfunction which is highly nonlinear. Here we propose a simple solution to this issue based on the application of Eq.(1) as a post-processing step to achieve mass-conservation.

Using the same discretization as for Eq.(1) and an a priori non-conservativeh-form as Eq.(12), we can yield two equations with the same right-hand side:

2.3 Explicit time discretization

Equations (21) and (22) are discretized in time with an explicit Euler scheme. Using a matrix notation, the overall solution procedure is the following:

where we have replaced the reference value “H～”used in Eq.(23) byHn+1.

In the DZ, this algorithm reduces to theθ-form of RE. This form is known to be numerically better suited than theh-form, especially in very dry cases. In the SZ, it reduces to an approximation of thehform which exponentially converges to Eq.(10) asm increases. The resulting approximation error leads to approximate fluxes in the SZ but the total mass is exactly conserved. Step 4 allows us to bringθtowards its correct value in both SZ and WZ.

It should be noted that when Eq.(36) has converged in the SZ, the result is equal to the one obtained with an implicit time integration scheme. Indeed, if Hn+1,m=Hn+1,m-1, we have reached the incompressible state which corresponds to the solution of the implicit equation forh , and also forθonce inserted in Eq.(38). To satisfy this property, the right-hand sides of Eqs.(21) and (22) have to be identical. Another point is that increasingmdecreases directly the strength of the approximation. Indeed, the division of Δt bym in Eq.(36) allows us to magnifyτby a factor m.

2.4 Selecting the values ofτandλ

Any explicit time integration scheme of Eq.(21) is fully mass-conservative and valid both in the UZ and SZ but requires the field hjto be known. Inverting the retention curve, i.e.,h=fθ-1(θ), is the ea siest way to obtain it. However there is an underlying Courant-Friedrich-Lewy (CFL) stability condition. For instance, the explicit Euler time discretization reads:

where the superscripts indicate the time step,Δtrepresents the time step andF(h)the right-hand side of Eq.(11).

One could see that this collapses exactly into the θ-form, and simple manipulations show that this is unstable whenK/Cis beyond the stability limit for diffusion, which is

where Δlis the smallest element length and 0＜p＜1depends on the type of explicit solver used. Although our algorithm is slightly more complex than that given by Eqs.(43) and (44), we hypothesize that its stability condition is similar and, with Eq.(13), it leads to the following condition onτ

where hmaxis soil-dependent.

The parameterτcontrols the position of the separation between the DZ and WZ. As the solution in the WZ is approximated, it is better to limit its extent as much as possible by increasingτ. However, to respect Eq.(46) for a constant mesh size Δl, increasing τleads to decreasing the time step. In the following, we take p =1/5in 1D case and p=1/15in 3D case.

The quasi-elliptic approximation produces fluxes in the SZ resulting in a local mass excess or deficit which has to be corrected with Eq.(41). This equation has the following stability condition

To optimize the correction, we take the explicit diffusive limit, i.e., we take the larger stable value of λ given the sizes of the adjacent elements, in a similar way as forτ. It has been observed that a smaller value pθ=1/20has to be taken for stability, due to complex interactions with the UZ/SZ interface.

The present model could develop unphysical fluxes in the WZ and SZ because of the approximation in Eq.(12). The parameterspandpθare a priori fixed once for all, andτandλcould then be fixed by the choice of Δt. The free parameters which act on the approximation are thenΔt,m andq. One could easily obtain a local error of the quasi-elliptical approximation

where overlines represent the mean over one element, the indiceseorj the affiliation of the variable to the elemente or to the node j. This expression is correct and mass conservative only if the element has the property that the sum of its nodes weighted by their number is the mean of the element, i.e.,∑jUej= neUe. All elements present in the following test cases respect this rule.

The limited valuesU?are unchanged (i.e.,L=eje 1) if the values inside the element (i.e.,Uej,?j∈e) are between the maximum and the minimum of the means of their neighbors. Otherwise, the value ofLeis less than one, and the reduction in the gradients applies in all directions. It would be possible to increase the accuracy of this limiter by considering separately each spatial direction and space derivative. However, the simplicity of this algorithm has been preferred to reduce the computational cost. A 1D example is shown in Fig.2.

A different behavior is expected on the boundaries of the domain. Only tangential components to the boundary have to be limited while the normal one should not. Indeed, we want to stabilize the scheme, without limiting the extrema on the boundaries as they appear frequently in physical cases. This is achieved by using a mirror image of the solution outside of the domain.Eventually, special care must be taken for the limiter and continuity constraint applied toθ. Indeed, both have to preserve the physical discontinuities ofθbetween different mediaAandB . The way to proceed is to translateθintoh by means of the inverse of the retention curve(θ), as h has the property to be physically continuous. Onceθhas been translated in terms of hwith the properties of the medium it belongs to i.e.AorB, it is translated back into θusing the properties of the medium where the computation is necessary (resp.BorA ). Then, if subscripts are used to represent the medium, one could define as

the continuity constraint and extrema per node. The differences for the continuity constraint and the means for the limiter are then coherent with physical discontinuities ofθ.

3. Numerical examples

In this section, we present three numerical examples demonstrating the ability of the model to produce results similar to the widely-used model Hydrus-1D[18], confirming its convergence as the number of iterations increase, and showing its scalability in a 3D application.

3.1 Unsaturated infiltration

This simulation is based on data numerically reproduced as a benchmark test in the Hydrus-1D code[18]. The experimental setup consists in a homogeneous column of soil of 0.6 m that is assumed to have an initial constant pressure head of -1.5 m. The characteristics of the soil are the same as in the Hydrus code and are represented by the modified van Genuchten-Mualem relations from Vogel and Cislerova (1988) summarized in Appendix, with the parameter values given in Table 1. The material properties are homogeneous and isotropic. At the beginning of the simulation, the pressure of a thin layer of water with the height assumed to be approximately zero is imposed at the top of the column. Mathematically, this is done by imposing the Dirichlet boundary condition h(z=0)=0. A zero-flux boundary condition is imposed on all the other boundaries. The spatial discretization is made with 30 equidistant layers in the vertical direction. Here, a time step of 1 s was used, with only one sub-iteration forhandθ(i.e.,m= q=1).

The mass balance error of Hydrus at the end of the simulation is 1.95×10-6m versus 2×10-16m for our code (with a maximum at 6×10-16m). Spatial or temporal convergence studies are difficult to fulfill with the proposed model. Indeed classical convergences orders are disrupted by the false transient approximation in the WZ and SZ, as will be shown in the next subsection. Those zones should therefore be limited to isolate proper convergences when keepingm andq constant as they are not the studied variables. Since the size of the WZ increases with the time step and decreases with the mesh size, we have to impose a maximum time step of 1 s and minimum mesh size of0.01 m. For the mesh convergence study, the top fifth of the domain is kept fixed with a discretization of 0.02 m to make sure the WZ is always discretized with the same resolution and hence prevent the false transient approximation from interfering with the convergence analysis. Figure 4 shows a temporal convergence rate of 1 and a spatial convergence rate of 1.5. The first order explicit time integration scheme behaves as expected but it is not the case of the spatial integration scheme which is theoretically second-order accurate. Such a discrepancy is likely due to the slope limiters.

3.2 Filling of groundwaters

This second test case highlights the effects of the“quasi-elliptic” approximation in the SZ given by Eq.(11). We consider a soil of one-meter depth described by four equally-thick layers of sand, loam, clay and loam whose properties are given in Table 2. Each layer is assumed homogeneous and isotropic. Initially, the groundwater fills half of the domain and the system is at physical equilibrium. This is achieved by taking a linear initial pressure field going from 0.5 m at the bottom to –0.5 m at the top. A Neumann boundary condition is imposed at the top to represent a strong infiltration of 10-5m/s that stops after 2 h. All the other boundaries are impervious. The soil column is discretized into 40 equidistant layers, the time step was equal to 1 s.

Figure 5 shows the evolution of the h and θ profiles during the simulation. The discontinuity between different types of soils in the θ-profile is in good agreement with the physics, while a classic continuous representation would have to rely on the mean of the soil properties at the interfaces. In this model each DOF belongs to a particular material and all properties are well defined, even at the interfaces between different soil layers.

The results are similar for both models until the infiltration front reaches the groundwater depth, at aroundt =2.15 h. Before this point, we set m=q =1.

The succession of the large conductivity of the sand and the conductivity of the loam inferior to the incoming flux causes an accumulation of water at this interface. However, some of water passes through this capillary barrier and eventually fills the loam layer. Once the loam layer has been filled, the depth of the incompressible water column instantaneously increases from half to nearly 3/4 of the computational domain.

After t =2.15, the h-profile is not instantaneously adapted as it should, and a curvature appears, as shown in Fig.6. This generates fluxes in the SZ leading to an increase in the water content, which is visible in the last panel of Fig.5. However, this excess is spread over time and the model converges towards the steady solution. Increasingmandqimproves the convergence of the solution, as shown in Figs.6 and 7. It can be seen that a value ofm≈25 is sufficient on this quick event to transport the information through 30 saturated elements.

The mass balance error of Hydrus-1D at the end of the simulation is 4.78×10-6m versus 5×10-16m for our code (with a maximum at 10-15m). For this result, the limiter onθhas been removed. It adds robustness to the algorithm but is not mandatory. The limiter applies some additional multiplications and additions to the mass variableθ, which increase the mass balance error. Despite that, with the limiter the mass balance error is not larger than 1.2×10-14m.

3.3 Capillary barrier in a simple landfill design

For this theoretical test case, a simulation of the water dynamics within the simple landfill represented in Fig.8 is considered. The waste is stocked over a capillary barrier made of sand and gravel. An impervious geotextile placed under the gravel diverts the water to the bottom where it is drained. The efficiency of the capillary barrier is assessed for an important flux of water into an initially dry system. We take advantage of the symmetries of the computational domain to model the problem only over one eighth of the domain. The waste layer is defined by soil-like properties given in Table 3 with the properties of the sand and gravel. The material properties are considered homogeneous and isotropic within each layer. The mesh is constituted of nine layers of 2 892 prisms for a total of 156 168 DOF, refined over corners as shown in Fig.8. Water input is taken as a homogeneous flux of 0.005 m/h on the top of the waste that is stopped after 12 h. The geotextile is assumed to be perfectly impervious while the drain continuously covers the bottom of the landfill. The drain boundary condition is defined as follows:

wherehcis the air-entry pressure of the gravel layer. It should be mentioned that such a boundary condition is quite difficult to apply in an implicit model. Indeed, conditional statements frequently produce oscillatory states into the convergence process, as a continuous and monotone system is theoretically required to ensure convergence. The mesh used is made of prisms obtained by extruding a 2D triangular mesh. Here, we setm=q=1and use an adaptive time step which varies from 50 s at the beginning of the simulation to 0.075 s when the gravel, which is very conductive, reaches saturation. The adaptive time step algorithm simply set the time step as the minimum between 50 s and that given by the CFL condition.

Water infiltration at the beginning of the simulation is visible on Fig.9. As expected, water crosses the capillary barrier first at the external lower corner of the domain and at the bottom of the slope, where it accumulates before being drained out, as shown on Fig.10. In the middle of the slope plane, the capillary barrier plays its role and diverts most of the water which slowly flows to the bottom of the sand layer, towards the drainage zone. The effect of the capillary barrier is also visible at the top of the corner, where a small area stays dry. On such a configuration, the most sensitive parts are the corners where the geotextile has to divert the strongest fluxes and sustains the largest water.

3.4 Parallel efficiency

The first reason for the suboptimal scaling observed with all these models is the reduction of the computational domain (CD) related to each node, compared to the extent of the interfaces (I) between different CDs which require communication between nodes. A second reason is the elliptic behavior of the equation in the SZ, which implies that any change impacts the entire SZ. The information has to pass through several CDs to cover the whole SZ. It therefore requires a lot of communication (often via additional linear solver iterations). The first case could be magnified by the so-called “strong” scaling test cases, which keep thesame mesh when increasing N. Then the ratio I/CD increases and the scaling is expected to decrease. The second case is magnified by “weak” scaling test cases, which keep the sizes of the CD and the I constant when increasing N and thus the domain size as well. This results in an increase in the number of CD, which in turn increases the number of linear solver iterations and reduces the scalability. But strong scaling test cases are also affected by this situation. That is why we choose to show a strong scaling test case.

In studying strong scaling, it is important to look at the number of DOF per node, as it is closely related to the ratio I/CD. Any model will see its performances degraded below a certain ratio I/CD for which the communication time begins to dominate the overall computation time. It is then more difficult to achieve a good scaling with a small number of DOF per node. Strong scaling test cases could benefit from additional fast-memory caches asN increases, as the memory load per node decreases accordingly. This can lead to super-linear scaling.

It is conceptually difficult to develop a model of an elliptic-type equation with a perfect scaling. Mathematically, any small local change in the model solution will have a global influence. This leads to a strong exchange of information between the sub-domains linked to each computational node, and thus a poor scaling. The multigrid method is specially designed to overcome this issue. The false transient method that is used here simply reduces this exchange and thus achieves a perfect scaling by transferring information at a finite speed. The consequence that the information in an explicit method is transferred only to neighboring elements is that the number of iterationsm could increase and hence reduce the efficiency as compared to implicit methods. The model efficiency would therefore be reduced in cases where the SZ occupies a large portion of the domain or is discretized with a large number of DOF’s. The minimum of m in this case is of the order of exp(Nc)where Ncis the number of saturated elements in a line.

The scaling of the model has been evaluated on a cluster of 64 nodes by running the model on a testcase similar to the last one, except for the 3D mesh which is now composed of 101 508 prisms, or 609 048 DOF. Speed-up results are shown in Fig.11 and exhibit an optimal scaling.

4. Conclusions

We have developed a 3D discontinuous-Galerkin model of the Richards equation using an explicit solver. As the mass variableθis updated with the mixed-form of the equation, the mass is conserved up to machine-error precision. The DGFEM is well-suited to the advection-dominated physics that occurs at sharp gradient fronts. The DGFEM also allows the use of stabilizing techniques such as slope limiters. These limiters remove the oscillations which appear with the classical FEM. A discontinuous approximation can also nicely capture the physical discontinuities of the water content. In our model, the UZ is modeled with the θ-form of RE and the SZ by an approximation of theh-form of this equation. Such an approximation, which results in a finite pressure propagation speed, allows us to use the same explicit time integration scheme for both the saturated and unsaturated zones.

Accuracy and efficiency could be further improved by using a better method to solve the purely elliptic part of the equation, such as the multigrid method. Indeed, if the SZ occupies a large part of the domain, the present model could require an important number of iterations to converge. That would dramatically reduce the overall efficiency of the method. As the elliptic part is linear, a linear solver alone is sufficient, without the burden of a non-linear solver. Speed improvements could be achieved by using an interpolation for the constitutive relationships, as the power function is very time-consuming. This could be alsoaccomplished by using a multi-rate method, which consists in reducing the time steps only on the part of the domain where the physical processes have the smallest time scale[21]. The latter seems very promising for further developments.

Acknowledgements

De Maet T. is a research fellow with the Fonds Spécial de la Recherche (FSR) of the Université catholique de Louvain. The present study was carried out in the framework of the project taking up the challenges of multi-scale marine modelling, which is funded by the Communauté Fran?aise de Belgique under contract ARC 10/15-028 (Actions de recherche concertées) with the aim of developing and using SLIM (www.climate.be/slim). Computational resources have been provided by the supercomputing facilities of the Université catholique de Louvain (CISM/UCL) and the Consortium des Equipements de Calcul Intensif en Fédération Wallonie-Bruxelles (CECI) funded by the Fond de la Recherche Scientifique de Belgique (FRSFNRS). The authors gratefully acknowledge Lambrechts J., De Brye B., Comblen R. and Javaux M. for their comments and help provided during the preparation of the paper.

[1] ZHANG Qian-fei, LAN Shou-qi and WANG Yan-ming et al. A new numerical method for groundwater flow and solute transport using velocity field[J]. Journal of Hydrodynamics, 2008, 20(3): 356-364.

[2] LIN Lin, YANG Jin-zhong and ZHANG Bin et al. A simplified numerical model of 3-D groundwater and solute transport at large scale area[J]. Journal of Hydrodynamics, 2010, 22(3): 319-328.

[3] LOTT P. A., WALKER H. F. and WOODWARD C. S. et al. An accelerated Picard method for nonlinear systems related to variably saturated flow[J]. Advances in Water Resources, 2012, 38: 92-101.

[4] KAVETSKI D., BINNING P. and SLOAN S. Adaptive time stepping and error control in a mass conservative numerical solution of the mixed form of Richards equation[J]. Advances in Water Resources, 2001, 24(6): 595-605.

[5] KOLLET S., MAXWELL R., Integrated surface-groundwater flow modeling: A free-surface overland flow boundary condition in a parallel groundwater flow model[J]. Advances in Water Resources, 2006, 29(7): 945-958.

[6] MAXWELL R. A terrain-following grid transform and preconditioner for parallel, large-scale, integrated hydrologic modeling[J]. Advances in Water Resources, 2012, 53: 109-117.

[7] HARDELAUF H., JAVAUX M. and HERBST M. et al. PARSWMS: A parallelized model for simulating threedimensional water flow and solute transport in variably saturated soils[J]. Vadose Zone Journal, 2007, 6(2): 255-259.

[8] HERBST M., GOTTSCHALK S. On preconditioning for a parallel solution of the Richards equation[J]. Computers and Geosciences, 2008, 34(12): 1958-1963.

[9] FLEMISCH B., DARCIS M. and ERBERTSEDER K. et al. DuMux: DUNE for multi-phase, component, scale, physics,... flow and transport in porous media[J]. Advances in Water Resources, 2011, 34(9): 1102-1112.

[10] HANERT E., LEGAT V. How to save a bad element with weak boundary conditions[J]. Computers and Fluids, 2006, 35(5): 477-484.

[11] BAZILEVS Y., HUGHES T. Weak imposition of Dirichlet boundary conditions in fluid mechanics[J]. Computers and Fluids, 2007, 36(1): 12-26.

[12] SHAHBAZI K. An explicit expression for the penalty parameter of the interior penalty method[J]. Journal of Computational Physics, 2005, 205(2): 401-407.

[13] CELIA M., ZARBA R. and BOULOUTAS E. A general mass-conservative numerical solution for the unsaturated flow equation[J]. Water Resources Research, 1990, 26(7): 1483-1496.

[14] RATHFELDER K., ABRIOLA L. Mass conservative numerical solutions of the head-based Richards equation[J]. Water Resources Research, 1994, 30(9): 2579-2586.

[15] KEES C., FARTHING M. and DAWSON C. Locally conservative, stabilized finite element methods for variably saturated flow[J]. Computer Methods in Applied Mechanics and Engineering, 2008, 197(51-52): 4610-4625.

[16] COCKBURN B., SHU C. Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[J]. Journal of Scientific Computing, 2001, 16(3): 173-261.

[17] AIZINGER V. A geometry independent slope limiter for the discontinuous Galerkin method[C]. The 4th Russian-German Advanced Research Workshop. Freiburg, Germany, 2009, 207-217.

[18] SIMUNEK J., SEJNA M. and Van GENUCHTEN M. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media[R]. University of California-Riverside Research Reports, 2005, 1-281.

[19] VANDERKWAAK J. Numerical simulation of flow and chemical transport in integrated surface-subsurface hydrologic systems[J]. Dissertation Abstracts International Part B: Science and Engineering, 2000, 60(7): 3170.

[20] BERNINGER H., KORNHUBER R. and SANDER O. Fast and robust numerical solution of the Richards equation in homogeneous soil[J]. SIAM Journal on Numerical Analysis, 2011, 49(6): 2576-2597.

[21] CONSTANTINESCU E., SANDU A. Multirate timestepping methods for hyperbolic conservation laws[J]. Journal of Scientific Computing, 2007, 33(3): 239-278.

Appendix: Modified Van Genuchten-Mualem relations

Following Vogel and Cislerova (1988) the retention function and the conductivity functions are described as follows:

10.1016/S1001-6058(14)60067-6

* Biography: De NATE T. (1986-), Male, Ph. D. Candidate