 

THESIS
I

 

URnAnv

Michigan State
University

 

 

 

This is to certify that the
thesis entitled

HIERARCHICAL GROUNDWATER MODELING AT THE
ZEPHYR SITE, MUSKEGON, MICHIGAN: EVALUATING THE

EFFECTIVENESS OF THE REMEDIATION CAPTURE
SYSTEM

presented by

Yang Li

has been accepted towards fulﬁllment
of the requirements for the

MS. degree in Civil Engineering

 

 

“gx LM

Major;Q Profe or’s Signature
(9/ /C><) C 07

//Date

 

MSU is an Alﬁnnative Action/Equal Opportunity Employer

A-.-o—------n-o-D-O-O-c-I-I--o-a-n-121-tgLQo-l-‘_n-|-l-O-t-l-I‘D-0-0-o-o-I— -

-<—-—A-l-‘-l-l-l-l-l

 

PLACE IN RETURN BOX to remove this checkout from your record.
To AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5/08 K./Proj/Acc&Pres/ClRC/DateDue.indd

 

HIERARCHICAL GROUNDWATER MODELING AT THE ZEPHYR SITE,
MUSKEGON, MICHIGAN: EVALUATING THE EFFECTIVENESS OF THE
REMEDIATION CAPTURE SYSTEM

By

Yang Li

A THESIS

Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of

MASTER OF SCIENCE
Civil Engineering

2009

UMI Number: 1485580

All rights reserved

INFORMATION TO ALL USERS
The quality of this reproduction is dependent upon the quality of the copy submitted.

In the unlikely event that the author did not send a complete manuscript
and there are missing pages, these will be noted. Also, if material had to be removed,
a note will indicate the deletion.

UMI

A
Dissertation Publishing
UMI 1485580
Copyright 2010 by ProQuest LLC.
All rights reserved. This edition of the work is protected against
unauthorized copying under Title 17, United States Code.

Pro uest'

 

ProQuest LLC
789 East Eisenhower Parkway
PO. Box 1346
Ann Arbor, MI 48106-1346

ABSTRACT

HIERARCHICAL GROUNDWATER MODELING AT THE ZEPHYR SITE,
MUSKEGON, MICHIGAN: EVALUATING THE EFFECTIVENESS OF THE
REMEDIATION CAPTURE SYSTEM

By
Yang Li
In this thesis, a hierarchical groundwater modeling system integrating the 10m Digital
Elevation Model and site-speciﬁc hydrogeologic data was developed for the Zephyr site,
Muskegon, Michigan. The site is currently under active remediation using pump and
treats through a network of 60 low capacity wells and 10 injection trenches. The
hierarchical system was calibrated to steady state water level data from 128 on-site

monitoring wells with and without remediation stresses.

Through a series of nested “patch dynamic” models, the hierarchical system enables
accurate simulation of detailed, complex well dynamics and capture zone distributions
associated with very low capacity wells. Systematic model simulations, including
forward and reverse particle tracking, as well as, integrated water budget analyses, were

performed to evaluate the effectiveness of the remediation.

The results show that current remediation system is unable to capture contaminated

groundwater and signiﬁcant gaps exist between purge wells and offsite migration exists.

DEDICATION

TO MY PARENTS

FOR THEIR SUPPORT AND UNDERSTANDING

iii

ACKNOWLEDGEMENTS

I grateﬁilly acknowledge the contribution of those who have made this thesis possible.
My deepest gratitude goes to Dr. Shu-guang Li, Professor, Department of Civil and
Environmental Engineering, Michigan State University, for his advice, patience, and

encouragement.

I also extend my appreciation to Dr. Hua-sheng Liao for his restless help during my

study.

iv

TABLE OF CONTENTS

LIST OF TABLE ................................................................................................................ vi
LIST OF FIGURES .......................................................................................................... vii
ABSTRACT ........................................................................................................................ 1
INTRODUCTION .............................................................................................................. 2
Local Grid Reﬁnement. ................................................................................................... 6
Local Analytical Correction ............................................................................................ 7
Local Numerical Correction ............................................................................................ 8
Hierarchical Groundwater Modeling ............................................................................... 9
ZEPHYR GROUNDWATER MODELING SYSTEM .................................................... 12
Conceptual representation ............................................................................................. 12
Hierarchical Flow Modeling ......................................................................................... 21
Model Calibration .......................................................................................................... 22
Model Sensitivity Analysis ............................................................................................ 28
ANALYSIS OF REMEDIATION SYSTEM PERFORMANCE ..................................... 31
Hierarchical Capture Zone Analysis .............................................................................. 32
Hierarchical Forward Particle Tracking ........................................................................ 35
Water Balance Analysis ................................................................................................. 40
CONCLUSION ................................................................................................................. 42
APEENDIX A ................................................................................................................... 43
REFERENCE .................................................................................................................... 45

LIST OF TABLE

Table 1 Flux Values for All Injection ................................................................................ 18
Table 2 Flux Values for all Extraction .............................................................................. 19
Table 3 Values of Parameters selected as model input .................................................. 25
Table 4 Values of Parameters selected as model input for extraction trench system ........ 25
Table 5 Results for Sensitivity Analysis ........................................................................... 29

vi

LIST OF FIGURES

Figure 1 Site Map ................................................................................................................ 3
Figure 2 Purge Wells, and Injection! Extraction Trench System ........................................ 4
Figure 3 Sampling Network on site and in the vicinity .................................................... 14

Figure 4 Conceptual Representation of key Sources, Sinks and Remediation Components

........................................................................................................................................... 15
Figure 5. Conceptual Cross-section for Zephyr Site ......................................................... 16
Figure 6 Divisions of Recharge Zone ............................................................................... 20
Figure 7 Divisions of Conductivity Zone ......................................................................... 20

Figure 8 Determination of Site Scale Model based on Simulation Results from Regional

Scale Model ...................................................................................................................... 23
Figure 9 Divisions of Recharge Zones for Site Scale Size Model .................................... 24
Figure 10 Divisions of Conductivity Zones for Site Scale Size Model ............................ 24

Figure 11 Comparison of observed and simulated heads at steady state for 176 monitoring
wells in the site scale size model ...................................................................................... 26

Figure 12 Comparison of observed and simulated ﬂux at steady state for extraction

trenches ............................................................................................................................. 27
Figure 13 Head Distribution with remediation system on under steady state .................. 30
Figure 14 Head distribution with remediation system off under steady state ................... 31
Figure 15 Hierarchical Reverse Particle Traching tracking .............................................. 34

Figure 16 hierarchical “family tree” to illustrate the relationships between the patch
models. .............................................................................................................................. 35

vii

Figure 17 2008 BTEX Analysis Result for Sampling ....................................................... 38

Figure 18 Hierarchical Forward Particle Tracking ........................................................... 39

Figure 19 Water Budget for Site Area ............................................................................... 41

viii

LIST OF ABBREVIATIONS

GPM: Gallon per Minute;

BETX: Benzene, Toluene, Ethylbenzene, Xylene;

ix

Hierarchical Groundwater Modeling of a Remediation Capture System in

Muskegon, Michigan
Hua-Sheng Liao, Yang Li, Shu-Guang Li
Department of Civil & Environmental Engineering
Michigan State University, East Lansing, MI 48824

S. Mathuram, Wilhelmus VanLeeuwen, Rimple

ABSTRACT

In this paper, we demonstrate a generalized hierarchical approach for modeling complex
groundwater systems. In particular, we demonstrate how the “hierarchical patch
dynamics modeling” approach developed by Li and his co-workers [Li et al., 2005, 2006,
2007; 2009; Afshari et al., 2009] enable simulating, in a highly ﬂexible and efﬁcient
manner, a complex plume capture system at one of the largest groundwater remediation
operations in Michigan. The groundwater flow system at the site, because of the complex
interaction between ambient hydrologic stresses and on-site remediation operations,
exhibits a unique multiscale pattern that proves to be difﬁcult to simulate using standard
modeling tools. The hierarchical modeling system was calibrated to water level
measurements collected from more than 120 monitoring wells located both on site and in '
the immediate proximity. Systematic hierarchical simulations, including forward and
reverse particle tracking, as well as, integrated water budget analyses, were performed to
evaluate the effectiveness of the on-going remediation. The results provided critical

information that is being used to improve the design of the pump and treat system.

INTRODUCTION

The site under study is located in Muskegon Township, a historically industrial area of
approximately 150 acres surrounded by wetlands and surface water bodies, including
Muskegon Lake, Bear Lake, Muskegon River, and Bear Creek (Figure 1). Decades of
industrial activities in the area including petroleum reﬁnery operations resulted in
signiﬁcant soil and groundwater contamination characterized by high concentrations in
BETXs, chlorinated organic compounds, semi-volatile organic compounds, and free
phases. Monitoring suggested that contamination had migrated offsite.

Since 2000, the State of Michigan has been operating and maintaining a pump and treat
remediation system on the site, with a goal to hydraulically contain the large plumes,
minimize offsite impact, recover free products, and ultimately cleanup the contamination.
The remediation system is composed of a large network of 59 low capacity extraction
wells (1 to 10 GPM), an extraction trench system, and 10 injection trenches, two 60,000
gallon above-ground bioreactors, and one 60,000 gallon clariﬁer tank (Figure 2). The
extracted groundwater is treated on-site in the bioreactors and clariﬁer tank prior to
re-injecting the treated groundwater into the ground through inﬁltration trenches. The

network was designed to operate at a total ﬂow rate of approximately 250 GPM.

as: an _ 2&5

 

 

 

 

 

.. .83st2 .._\.H\

 

      

 

 

Ly

  

 

 

m2zanHHHHIlllll
v; 5.0 o

2.233 x...

e853. llr

 

gem

$3M

 

 

vcowoq

 

 

82 ea _ semi

 

 

 

 

 

 

Sim
530352 \,

   

 

   

Il\tiioo \\l.\\\k .

snag:
8:5
834
cowoxmsz
8:2 ”I

 

255% lit
33.

$3M

 

 

vcowoq

 

 

 

883m 55C. :ouombxm E2803 98 n2:“. >9 ems; N 98mm.

////
:2 E4 8892 \ \x
2‘: 1.22292 e/K
rhlPlbiei 20.3. e , / e e e \\
. \

l 1 3.
23:2 arm: v 8% _ / ti .3 e . 1
0mm ir—tH V NEH \\ i//_ 9 \

.
m E / _,_ V \\
Ilo E N E ,, ,,,../I\w \\ \.
. ///% \ mm \\
m E E E . o x \ . \.
ee .\ , e e \.

 

 

 

 

 

   

 

 

._ a
e
/b ﬂ \il ll/ 99/ \\ / $ 9 2 ﬂ
,Illii / a \\ e E ., a / 3:2 U
\\ , a / S o 36 o
..\ / 9 e .\\,

«was. a e a «aw/ix .\ / 5%.; 883E 302 E2

4 I. ii I \v\. / / account aouoabxm BO Fm

I. , mm”. c, / some“... cocoon: w:

\w/ Nil... \x: e / /
\ //P/. x Y e / / :03 85 e
// ii, \wmet I/ / \K._ / UCQWDJ
\ \ Ar K \y/\

 

 

Systematic sampling over years shows that concentrations of key contaminants in

many monitoring wells, especially those in the source areas, have signiﬁcantly declined.
Many questions, however, still remain:
o How are the site wide plumes responding to the complex pumping and injection
stresses?
0 Is the remediation system containing the contamination?
0 Can the remediation system be optimized to signiﬁcantly improve cleanup

efficiency and reduce operational costs?

One earlier investigation tackled these issues through a systematic groundwater modeling
study, but the effort was not particularly effective. One major limitation of the earlier
effort lies in that the model developed was unable to capture critical details and
multiscale interactions essential to remediation system evaluation and optimization.

The groundwater ﬂow in the area is characterized by complex interplay of signiﬁcant
variability across disparate length scales. These include variations at “well scale”, “site
scale”, and “regional scale”. At the well scale, pumping by extraction wells across the site
creates a network of small drawdown cones, each with a characteristic length scale of 1 to
Sit. At the site scale, the combined effects of aggregated pumping, injection, natural
recharge, and drainage to extraction trenches, surface seeps, local wetlands, and surface

water bodies create a complex, large groundwater mound spanning across thousands of

feet. Regionally, groundwater ﬂow is controlled by Muskegon Lake and Muskegon River,

the global sink in the Muskegon watershed. These different scales of variations must all

be properly taken into account in order to simulate properly the groundwater ﬂow system

at the site and to quantify and optimize remediation system performance.

Several methods can in principle be used to model the multiscale groundwater ﬂow
system, but they all have limitations. The most commonly used methods include: 1) local
grid reﬁnement (e.g., Fung, L.S.-I(, 1992; Gable et al., 1996; Heinemann, 1983), 2) local
nested analytical correction (Prickett, 1967; Peaceman, 1978; Pritchett and Garg, 1980),
and 3) local nested numerical correction (e.g., Ward et al., 1987; Efendiev et al., 2000;
Mehl and Hill, 2002). In the following paragraphs, we brieﬂy review the advantages and

disadvantages of these methods.

Local Grid Reﬁnement.

Local grid reﬁnement represents the most commonly used method to predict detailed well
dynamics in a numerical model. Relatively-large-size cells in a numerical model grid are
subdivided or reﬁned with progressively smaller spatial dimensions in the areas of
interest (e.g., Fung, L.S.-K, 1992; Gable et al., 1996; Heinemann, 1983). This results in a
more accurate estimation of hydraulic head or drawdown at the well scale. In other words,
there will be less approximation error involved with averaging the predicted head for that

particular well node. In using local reﬁnement for simple or small-scale problems, the

solution is obtained quickly, and the consistency between the regional and local area

 

around the wells is maintained. However, for large-scale, regional groundwater models

where there may be a signiﬁcant increase in the number of nodes (e.g., millions rather
than thousands), the cost of computation increases exponentially and the process can

become very expensive in terms of time and computer resources required.

Local Analytical Correction

An alternative approach to represent detailed well dynamics in a regional model is to use
a local analytical correction within the cell containing the well based on the steady-state
Thiem equation (Thiem, 1906; Anderson and Woessner, 1992). Trescott et a1. (1976)
systematically evaluated the performance of local analytical correction and showed that
corrected drawdown is approximately applicable if the following assumptions are
satisﬁed: (1) flow to the well is within a square ﬁnite-difference cell and can be described
by a steady-state equation with no source term except for the well discharge, (2) the
aquifer is isotropic and homogenous in the well cell, (3) only one well, located at the cell
center, is in the well cell, and (4) the well fully penetrates the aquifer. Prichett and Garg
(1980) provided formulas for analytical correction for model grid cell geometries other
than square. Our own empirical experience shows that local analytical correction can lead
to problematic flow pattern around wells although the corrected drawdown is reasonably

accurate .

Local Numerical Correction

A more general approach for modeling detailed local well dynamics is through “local
numerical correction” or nested grid modeling. In this approach, grid-dependent
information from the regional model is used to construct a separate model with ﬁner grid
spacing to obtain more information around the area of interest (Ward et al., 1987;
Efendiev et al., 2000; Mehl and Hill, 2002). The model with the ﬁner grid is called a
“submodel”, “patch model”, or “local model” which derives its initial and boundary
conditions from the parent model (Townley and Wilson 1980; Ward et a1. 1987; Buxton
and Reilly 1986). The boundary conditions can be either interpolated heads or ﬂuxes.
Once the local model is created, it performs as an independent model. A nested grid
approach avoids the potential difﬁculties encountered with local grid reﬁnement by
converting the original problem of solving a very large matrix system into one that solves

multiple, much smaller matrix systems.

A major drawback in implementing the nested model approach is that the interaction
between the parent and local models (of which there can be many) depends on the offline
analysis and processing of model modiﬁcations or simulation results from the parent
model to obtain the boundary and starting conditions for the local model. For example,
once a new simulation is completed using the regional model, new boundary conditions
and starting ﬂow and solute transport conditions, if applicable, are determined for each

local model at the local model grid spacing. Making modiﬁcations to models or

processing simulation results for use in different scales of models can be very time

consuming; especially when the problem is a transient-ﬂow or coupled with solute
transport simulation, or a feedback loop is needed to account for potential signiﬁcant
two-way interaction between parent-nested submodels. The effort involved may become
impractical when the ofﬂine conceptual changes must be made iteratively or in more than
one model. Because of this, applications of the nested grid approach are limited, in most
cases, to a very small number of submodels (e.g., l or 2), and are implemented with little
ﬂexibility. This practical implementation difficulty has severely limited our ability to take
full advantage of the nested grid approach for solving complex groundwater problems,

especially those that span across multiple spatial scales.

Hierarchical Groundwater Modeling

Li et al. (2006) recently developed a new method of implementing nested grid modeling.
In particular, Li et al. (2006) took advantage of two emerging computing strategies —
“dynamic fusion” and “interactive steering” to develop a dynamically-integrated,
human-centered “hierarchical patch dynamics paradigm (HPDP)” for generalized
hierarchical groundwater modeling. The HPDP was originally developed by Wu and his
co-workers to address complexities and scale interactions in the context of landscape
ecology and ecological modeling (Wu and Loucks 1995, Wu 1999, Wu and David 2002).
The two modern computing strategies make it possible to employ arbitrary numbers of

nested submodels (patch models) for solving much more complex groundwater

problems. Dynamic fusion enables realtime integration of modeling, subscale modeling,

and visualization, eliminating a major bottleneck in nested grid modeling. In particular, it
allows dynamically integrate all patch models at the end of each time-step during the
simulation so that results from the parent-models (e.g., boundary-conditions), and any
changes to the parent-models propagate dynamically and automatically to all patch

models, without the need for ofﬂine post-processing.

The basic concept behind this strategy is straightforward. Instead of treating the
simulation of groundwater ﬂow and solute transport separately, Li et al. (2006) model
both processes concurrently (sequential within a time step). Instead of treating the various
scales of patch dynamic modeling (regional, subregional, local, site, hot spots) as
different phases in a long sequential batch-process, the multi-scaled processes are
coupled and modeled simultaneously. Instead of relegating the presentation and analysis
of graphical results for each patch of interest to the “post-processing” phase, at the end of
a time-consuming sequence of many disjointed steps, all results are presented at the end
of each simulation step to permit the interpretation of results as soon as they become

available.

Computational steering enables a new way of model-user interaction that is particularly
beneﬁcial for hierarchical modeling. It provides a mechanism to incorporate human

intelligence in the incremental simulation process. It allows a modeler to interactively

10

steer the hierarchical computation, to control the program execution-sequence, to guide

the evolution of the aquifer dynamics, to control the visual data representation during
processing, and to dynamically modify the hierarchical computational process during its
execution. At anytime, the modeler can interrupt the program to assess the solution, and
based on this assessment, to insert patches in selected areas and layers of interest. As the
modeling-hierarchies are dynamically-coupled, the modeler can work with large numbers
of modeling patches, as if they are working with a single modeling system. The modeler
can perform sensitivity analysis of the integrated system, making changes to patch
boundaries, resolutions, conceptual representations, solvers, and iteratively assessing the

inﬂuence of these changes on the ﬁnal solutions.

11

ZEPHYR GROUNDWATER MODELING SYSTEM

We apply in this paper the dynamically integrated HPDP to model the multiscale

groundwater ﬂow system at the Zephyr site. Images in this thesis are presented in color.

Conceptual representation

Our modeling effort was systematically integrated with ﬁeld investigations. Data were
collected on-site to speciﬁcally support the modeling study. Additional data from the
Michigan statewide geospatial databases and the statewide groundwater database were
used to construct the modeling system. Figure 3 shows the sampling network that consists
of 128 monitoring wells onsite and in the vicinity. Figure 3 shows conceptual aquifer
cross section for Zephyr Site. Figure 4 shows a conceptual site representation created
using the data available, depicting key sources, and sinks. Details of remediation
components are shown by Figure 2. Some of the most salient hydrological features are

summarized below:

> The site is located on a topographic plateau bordered immediately by wetlands
and surface waters. A cliff exists between the site in the highland area and
surrounding wetlands at lower elevations. Surface seeps can be detected at the

cliff at some locations.

12

The aquifer, formed of glacioﬂuvial deposits, outwash and till, has a relatively

low conductivity, especially at lower elevations in the wetland area.

The aquifer is underlain by a continuous layer of clay, providing an effective
barrier to vertical migration of contaminants. The clay layer was encountered at
approximately 115 boring locations and was not completely penetrated.

The area within the Zephyr site, being relatively ﬂat and covered by permeable
top soil, has high potential for inﬁltration.

Extracted and treated water is injected entirely back into the aquifer in and around

the contamination source areas.

13

 

386$ 08. E 98 85 no #52502 miiﬁwm m BEE

 

 

e eew, . .
Illll . e 3.. / me5):.HHI
al\ mod 35 o

 

 

l/ x .... .
\\ ﬂ»... i. / e 8 ea .3 / :03 33:82 8
xi/l- \ av: as as ./ /
/./ at & 3, / :owo
\\\\ all! cw ./ U 1H
\\ L, / a. ,.
xx. / // I E \)/ \ /

 

l4

cocamcoﬁom can 835 .8058 m8. mo cosﬁaowuaom 339800 v BEE

 

 

 

3:: HI

v.0 Nd o

833m :25;

bavczom .252 28m Ezewom D
23>» owbi .

 

enema;

 

 

15

       
 

52M cowoxmsE

use—$3

“Em 53qu com souoomémEU 3:38:00 .m Emmi

 
 

 

H93..— >20

shag. Suﬁsm:

    

( x020 Rom

 

\

  

some;
cocoa .m:

:03
cocoabxm
=03

accustom
=03 . m
wstoxcoz ES.
Eogmoﬁt

              

l6

All of this contributes to the signiﬁcant, site-wide groundwater mound, with a network of

“embedded”, well-scale drawdown cones, resulting in a multiscale groundwater ﬂow

system.

To quantify such a complex ﬂow system, the following sources and sinks are explicitly
modeled: natural recharge, extraction wells, extraction trenches, injection trenches, Bear

Creek, Bear Lake, Muskegon River, Muskegon Lake, surface seeps, and wetlands.

The highly variable terrain in the model area is represented using the 10m high resolution
digital elevation model (DEM). The land surface is treated as a drain boundary, allowing
groundwater to discharge to surface where groundwater level intercepts the land surface.
The drain elevation is set to be equal to the detailed, DEM-based groundwater surface
elevation. All major surface water bodies, including Muskegon River, Muskegon Lake,
Bear Creek, and Bear Lake are represented as constant head boundaries, with the water
levels assumed to be approximately equal to the DEM elevations. The Muskegon
wetlands, Bear Creek Wetlands, and other surface seeps were represented as part of the

land surface drains.

A total of 55 wells out of the 59 extraction wells are currently in operation and are
included in the model. Although pumping rates of extraction wells vary at times, the

slight changes were not represented as we focus on long term mean conditions. Pumping

l7

rates used in the model are summarized in Table A1 in the Appendix. The extraction

trench system, including both the new and old trenches, is modeled as line drains or head
dependent ﬂuxes. The new trenches are relatively deep (14 it) but narrow (1ft) while the
old trench is shallow (4 ft) but wide (311). The injection trenches, 8 in the south part of the
site, and 2 in the north, are represented as line sources with prescribed ﬂuxes. The
injection trenches are ﬁlled with highly permeable peat gravels and treated water is
continuously injected to the trench through a well. The prescribed line ﬂux per length is
set to be equal to the rate of injection well divided by the length of the trench. The ﬂux
values for all injection and extraction under the current design condition are presented in
Table 1 and Table 2. The total injection rate for all injection trenches is equal to the total

pumping rate of all extraction wells and trenches.

Table l Flux Values for All Injection

 

 

 

 

 

 

 

 

 

 

 

Injection Trench Injection Rate in GPM
IT-l 5.28
IT-2 20.52
IT-3 32.25
IT-4 26.97
IT-S , 29.32
IT-6 32.83
IT-7 31.66
IT-8 I 24.04
IT-9 25.80

IT-lO 19.93

 

 

 

 

18

Table 2 Flux Values for all Extraction

 

 

 

 

 

 

 

 

 

 

 

Extraction Trench Extraction Rate in GPM

NT0809 0

NT1011 -2.6

NT1213 -6.1

NT1415 -2.7

NT1617 -7.6

NT1819 -6.5

Old Extraction Trench -64

Test Trench for Free Product -2.4

 

Natural recharge is modeled as spatially variable and divided intoiﬁve zones based on
' land use (Figure 6). The ﬁve recharge zones are: site proper, two residential districts on
the plateau (NE and SW side of site), the Muskegon wetland area, and the Bear Creek
Wetland area. Hydraulic conductivity is divided into three zones based on how the
aquifer materials were generated geologically and the Michigan land system information

available from the statewide groundwater database. The different conductivity zones are

shown in Figure 7.

19

 

 

 

  
  
  

Legend

|: Recharge Zones

0 0.25 0.5
-—_-:3Miles

 

 

Figure 6 Divisions of Recharge Zone

 

 

Legend

I:I Conductivity Zones

0 0.25 0.5
_:':_1Miles

 
   

 

 

Figure 7 Divisions of Conductivity Zone

20

-. n F; wry-.-

 

The clay layer underneath the aquifer is assumed to be impervious. The clay elevation

was interpolated using universal Kriging with a linear drift. The northeastern model
boundary is selected such that it coincides with a streamline and can thus be represented
as an approximately no ﬂow hydraulic boundary. The “remote” streamline is delineated

based on regional static water levels from the Michigan statewide groundwater database.

Hierarchical Flow Modeling

Under the HPDP, we model the complex ﬂow system incrementally, visualize the results
in realtime, and “zoom” into subareas when and where we feel there is a need to.
Groundwater ﬂow is modeled as unconﬁned and two-dimensional. Slight vertical

variation in the center of the mound and in the discharge area is ignored.

We begin with modeling the entire area using a coarse-grid and then make localized
corrections by adding patches or patches-in-a—patch (submodels nested within submodels)
where the solution is judged to be inaccurate. The patch boundaries are interactively
located where the parent model dynamics are deemed to be adequately resolved. The
boundary conditions for each new submodel are represented as prescribed heads based on
linear interpolation of the nearest 4 nodal heads from the preceding parent model under a
shared node based grid system. For this example, no feedback loop is considered between

parent models and their nested submodels (only one-way interaction is implemented). A

21

 

visual sensitivity analysis is used to evaluate the sufficiency of the patch boundary

location.

Model results, including calibration, sensitivity analysis, and application to remediation

performance evaluation are presented in the sections that follow.

Model Calibration

The groundwater ﬂow model was calibrated under a steady state condition. Considering
computer computation compatibility and to increase calculation efﬁciency, one site scale
model is created with boundary perpendicular to head contours generated by regional
scale size model (Figure 8). The calibration is conducted based on Site Scale Model. The
calibration parameters are conductivities and recharge in different zones (Figure 9 and
Figure 10), and bottom elevation and leakancy of new extraction trenches, old extraction
trench and test trench of free product. The calibration targets are static water levels
collected at 128 monitoring wells throughout the site on October 28, 2008 with
remediation system on and seepage ﬂux from new trenches, test trench and old extraction
trench. Calibration is achieved by minimizing the Sum of Squared Differences between
the simulated and observed heads, and simulated and observed ﬂux. The ﬁnal calibrated
conductivities and recharge are presented in Table 3. The bullets below provide several

comments regarding the calibration.

22

 

 

 

 

 

 

 

 

   

Legend

|:l Recharge Zones

0 0.2 0.4
_= Miles

 

Figure 9 Divisions of Recharge Zones for Site Scale Size Model

 

 

   

Legend

I:I Conductivity Zones

0 0.2 0.4 ,
_::r Miles

 

Figure 10 Divisions of Conductivity Zones for Site Scale Size Model

24

 

 

 

Table 3 Values of Parameters selected as model input

 

 

 

 

 

 

 

Conductivity Recharge
(ﬁ/ day) (inch/ year)
Zephyr Site 12.39 18.29
Muskegon Wetland 11.41 -1.00
Bear Creek Wetland 2.41 1.00
Upper Community 12.39 1.00><10'5
Lower Community 12.39 10.37

 

 

 

 

Table 4 Values of Parameters selected as model input for extraction trench system

 

 

 

 

 

 

 

 

 

 

Bottom Elevation in m Leakancy in m/ day
NT1011 178.43 0.51
NT1213 175.82 0.92
NT1617 178.46 1.40
NT1819 179.04 92.48
Test Trench 178.96 0.32
Old Extraction Trench 177.12 2.09

 

 

> Although same uniform values in conductivity and recharge are used as initial
guesses for all zones, automatic calibration leads to spatially variable values in
different zones.

> Different initial guesses for the calibration parameters lead to essentially the same
estimated values for recharge and conductivity in the site polygon. The large

numbers of water level measurements enable unique calibration, given the

25

relatively physically-based parameter zonation.

> The distribution of the calibrated values for recharge, conductivity, and bottom
elevation for extraction trenches makes sense and is consistent with data and our
conceptual understanding

> Estimated recharge at the site polygon is signiﬁcantly higher than the estimated
values in the residential areas and in the wetland areas.

> Estimated conductivity at the site polygon is relatively low, but higher than that in

the wetland areas

 

 

 

187 ""‘l’_"1" 1 1 'l r"'r‘ T ' l—'T‘ITT

I I I I I I I I I xI’x".

186 . X":
185 ,. ,1

e 184 W”
E I I I | I I / I I J
a 18311:
8 I I I I I .g I I I I
m 182 ": ‘ or“: I '1
B 181 ,. fin-Lu: --I
2 ' ,yI /9'I'
_§18°ii.'//::r:II:
m 179 -~--:————I;.'- .
.78 ,

_ ./_I_--I_-_I --_I---_L_--I-_-I-_-I_--I---I

177 XX : : : : : I : : : I
176 II I I I I I I I I i i
176177178179180181182183184185186187

Observed Head in m

 

 

 

Figure 11 Comparison of observed and simulated heads at steady state for 176 monitoring
wells in the site scale size model

26

 

 

 

O “c"'”"r “““ r “““ I ““““ T """" I """ 'I “““ r’“"

.s : : : : : : : :
a -50 ------ -: ----- I ----- : ----- I ----- :-----: -----
£1 I I I I I I I I
0 I I I I I I I I
E l I I I I I I I
a 400%": “““ r """ : """ I """ : “““ . """ rm":
I: I I I I I I I I
9 I I I I I I I I
E. -150 - ----- 1 ----- :- ----- : ----- I ----- . ----- r-----I
>< 5‘ : : : : : : : :
In"? -200 ‘l'"“'I """ r """ I --------- I ---- 'l ''''' r-"""1
E e : : : : : : : :
9% I I I I I I I I
I. -250« ----- -: ----- . ----- r----I ----- r----:
2 I I I I I I I I
e : : : : : : : :
B ‘300‘ “““ 1 """ . “““ : """ T """ : """ 1 """ IN“.
.53 I I I I I I I I
E ~350- ----- . ----- : ----- I ----- :-----1 ----- :-----1
m I I I I I I I I
-400 F I I I I I I I

-400 -350 -300 -250 -200 -150 ~100 -50 0

Observed Flux from Extraction Trenches in m3/day

 

 

 

Figure 12 Comparison of observed and simulated ﬂux at steady state for extraction
trenches

Figure 11 presents the comparison of observed and simulated heads at steady state for
176 monitoring wells in the site model. The predicted head distribution matches well with
the observed. The root mean square error is approximately 0.40 m or 5% of the maximum
observed head difference across the site. The arithmetic mean error is almost zero at 0.1m
or 0.5% of the maximum observed head difference. The errors at the vast majority of the
128 monitoring wells are within one standard deviation. These results clearly show that
the model is accurate, unbiased and able to capture the dominant processes. Figure 12

presents the comparison of observed and simulated ﬂux from extraction trench at steady

27

 

states. The calibration result shows the simulated extraction trench system in the model is

in a working condition close to reality.

Model Sensitivity Analysis

A sensitivity analysis of the steady-state model was conducted to assess the uncertainty in
the input values. The analysis was used to determine if the differences between simulated
and observed data values could be accounted for by the range of uncertainty in the values
of input parameters. The analysis provided a measure of the sensitivity of the model
results to changes in the values of key parameters and, thus, provided a check on the
calibrated model. Throughout the model area, the principal input parameters were
independently decreased by a constant factor of 2%, while other parameters were left
unchanged. The gradient of the objective function based on parameters was used as a
measure of sensitivity. A summary of the gradient for all parameters is included in Table
5. The analysis indicated that model simulations are most sensitive, in decreasing order of
importance, to (a) the conductivity of the plateau and recharge of the site (b) the recharge
in the Muskegon wetland area and (c) conductivity in Muskegon wetland and Bear Creek
wetland respectively. Other parameters including grid size, and recharge of lower

community had a relatively small inﬂuence.

28

Table 5 Results for Sensitivity Analysis

 

 

 

 

 

 

 

 

 

 

 

 

Parameters Gradient
Recharge of the site 0.0407
Conductivity of Plateau 0.0360
Recharge of Muskegon Wetland 0.0170
Conductivity of Bear Creek Wetland 0.0128
Conductivity of Muskegon Wetland 0.0103
Grid Size 0.0018
Recharge of Lower Community 0.0006
Recharge of Bear Creek Wetland 0.0000
Recharge of Upper Community 0.0000

 

29

 

88m 38% 595 so Samba. couﬂcoEB HEB corn—£55 exam 2 SEE

 

 

 

 

 

(1.x

 

on: e

 

 

w—

./ mozznl
., N_ .o cod 0
.552 5 ﬁne: ow—

as:

 

“snowed

 

 

30

88m 38% cops: to 889$ 533382 5;» cousatummv Ron 3 oSmE

 

 

 

VI“.
1
,

wt:

/ s

’

wwto:

8:2

\/\\

         

Wot

N_.o 0V6 o

 

L808 5 can: 2:

use:

 

ncowoq

 

 

31

ANALYSIS OF REMEDIATION SYSTEM PERFORMANCE

In this section, we investigate systematically the remediation system performance at the
Zephyr site, taking advantage of the model’s ability to simulate not only large scale

dynamics but also detailed near-well dynamics.

Hierarchical Capture Zone Analysis

We begin with an analysis to quantify the capture zones for the extraction system. We
achieve this by performing reverse particle tracking for all capture wells based on

hierarchically—modeled velocity ﬁelds.

Figure 15 presents results from hierarchical reverse particle tracking. Up to three levels of
sub-models are created with spatial step size of 30 ft, 10ft, and 3.3 ft, respectively. At
each level, multiple patch models are created. Level I model represents the regional
model described above. Level II model is the ﬁrst level of sub-model that uses the
prescribed heads simulated in the level I model. Similarly, level 111 models were created,
using prescribed heads at the boundaries simulated in the previous level of the model.
Figure 16 presents a hierarchical “family tree” illustrating the relationships between the
patch models.

Although the models were developed hierarchically at multiple resolutions, only the ones

32

ﬁnest resolution properly represent rapidly varying well dynamics and their capture

zones. Level I model, at a resolution of 30ft, is only sufficient in delineating the general
orientation of the capture zones. Level 2 model, at an improved resolution of 10 ft, can be
used to delineate the general outline of the capture zones, but signiﬁcantly underestimates
the capture widths. Level 111 models, zooming into ﬁve focused areas at a substantially
reﬁned resolution at 3ft, enable resolving detailed capture zones for wells that pump at
relatively high rates (closer to 10 GPM). making it possible to delineate in detail
small-scale drawdown cones of depression and the associated capture zones even for

wells with pumping rates less than 5gpm in a large, complex model.

The hierarchical modeling results clearly show the challenges in containing the
contamination at the Zephyr site, given the large plume size, small sustainable extraction
rates, and diverging ﬂow pattern. The models show that signiﬁcant gaps exist in the
current capture system despite the large numbers of extraction wells used. The area that is
most problematic is north side of the site where the well distribution is relatively sparse
and hydraulic gradient is strong. The widths of the predicted gaps between individual

capture zones range from 100ft to 500 it.

The models predict even wider gaps in the well capture zones on the south side of the site
(north of the cliff) where the hydraulic gradient present is strongest. But this area, unlike

the north site, is also protected a series of “new” and “old” extraction trenches further

33

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 15 Hierarchical Reverse Particle Traching

34

 

 

 

 

Level [-1 Coarse
Level I

Fine
Level 1 +1

 

 

Figure 16 Hierarchical “family tree” to illustrate the relationships between the patch
models.
downstream. Contamination that escapes the wells can be potentially captured by the
extraction trenches. The effectiveness of the trench system is further evaluated in the next

section.

The area that is best protected, based on the model prediction, is in the southwest portion
of the site where the well distribution is denser and pumping rates are higher. Individual

capture zones for the different wells overlap to form a contiguous, larger capture zone.

Hierarchical Forward Particle Tracking

In this section, we investigate how contamination at the Zephyr site can potentially

migrate offsite, given the predicted gaps in the capture system.

To obtain a conservative estimate about the plume migration, we model only advective

transport, ignoring degradation and sorption. We release particles throughout the site

35

where data shows elevated contaminant concentrations and use forward particle tracking

to simulate plume evolution under steady ﬂow condition.

Fig. 18 presents the ﬁnal hierarchical forward particle tracking solutions presented in
7patch models across 3 levels. Figure 17 presents data showing total BTEX distribution
at the site. Fig. 18 MI shows the regional solution at the top of the hierarchical-tree. Fig.
18 Mi? describes in more detail the subregional dynamics and wellﬁeld interactions. Fig.
17, from M? to M; , presents, respectively, a subsubmodel for each wellﬁeld capturing
more detailed well interferences. The rest of the subplots present successions of
patch-models zooming into 5 focused-areas that play critical roles in the overall scheme

of integrated management.

The models clearly show that the particles expand in all directions. The most dominant
migration pathways are towards west, south, and north. The west movement eventually
branches into two directions - northwest and southwest. The northwest branch moves
towards Bear Creek but stopped expanding in the wetland area near Bear Creek. The
southwest branch moves toward the Muskegon River but stopped spreading in the

Muskegon wetlands.

Particle transport towards the south boundary is probably most intense because eight of

the ten injection trenches are located on the south site and are close to the site boundary.

36

Models xxx clearly show that particles migrate past the capture wells through the gaps

between the well capture zones. The models also show that particles that escape the
capture wells also passes through the new trenches and old capture trench before being

captured by the Muskegon wetlands.

Particle transport towards north boundary is also signiﬁcant because Bear Creek bends
towards the site and is closest to contamination at the junction where Little Bear Creek
and Bear Creek meet. At this location Bear Creek almost directly borders the site with
virtually no buffer zone in between. Forward Particle Tracking shows that particles
eventually expands into Little Bear Creek, which is connected with Bear Creek, and will

ﬁnally move towards the Bear Lake at the downstream end.

The models also show that particles released at the site migrate east in a direction
opposite to the natural regional gradient. The cast movement also bifurcates in two
directions - northeast toward Bear Creek and southeast toward Muskegon River. But the
movement in both directions stopped in the wetland areas near Bear Creek and Muskegon

River.

It becomes clear that the wetlands on both side of the site act as buffer zones that absorb

contamination and arrest further migration of the contamination.

37

335% Le ease warez xmem 88 E 2sz

 

 

 

         

/ :/ -\\
anon, .
_ $295 Nam-11" // E .-
._ 0 01-10II‘ , /,-..- ..\.\
3:. we? ., --.11 Rm... -61.... o _
u....-.1wwm_.1-..Ia . R: moo.» ... a.» a .
. ea
-.......,-1. . ,
x O o d ,,
vmm ..0. 1 N
3.3 m: .-
- I..- -4.
36mm SEEHI
3mm 3 3o o

ammo—om 2229* E8 guczom
83:4 - 8.88 .
8.38 - :32 .
8.2K - :82 .
8.8. - :2 .
2: - 85 o
z a: 5. $588
@834

 

 

38

 

 

 

 

 

 

 

 

 

Figure 18 Hierarchical Forward Particle Trackil

39

 

 

 

 

 

 

 

 

 

 

 

 

1PM

. s
1.1.... w. -

.- u-«I. Juana/.3; . .

- I .s11
M1,? J 2}. - 1.,
. s. I.

. .. \

   

. a-

.f I

I
s
v .

 

.

 

 

.\.
0 ’ I\
..v. f. . . .
o _\
u r

 

 

Water Balance Analysis

In this section, we quantify the total seepage ﬂux off-site. We achieve this by creating a
water budget for the aquifer in the site area. The budget analysis is performed based on
modeling results, both for natural condition and under active remediation. The analysis

area and the results are presented in Figure 19.

One signiﬁcant message from the budget analysis is: with or without remediation the
seepage ﬂux offsite is about the same and is approximately equal to natural recharge in
the site area. Speciﬁcally, the predicted seepage ﬂux out of site boundaries was
approximately 349m3 /day without remediation, and approximately 356 m3 /day with

remediation.

We were initially puzzled by the ﬁnding, but it becomes immediately obvious upon a
closer examination of the water budget. The ﬂux offsite has not change much is simply
because what is pumped out of the site is injected back into the ground within the same
area. Under natural condition, ﬂow out of the site at steady state is essentially balanced
by natural recharge, as other sources and sinks (e.g., inﬂow from the boundaries and
surface seeps) are relatively small. Under active remediation, the ﬂow offsite is still

balanced by natural recharge, since total injection on-site is designed to be equal to total

40

 

 

1500

 

 

1000 -

 

 

 

 

         

 

 

 

é“
-o
33'
‘3- 500
16..
‘5’
,2 Extraction
g 0 _ Flux out Trench
~53 Flux in Injection
.3. Trench
L1.-
-500
-1000

 

Sources and Sinks

I System Off (Natural Condition) I-i System on - Steady State

 

 

 

Figure 19 Water Budget for Site Area

pumping, and boundary inﬂows and surface seeps are again relatively small. Of course,
the quality of the water moving offsite with and without remediation should be

signiﬁcantly different, and contaminant ﬂux offsite should decrease with time.

41

CONCLUSION

In this thesis, we investigated the impact of the current remediation system on the
groundwater ﬂow and contaminant migration at the Zephyr site in Muskegon, Michigan.
A two-dimensional hierarchical groundwater modeling system is created to evaluate in
detail the efficiency and effectiveness of the remediation system. Water balance analysis
based on the modeling results show that it is not possible to capture all the groundwater
moving offsite unless part of the extracted water is disposed offsite. Hierarchical reverse
particle tracking reveals that signiﬁcant gaps exist among purge well capture zones,
contributing to leakage of contaminant. Hierarchical forward particle tracking shows that
contamination that escaped capture is eventually “arrested” by the Muskegon and Bear
Creek wetland system. The hierarchical modeling system is being applied to redesign
remediation injection, extraction, and ‘offsite disposal needed for maximal cleanup

efficiency and effectiveness.

42

APEENDIX A

Table A l Pumping Rate for Purge Wells

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Well Longitude in Meter Latitude in Meter Pumping Rate in GPM
PW—l 480465.95 302234.15 -2.40
PW-2 480419.61 302194.20 -6.19
RWl 480231.61 302163.02 -l.60
RW2 480253.30 302169.86 -0.80
RW3 480275.18 302161.82 -2.80
RW4 480295.92 302162.52 -4.19
RW5 480317.82 302162.41 -3.69
RW6 480337.97 302126.56 -2.70
RW7 480364.57 302126.92 -l.90
RW-8 480453.37 302113.75 0.00
RW-9 480477.47 302117.67 0.00
RW-10 480509.74 302115.80 -1.36
RW-ll 480535.92 302117.83 -1.20
RW-12 480564.73 302124.77 -3.52
RW-13 480587.74 302133.05 -2.64
RW- 14 480620.19 302146.54 0.00
RW-lS 480643.74 302163.01 -2.72
RW-16 480688.74 302171.46 -3.28
RW-l7 480714.37 302176.93 -4.32
RW-18 480754.75 302181.56 -3.04
RW-19 480782.44 302184.78 -3.52
RW—20 480775.28 302238.36 -2.40
RW-21 480730.23 302237.75 -1.90
RW-22 480663.90 302234.25 -1.60
RW-23 480587.58 302235.70 -9.99
RW-24 480481.17 302168.51 -6.39
RW-25 480210.44 302215.59 -4.59
RW-26 480210.17 302297.20 -3.99
RW-27 480424.57 302432.15 -9.59
RW-28 480487.69 302445.71 -2.60
RW-29 480380.98 302486.25 -4.99
RW-30 480419.73 302518.00 -6.19
NW—l 480474.43 302501.45 -3 .79
NW—2 480507.12 302526.98 -3 .40
NW-3 480499.89 302553.87 -3.99
NW-4 480386.86 302766.01 -2.00

 

43

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

NW-5 480382.63 302624.29 -1.70
NW—6 480316.46 302637.67 .3.79
NW—7 480273.32 302630.13 -1.60
NW—8 480228.50 302604.85 2.00
NW—9 480177.98 302541.02 -1.00
NW—lO 480179.30 302513.72 -1.60
NW—ll 480179.34 302487.45 -1.60
NW—12 480105.64 302469.06 -0.80
NW—13 480090.74 302447.85 -O.80
NW-14 480095.59 302395.75 .1.20
NW-lS 480299.89 302403.25 .3.20
NW-16 480601.20 302438.43 -O.80
NW-17 480202.76 302581.23 -0.80
NW-18 480300.88 302514.57 -2.80
NW-19 480288.58 302555.75 -l.60
NW-20 480293.06 302426.92 -120
NW—21 480304.38 302379.31 .2.40
NW—22 480344.45 302736.74 -1.80
NW—23 480809.61 302247.85 .4.39
NW—24 480061.82 302348.25 -1.60

SL-l 480248.41 302206.09 .5.59

SL-2 480323.22 302212.44 .5.59

SL-3 480391.17 302151.56 -6.39

SL-4 480435.43 302167.97 .4.79

SL-S 480450.21 302169.26 .3.79

 

 

 

 

 

REFERENCE

Li, S.G, and Q. Liu. 2003. Interactive Ground Water (IGW): An Innovative Digital
Laboratory For Groundwater Education and Research. Computer Applications in
Engineering Education 11 no. 4: 179-202.

Li, 8.0, Q. Liu, and S. Afshari. 2006. An object oriented hierarchical patch dynamics

paradigm (HPDP) for groundwater modeling. Environmental Modeling and
Software 21 no. 5: 601-758.

45