This repository contains the unconfined (moving water table - aka
Neuman 1972) slug test solution given in

Malama, B., K.L. Kuhlman, W. Barrash, M. Cardiff & M. Thoma,
2011. "Modeling slug tests in unconfined aquifers taking into account
water table kinematics, wellbore skin and inertial effects", Journal
of Hydrology
, 408(1–2), 113–126.

which is updated and applied to cross-hole slug tests in

Malama, B., K.L. Kuhlman, R. Brauchler & P. Bayer, 2016. "Modeling
cross-hole slug tests in an unconfined aquifer", Journal of
, 540:784–796.

2015 Update

The model has been modified since publication in 2011, to
compute the response at an observation location outside the source
well (i.e., cross-hole slug tests), to use the finite Hankel transform
for inversion, include effective source/observation well skin, and to
include observation well storage and inertial effects. The current
version of the code still works for its original purpose as a
traditional slug-test analyzer.

Cross-hole slug test data from Widen, Switzerland and the associated
parameter estimation results (PEST and DREAM) are located in the
"swiss-data" directory of the repository.

Code License

The code is available under the MIT open-source license, which allows
free use and redistribution; please acknowledge the original authors
by citing our papers.

Compilation Directions

The code requires gfortran (>= 4.6) compiler to properly compile (one
that includes libquadmath), but I have gotten it to work with some
minor tweaks using the latest version of the Intel Fortran compiler
(which supports quad precision, but not for the complex hyperbolic
trig functions -- so they must be implemented in Fortran, rather than
in the compiler library).

The code can be compiled, given you have the correct compiler set up
already, by calling the compiler-specific makefile:

make -f Makefile.gfortran debug_driver
make -f Makefile.gfortran driver

these two makefile targets build the debug target (slower, checks for
bounds overflows, etc) and the standard target (builds using OpenMP
for a parallel build which is optimized and much faster.

Input Instructions

The input.dat file has all the input parameters read by the 2011
program. The finite-Hankel cross-hole slug test program (2015) uses
input-finite.dat. Parameters include physical properties of the
formation, well geometry parameters, numerical convergence parameters
and the locations/times to compute the solution at. Anything to the
right is a comment and is ignored by the program when reading in the
input. Units are arbitrary but consistent.

For the cross-hole version of the program, consult
"cross-hole-slugtest-input-explanation.txt" for discussion about each
variable appearing in the input file (input-finite.dat).

Numerical Approach

The numerical approach used to invert the Laplace- and infinite Hankel
numerically is a slightly modified and updated version of
the approach documented in Appendix B (repeated below citation) of:

Malama, B., K.L. Kuhlman, and A. Revil, 2009. "Theory of transient
streaming potentials associated with axial-symmetric flow in unconfined
aquifers", Geophysical Journal International, 179(2), 990-1003.

The numerical approach used to invert the double Laplace-Hankel
inversion is an improvement on the method in Malama et al. (2009). In
this study the Laplace transform is now inverted using the Stehfest
algorithm (Stehfest 1970) and a portion of the Hankel integral is
integrated using an accelerated form of tanh-sinh quadrature (Takahasi
& Mori 1974). The infinite integral in eq. (A2) for inverting the
Hankel transform is evaluated in two steps (Wieder 1999). One finite
(0 <= a <= j_{0,1}), the other approaching infinite length (j_{0,1} <= a <=
j_{0,m}), where j_{0,n} is the nth zero of J_0(ar_D) and m is typically
10-12. The finite portion is evaluated using tanh-sinh quadrature
Takahasi & Mori (1974) where each level k of abscissa-weight pairs uses
a step size h = 4/2k and a number of abscissa N = 2k , similar to the
approach of Bailey et al. (2000). Evaluating the function at Carissa
related to level k, the function evaluations at lower levels k-j
are simply the 2^{k-j} th terms in the series (only requiring the
weights to be computed for each level, reusing the function
evaluations from the largest k). Richardson’s deferred approach to the
limit (Press et al. 2007) is used to extrapolate the approximate
solutions to the limit as h -> 0 using a polynomial extrapolation. The
infinite portion of the Hankel integral is integrated between
successive j_{0,n} using Gauss-Lobatto quadrature and these alternating
sums are accelerated using Wynn’s epsilon algorithm.

The finite Hankel transform used in the 2016 paper is a simple sum, and
is therefore much simpler than the above approach.

Contact Code Author

If you have any questions or require assistance using the code, please
contact me (Kris Kuhlman) at klkuhlm at sandia dot gov. I will
try to help you. I mostly work in Linux and on a Mac, but I have
gotten codes to work in Windows before.

July 2016