COMSOL 6.4 - Aquifer Water Table Calculation

This model demonstrates the application of COMSOL Multiphysics to a benchmark case of steady-state subsurface fluid flow and transient solute transport along a vertical cross section in an unconfined aquifer. Because of profound geologic heterogeneity, the model must estimate solute transport subject to highly irregular flow conditions with strong anisotropic dispersion. Van der Heijde (Ref. 1) classifies this case as “Level 2,” with enough potentially difficult parameter combinations to test a code’s ability to tackle realistic hydrological situations. Sudicky (Ref. 2) developed this problem to demonstrate a Laplace transform Galerkin code.

Model Definition

The hydrological setting for this problem is described in Figure 1, for groundwater flow at steady state. The aquifer is composed largely of fine-grained silty sand of hydraulic conductivity K1 = 5·10−4 cm/s with lenses of relatively course material of hydraulic conductivity K2 = 1·10−2 cm/s. Generally, groundwater moves from the upper surface of the saturated zone, the water table, to the outlet at x = 250 m. The water table is a free surface, that is, fluid pressure equals zero, across which there is vertical recharge, denoted by R, of 10 cm/a. The groundwater divide, a line of symmetry, occurs at x = 0. The base of the aquifer is impermeable.

Figure 2 shows conditions related to solute transport. The aquifer initially is pristine, and concentrations equal zero. For the first five years, a relative concentration of 1 is loaded over the interval 40 m < x < 80 m at the water table. The solute source is removed in year 5, and the concentration along this segment immediately drops to zero. The contaminant migrates within the aquifer via advection and dispersion. Throughout the domain, porosity, denoted by εp, is 0.35 as given by the benchmark, the longitudinal dispersivity, αL, and transverse vertical dispersivity, αT, are 0.5 m and 0.005 m, respectively. The effective molecular diffusion coefficient, Dm, is 1.34·10−5 cm2/s.

Figure 1: Settings for Darcy’s Law.

Figure 2: Definition of the transport problem.

Fluid Flow

Steady groundwater flow generally is expressed with a Darcy’s law for hydraulic head

(1)

where K (m/s) is the hydraulic conductivity, Qm (m3/s) is the rate of recharge to water table per unit volume of aquifer, and H is the hydraulic head (m). Darcy’s Law in COMSOL is formulated with pressure p as the dependent variable. Hydraulic head and pressure are related by

where D (m) is the elevation.

The equations for groundwater flow and solute transport are linked by the average linear velocity:

(2)

where the porosity εp or the fraction of the aquifer containing water accounts for the fact that only a portion of a given aquifer block is available for flow.

The boundary conditions for the groundwater flow problem are shown in Figure 1 and stated below. No flow condition and symmetry condition are applied at x = 0 and y = 0. Both are mathematically the same:

Hydraulic head is specified at x = 250 m. To model the recharge of 10 cm/a this velocity is specified at the upper boundary.

Solute Transport

Species transport typically is time dependent for geologic problems. The transport of a dissolved concentration c (mol/m3) is described with the advection-dispersion equation:

(3)

where DD (m2/s) is the dispersion tensor and De (m2/s) the effective diffusion; u is the velocity field; and t is time.

The dispersion tensor defines solute spreading by mechanical mixing and molecular diffusion. Typically, in porous media the spreading parallel to the flow (longitudinal dispersivity) exceeds the spreading perpendicular to the direction of flow (transverse dispersivity) by a multiple. This is an important aspect of the benchmark problem.

At the top boundary a space- and time-dependent function for the concentration is defined, according to

(4)

At the left boundary the concentration is set to zero and a no flux condition is applied to the right and bottom boundary.

Results and Discussion

Fluid Flow

Figure 3 provides hydraulic heads estimated with the COMSOL Multiphysics steady-state groundwater flow simulation. The hydraulic heads and streamlines correspond nicely to the benchmark results provided by Sudicky (Ref. 2). The water table geometry determined for the simulation nearly duplicates the benchmark geometry. The good match between the flow fields is expected because the initial water table geometry used with COMSOL Multiphysics was designed to closely resemble the benchmark geometry.

Figure 3: Estimated hydraulic head and flow lines.

Solute Transport

Solute transport solutions from the model are almost identical to the ones presented in Ref. 2. This is shown in the contour intervals for three times in Figure 4. In 1989, Sudicky concluded that the results illustrated in Ref. 2 are relatively free from numerical dispersion, as the low concentration contours closely follow the flow pattern. The surface plot for COMSOL also displays this property, in that even the lowest concentrations in Figure 4 still follow the irregular flow lines of Figure 3.

Figure 4: Solute concentration after 8, 12, and 20 years.

Figure 5: Solute concentration after 8, 12, and 20 years, displayed in one array.

To show how the contaminant spreads within the soil while the concentration maximum decreases with time, it is possible to plot all images of Figure 4 in one array. The result is displayed in Figure 5.

References

1. P.K.M. van der Heijde, “Model Testing: A Functionality Analysis, Performance Evaluation, and Applicability Assessment Protocol,” Groundwater Models for Resources Analysis and Management, A.I. El-Kadi (ed.), CRC Press, Lewis Publishers, Boca Raton, FL, pp. 39–58, 1995.

2. E.A.Sudicky, “The Laplace Transform Galerkin Technique: A Time-Continuous Finite Element Theory and Application to Mass Transport in Groundwater,” Water Resour. Res., vol. 25, no. 8, 1989.

Subsurface_Flow_Module/Solute_Transport/aquifer_water_table

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

In the tree, select > .

Click .

Click

In the tree, select > .

Click

Geometry 1

In the toolbar, click

Polygon 1 (pol1)

In the toolbar, click

In the window for , locate the section.

In the table, enter the following settings:


x (m)	y (m)
0	0
250	0
250	5.35
180	5.575
155	6.045
127	6.455
80	6.603
0	6.645

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 120.

In the text field, type 2.

Locate the section. In the text field, type 2.

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 70.

In the text field, type 2.

Locate the section. In the text field, type 180.

In the text field, type 2.

Click

The geometry is elongated. Define another which scales the view in the window for easier setup.

Definitions

In the window, under right-click and choose .

Axis

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

Click

Click the

button in the toolbar.

Darcy’s Law (dl)

Porous Matrix 1

In the window, under > > click .

In the window for , locate the section.

From the list, choose .

In the K text field, type 5e-4[cm/s].

The rectangles have a higher hydraulic conductivity.

Porous Medium 2

In the window, under > right-click and choose .

Select Domains 2 and 3 only.

Porous Matrix 1

In the window, expand the node, then click .

In the window for , locate the section.

In the K text field, type 1e-2[cm/s].

Next, set up the boundary conditions.

Inlet 1

In the toolbar, click

and choose .

Select Boundaries 7, 8, 10, 11, and 15 only. This corresponds to all upper boundaries.

In the window for , locate the section.

In the U0 text field, type 10[cm/a].

Hydraulic Head 1

In the toolbar, click

and choose .

Select Boundaries 16–18 only.

In the window for , locate the section.

In the H0 text field, type 5.3486.

Symmetry 1

In the toolbar, click

and choose .

Select Boundaries 1, 3, and 5 only.

Next, set up the transport equations. One important aspect of this benchmark problem is anisotropic dispersivity with a ratio of 100 between longitudinal and transverse dispersivity.

Transport of Diluted Species in Porous Media (tds)

Fluid 1

In the window, under > > click .

In the window for , locate the section.

From the u list, choose .

Locate the section. In the DF,c text field, type 1.34e-5[cm^2/s].

Porous Medium 1

In the window, click .

Dispersion 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the αL text field, type 0.5.

In the αT text field, type 0.005.

The next step is to set up the boundary conditions. At the top boundary the concentration is 0 except for the first five years within 40 m and 80 m. Define a function that will help you setting up this boundary condition.

Definitions

Rectangle 1 (rect1)

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type 40.

In the text field, type 80.

Click

This defines the local position of the concentration source.

Transport of Diluted Species in Porous Media (tds)

Concentration 1

In the toolbar, click

and choose .

Select Boundaries 7, 8, 10, 11, and 15 only.

In the window for , locate the section.

Select the checkbox.

In the c0,c text field, type rect1(x[1/m])*(t<=5[a]).

This expression calls the rectangle function with the argument x and gives the location of the species source. It is followed by a logic expression that applies the concentration for the first five years only. After that it remains 0.

Concentration 2

In the toolbar, click

and choose .

Select Boundaries 1, 3, and 5 only.

In the window for , locate the section.

Select the checkbox.

Outflow 1

In the toolbar, click

and choose .

Select Boundaries 2 and 16–18 only.

Add a material and fill out the missing values for the porosity and density.

Add Material

In the toolbar, click

to open the window.

Go to the window.

Material 1 (mat1)

In the window, under right-click and choose .

In the toolbar, click

to close the window.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	1000	kg/m³	Basic
Porosity	epsilon	0.35	1	Basic

Mesh 1

A proper mesh resolution ensures smooth and accurate results. Restrict the maximum element size to prevent too large mesh elements.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Size

In the window, under > click .

In the window for , locate the section.

Click the button.

Locate the section. In the text field, type 0.5.

Click

Study 1

Now, set up the solver sequence. First calculate a stationary Darcy velocity field. Use this for the second time-dependent step which computes the transport of the species over 20 years.

Step 1: Stationary

In the window, under click .

In the window for , locate the section.

In the column of the table, under , clear the checkbox for .

Step 2: Time Dependent

In the toolbar, click

In the window for , locate the section.

In the column of the table, under , clear the checkbox for .

Locate the section. From the list, choose .

In the text field, type range(0,1,20).

These are the output time steps. The solver will choose the computational time steps according to a convergence criteria. When using time dependent functions it is often recommended to restrict the computational time steps. Otherwise, if the convergence of the time-dependent solver is good, it can happen that the solver overestimates the required time step size and hence the user-defined time dependent function is not resolved properly.

Solution 1 (sol1)

In the toolbar, click

In the window, expand the node, then click .

In the window for , click to expand the section.

From the list, choose .

In the window, click .

In the window for , locate the section.

Clear the checkbox, to disable the automatic generation of default plots. Instead, choose the relevant plots from the Result Templates after computation.

In the toolbar, click