Hydrology [H]

H13B  MS:Exh Hall B   Monday
Analytical and Semianalytical Models of Subsurface Flow and Transport I Posters
Presiding: J R Craig, University of Waterloo; M Bakker, Delft University of Technology

H13B-1241 

Analyzing non-Darcian flow in a confined aquifer toward a well with a linearization method

Wen, Z (sl002wenzhang@126.com), China Agricultural University, Department of Irrigation and Drainage, College of Water Conservancy and Civil Engineering, Beijing, 100083, China Huang, G (ghuang@cau.edu.cn), China Agricultural University, Department of Irrigation and Drainage, College of Water Conservancy and Civil Engineering, Beijing, 100083, China Huang, G (ghuang@cau.edu.cn), China Agricultural University, Chinese-Israeli International Center for Research in Agriculture, Beijing, 100083, China * Zhan, H (zhan@geo.tamu.edu), Texas A & M University, Department of Geology and Geophysics, College Station, 77843- 3115, United States

In this study, we have developed a new method to analyze non-Darcian flow toward a well in a confined aquifer with and without wellbore storage. The power law has been used to describe the relationship of the specific discharge and hydraulic gradient for non-Darcian flow. This new method is based on a combination of the linearization approximation of the non-Darcian flow equation and the Laplace transform. Approximate analytical solutions of steady-state and late-time drawdowns are also obtained. The drawdowns at any distance and time are computed by using the Stehfest numerical inverse Laplace transform with MATLAB programs. The results of this study agree perfectly with previous Theis solution for an infinitesimal well and with the Papadopulos and Cooper's solution for a finite-diameter well under the special case of Darcian flow. The Boltzmann transform, which is commonly employed for solving non-Darcian flow problems before, is problematic for studying radial non-Darcian flow. Comparison of drawdowns obtained by our proposed method and the Boltzmann transform method suggests that the Boltzmann transform method differs from the linearization method at early and moderate times, and it yields similar results as the linearization method at late times. The drawdowns decrease at late times as the power index n or the quasi hydraulic conductivity k increases, regardless of the wellbore storage. It has also been found when n is larger, flow approaches steady state earlier. The approximate analytical solutions indicate that the drawdown at steady state is approximately proportional to r over a power of (1-n), where r is the radial distance from the pumping well; the late time drawdown is a superposition of the steady- state solution and a negative time-dependent term that is proportional to time t over a power of (1-n)/(3-n).

H13B-1242 

Three-Dimensional Flow Generated by a Partially Penetrating Well in a Two-Aquifer System

* Sepulveda, N (nsepul@usgs.gov), U.S. Geological Survey - FISC, 12703 Research Parkway, Orlando, FL 32826,

An analytical solution is presented for three-dimensional (3D) flow in a confined aquifer and the overlying storative semiconfining layer and unconfined aquifer. The equation describing flow caused by a partially penetrating production well is solved analytically to provide a method to accurately determine the hydraulic parameters in the confined aquifer, semiconfining layer, and unconfined aquifer from aquifer-test data. Previous solutions for a partially penetrating well did not account for 3D flow or storativity in the semiconfining unit. The 3D and two- dimensional (2D) flow solutions in the semiconfining layer are compared for various hydraulic conductivity ratios between the aquifer and the semiconfining layer. Analysis of the drawdown data from an aquifer test in central Florida showed that the 3D solution in the semiconfining layer provides a more unique identification of the hydraulic parameters than the 2D solution. The analytical solution could be used to analyze, with higher accuracy, the effect that pumping water from the lower aquifer in a two-aquifer system has on wetlands.

H13B-1243 

Using Analytical Solution Methods to Analyze and Combat Alarming Growth of Errors in Traditional Unsaturated Flow Numerical Computations

* Tracy, F T (Fred.T.Tracy@erdc.usace.army.mil), Engineer Research and Development Center, 3909 Halls Ferry Road, Vicksburg, MS 39180, United States

Recently (WRR, 2006), analytical solutions for Richards' equation, \begin{equation} \nabla · \left (k \nabla h \right ) + \frac{∂ k}{∂ z} = \frac{∂ θ}{∂ t} \end{equation} where h is pressure head, k is the hydraulic conductivity, θ is moisture content, z is the z coordinate, and t is time, were derived for unsaturated flow in a homogeneous three-dimensional (3-D) box region. Also, \begin{equation} k = kr ks \end{equation} where kr is the relative hydraulic conductivity, and ks is the saturated hydraulic conductivity. Here, kr and θ are both functions of h, thus creating a strong nonlinearity in the partial differential equation (PDE). Tools such as the quasi-linear approximation, \begin{equation} \ln k = \ln ks + α h \end{equation} separation of variables, Fourier series, and the change of variable \begin{equation} \bar h = eα h - eα hd \end{equation} were all used. α is a constant, and hd is the pressure head when the soil is very dry. The bigger α is, the more nonlinear the problem is. This 3-D problem and a two-dimensional (2-D) version of the problem in the x-z plane were then solved by traditional finite element and finite volume methods, respectively. It was startling to observe how much the error increased with increasing nonlinearity. This presentation will first summarize the box-shaped problem and its solution. It will then highlight the rapid growth of errors as α is increased (a phenomenon not readily understood by many practicing engineers), even when high performance computing and parallel processing are used to solve for 501 × 501 grids in 2-D and 161 × 161 × 161 grids in 3-D. The presentation will then show great improvement in these errors by using the analytical/numerical combination from first doing the following change of variables: \begin{equation} \bar h = eα h - 1, ~~~~~ \bar kr = kr e-α h, ~~~~~ h ≤ 0 \end{equation} \begin{equation} \bar h = h, ~~~~~ \bar kr = kr = 1,~~~~~ h ≥ 0 \end{equation} before doing the numerical solution. To illustrate, Eq. 1 for steady state becomes \begin{equation} \nabla · \left (\bar kr \nabla \bar h \right ) + α \frac {∂}{∂ z} \left [\left (\bar h +1 \right ) \bar kr \right ] = 0, ~~~~~ h ≤ 0 \end{equation} \begin{equation} \nabla2 \bar h = 0, ~~~~~ h > 0 \end{equation} Using the quasi-linear approximation, Eq. 7 becomes the linear PDE \begin{equation} \nabla2 \bar h + α \frac {∂ \bar h}{∂ z} = 0 \end{equation} Eq. 9 is similar to the constant diffusivity moisture-content-based equation often used in the analytic finite element. Finally, this presentation will show great improvement in convergence of the nonlinear system of equations generated from Eq. 7 as compared to the original nonlinear system. ERDC's Adaptive Hydrology (ADH) 3-D finite element program will be improved based on this work.

H13B-1244 

A Study for Residual Drawdown Solution with Considering the Wellbore Storage

* Wang, C (ctwang.ev90g@nctu.edu.tw), Institute of Environmental Engineering, National Chiao Tung University, No. 75, Po-Ai Street, Hsinchu, 300, Taiwan Yeh, H (hdyeh@mail.nctu.edu.tw), Institute of Environmental Engineering, National Chiao Tung University, No. 75, Po-Ai Street, Hsinchu, 300, Taiwan

A recovery test is to measure the residual drawdown, which is the well water level after the shutdown of pumping test, and analyze the residual drawdown data for the estimation of the aquifer parameters such as transmissivity and storativity. In the past, the methodology of the recovery data analysis is based on the Theis equation associated with superposition principle and thus neglects the effect of wellbore storage even the test well has a finite-diameter. In this paper, we develop a new residual drawdown solution for describing the recovery water level in the pumping well with considering the wellbore storage. The solution shows that the decrease of the residual drawdown depends on the boundary condition related to the well drawdown and the initial condition related to the aquifer drawdown. In addition, the well residual drawdown can be approximated by Theis equation associated with superposition principle with a good result when recovery time is large.

H13B-1245 

Coordinate Mapping of Analytical Transport Solutions to Non-Uniform Flow Fields

* Craig, J R (jrcraig@uwaterloo.ca), University of Waterloo, 200 University Ave West, Waterloo, ON N2L 3G1, Canada

While analytical solutions for solute transport continue to be used as screening and predictive tools for contaminated site analysis, there remain many limitations to their applicability. Despite the recent advances in their ability to simulate complex chemistry (e.g., parent-daughter decay and mixing-driven chemical equilibrium), most existing analytical transport models depend upon the critical assumptions of uniform flow and uniform saturated thickness. A time-of-flight-based mapping procedure has been developed to take existing analytical solutions for transport in uniform flow and map them to simple non-uniform flow fields. This presentation will discuss some of the inherent assumptions in the mapping method, as well as its strengths and its weaknesses. It was found that the mapping approach, exact for purely advective transport problems, can be reliably used in steady-state flow systems with wells, surface water features, and (to a limited degree) mild recharge and heterogeneity. The mapping errors are on the same magnitude as those due to discretization error or upstream weighting in numerical simulation methods.

H13B-1246 

Analytical Solutions for Examining Spatial Variations in Evapotranspiration-Driven Fluctuations in the Water Table in Vegetated Riparian Zones

* Jin, W (jacking@ksu.edu), Kansas Geological Survey, 1930 Constant Avenue, Lawrence, KS 66047-3726, United States * Jin, W (jacking@ksu.edu), Department of Agronomy, 2004 Throckmorton Plant Sciences Center, Kansas State University, Manhattan, KS 66506, United States Kluitenberg, G J (gjk@ksu.edu), Department of Agronomy, 2004 Throckmorton Plant Sciences Center, Kansas State University, Manhattan, KS 66506, United States Butler, J J (jbutler@kgs.ku.edu), Kansas Geological Survey, 1930 Constant Avenue, Lawrence, KS 66047-3726, United States

Diurnal water-table fluctuations observed in shallow wells in vegetated riparian zones are a diagnostic indicator of groundwater consumption by evapotranspiration (ETG). In order to better understand how these fluctuations vary across a riparian zone, analytical solutions for strip-sinks with periodic forcing functions were developed. The periodic forcing functions, which represent the evapotranspirative consumption of groundwater, are finite within the strips and zero outside. Linearity allows superposition of the strip solutions to simulate spatial variations in ETG due to variations in vegetation density within the riparian zone. The solution is applied to data from a field site in a riparian zone along the Arkansas River in west-central Kansas where groundwater levels in shallow wells have been monitored at 15-minute intervals for up to five years. A series of strips were used to represent the largely dry unvegetated river channel, two zones within the vegetated riparian zone, and adjacent pastures and cultivated fields. The spatial variations in the phase and amplitude of the simulated fluctuations were consistent with general patterns observed in the field data during both wet and dry years. Several empirical methods have been developed for estimating ET G using characteristics of the diurnal water-table fluctuations. This solution reveals where best to place wells for use with those methods. Similar solutions have been developed for circular sinks to assess the impact of phreatophyte-control activities on water-table fluctuations at wells in areas of cleared vegetation.

H13B-1247 

Flux-based semi-analytic prediction contaminant transport in a GIS environment

* Becker, M W (mwbecker@geology.buffalo.edu), Department of Geology, University at Buffalo, SUNY, 876 NSC, Buffalo, NY 14260, United States Jiang, Z (zjiang4@buffalo.edu), Dept. of Civil, Structural, and Environmental Engineering, University at Buffalo, SUNY, 212 Ketter Hall, Buffalo, NY 14260, United States

A computationally efficient method for predicting contaminant mass flux to a specified boundary is presented. The method combines streamlines calculated using the analytic element method with a first-passage-time (residence time) approach to calculating mass flux. Computations exploit an analytic forward solution of convolution in Laplace space followed by numerical inversion. The technique is carried out in a geographic information system (GIS) that takes full advantage of widely available digital hydrologic data. The flux-based estimations allow efficient coupling of contaminant mass across hydrologic interfaces. This combined approach is useful for rapid but approximate transport computations either at the local or regional scale. We demonstrate the approach by considering nitrogen-loading to streams due to surface application of animal feedlot waste. Nitrogen is predicted to migrate through ground-water, resulting in a long-term loading to nearby streams. The object-oriented implementation of this modeling approach improves integration between the GIS and the transport modeling, allowing transport results to be returned to the database.

H13B-1248 

An Open-source Community Web Site To Support Ground-Water Model Testing

* Kraemer, S R (kraemer.stephen@epa.gov), US Environmental Protection Agency, Ecosystems Research Division 960 College Station Road, Athens, GA 30605-2700, United States Bakker, M (Mark.Bakker@tudelft.nl), Delft University of Technology, Faculty of Civil Engineering and Geosciences, Delft, 2628 CN, Netherlands Craig, J R (jrcraig@uwaterloo.ca), University of Waterloo, Department of Civil and Environmental Engineering, Waterloo, ON N2L 3G1, Canada

A community wiki wiki web site has been created as a resource to support ground-water model development and testing. The Groundwater Gourmet wiki is a repository for user supplied analytical and numerical recipes, howtos, and examples. Members are encouraged to submit analytical solutions, including source code and documentation. A diversity of code snippets are sought in a variety of languages, including Fortran, C, C++, Matlab, Python. In the spirit of a wiki, all contributions may be edited and altered by other users, and open source licensing is promoted. Community accepted contributions are graduated into the library of analytic solutions and organized into either a Strack (Groundwater Mechanics, 1989) or Bruggeman (Analytical Solutions of Geohydrological Problems, 1999) classification. The examples section of the wiki are meant to include laboratory experiments (e.g., Hele Shaw), classical benchmark problems (e.g., Henry Problem), and controlled field experiments (e.g., Borden landfill and Cape Cod tracer tests). Although this work was reviewed by EPA and approved for publication, it may not necessarily reflect official Agency policy. Mention of trade names or commercial products does not constitute endorsement or recommendation for use. http://www.analyticelements.org

H13B-1249 

Approximate modeling of transient wells and line-sinks with analytic response functions

* Bakker, M (Mark.Bakker@tudelft.nl), Water Resources Section, Faculty of Civil Engineering and Geosciences, Delft University of Technology, Stevinweg 1, Delft, 2628CN, Netherlands * Bakker, M (Mark.Bakker@tudelft.nl), Kiwa Water Research, Groningenhaven 7, 3433PE, Nieuwegein, Netherlands Maas, K (Kees.Maas@kiwa.nl), Water Resources Section, Faculty of Civil Engineering and Geosciences, Delft University of Technology, Stevinweg 1, Delft, 2628CN, Netherlands Maas, K (Kees.Maas@kiwa.nl), Kiwa Water Research, Groningenhaven 7, 3433PE, Nieuwegein, Netherlands Veling, E (E.J.M.Veling@tudelft.nl), Water Resources Section, Faculty of Civil Engineering and Geosciences, Delft University of Technology, Stevinweg 1, Delft, 2628CN, Netherlands

We propose a new approach for the modeling of groundwater head fluctuations due to wells and line-sinks with discharges that are highly variable in time. The approach is based on the use of impulse response functions, the head response due an impulse of discharge. Once the impulse response function is known, the response to any stress that varies with time may be obtained through convolution. Superposition and convolution of the classic impulse response function for a well allows for the simulation of well fields with complicated discharge distributions. This classic approach becomes cumbersome or may break down when other aquifer features are present. When the well field is located near a stream, for example, it is necessary to include the effect of the stream on the impulse response function, and thus on the heads in the aquifer. We have developed an approach where the impulse response function of a well or line-sink may be approximated by a parametric, analytic function. Impulse response functions may be viewed as scaled probability density functions and may be characterized by certain integral characteristics such as the area, mean, and standard deviation. Alternative integral characteristics are the temporal moments and the exponentially-scaled temporal moments (the mean and standard deviation may be expressed in terms of moments). Temporal moments are specifically useful, as they fulfill steady, Poisson-type differential equations, with known values along the boundaries. Hence, the moments of the impulse response function may be computed exactly by constructing (relatively simple) steady models. The exact impulse response function may then be approximated by predefined analytical functions through moment matching. We propose an approximate impulse response function that has four parameters, so we need to determine four moments at a point to compute the approximate impulse response function there. This means that we have to build four steady models for each transient stress. At any point in space we can then compute four moments, and thus the four parameters of the approximate impulse response function. It will be shown that this approach gives accurate transient results for a number of cases, while it requires the solution of a small number of steady models only.

H13B-1250 

Analytic element formulations for flow in leaky aquifers and for transient flow, with discontinuous aquifer parameters.

* Strack, O D (strac001@umn.edu

Analytic elements are functions that satisfy a given differential equation, mostly applicable to groundwater flow, and have properties that make it possible to meet a variety of boundary conditions; the solution to the flow problem is obtained by superposition of suitable analytic elements. Analytic elements are similar to boundary elements, but differ from these in two respects. First, they need not be obtained by integration, and second, and more important, analytic elements are defined throughout the infinite domain, even for cases where properties are discontinuous. In the analytic element formulation presented here all analytic elements are defined throughout the infinite domain, and special analytic elements are added to meet the conditions along the boundaries of inhomogeneities in properties. This is in contrast with boundary element and boundary integral formulations, where the domain is subdivided, with each element defined in its own sub domain, and with boundary elements used to stitch the sub domains together. For the case of problems governed by the Poisson equation, discontinuities in properties imply that the potential function representing the mathematical solution to the problem is discontinuous along with the component of flow tangential to the discontinuity boundary, but the differential equation itself is not affected. Such discontinuities can be dealt with conveniently by the use of the appropriate discontinuous function. If we consider, however, more general cases of flow, such as leaky aquifer flow and transient flow, then the coefficients that appear in the differential equation itself are discontinuous across the discontinuity line, which makes not obvious that the analytic elements can be defined throughout the domain, the discontinuities notwithstanding, and that correcting discontinuous functions can be constructed to include the discontinuities correctly in the solution obtained by superposition. We demonstrate in this presentation, however, that this is possible both for the case of leaky aquifer flow, governed by the modified Helmholz equation, and for transient flow, governed by the heat equation.

H13B-1251 

Contaminants in Coastal Aquifers - An Analytical Approach

* Bolster, D (diogobolster@gmail.com), Universitat Polytecnica de Catalunya, Campus Nord, edifici D2, C. Jordi Girona, 1-3., Barcelona, 08034, Spain * Bolster, D (diogobolster@gmail.com), University of California - San Diego, 9500 Gilman Drive, San Diego, CA 92093, Tartakovsky, D (dmt@ucsd.edu), University of California - San Diego, 9500 Gilman Drive, San Diego, CA 92093, Dentz, M (marco.dentz@upc.edu), Universitat Polytecnica de Catalunya, Campus Nord, edifici D2, C. Jordi Girona, 1-3., Barcelona, 08034, Spain

The Henry formulation, which couples subsurface flow and salt transport via a variable-density flow formulation, can be used to evaluate the extent of sea water intrusion into coastal aquifers. The coupling gives rise to nontrivial flow patterns that are very different from those observed in inland aquifers. We investigate the influence of these flow patterns on the transport of conservative contaminants in a coastal aquifer. The flow is characterized by two dimensionless parameters: the Péclet number, which compares the relative effects of advective and dispersive transport mechanisms, and a coupling parameter, which describes the importance of the salt water boundary on the flow. We focus our attention on two regimes – low and intermediate Péclet number flows. Two transport scenarios are solved analytically by means of a perturbation analysis. The first, a natural attenuation scenario, describes the flushing of a contaminant from a coastal aquifer by clean fresh water, while the second, a contaminant spill scenario, considers an isolated point source.

H13B-1252 

Contaminant transport modeling using the Analytic Element Method and the deterministic Streamline Method

* Bandilla, K w (bandilla@eng.buffalo.edu), University at Buffalo, Dept. Civil, Structural, and Environmental Engineering 207 Jarvis Hall, Buffalo, NY 14260, United States Jankovic, I (ijankovi@eng.buffalo.edu), University at Buffalo, Dept. Civil, Structural, and Environmental Engineering 207 Jarvis Hall, Buffalo, NY 14260, United States Rabideau, A (rabideau@eng.buffalo.edu), University at Buffalo, Dept. Civil, Structural, and Environmental Engineering 207 Jarvis Hall, Buffalo, NY 14260, United States

The Analytic Element Method (AEM) formulation for steady 2D groundwater flow is combined with deterministic Streamline Method (SM) to model large-scale transport of reactive contaminants. AEM is an alternative to the Finite Element (FEM) and Finite Difference Methods (FDM) for solving subsurface flow problems on large scales. The domain is discretized along the hydrogeologic elements (e.g. surface water features, zones where conductivity differs from the surrounding conductivity, etc.) instead of using a grid discretization as in FEM and FDM. Two features that make AEM well suited for a basis for the SM are particle tracking without interpolation and the strong parallel processing capabilities. In the implementation presented here a 2D steady-state groundwater flow simulator is used to solve the flow problem in the horizontal plane. Vertical velocities are computed based on mass balance considerations, thus leading to a quasi 3D flow field. The 3D particle tracks are then used in the SM for contaminant transport. The SM discretizes the transport domain by converting curvy 3D streamlines into straight 1D streamlines. The conversion is achieved by transforming the Cartesian coordinates of the transport domain into a 1D coordinate system based on the `time-of-flight' coordinate, which describes the time for a particle to travel a distance along the streamline. The transport along the streamline can then be solved using FEM or FDM. Two beneficial features of the Streamline Method are the comparative ease of solving a set of uncoupled 1D models instead of a single fully-coupled 3D model and the independence of streamlines which leads to efficient parallel processing. The capability of this approach to model large scale reactive contaminant transport is shown based on a test case. The influence of reaction complexity and parallel processing will be shown.

H13B-1253 

Analytic Element Solution For Inhomogeneities In Hydraulic Conductivity Placed In Anisotropic Porous Background

Jankovic, I (ijankovi@buffalo.edu), Department of Civil, Structural Engineering, University at Buffalo, 207 Jarvis Hall, Buffalo, NY 14260, United States * Suribhatla, R (rms29@buffalo.edu), Department of Civil, Structural Engineering, University at Buffalo, 207 Jarvis Hall, Buffalo, NY 14260, United States

The Analytic Element Method (AEM) was originally developed (e.g. Strack, 1989) for solving steady-state groundwater flow problems in two-dimensional (2D) and three-dimensional (3D) isotropic domains. In the present study the AEM approach is extended to solve for head and specific discharge due to inhomogeneities (inclusions) of isotropic conductivity that are embedded in anisotropic background (the 2D/3D conductivity tensor has two different principal components). The original AEM approach uses a single coordinate system (e.g. radial for circular inclusions and elliptical for elliptical inclusions) to construct AEM solutions from general solutions to the Laplace equation. For the anisotropic case presented here, the solution procedure involves scaling the original Cartesian coordinate system (and inclusions' boundaries) to transform the governing equation in the anisotropic background to the Laplace equation; no scaling is used for inclusions' interiors. The complete solution is hence based on two separate coordinate systems. Boundary conditions (continuity of head and normal component of specific discharge) are rewritten in terms of the two corresponding coordinate systems. Scaling relations between the two coordinate systems are established in order to derive explicit linear relations between various solution coefficients. The infinite series solutions are truncated for numerical implementation. During implementation, both boundary conditions are satisfied approximately but with high precision using a moderate number of coefficients. To demonstrate accuracy of the AEM solutions, comparison with highly-resolved Finite-Difference based solution is presented.

H13B-1254 

An exact solution to a line-sink in a leaky aquifer

* Gusyev, M A (mgusyev@indiana.edu), School of Public and Environmental Affairs, 1315 East Tenth Street Room 439, Bloomington, IN 47405, United States Haitjema, H M (haitjema@indiana.edu), School of Public and Environmental Affairs, 1315 East Tenth Street Room 439, Bloomington, IN 47405, United States

By use of Wirtinger calculus we obtained an exact solution for a line-sink in a leaky aquifer by integrating the potential for a well in a leaky aquifer. The potential for a well in a leaky aquifer is the modified Bessel function of the second kind and zero order K0, which can be represented by an infinite series. Theoretically, this series expansion for the well is exact, although numerical evaluation will only give exact results within some finite distance from the well, depending on machine accuracy. For a double precision machine this distance is about to 18λ, whereby λ is the "leakage factor" or "characteristic leakage length" which depends on the aquifer properties. Earlier solutions that were based on an approximation to the function K0 limited the domain of validity even more; from 2λ to 8λ away from the well. As a result, these earlier (approximate) solutions for a well in a leaky aquifer limited the length of the line-sink along which it could be integrated to approximately λ. It appears that our use of the infinite series (exact representation of K0), makes it possible to formulate a solution for a line-sink of any length, thus avoiding to need to break up line-sinks into smaller sections as has been done to date. Formulating our solution in terms of the complex variable z and its conjugate \bar{z}, using Wirtinger calculus, also allows us to calculate the exact integrated steady-state flow induced by the line-sink across an arbitrarily placed line element. This feature is often necessary in the context of the analytic element method in order to satisfy boundary conditions in terms of integrated fluxes, such as no-flow boundaries or leaky walls. The capability to accurately calculate such integrated fluxes across line elements is also important in order to obtain the integrated leakage over a domain by applying water balance rather than (numerically) integrating the leakage directly. The new solution is particularly suitable for use in analytic element models of leaky aquifer systems.