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.

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 Resources Research, 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

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 set-up.

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

Darcy’s Law (dl)

Fluid and Matrix Properties 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.

Fluid and Matrix Properties 2

Right-click and choose .

Select Domains 2 and 3 only.

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 check box.

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

This expression calls the rectangle function width 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