1. Kris Kuhlman
  2. slug-osc



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 Hydrology, 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 or 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 transform 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