Free Surface
|
White Paper Vol. I, Chapter 2 |
Two-dimensional free surface definition
For 2D models, the settings for unconfined/confined aquifers are defined in Problem Class. In case of an unconfined 2D model, the Head limits for unconfined models and the Residual water depth are defined on this page.
Three-dimensional free surface definition
If the setting Unconfined aquifer (s) is chosen, two basic modes to consider the phreatic surface are available for 3D models:
As also combinations of the two modes are possible in very specific cases, the settings for each slice has to be defined separately. The setting for a slice determines the behaviour of the layer below. The last model slice is always set to fixed. Four different possible states are available:
- Free:
Only available for the top slice of the model. A movable slice tops an unconfined layer. The surface slice follows the water table according to the free & movable method. - Phreatic:
Applies phreatic mode to an unconfined layer. Typically only applicable to the top slice. Exceptions are discussed in the detailed description of the phreatic method. - Dependent:
Setting is defined by the nearest non-dependent slice above. Slices below a free slice may move, too, layers below a phreatic layer may also be partially saturated. - Confined:
A confined slice cannot move. The layer below a confined slice is always calculated as fully saturated. This setting can be used to limit slice movement in free & movable mode to slices above a slice that is known to be always saturated.
|
The phreatic and free & movable mode which consider a phreatic surface in a basically saturated model based on Darcy's equation are mostly applicable to typical groundwater modeling cases with a regional water table. For more specific cases involving perched water tables or unsaturated areas below saturated aquifers typically an unsaturated simulation applying Richards' equation is required. |
|
Models with a partially or fully unstructured mesh cannot use any of the options for unconfined aquifer (s). They need to be run with the unsaturated/variably saturated flow option using Richards' equation instead. |
Free & movable
The phreatic surface is considered by moving the top boundary of the model in a way that the elevation of the first model slice always corresponds to the water table elevation.
In case that the water table drops below the first model layer, not only the first slice is moved, but underlying dependent slices may also move. At each time step, material properties for each element are determined by the actual location of each computational layer. For example, if the entire first stratigraphic layer is dry, the first computational layer will be in the second stratigraphic layer, inheriting all material properties of this. As the technique for moving the slices (BASD - Best Adaptation to Stratigraphic Data) focuses on original slice elevations, for most of the elements this inheritance is unambiguous. In elements, however, where the computational slice crosses an original slice location, an averaging of parameters from upper and lower layers has to be done. As a simple average is used, this can lead to making an aquitard conductive by artificially increasing its low conductivity values.
Here are the advantages and disadvantages of this approach:
![]() |
The model only contains the saturated zone. No additional assumptions are required for the capillary fringe and the unsaturated zone. |
![]() |
In cases where the water level does not stay within one stratigraphic layer, parameter interpolation can occur during the simulation in elements where the water table cuts through slices of the initial stratigraphy. |
|
Free & movable is a particularly good method for cases where water level movements are expected within one layer. In complex aquifer-aquitard systems with steep gradients the parameter interpolation resulting from slice movement may not be acceptable. |
Phreatic
In phreatic mode, the model stratigraphy is not changed if the water table changes. As a result, partially saturated or unsaturated elements may occur in the model domain. They are taken care of by applying a partial-saturation approach: In partially dry elements, the conductivity in the element (in all directions) is reduced by multiplying the saturated conductivity with the saturation of the element. Saturation is calculated from saturated thickness divided by element height. Completely dry elements (hydraulic head is lower than the elevation of the bottom nodes) a minimum (residual) saturated thickness is used for the saturation calculation, the Residual water depth. The Residual water depth (by default approx. 1 mm) is constant over the entire model domain and can be defined by the user.
The advantages and disadvantages of the phreatic approach are:
![]() |
This method does not suffer from parameter interpolation problems in case of steep hydraulic gradients and a water table cutting through slices. |
![]() |
Linear conductivity reduction for partially saturated elements may lead to high contrasts in actual hydraulic conductivity between dry and saturated layers. This complicates the numerical solution. In cases where recharge is applied to a dry top layer, oscillations may occur. |
|
It is not possible to set subsequent slices to phreatic or free as this would lead to an overestimation of the unconfined storage as FEFLOW assumes a separate water table for each slice that is set to phreatic or free. For a correct slice definition for an unconfined aquifer spanning several layers, see the example below. |
The correct slice definition of a model consisting of an unconfined aquifer spanning 5 layers looks as follows:
Configuration of slices in the Phreatic section.
An unconfined aquifer spanning several layers is correctly defined by setting the aquifer top slice to phreatic or free and all other slices of that aquifer to dependent. This means that the location of the free surface can change during the simulation run and that it will be located in one of the layers belonging to the aquifer.
The advantage of this free surface definition is that it is not necessary to know in advance which layer the free surface will be located in during the simulation run. FEFLOW will automatically derive the values for drain-/fillable porosity from the layer where the water table is located at a given time to correctly account for storage changes.
|
The phreatic mode can also be used to model different water tables, one over another. Use, for example 'free & movable' at the top and 'phreatic' for the slice covering the layer where the second water table is expected. The top of the aquitard between the two aquifers has to be set to 'fixed'. Multiple water tables can only be modeled correctly using the phreatic method if the upper aquifer holds water over the entire model domain. The lower aquifer can be fully unconfined or partially confined. Perched water tables in parts of the model domain can only be simulated by doing an unsaturated simulation! |
Storage change in phreatic top layer where water table exceeds surface
Two different settings are available:
- Extend storage of unconfined layer to the water table: The local drain-/fillable porosity of the top layer is assumed for the virtual above-ground part as well.
- Treat as confined: In the confined sections of the top layer, storage changes only due to compressibility effects.
Sources/sinks in phreatic layers
Source and sink parameters assigned within phreatic layers can be scaled by the local pseudo-saturation. This will result in reduced source/sink values wherever the hydraulic head falls below the element top elevation (elements are not fully saturated).
By default, the scaling of source/sink values is enabled.
Head limits for unconfined conditions
For unconfined models in 2D and 3D, the model behaviour has to be defined for cases where the the model bottom is completely dry at a node, and where the water level exceeds the model top.
Head limits at model bottom
Option | Description |
---|---|
Unconstrained head (default) | The water level is allowed to drop below the model bottom. The remaining transmissivity for otherwise dry elements of the model is calculated from the conductivity multiplied with the residual water depth. |
Constrained head | The water level is not allowed to drop below the bottom elevation. In otherwise dry mesh nodes on the model bottom, a hydraulic head boundary condition with a value of bottom elevation plus residual water depth is set automatically . |
Head limits at the model top
Option | Description |
---|---|
Unconstrained head (default) | The water level is allowed to exceed the model top by extending the aquifer up to the actual water level. |
Constrained head | The water level is not allowed to exceed the model top. At the nodes at which the water level would otherwise rise above the top, a hydraulic-head boundary condition with a value of the top elevation is set automatically. Thus, all water exceeding the top is removed via these boundary conditions. |
|
Applying the Constrained head
option influences the water balance of the model: Internally,
hydraulic-head boundary conditions are set. The flow generated
by these additional boundary conditions is displayed in the Neumann
BCs section in the ![]() |
Optionally, the setting Prevent inflow can be enabled for both head limits. With this setting, the fixed-head boundary conditions act as Seepage Faces, i.e., FEFLOW applies a constraint condition that only allows outflow. This causes an iterative setting of the condition per time step. Only up to 30 iterations are performed. If there are still changes in the location of these additional head boundary conditions, FEFLOW proceeds to the next time step.
The checking for nodes with a water level higher than the top/lower than the model bottom is only done once per time step (option Prevent inflow disabled), setting the head boundary conditions for the next time step.
The Prevent inflow option is more accurate, but requires more computational effort as up to 30 iterations per time step might be performed.
Residual water depth for unconfined layers
Unit: [m]
Default value: 0.00101 m
The residual water depth is applied for three different tasks:
- Calculation of residual saturation for dry layers in phreatic mode.
- Minimum water level that is always kept at the model bottom if Head limits for unconfined layers are defined.
- Calculation of a minimum transmissivity in the bottom layer if Head limits for unconfined layers are defined for the model bottom and the bottom falls dry.