Introduction

The objective of this article is to give a brief overview of core groundwater flow concepts implemented in PLAXIS and also to provide helpful tips for carrying out a groundwater flow analysis in saturated and unsaturated conditions.

Figure 1. Groundwater flow in an embankment dam

A groundwater flow analysis involves a flow of water in saturated as well as in unsaturated soils. The Finite Element method in PLAXIS can also be used for groundwater flow modeling, such as dewatering, flow through an embankment, water flow in a tunnel and many other cases. These analyses can be performed in a steady-state model or in a time-dependent model.

PLAXIS can also be used to calculate leakage through structures such as sheet pile walls or cut-off walls. Such groundwater flow analyses require the additional parameters for the constitutive materials, such as permeability values and its variation with suction, porosity and unsaturated zone behavior along with specific flow boundary conditions such as prescribed head of water flow, closed flow, seepage flow boundary, drains, wells (both infiltration and extraction), and precipitations.

Soils can be considered as dry (all voids are filled with air), fully saturated (all voids are filled with water) or unsaturated. In unsaturated soil, voids are considered to be filled with both liquid (water) and gas (air). Generally, when the soil is above the water table it is considered as unsaturated and it has a degree of saturation smaller than 1.0.

Suction

Suction is caused by the capillary flow against gravity, so the suction at a point depends on how deep the saturated flow is below that point. Pore water pressures are considered to be negative in compression. In PLAXIS suction pore water pressures are positive with respect to the atmospheric pressure. Below the phreatic level, the pore water pressures are negative and the soil is usually saturated.

Total suction is the summation of matric and osmotic suction of which osmotic suction for most practical applications is considered zero. So, total suction is considered equal to the matric suction which is related to the soil matrix (absorption and capillarity due to soil matrix) and the difference between soil water pressure and gas pressure:
S = p_a – p_w

Where p_a and p_w are the pore air pressure and the pore water pressure. In most cases, pore air pressure is constant and small enough to be neglected. Therefore, the matric suction is negative of the pore water pressure:
S =  – p_w

Groundwater flow

The flow in saturated soils is commonly described using Darcy’s law in which the rate of water flow through a soil mass is proportional to the hydraulic head gradient. In an unsaturated state, the hydraulic conductivity is dependent on the soil saturation. The relative permeability k_rel(S) is defined as the ratio of the permeability at a given saturation to the permeability in a saturated state:

k = k_rel k_sat

As the water flow in the unsaturated zone is dependent on soil saturation, hence groundwater flow modelling in this zone requires a function that relates hydraulic conductivity to the saturation level. PLAXIS has five types of hydraulic models that relate hydraulic conductivity to soil saturation. These are Van Genuchten, Mualem (simplified Van Genuchten), linearized Van Genuchten, spline and fully saturated.

For a fully coupled flow and deformation analysis, a more general formulation according to Biot’s theory of consolidation is used enabling users to simultaneously calculate deformation and groundwater flow with time-dependent boundary conditions in saturated and unsaturated soil.

The equations for hydraulic models and fully coupled flow deformation analysis can be found in the following article: Groundwater flow, fully coupled flow deformation and undrained analyses in PLAXIS 2D and 3D.

Boundary conditions

There are numerous boundary conditions that are applied either on the outer boundaries of the model or applied internally inside the model. A phreatic surface can be specified in the model or it can be determined as a result of the specified boundary conditions. These boundary conditions are summarized below:

Boundary condition	Information
Closed	Zero Darcy flux over the boundary.
Inflow	Non-zero Darcy flux over the boundary is set by a prescribed recharge value.
Outflow	Non-zero Darcy flux over the boundary is set by a prescribed Darcy flux.
Head	Groundwater head can be specified.
Infiltration/Evaporation	An inflow value may depend on time and as in nature the amount of inflow is limited by the capacity of the soil. If the precipitation rate exceeds this capacity, ponding takes place and the boundary condition switches from inflow to prescribed head. As soon as the soil capacity meets the infiltration rate the condition switches back.
Seepage	A seepage surface is used for free water level for example the downstream toe of an embankment dam where water exits the downstream slope surface. At the seepage surface, the groundwater head, h, equals the elevation head and so the water pressure is zero which is the same condition that exists at the water level. Seepage surface can be assigned to the whole boundary where seepage is expected to occur. Seepage boundaries may therefore be specified at all boundaries where the hydraulic head is unknown.
Infiltration well	Inflowing water flux per meter is specified.
Extraction well	A discharge simulating amount of water leaving the domain is specified.
Drain	Drains are like seepage boundaries. However, they are located inside the domain and cannot permit water leaving the domain at atmospheric pressure. Therefore, a prescribed head should be considered for the part of the drain below the water level.
Discontinuity element	Discontinuities can be drained or non-porous. When set to drained they always have an infinite cross permeability.
Interface element	Interfaces can be permeable, impermeable or semi-permeable. More information on interfaces and their use in groundwater flow can be found in the following article: Permeability in interfaces.
Time-dependent	PLAXIS provides several features for analysis of transient groundwater flow and fully coupled flow-deformation problems with varying conditions in time. Seasonal or irregular variations in water levels can be modelled using linear, harmonic, or user-defined time distribution. Four different functions can be assigned for this purpose, namely constant, linear, harmonic, and user-defined functions.

Analysis of groundwater flow

There are two ways of proceeding with a groundwater flow analysis. The first option is to start the analysis with a ‘flow only’ analysis as the initial phase. This type of analysis precludes the soil stiffness and so deformations are not calculated in any subsequent phases. This is useful when only uncoupled groundwater flow is required.

The second option is when soil deformation is also considered and for that, the initial phase is set to either a gravity/ko/field stress to generate the initial stresses. All subsequent phases will then consider soil deformation. In this type of analysis, a fully coupled flow deformation or a consolidation analysis can also be performed.

Both of these options are briefly described below.

• Flow only analysis: In the initial phase, there is no option to include or exclude suction as it is included by default. Pore pressures for the initial phase can be calculated by assigning Phreatic level or a steady-state groundwater flow. The subsequent phases can have a Phreatic level or use pressures from the previous phase or it can be a steady-state groundwater flow in which the hydraulic head and the hydraulic conductivity at any point in the soil mass remain constant with respect to time or it can be a transient groundwater flow analysis in which the hydraulic head (and possibly the hydraulic conductivity) changes with respect to time.
• Gravity/K0/Field stress: In the initial phase there is an option to include or exclude suction which by default is ignored. Pore pressures for the initial phase can be phreatic or steady-state groundwater flow. The subsequent phases can be a Plastic analysis or fully coupled flow deformation analysis or a consolidation analysis. In the plastic analysis, suction is included by default and there is an option to include or exclude suction. Pore pressures can be from a phreatic level or pressures from the previous phase or steady-state groundwater flow. In the fully coupled flow deformation analysis suction is included by default.

Water flow in the unsaturated zone when only saturated hydraulic conductivity is defined

If no unsaturated hydraulic conductivity is defined in the model, PLAXIS still can calculate small unsaturated flows that are present above the water table. The PLAXIS Groundwater tab in the materials dataset contains some "default" parameters for the Van Genuchten model. The sets contain parameters for different soils (coarse/fine/....), see Reference Manual Section 6.1.5 Groundwater tabsheet. For the calculation of the relative permeabilities according to Van Genuchten, see Material Models Manual Section 19.1 Van Genuchten model. In the background, these parameters are generated which results in a reduction of the permeabilities above the phreatic line. This means that (small) flow above the phreatic line is possible.

Modelling a wetting front in a transient analysis

Models involving a moving wetting front in a transient analysis require a dense mesh to give accurate results.

In the figure on the left, a coarse mesh was used and we can see that the wetting front is much more turbulent than the one in the figure on the right which has a much denser mesh.

Relationship between element size and time step

The accuracy of results in a transient analysis and fully coupled flow deformation analysis is influenced by the element size as well as the time step size. The calculation automatically detects a minimum (critical) time step and to some extent, a maximum time step based on the consolidation theory stated in Section 4.4 of the Scientific manual. It depends on element size, hydraulic conductivity, soil and water compressibilities.

When performing a Fully coupled-flow deformation analysis, by default the Time step determination in the phase’s Numerical control parameters is set to Automatic. The time step considered for first, min and max is used to start the calculation are automatically determined by PLAXIS. However, a user is allowed to switch to Manual by deselecting the Use default iter parameters option. In those cases, different results can be observed, especially for the 1st time step as PLAXIS will directly use what is provided as input.

The default time step size is calculated based on the average element size in the model. This means that when manually changing the default time step size to a larger time step size, the element size would need to be made small or in other words, the mesh needs to be finer, to prevent inaccuracy introduced.
When switching to Manual Time step determination you can also click on one of the time step fields (first, min or max) and a side panel will appear with the option to retrieve the suggested parameters for the analysis.