We were fascinated to learn about how Dr. Fitts became
interested in learning, coding, and applying AEM models. But we are especially
grateful for the in-depth responses that he provided on some of AnAqSim’s
features, complete with references for those who want to explore this topic
further. We want to offer our thanks to Dr. Fitts, and we hope you
enjoy his responses below.
Q: How did you first become interested in groundwater flow, mathematical modeling, and the AEM?
A: I always
liked mathematics, maps, water, and geology. As a teen, my parents called me
“Chuck Atlas” because I spent so much time perusing the Rand McNally atlas. As an undergrad I majored in Geology/Biology
at Colby and then I studied engineering geology and soil mechanics at Cornell
for an MS. Then in the early 1980s I worked
for four years doing geotechnical consulting in the Boston area. During those years I worked on many projects
with seepage and groundwater issues. I took
a semester off work to take engineering classes at U. Minnesota (my home state)
to prepare for PE exams, and one of the courses I took was Groundwater
Mechanics with Dr. Otto Strack, the founder of the analytic element method
(AEM). He is a gifted teacher and I really
enjoyed learning analytic techniques and applying them to groundwater flow problems. I decided to pursue more graduate study under
Otto, and eventually got a PhD in civil engineering. My research there focused on AEM methods for
vertical plane flow, 2D heterogeneities, and 3D heterogeneities. I was drawn to the AEM because of the beauty of
the mathematics and the promise of the method with its high accuracy, minimal
inputs, and geometric flexibility.
Q: What led you to develop AnAqSim with an approach that is quite different from previous AEM programs?
A: For my
research at U. Minnesota, I wrote Fortran programs which I kept improving and
using in my consulting work after leaving school. Ultimately this code became TWODAN
(two-dimensional analytic groundwater model), which I distributed in the 1990s
and 2000s. TWODAN was a good tool with a
decent interface, but like most AEM programs at that time, it was limited to a
single layer (2D) and had almost no transient capabilities (merely adding Theis
solutions to the steady solution). In
the 2000s, the AEM declined in popularity while numerical methods like FDM
(MODFLOW) and FEM were on the rise. This
trend was probably due to the robust multilayer (3D) and transient capabilities
of numerical models, whereas AEM programs were generally limited to 2D steady
flow. Rather than keep up with multiple
modeling techniques and programs, many modelers just used one numerical modeling
program. I felt that the key to growing
AEM use was to improve its 3D and transient capabilities.
I saw a
significant break when Strack and Jankovic (1999) published their paper on solutions
for spatially variable area sinks (SVAS). These are 2D functions with a sink term whose strength (rate) varies
with position in the plane using radial basis functions. SVAS functions create a smooth, irregular
surface of sink strength that passes through any number of basis points where
the strength is set. The radial basis
functions used by Strack and Jankovic (1999) were first employed by Hardy
(1971) to interpolate smooth topographic surfaces from point elevation
data. Strack and Jankovic demonstrated the
SVAS functions in a model simulating leakage from a single-layer aquifer to an
overlying surface water using relatively few basis points. I saw that these SVAS functions could be used
in a different way with larger numbers of basis points to dramatically improve multi-layer
and transient modeling with the AEM. The
distributions of vertical leakage between layers in multi-layer models and the storage
fluxes in transient models may be combined and viewed as spatially variable
distributed sinks. The functions
presented by Strack and Jankovic could be used to approximate the required sink
distribution.
On another
front, I worked on methods for modeling anisotropy in the plane of the domain,
using some old (Muskat, 1937) and new coordinate transformations (Fitts,
2006). At this time, I also worked on a
subdomain approach in AEM where you break up a model layer into separate finite
subdomains instead of one infinite domain (Fitts, 2004; Fitts, 2010). Each subdomain has its own mathematical model
and can have its own direction and ratio of anisotropy, making for very robust
anisotropy capabilities. A subdomain
approach has been used with boundary element methods, which are like AEM but
use numerical integrations along boundary elements instead of analytic
solutions (Liggett and Liu, 1983; Bruch and Lejeune, 1989). The subdomain approach has other nice
side-effects: (1) it makes the equations for potential or discharge shorter,
(2) model layering can vary from one region of a model to another, and (3)
heterogeneity boundaries can be a mix of different boundary types such as
specified head, specified normal flux, and interdomain.
It was exciting
to see how all these pieces could come together. I began writing code that used subdomains,
SVAS, and coordinate transformations to make AEM models that had multi-layer
and general anisotropy capabilities. I
also incorporated high-order line elements which make it easier to tune
accuracy and use longer line segments (Jankovic and Barnes, 1999). I later added storage terms and transience using
the same SVAS functions as are used for vertical leakage between layers. The result is a program that has the
accuracy, minimal inputs and other benefits of AEM, but with comprehensive
capabilities for anisotropy, multi-layer systems, and transient flow.
The vertical
leakages are computed using finite difference approximations of dh/dz
and storage fluxes are computed using finite difference approximations of dh/dt. The flow models are based on the AEM while
these derivatives are approximated with finite differences. Such approximations are employed by MODFLOW
and other numerical modeling programs as well and afford a wide range of
application to both AnAqSim and numerical models.
Q: What are the main differences between AnAqSim’s methods and methods used in other AEM modeling codes?
A: This topic has many facets, and for a quick
overview you may want to see this web page comparing AEM program capabilities: https://www.fittsgeosolutions.com/analytic-element-method/. I will discuss the major differences below.
Most AEM programs
have elements that are all internal to a single infinite domain in each
layer. In contrast, AnAqSim uses finite
subdomains that extend laterally only as far as their polygonal external line
boundaries. Each subdomain in AnAqSim
has its own mathematical model, whereas the infinite domain approach has a
single model per layer. This means the
equations in AnAqSim are shorter (only elements within or on the boundary of a
subdomain contribute to the potential or discharge functions of that
subdomain), which leads to computational savings while solving the system of
equations and while generating graphic outputs like contours and
pathlines. Subdomains also allow complex
layering schemes where, for example, three subdomains (layers) on one side of
an interdomain boundary abut a single subdomain (layer) on the other. This lets you easily have more layers in the
areas of interest and fewer layers farther out. No other AEM program has this feature.
Other AEM
programs typically specify a reference head at a point outside the area of
interest, which determines the flow to or from infinity. Since AnAqSim uses finite subdomains, this is
not needed. The external boundary
conditions in AnAqSim’s bounded models are more clearly defined than the
infinity boundary condition generated by the reference head in infinite-domain
AEM models. I know from my years of
supporting TWODAN that the reference head was a source of confusion.
Most AEM
programs are single layer (2D), but AnAqSim, MLAEM, and TimML all offer
multi-level aquifer simulation (3D). TimML (and MLU, which is pumping test analysis
software based on the same approach) is a multi-level simulator that provides
exact distributions of the sink terms simulating vertical leakage and storage,
whereas AnAqSim approximates these distributions between basis points. TimML has some limitations at this point in
its development. For example all domains
must be fixed transmissivity (confined), layers must be infinite, fewer kinds
of line and area element boundary conditions are possible, and it lacks a
graphic user interface. MLAEM appears to
allow simulation of leakage in multi-level aquifer systems but after
considerable online searching, I could not find documentation of how that is
accomplished with area sinks, and the models are only for steady-state flow
(not transient).
AnAqSim has
comprehensive transient capabilities meaning that just about all key inputs
(discharges, heads, recharge rates, normal fluxes, river stages) can be
transient, varying from one time period to the next. Head-specified wells and line boundaries can
have transient heads and they can turn off during specified periods. This on/off feature is handy for dewatering
simulations. TTim is a transient AEM
program like TimML that can also do transient simulations with transient
discharges and heads, but TTim is subject to the same kinds of limitations
listed above for TimML.
AnAqSim has
general anisotropy capabilities, allowing unique directions and ratios in each
subdomain. There are no other AEM
programs that allow anisotropy in the plane of domains.
Q: Describe the process of creating and updating AnAqSim. What were (are) some of the big hurdles?
A: AnAqSim
is written in C#, an object-oriented language with many nice built-in
components supplied by Microsoft. In the
2000s I switched from FORTRAN to C# because AnAqSim’s subdomains and transient
time-stepping are complex and require a lot of bookkeeping that C#’s classes
and objects are good at. I chose C#
because it had the right language features, Microsoft is likely to be around
for a while, it comes with a great development environment (Visual Studio), and
there were good libraries for math and graphics. AnAqSim has tens of thousands of lines of
code, so I had to build it incrementally and slowly from the ground up. I started with fundamental classes to
represent spots (x, y, level, subdomain, potential and discharge functions),
subdomains (elevations, properties, coordinate transformations, conversions
between head and potential), and elements (wells, line elements, area
sinks). I tested the daylights out of
these fundamental classes before building higher-level classes (e.g. boundary
conditions, equation-building) that use these lower-level classes.
It’s
important to test the code thoroughly every step of the way and whenever
anything is modified or added. Bugs are a
fact of life with complex software and it takes time to root them out. Fortunately, AnAqSim has now been in use for
almost a decade with many hundreds of users, so it is mature in the sense that
serious bugs have been found and cured. I built lots of analysis features into
AnAqSim, in part because the outputs are useful to users, but also because it
makes thorough testing of the code much easier for me.
The biggest
hurdle was the sheer magnitude of the project. It is a lot of work to write the computational parts of the code, the
user interface, the documentation, the training tools, and the web site. I like the challenge of big projects and at this
point I have invested probably a decade of work on this. I use bought components for the plotting,
graphics, and math (complex variables, solvers), and these have saved time and
made it possible. A difficult, but
unseen part is creating an installation package that a user runs to install the
software. This is a whole field unto itself,
with specialized software that creates the install file and numerous arcane
considerations such as .net frameworks, C++ libraries, dynamic link libraries, digital
signatures, target platforms, etc.
Q: How has your use of groundwater modeling software evolved over time?
A: In the
1990s I used TWODAN for 2-D steady flow analyses and MODFLOW for anything that
was 3D and/or transient. I also used
Aqtesolv, which is good software for analyzing pumping tests using various
radial single-layer analytic solutions. I still use Aqtesolv for some pumping test analyses, but it has been a
few years since I’ve used MODFLOW. I find
that I can reasonably model most any situation with AnAqSim. I have been using AnAqSim often to analyze
pumping tests in a more sophisticated way when one or more of these conditions
are present: partial penetration of pumping or observation wells,
heterogeneity, anisotropy, or irregular boundary conditions. With AnAqSim’s hydrograph features, it is
easy to plot simulated vs. observed heads at all observation wells and tweak aquifer
parameters to get a good match of those hydrographs.
Q: What types of groundwater projects have you or others applied AnAqSim to?
A: AnAqSim
has been used on many projects for construction or mining dewatering
simulations, which appears to be the most common application. Dewatering
simulations have involved wells, sumps, trenches, sheet-pile barriers, and
recharge trenches. The ability to turn
specified head wells and line boundaries off is handy for dewatering projects
that proceed in stages (e.g. an expanding and deepening excavation). I have also done transient analysis of the
impact of constructed ponds on local groundwater flows and heads.
Another
common application is capture zone analysis for remediation design and for
water supply zone of contribution studies. These often involve multilayer models in the area of interest to simulate
partially penetrating wells, drains, and surface waters; and they can be
performed as part of groundwater protection zone delineation or groundwater
remediation projects. AnAqSim’s area
pathlines and capture constrain features give it unique capabilities for
capture zone analysis.
Several of
my projects have involved simulating vertical plane flow. One was analyzing discharge to a ship channel
with bulkheads and heterogeneous subsurface conditions, and another was looking
at patterns of discharge from a sand and gravel aquifer to a lake. The anisotropy features are helpful with
vertical plane models, where stratification typically makes subdomains
anisotropic.
I often use
AnAqSim to analyze pumping tests with partially penetrating wells, heterogeneous
and layered systems, and surface water connections near the tested well. An analysis done this way can be
accomplished in just a few hours. You
get a much more sophisticated and realistic analysis compared to the typical pump
test solution for a single layer, no boundaries, and radial symmetry.
I’ve
analyzed fresh/salt interface aquifers a couple of times and I know of other
consultants who have done many projects analyzing interface positions in sandy
coastal aquifers.
Q: What materials are available to help users learn to use AnAqSim?
A: The more
capable a modeling software is, the more there is to learn. Since AnAqSim can do so many things, there is
a lot to learn. So, I have worked hard to
provide detailed documentation and training for users. There is a 120+ page user guide that is
available online and with the software, which explains all inputs and
outputs. There are also instructional
videos and tutorials at this web page, which are a good way to get up and
running: https://www.fittsgeosolutions.com/learning-anaqsim/. McLane Environmental’s flexAEM training materials
http://www.flexaem.com/ are well-designed
and are a great way to learn AEM modeling and AnAqSim in particular.
The https://www.fittsgeosolutions.com/
website has some other useful learning resources, including a detailed
explanation of a transient construction dewatering model, a page explaining how
to use AnAqSim with PEST (parameter estimation software), 19 different checks
of AnAqSim complete with model files, and a reference list of helpful books and
journal articles.
Recently, we
have been offering a six-session webinar series called AnAqSim in Depth
at least once a year. The series aims to
help a user learn how AnAqSim works and how to apply it. Each session is one hour and there is one
session per week. This gives time for
optional home exercises each week that walk you through building or modifying
AnAqSim models that use the content of that week’s session. The next AnAqSim in Depth series starts March 10, 2020 and it is
described at this web page: https://www.fittsgeosolutions.com/webinar/.
Q: What are some useful, but possibly obscure, modeling techniques that you could share with AnAqSim users?
A: My first
tips are about analyzing pumping tests. I often do this using multi-level models with an initial condition of
head = zero everywhere, and just simulate the deviation from zero with the
transient model, like is done in typical pumping test analysis. To do this, set the average head in all
domains to zero, and have your elevation datum be the static water level if you
have unconfined domains at the top. I
select the constant head option for Starting Heads in the solve settings,
which uses a domain’s average head as the starting head. I typically use Head-Specified External
with Gradient line boundaries for the far-field boundary condition in all
levels of the model, and I set the head and gradient at these to zero. Place these external boundaries far enough
out so the drawdown doesn’t reach there by the end of the test. If there are constant head boundaries such as
a surface water nearby, simulate them with a specified head = 0, so there is no
drawdown there. These sorts of models set
up quickly in an hour or so. I calibrate
them manually using the Analysis/Hydrograph features in AnAqSim, but for PEST
enthusiasts calibration can be automated that way.
Another feature that few people seem to use are the SVAS
top/bottom condition surfaces. Most
people apply SVAS top conditions of a constant recharge rate or sometimes of a
constant head for leakage to an overlying surface water. But it is possible to use these surfaces to
model recharge that varies with location, or leakage to an overlying surface
water with a sloping surface (e.g. a wetland or stream). Likewise, at the bottom of a model, it is
most common to specify zero flux (impermeable) everywhere. If there is leakage to/from a deeper layer
that is not in your model, you can use a surface to represent the heads in that
deeper layer below your lowest model layer, and simulate leakage to/from that
deeper layer using a head-dependent flux for the bottom condition.
Q: Any closing remarks?
A: Thank you for the opportunity to do this Q&A. I want to encourage feedback from users, who
help me track down bugs but also come up with suggestions for capabilities
they’d like to see added. I keep a
database “wish list” of these and chip away at them over time.
References
Bruch, E.,
and A. Lejeune (1989), An effective solution of the numerical problems at
multi‐domain points for anisotropic Laplace problems, in
Advances
in Boundary Elements, vol. 2, Field and Flow Solutions, edited by C. A.
Brebbia and J. J. Connor, pp. 3–14, Comput. Mech., Boston.
Fitts, C.R.
(2010), Modeling Aquifer Systems with Analytic Elements and Subdomains, Water
Resources Research, 46, W07521, doi:10.1029/2009WR008331.
Fitts, C.R.
(2006), Exact Solution for Two-Dimensional Flow to a Well in an Anisotropic
Domain, Groundwater, 44(1), 99-101.
Fitts, C.R.
(2004) Discrete Analytic Domains: A New Technique for Groundwater Flow Modeling
in Layered, Anisotropic, and Heterogeneous Aquifer Systems, American
Geophysical Union Fall Meeting, San Francisco, CA.
Hardy, R.L.
(1971), Multiquadric equations of topography and other irregular surfaces. Journal
of Geophysical Research 76, 1905–
1915.
Jankovic,
I., and R. Barnes (1999), High‐order line elements in modeling two‐dimensional
groundwater flow, Journal of Hydrology, 226, 211–223.
Liggett, J.
A., and P. L‐F. Liu (1983), The Boundary Integral Equation Method for Porous
Media Flow, Allen and Unwin, London.
Muskat, M.
(1937) The Flow of Homogeneous Fluids through Porous Media. New York:
McGraw-Hill.
