Groundwater flow analysis in PLAXIS – Brief overview

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:

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.