Lesson 2: Contaminant Sources

Sections in this lesson:

• 2.1. Lesson Introduction
• 2.2. Instantaneous Sources
  • 1D and 2D Models
  • Sensitivity Analysis
• 2.3. Continuous Sources
  • 1D Model and 2D Model
  • Sensitivity Analysis
• 2.4. Lesson Summary
• 2.5. Step-by-Step Tutorials
• 2.6. References

Acknowledgement: This lesson and the accompanying tutorials were originally written by Carmen Nale in 2004 as part of her Master's Project in Environmental Engineering.

2.1. Lesson Introduction

The main transport mechanism for aqueous chemical species is advection in the dominant flow direction -- the solute moves with groundwater seepage. If this were the only mechanism governing the contaminant's movement, it would move through the porous medium along streamlines. This does not occur because two other mechanisms, diffusion and dispersion, influence the solute's movement. Diffusion is a molecular mass transport process in which solutes travel from areas of higher concentration to areas of lower concentration. Its importance increases with the very low velocities found in clay and other nonporous media. Dispersion is a mixing process due to velocity variations and pathways within a porous medium that cause the front edge of a solute to spread out and lower in concentration. As the heterogeneities in a porous media increase, dispersion becomes more important. The incorporation of these processes into the two-dimensional (2D) advection-dispersion equation is shown below in the governing nonsteady-state two-dimensional contaminant transport equation (Bedient et al., 1999).

Where and with representing dispersivity [m] and D representing the dispersion coefficient (m²/day) in the x- and y-direction. Several different solutions can be derived from this conservation equation depending on the initial and boundary conditions (Bedient et al., 1999). For each type of source the following assumptions are made: the contaminant is conservative, the fluid is incompressible, the medium is homogeneous and isotropic, and saturated-flow conditions occur.

The types of contaminate sources discussed in this lesson are: a continuous source with constant input of contaminant and an instantaneous source, a pulse source or a spill. To simulate the transport of both types of plumes, a line of constant head on each side of the source was used to create a uniform flow field so advection is only in the x-direction. The change in head creates a uniform velocity gradient that can be used to simulate the transport of a plume. Figure 2.1 displays the breakthrough curves for a continuous versus an instantaneous source using the analytical solutions with the corresponding data shown in the tables below.

Figure 2.1. (A) 1D continuous source and (B) 1D instantaneous source

Figure 2.1(A) shows the typical breakthrough curve of a continuous contaminated source at x equals 60 m. Notice that as the contamination approaches a certain point, x equals 60 m, there is an initial rapid increase in concentration. The concentration then levels off when the initial concentration of the contaminant is reached at the point. This result is predicted by the 1D instantaneous source solution to the advection-dispersion equation. Figure 2.1(B) shows the breakthrough curve of an instantaneous contaminant source at x equals 60 m and 120 m, again with the models formulated to simulate 1D transport. As the pulse source moves through the subsurface, the concentration decreases due to the dilution effects caused by dispersion and diffusion. The area under the two breakthrough curves is equal as predicted by transport theory. When comparing Figure 2.1(A) to 2.1(B), the continuous source produces a higher concentration of contaminant within the aquifer versus the instantaneous source that is increasingly effected by dilution as it moves through the aquifer. The following sections will further discuss both types of sources in detail.

2.2. Instantaneous Sources

The concentration of an instantaneous source of contaminant is best described using a Gaussian distribution curve. This variation of the advection-dispersion equation assumes a homogeneous, isotropic, and saturated porous medium; steady state flow; and conditions where Darcy’s Law applies (Bedient et al., 1999). Equations 2.2 and 2.3, predict the change in a pulse source concentration with time in one or two dimensions, respectively. The variables used in Equation 2.2 and 2.3 are described in Table 2.1. These equations assume that the instantaneous source of contaminant is initially an infinitesimal line and point instantaneously distributed over the entire thickness of the aquifer into a one- and two-dimensional flow field, respectively.

1D and 2D Models

We can use IGW to visualize transport of an instantaneous source in a one- and two-dimensional flow field. Table 2.2 shows input parameters for the series of simulations described below. For comparison, the same parameters were used to derive analytical solutions to Equation 2.2 and 2.3 in an excel spreadsheet.

Table 2.2. 1D and 2D instantaneous source parameters used in the numerical and analytical solutions.
one-dimensional: Constants: Calculations: xo= 134 m *M = 5.0E+03 g/m2 x length = 4.7 m *Dx = 3.54 m2/d y length = 501 m Stability: B = 10 m x = 5 m Co = 1000 g/m3 t = 2.5 d K = 10 m/day Cr = 0.2   vx = 0.354 m/d Pe = 0.5   x = 10 m	two-dimensional: Constants: Calculations: xo= 230 m Dx = 3.46 m2/d yo= 400 m Dy = 3.46 m2/d Area = 2311 m2 Stability: B = 10 m x = 5 m Co = 1000 g/m3 t = 2.5 d K = 10 m/day Cr = 0.2   vx = 0.346 m/d Pe = 0.5   x = 10 m       y 10 m
Where: B = Aquifer thickness K = hydraulic conductivity *M = Co* x length = injected mass per unit cross-sectional area *Dx = x* vx

The basic setup of the model in IGW includes a parent zone with the constants described in Table 2.2, the creation of a velocity field, and then the addition of a zone with an initial instantaneous contaminant concentration of 1 ppm. Figure 2.2 shows the IGW simulations and the breakthrough curves for a one- and two-dimensional instantaneous source. The differences between the one- and two-dimensional analytical solutions are the added dispersivity constant in the y-direction for the two-dimensional solution and boundary conditions (i.e., line vs. point).

Figure 2.2. (A) IGW & graph for 1D Instant Source and (B) IGW & graph for 2D Instant Source

The figure shows the model comparing well with the analytical solution. The model curves for both graphs have a slightly lower maximum concentration than obtained using the analytical solution, which could result from the boundary conditions used in the model. As previously stated, the analytical solution assumes that an instantaneous infinitesimal line (1D) or point (2D) is added to the entire length of the aquifer. In the model, both are represented as having an area. This difference introduces some error.

Sensitivity Analysis

A sensitivity analysis on dispersivity helps evaluate the ability of the model to simulate sharp fronts by changing the dispersivity value by a factor of 10. Table 2.3 shows the model values used in the sensitivity analysis on dispersivity. Figure 2.3 shows the effects of the dispersivity () values of 10 m, 1 m, and 0.1 m, has on the maximum concentration determined from the numerical and analytical solutions.

Table 2.3. Dispersivity parameters used in the numerical and analytical solutions.

Constants: xo= 108 m yo= 442 m Area = 67 m2 B = 10 m Co = 1000 g/m3 vx = 0.091 m/d x = y = 10, 1, 0.1 m Calculations:     Dx = Dy = 0.91, 0.091, 0.0091 m2/d

Figure 2.3 shows that the best match between the model and the analytical solution is obtained when the highest dispersivity value (10 m) is used. Lowering the dispersivity value from 10 m (Figure 2.3(A)) to 1m (Figure 2.3B) decreases the ability of the model to accurately reproduce the concentration profile. Figure 2.3(C) shows the greatest difference between the model and analytical solution. It illustrates the difficulty of modeling a small dispersivity value with the spatial resolution of these simulation. It is difficult for the model to produce the sharp front necessary to simulate the analytical solution and a decreasing dispersivity can have a dramatic effect on the accuracy of the model compared to the analytical solution. Also notice that as the dispersivity value decreases, the maximum concentration of the contaminant increases.

Figure 2.3. Effects of the dispersivity values of (A) 10 m, (B) 1 m, and (C) 0.1 m.

An explanation for the dispersivity results is the basic stability of the numerical program during its computation of iterations. The stability criterion for an explicit numerical scheme with a mesh oriented in the direction of groundwater flow (i.e., x-direction) is determined using Equation 2.4. A value of the Peclet number (Pe) less than 2 ensures that oscillation will not occur at the front of the plume. When away from the zone of displacement of a sharp front, the value criteria for the Pe can be relaxed up to a value of 20 without major instabilities (de Marsily, 1986).

Pe = x / < 2

Equation 2.4. Peclet Number

Stability of explicit schemes also requires that the time step is appropriately set according so that the Courant number (Cr) equals about 1.

Equation 2.5. Courant Number

Where v is Darcy’s velocity and v_x is groundwater velocity, which is the velocity at the spatial node of interest in the model. The highest velocity value should be used to calculate the Cr. Table 2.4 provides the stability of each dispersion model shown in Figure 2.3.

Table 2.4. Stability criteria for dispersion simulations (v = 0.091m/day, n = 0.3)

Dispersion Model x t Cr Pe 10 m 8 m 10 days 0.1 0.8 1 m 8 m 10 days 0.1 8 0.1 8 m 10 days 0.1 84

The dispersion model with a dispersivity of 10 m is the only model that meets both the Cr and Pe criteria. The other two dispersion simulations have a large x and small dispersion values resulting in a high Pe. The Pe criteria could be met if the value of x was decreased to 2 m and 0.2 m for the 1 m and 0.1 m dispersion models, respectively. Interestingly, the model does not seem to be unstable when looking at Figure 2.3. This suggests, that the difference between the numerical solution and the analytical solution could be from numerical dispersion occurring in the model (Bedient et. al, 1999).

2.3. Continuous Sources

Equation 2.6 and 2.7 display the analytical solutions to the advection-dispersion equation for a transient one-dimensional and a steady-state two-dimensional continuous source model. Equation 2.6 was derived from Equation 2.1 with the following conditions: initial conditions are set equal to 0, i.e., C(x, t=0)=0; boundary conditions are set to the length of plume traveled, i.e., x; and continuous source load at x=0 and x=L are set at zero, i.e., C(x=0, t)=C_o for t>0 and C(x,)=0 for t>0, respectively (Bedient et al., 1999). Equation 2.7 is a two-dimensional steady-state solution with groundwater flowing in the x-direction. Notice that the exponential term decreases the concentration of the contaminant as the plume moves away from the source and disperses in the y-direction.

	Equation 2.6. 1D continous source step function
	Equation 2.7. 2D steady state continuous source
Table 2.5. Variables definitions for Equations 2.6 and 2.7.
Variable Description Units B Aquifer thickness m C Concentration of contaminant g/m3 Co Initial concentration of contaminant g/m3 Dx , Dy dispersion coefficient in x, y direction, respectively m2/day m Mass flux of contaminant g/day n Porosity of the medium - vx Groundwater velocity m/day x, y Distance in x, y direction, respectively m xo , yo Initial position of the plume in the x, y direction, respectively m

1D Model and 2D Model

IGW can be used to model the advection-dispersion equation for a transient one-dimensional and a steady-state two-dimensional continuous source. Table 2.5 and 2.6 describe additional variables and parameter values used in these analytical and numerical solutions.

Table 2.6. 1D and 2D steady state continuous source parameters used in the numerical and analytical solutions.
One Dimension Constants: xo = 241 m x' = 302 m Co = 1000 g/m3 vx = 0.34 m/day x = 10 m K = 10 m/day Calculations: Dx = 3.4 m2/day Stability: Dx= 9.2 m Dt = 5 d Cr = 0.2   Pe = 0.9   xo equals plume location (right side of area) and x' equals model distance value	Two Dimensions, Steady State Constants: xo = 113 m yo = 454 m Area = 67 m y plume = 8.2 m x plume = 8.2 m K = 10 m/day B = 10 m Co = 2400 g/m3 Cmodel = 1000 g/m3 vx = 0.34 m/day n = 0.3   y = 10 m Calculations: Dy = 3.4 m2/day m = 19861 g/day Stability: x = 8.4 m t = 10 day Cr = 0.4   Pe = 0.9

The basic setup of the one- and two-dimensional models includes a parent zone with the constants described in Table 2.6, the creation of a linear velocity field in the x-direction, and then the addition of a zone with the continuous contaminant concentration. The two-dimensional model includes the addition of the dispersivity constant in the y-direction and the requirement for the plume to reach steady state and stabilize (i.e., long simulation time). The breakthrough curve and IGW simulation for a 1D continuous source is presented in Figure 2.4 (A). As the contaminant concentration from the continuous source spreads through the subsurface heading towards a point x, the concentration initially increases rapidly then levels off at the concentration of the contaminant source.

Figure 2.4. (A) IGW image & graph for 1D breakthrough curve and (B) IGW image & graph for SS 2D concentration curve

Figure 2.4 (B) illustrates the IGW simulation and change in concentration with distance from the contaminant source for a two-dimensional, steady-state continuous source. When at the source, the analytical solution determines the concentration to be infinite, and then rapidly decreases and levels off to some value, a characteristic of the exponential function. This decrease in concentration is a result of transverse dispersion. When comparing the analytical solution to the numerical solution for both Figures 2.4 (A) and (B), they are in good agreement. This is due to the completion of an extensive sensitivity analysis (discussed next) in order to use the best area, grid size, and dispersion coefficient within the model to best simulate the analytical solution.

Sensitivity Analysis

Grid Resolution

A sensitivity analysis for different grid resolutions helps determine when good matches between the model and the analytical solution for a 1D continuous source are achieved. The results are presented in Figure 2.5. The stability criterion is then examined for each grid resolution in Table 2.7.

Table 2.7. Stability parameters for grid resolution (v = 0.34 m/day, y = 10 m, n = 0.3)

Grid Model (NX) x t Cr Pe 40 26 m 5 days 0.1 3 90 11 m 5 days 0.2 1 110 9 m 5 days 0.2 0.9

The input data needed for IGW to calculate x and y automatically is the NX parameter, which is used to represent each grid model shown in the first column in Table 2.7. Grid model NX = 40 is the only model where the Pe predicts instability, which is due the use of a coarse grid. The influence of using different grid resolutions is apparent in Figure 2.5.

Figure 2.5. (A) Breakthrough curve grid sensitivity and (B) C/Co vs. x-grid sensitivity.

As the grid resolution increases and the x decreases, a better match to the analytical solution is achieved. In Figure 2.5 (B), the grid model 40 seems to be showing a slight instability, since the line is not smooth. This instability may be due to the fact that grid model 40 has a Pe that predicts instability. There is a significant improvement in model accuracy as the grid resolution increases from 40 to 90, as indicated by both figures. However, there seems to be a point in diminishing returns when comparing grid model 90 to 110. When the model’s grid resolution increases, so does the accuracy of the model and time to run the model simulation. The ability to use a very fine grid resolution for simulations is limited by the hardware’s processor and memory capabilities.

Area Constraints

The ability of the model to simulate an instantaneous and continuous source depends on defining the source boundary condition. The analytical solution for the instantaneous and continuous source assumes that the contaminant source is an infinitesimal point while the model requires a grid area to be specified. Therefore, a sensitivity analysis on area was performed to determine an area for the model simulations that represent the analytical solution well. The steady-state two-dimensional continuous source model was the most sensitive to area, therefore, the model was used for the sensitivity analysis on area. The values used for the models and analytical solution (Equation 2.7) are shown in Table 2.6, with the following changes shown in Table 2.8.

Table 2.8. Parameters used in the sensitivity analysis on area for the numerical and analytical solutions

Area Model   2616 m2 1116 m2 67 m2 Constants: xo (m) = 138 128 113 yo (m) = 440 436 454 Area (m2) = 2616 1116 67 y plume (m) = 51 33 8.2 x plume (m) = 51 33 8.2 Co (g/m3) = 1700 1700 2400 Calculations: m (g/day) = 87906 57416 19861 Stability: x (m)= 8.4 8.4 8.4 t (d)= 10 10 10 Cr = 0.4 0.4 0.4 Pe = 0.9 0.9 0.9 B=10 m, Cmodel=1000 g/m3, n=0.3, vx=0.34 m/d, y=10 m, Dy=3.4 m2/d

Figure 2.6 shows the influence of the contaminant area, or boundary condition, on the ability of the model to simulate the concentration profiles as predicted by the analytical solution. The figure shows the log maximum contaminant concentration that occurs along the centerline (i.e., y equals 0) with distance at steady-state conditions. In Figure 2.6(A), an area of 2616 m² causes the model to deviate significantly from the analytical solution close to the source. As the area of the continuous contaminated source was minimized, the model was better able to simulate the analytical solutions. Figure 2.6(C) shows the best fit between the numerical and analytical solution. The figure illustrates the effect of the boundary condition has on the models ability to accurately simulate the analytical solution.

Figure 2.6. These graphs show the log centerline concentration of the plume with distance in the x-direction for a 2D steady-state continuous source plume. Each figure compares the analytical solution with the numerical model results.

Figure 2.7 shows the relative size of the contaminant source areas and the resulting steady-state plumes. As a result of both advection and transverse dispersion, the area of the source decreases along with the resulting size of the plume. Thus, the smaller size source yields a more accurate simulation along the centerline of the plume and in the concentration predicted in the transverse direction as well.

Figure 2.7. These images show the relative size of the area of the continuous source (red box) and the resulting 2D steady-state plume that develops for each source area given. (A) Model plane view of 2616 m2, model plane view of 1116 m2, model plane view of 67 m2.

Dispersion Effects

The two-dimensional steady-state continuous source model presented previously was used to conduct a sensitivity analysis of dispersivity, see Table 2.9. The illustration of the model and analytical solutions’ sensitivity to dispersivity is shown in Figure 2.8. Based on the earlier area analysis, a value of 67 m² was used for the initial contaminant area in order to get the best results for the dispersion simulations.

Figure 2.8. Sensitivity analysis on dispersivity.

The figure shows the general trend of contaminant concentration decreasing as dispersivity increases, which is caused by an increase mixing and dilution. As observed previously, the initial boundary conditions of the plume has a significant impact on the models ability to match the analytical solution close to the source. Further away from source, the model fits the analytical solution very well for each dispersivity simulation. When comparing the different dispersivities in Figure 2.8, there is a noticeable increase in variation from the solution for the 0.1 m dispersivity model compared to the other simulations. This is explained by viewing the stability criteria for each dispersivity value shown in Table 2.9.

Table 2.9. Parameters used in the sensitivity analysis on dispersivity for the numerical & analytical solutions.
Dispersivity Model   0.1 m 1 m 10 m Calculations: Dy= (m2/day) = 0.02 0.2 2 * m (g/day) (input value) = 4800 6800 12500 Stability: x (m)= 4 4 8 t (d)= 5 5 10 Cr = 0.2 0.2 0.2 Pe = 40 4 0.8	Constants: xo (m) = 113 yo (m) = 454 Area (m2) = 67 y plume (m) = 8.2 x plume (m) = 8.2 B (m) = 5 Co (g/m3) = 2037 vx(m/d) = 0.2 n = 0.3 y (m) = 10

The high Pe in the model causes a numerical instability resulting in deviations. This is evident in the figure due to the difficulty of the model to accurately simulate near source behavior.

Velocity Analysis

A sensitivity analysis of velocity for a two-dimensional steady-state continuous source is illustrated in Figure 2.9 with the models input data provided for in Table 2.10.

Figure 2.9. Sensitivity analysis on velocity.

Since a constant mass of contaminant introduced into the aquifer per day could not be input into IGW model, a constant concentration of 2037 g/m³ was used. This explains the visually undesirable graph given in Figure 2.9, but the graph does confirm the results of the velocity analysis.

Table 2.10. Parameters used in the sensitivity analysis on velocity for the numerical and analytical solutions.

Velocity Model   0.02 m 0.2 m 2 m Calculations: y= (m) 1 1 1 Dy= (m2/day) = 0.02 0.2 2 * m (g/day) (input value) = 670 6700 67000 Stability: x (m)= 4 4 4 t (d)= 20 5 0.625 Cr = 0.1 0.3 0.3 Pe = 4 4 4 n = 0.3

As groundwater velocity increases, the concentration of the contaminant decreases. This is not apparent in Figure 2.9, but can be explained by noticing the value of the contaminant mass rate and corresponding velocity value shown in Table 2.10. Notice that as velocity increases by a factor of 10, the mass rate of contaminant added to the aquifer also increases by a factor of 10. Since a single line is shown in Figure 2.9, it can be concluded that as the velocity increases, more mass per day has to be added to the aquifer to keep the same concentration profile. As the velocity increases, the effects of dilution also increases.

When comparing the analytical solution to the model results, the initial boundary conditions still have a significant effect on model accuracy. Near the boundary, the model deviates from the analytical solution; the model fits better when x is greater than 40 m from the source.

2.4. Lesson Summary

This lesson compared analytical solutions to IGW model simulations for 1D and 2D contaminant transport equations. It provides students with a way to visualize groundwater flow and contaminant transport within an aquifer system and reinforces important classroom concepts.

The comparisons between the analytical solutions and the model simulations along with sensitivity analyses, also allow students develop a better understanding the capabilities and limitations of the IGW model. Below are some key points that will be applied in Lesson 3, a case study at the Superfund Site at St. Joseph, MI, USA:

1. The analytical contaminant transport solutions have assumptions that must be incorporated into the model to accurately simulate the analytical solutions. This was evident for the boundary condition used in the IGW model (i.e., initial area of source) of both types of contaminant sources. The analytical solutions assume that the contaminant source is an infinitesimal point instantaneously distributed over the depth of the aquifer, meanwhile, IGW models it as an area. This causes the maximum concentration given in IGW to be less than predicted by the analytical solutions near the source. The result proves that the boundary conditions have a significant effect on model results. Therefore, caution should be utilized when using output data near a boundary condition.

2. The sensitivity analysis on dispersivity for an instantaneous source illustrates that the model is unable to accurately model sharp fronts (i.e., a very low dispersivity value) when using a large time step and grid size. In order for the simulation to produce the best results, the model must be given an appropriate x and t value, determined by the Pe and Cr, required for the stability of the model. However, the sensitivity analysis also showed that a slightly larger Pe and/or Cr value (one that should predict instability) may still provide reasonable accuracy. The sensitivity analysis also shows that the grid analysis is closely related to the stability parameters. As the grid resolution increases, ability of the model to simulate the analytical solutions describing groundwater flow and contaminant transport also increases.

3. The sensitivity analysis on the variation of dispersivity and velocity for a continuous source demonstrates other model limitations. When the dispersivity value increases, the concentration of the contaminant decreases due to dilution effects. When velocity decreases there is an increase in the concentration profile; this is also due to the dilution effect.

4. In this lesson, the model was unable to accurately simulate the analytical solutions when the dispersivity value was low and near boundary conditions. Reducing the grid size and time step to achieve a Pe less than 2 and reducing the initial area of the contaminant was required to improve the accuracy of the model.

Despite these limitations, the model does a good job of modeling a groundwater divide (Lesson 1) and the instantaneous and continuous source analytical solutions (Lesson 2). The skills learned in these first two lessons will be applied in Lesson 3.

2.5. Step-by-Step Tutorials

Click below to link to pdf documents with step-by-step instructions to complete these simulations in IGW. You can also download the IGW file that will result from successful completion of the tutorial. The IGW file is compressed as a ZIP file. To open it, you will need to save the ZIP file to your local computer and uncompress it with a program such as Winzip or Windows Explorer.

• Tutorial 2A: 1D Instantaneous Source
  • 1D_instanteous_source.zip
• Tutorial 2B: 2D Instantaneous Source
  • 2D_instantaneous_source.zip
• Tutorial 2C: 1D Continuous Source
  • 1D_continuous_source.zip
• Tutorial 2D: 2D Continuous Source
  • 2D_continuous_source.zip

Please note that the lessons and tutorials on our website were designed for IGW version 3.5.6 -- we recommend using this version of the software to ensure compatability between the step-by-step instructions and what you see on your computer screen.

2.6. References

Bear, Jacob and Arnold Verruijt, 1987, Modeling Groundwater Flow and Pollution, D. Reidel Publishing Company: Boston, p. 319-326.

Bedient, Philip B and Wayne C. Huber, 1992, Hydrology and Floodplain Analysis, 2nd ed, Addison-Wesley Publishing Company: Reading, Massachusetts, p. 7-37.

de Marsily, Ghislain, 1986, Quantitative Hydrogeology for Engineers, Academic Press: San Diego, pp. 394-400.

Back to top