Basic HTML Version
34
Geotechnical News • December 2012
www.geotechnicalnews.com
GROUNDWATER
Influence of element size in numerical studies of seepage:
Unsaturated zones, transient conditions
Robert P. Chapuis
In an unsaturated zone above the
water table, the hydraulic conductiv-
ity (
K
) can vary by several orders of
magnitude over a small suction (
u
)
range, or over short vertical distances.
In numerical models of unsaturated
seepage, small height elements are
needed: a previous paper examined the
convergence issues under steady–state
conditions (Chapuis 2012a). This
paper examines how the size of time
steps may influence the numerical
solution for a simple one–dimensional
(1D) transient example. It suggests a
procedure to select the time steps so as
to converge not only numerically, but
also towards the correct solution.
Background
A first paper (Chapuis 2010) explained
the problems with large elements,
especially when trying to study a
regional groundwater problem. A sec-
ond paper (Chapuis 2012a) provided
basic rules for assessing the influence
of small details especially in engineer-
ing projects. A third paper (Chapuis
2012b) explained how to select the
height of finite elements for unsatu-
rated zones under steady state condi-
tions. This paper is thus the fourth
dealing with the influence of element
size. General concluding remarks are
provided for the four short papers.
The example here is that of a verti-
cal column under transient condi-
tions. This 1D problem can be solved
numerically using many codes, but
correctly only by those using the
real conservation equation (Richards
1931), and not a simplified, linear-
ized, or modified equation. However,
each code has its own ways of treating
the Darcy and conservation equa-
tions, to select element sizes and time
steps, using more or less automatic
procedures. These internal features of
each numerical code are not discussed
hereafter. The results of this paper
were obtained using the same finite
element code as in the steady–state
paper (Chapuis 2012b). Similar results
can be obtained using any unsaturated
seepage code, although built–in spe-
cific procedures of a code may obscure
our understanding of convergence and
oscillation issues for transient condi-
tions.
The numerical study of unsaturated
seepage is more difficult than that
of saturated seepage. This happens
because seepage equations involve
the volumetric water content
θ(u
)
and unsaturated hydraulic conductiv-
ity
K
(
u
) versus pore water pressure
u, which are described by highly
non–linear functions. Of course, a few
numerical codes still make simplifi-
cations, and either use constants for
these functions, or strongly linearized
equations for highly non–linear physi-
cal phenomena that require highly
non–linear equations for their correct
description. To obtain correct solu-
tions (e.g., Chapuis and Dénes 2008;
Chapuis et al. 2005), the mesh must be
refined vertically in the vadose zone
to avoid convergence problems and to
overcome the difficulties associated
with the non–linearity of the equations
and hydraulic parameters.
Start with the mesh refinement
before refining the mesh step
The first action, when using any code
for unsaturated seepage, is to reduce
the element size in the unsaturated
zone to avoid potential numerical
oscillations, and ease the convergence
process. This is verified first using
steady–state (time independent) condi-
tions (Chapuis 2012b). As a rule of
thumb for meshing unsaturated zones,
the element height must not exceed the
value giving a maximum change of
one order of magnitude for
K
under a
no–flow condition.
The second action is to reduce the
time step to ensure proper conver-
gence towards the correct solution,
and not towards an incorrect solution.
This paper examines only the second
action, reducing the time step for tran-
sient conditions.
Example and numerical results
A 2 m high vertical column contains a
single soil with the functions of Fig.1,
the same functions as in Chapuis
(2012b) where a convergence study
under steady–state conditions has
shown that the element height must
be 10 cm or smaller to have small
numerical errors. For the transient
study, the initial conditions are steady–
state with the following boundary
conditions (BCs): the BC at
z
= 0 m
is
h
= 0 m; the BC at
z
= 2 m is
h
= 1
m. For the following transient–state,
the BCs are as follows: the BC at
z
=
0 m is
h
= 0 m; the BC at
z
= 2 m is a
periodic sinus–like imposed hydraulic
head with a period of 2 hours (Fig. 2).