- Email: [email protected]
Principles of the analytic element method
HYDROL 3853
Journal of Hydrology 226 (1999) 128–138 www.elsevier.com/locate/jhydrol
Principles of the analytic element method O.D.L. Strack* University of Minnesota, Department of Civil Engineering, Minneapolis, MN 55455-0220, USA Received 5 September 1997; accepted 7 October 1998
Abstract An overview of the development and current state of analytic element method, as well as a brief presentation of the principles on which the technique is based is presented in this article. Illustrations of a few representative analytic elements are also presented. A brief discussion is added of both rotational flow and transient flow. The latter two aspects of the method require further development before they are ready for implementation in computer programs and can be applied to the modeling of groundwater flow. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Analytic element method; Groundwater modeling; Mathematical modeling
1. Introduction The analytic element method was developed for the mathematical modeling of the flow of groundwater and was originally intended for problems of regional flow. The analytic element method is a numerical method that gives an approximate solution to a problem. Analytic elements are mathematical functions that differ from the many existing classical solutions in that they are not restricted to a single boundary value problem, but rather possess degrees of freedom that makes it possible to combine them. The method is based upon the superposition of analytic functions and might be considered as an extension of classical models of regional flow based on the superposition of elementary functions, for example the Thiem and Theiss solutions for wells, and the potential for uniform flow. Characteristic for such classical models is that each function is unaffected by the other functions in the solution. The discharge * Fax: 11-612-626-7750. E-mail address: [email protected] (O.D.L. Strack).
of each well in a field of wells, for example, is known a priori and is not affected by the discharges of the other ones. The analytic element method presented as a modeling approach by Strack (1989), differs from the classical analytical models in that the elementary solutions possess degrees of freedom that may be chosen in such a way that boundary conditions are met along certain internal boundaries. The technique resembles classical boundary element methods, and differs from these in the nature of the analytic solutions chosen. In boundary element methods, these solutions are boundary integrals of the Cauchy type (double or single layers), the boundaries are usually closed, the integrals are often evaluated numerically, and the formulation used to compute the unknown coefficients is usually different from that used to compute values in the field once the coefficients have been determined. The solutions used in the analytic element method are not necessarily Cauchy integrals; elements may be constructed, for example, by generalizing closed-form analytic solutions obtained by the use of conformal mapping techniques. Examples of such elements are the so-called double-root elements
0022-1694/99/$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S0022-1694(99 )00 144-4
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
and the elements for circular and elliptical inhomogeneities presented by Strack (1989), the analytic elements useful for modeling flow near corners formed by internal boundaries presented by Detournay (1990), the three-dimensional analytic elements for partially penetrating wells presented by Haitjema and Kraemer (1988), and the three-dimensional solutions for inhomogeneities bounded by ellipsoids presented by Fitts (1990). Haitjema (1982) introduced an approach for the regional modeling of the threedimensional flow in confined aquifers based on a combination of two-dimensional and three-dimensional solutions. A major difference between the analytic element method and classical boundary element methods is the use of area-sinks. Area-sinks are functions that represent extraction (or infiltration) over areas bounded by polygons, and are useful for modeling quasi three-dimensional flow. Such flows occur in shallow confined and unconfined aquifers where, at least in part of the domain, the Dupuit–Forchheimer assumption has been adopted. The area-sink, as presented by Strack (1989), is an analytic element that simulates uniform infiltration or extraction over an area that is bounded by a polygon. The use of area-sinks makes it possible to simulate leakage between aquifers, for example, by subdividing the domain into polygons, each with its own uniform rate of leakage. Area-sinks of this kind form the mainstay of the Dutch National Groundwater Model as presented by De Lange (1991) that was constructed by the use of the Multi-Layer Analytic Element Model (MLAEM). The individual leakages of the area-sinks are computed from the difference in head between the aquifers above and below the element. A slightly more sophisticated area-sink, bounded by a triangle and with a leakage that varies linearly, was also presented by Strack (1989) and was used by Haitjema and Strack (1985) to produce a crude model of transient flow. The triangular areasinks were used in the latter model to simulate the release of water from storage. Strack and Jankovic´ (1999) developed the potential for an area-sink that simulates extraction or infiltration that varies over the domain enclosed by a polygon as a radial-basis interpolator (Hardy, 1971). Only available resources, such as computer memory, limit the number of degrees of freedom of this interpolator. The development of this
129
new area-sink makes possible a broader definition of the analytic element method, presented here for the first time. This enables us, in theory, to determine approximate closed-form expressions for the components of general vector fields, i.e. vector fields that have both a non-zero divergence and a non-zero curl. We may use these expressions to apply the analytic element method to aquifers with continuously variable properties, such as hydraulic conductivity, base elevation, and thickness. Another development, based on the combination of the area-sinks with the technique of separation of variables, makes it feasible to create accurate and efficient analytic element models of transient flow. The combination of these advances with the technique of overspecification introduced by Jankovic´ and Barnes (1999a), and the superblock approach discussed briefly at the end of this paper and presented in Strack et al. (1999), will make it possible to develop the analytic element method into an accurate and efficient method for the simulation of general cases of groundwater flow. We describe the analytic element method briefly in what follows. We do this for the two-dimensional steady flow, for the sake of brevity, but the technique is not restricted to two dimensions. A few comments will also be made with regard to transient flow, overspecification, and the superblock approach. It is important to note that some of the recent advances in the analytic element method exist at present only in theory or have been implemented only in part.
2. General approach We formulate the analytic element method in terms of the discharge vector which represents the total flow per unit width that occurs over the saturated thickness of the aquifer. If the flow is irrotational, then the discharge vector may be expressed as the gradient of a scalar potential. The piezometric head may be determined directly from this potential. This is not possible if the flow is rotational; in that case a potential does not exist and the piezometric head must be determined after the discharge vector has been obtained as a function of position. We restrict this discussion to problems of regional
130
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
flow in systems of confined, unconfined, or semiconfined aquifers. The flow in each aquifer is described in terms of its own two-dimensional discharge vector. 2.1. Irrotational flow The discharge vector field is defined by suitable boundary conditions and expressions for the divergence and the curl of the field. We first consider the case of irrotational flow. In this case there exists a potential, F ; the discharge vector may be represented as minus the gradient of this potential, i.e. Qi 22i F
1
where 2i represents the derivatives with respect to the global Cartesian coordinates x1 and x2, or xi, i 1; 2; … i.e. 2i
2 2xi
2
We adopt the Dupuit–Forchheimer assumption, so that the potential can always be expressed in terms of the piezometric head f . The relation between head and potential may be either linear (e.g. for the cases of semi-confined and confined flow) or nonlinear (e.g. for the case of unconfined flow.) We list as examples the expression for the potential for shallow confined flow,
F kH f 2 12 kH 2
3
F
kf
2
1ij 2i Qj 0
5
and 2i Qi g
and the one for shallow unconfined flow, 1 2
layer. This expression for the divergence may be nonlinear in terms of the potential, for example in the case of leakage between an unconfined and a confined aquifer. The discharge vector is expressed in terms of the gradient of the potential rather than head, and therefore an iterative method of solution is required when the expression for the divergence of the discharge vector is non-linear in terms of potential. We approach the problem of determining the head and discharge vector throughout the aquifer in the following manner. We first estimate the divergence throughout the flow domain and represent it as a function of position. This function is taken to be piecewise smooth; the domain is subdivided into areas bounded by polygons and the divergence is represented as either a constant or an interpolator inside each polygon. The divergence jumps across the boundaries of the polygons and is zero in the area not covered by a polygon. The potential function for an area-sink where the leakage is represented by a radial basis interpolator (the multi-quadric interpolator, see Hardy, 1971) is presented by Strack and Jankovic´ (1999). Next we write the basic equations to be solved, adopting the Einstein summation convention. We represent the divergence of the discharge vector as g and the governing equations for Qi become
4
where k represents the hydraulic conductivity, H the thickness of the confined aquifer and f the piezometric head. The divergence of the discharge vector for any given aquifer equals the net extraction rate per unit area for each aquifer. This divergence may be estimated a priori in the cases of rainfall or leakage through the bottoms of rivers or lakes that are above the groundwater table. In many cases of practical interest, however, the divergence is not known a priori but may be expressed in terms of parameters of the problem. For the case of leakage between aquifers, for example, the divergence equals the difference in head across the leaky layer divided by the resistivity of that
6
where 1 ij is the alternating tensor of the second rank, and 1ij 2i Qj represents the curl of the discharge vector. Note that Eq. (5) is satisfied by representing the discharge vector Qi as the gradient of a potential as in Eq. (1). The analytic element method makes use of superposition of solutions; the discharge vector, Qi, is composed of a sum of functions. These functions may be combined in two groups, the first contains g all area-sinks, and may be represented as Q i ; the h
second is written as Q i ; and represents the sum of all remaining functions, which can be expressed as a sum of the gradients of harmonic potentials, whence the superscript h. The extraction rates of the areasinks have been chosen to match the known function
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
131
Fig. 1. Piezometric contours for an area-sink with infiltration rate shown in Fig. 2; maximum and minimum values are 49.970 and 50.085 m, respectively. The contour interval is 0.5 cm. g
g and thus the discharge vector Q i satisfies Eq. (6), i.e. g
2i Q i g
7
The total discharge is represented as the sum of two parts, g
h
Qi Q i 1 Q i
8
and in view of (7) the second part is the gradient of a harmonic function, h
h
Q i 22i F
9
where h
72 F 0
10
Thus, all functions other than the functions that represent the area-sinks are harmonic. The complete solution to the problem thus may be obtained by
superposition of the potentials for area-sinks and additional harmonic potentials that simulate the effect of a variety of internal boundaries upon the flow. Examples of the latter functions are the potentials for linesinks (sources or sinks distributed along lines or curves), line-doublets (double layers, representing a jump in potential along lines or curves) and line dipoles (distributed dipoles with their axes tangent to lines or curves, representing a jump in the stream function along the lines or curves). An efficient manner of representing such functions is by the use of complex variable techniques. The parameters in an analytic element model must be determined via an iterative process when the problem is non-linear. Many of the boundary conditions applied in regional models involve the piezometric head, and the expression for the divergence (often equal to leakage between aquifers) also is a function of the piezometric head in the aquifer. The discharge vector is the gradient of the discharge
132
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
Fig. 2. Contours of constant infiltration rate inside the polygon that bounds an area-sink; the maximum and minimum values plotted are 0:9 × 103 and 20:9 × 1023 m=day; respectively. The contour interval is 0.001 m/day.
potential. If one of the aquifers in the system is unconfined, the relation between the potential and head of that aquifer is non-linear and an iterative procedure must be followed in order to determine all parameters in the solution. The expression for the divergence g is updated at each step in the iteration until sufficient accuracy has been achieved, see Strack (1989). The theory as outlined above has been implemented and is used in computer models of groundwater flow, particularly in large-scale computer models such as those of The Netherlands (De Lange, 1991), of a future retention basin in The Netherlands (Moorman, 1997), and the Twin Cities Metropolitan Area (Seaberg et al., 1997). As an illustration, we present plots of a few representative analytic elements. The effect of an area-sink by itself is illustrated in Fig. 1, which shows the pattern of piezometric contours, produced by an area-sink with a given rate of extraction that varies over the area. The element is bounded by a polygon and resides in an infinite aquifer; there are no other features present to affect the flow. The aquifer is unconfined with a hydraulic conductivity of 10 m/day. The width of the area shown is 400 m and the head at a distance of 1000 m from the center of the area-sink is fixed at
50 m. The values of heads vary between a maximum of 50.085 and a minimum of 49.97 m. The contour interval is 0.5 cm. The maximum value occurs to the left of the center of the figure and the minimum value to the right. The distribution of the extraction rate over the area-sink is shown in Fig. 2. The extraction varies nearly linearly from a value of 210 22 m/day at the extreme left of the area-sink to a value of 10 22 m/day at the extreme right. Note that a rate of infiltration is modeled by means of a negative rate of extraction. The variation in infiltration rate from positive to negative values causes the maximum and the minimum values to be visible in the plot of piezometric contours as shown in Fig. 1. Doublet elements are used to produce a jump in the potential (see Strack and Haitjema,1981; Strack, 1989) that occurs across the boundary of an inhomogeneity in the aquifer properties. An example of the use of doublets to simulate a jump in aquifer properties is shown in Fig. 3. There is uniform infiltration of 0.1 m/day over an area that includes the inhomogeneity, which is bounded by a polygon. The width of the area shown is 400 m. The flow is unconfined everywhere; the base of the aquifer is at elevation 0 m, the hydraulic conductivity is 10 m/day, and the head at a
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
133
Fig. 3. Piezometric contours for flow in an aquifer with an inhomogeneity in aquifer properties; maximum and minimum values are 65.3 and 59.6 m, respectively, and the contour interval is 0.1 m.
reference point 10,000 m away from the inhomogeneity is 50 m. The aquifer is unconfined and the base inside the inhomogeneity is at an elevation of 50 m, i.e. well above the base outside the inhomogeneity. The hydraulic conductivity inside the inhomogeneity is 0.1 m/day, i.e. hundred times lower than that outside. The pattern of piezometric contours is shown in Fig. 3; the heads vary between a minimum of 59.5 m occurring outside the inhomogeneity and a maximum of 65.3 m, visible just to the left of the center of the polygon that binds the inhomogeneity. The contour interval is 0.1 m. Note that the contours are much closer together inside the inhomogeneity than outside as a result of the contrast in the hydraulic conductivity. The jump across the doublet elements is approximated by a polynomial of the order of 20; the number of degrees of freedom associated with each side is 20. The solution has been obtained using the
method of over-specification proposed by Jankovic´ and Barnes (1999b). The technique of over-specification relies on allocating a larger number of control points along the boundary segment than there are degrees of freedom in the doublet function used to model the jump. The unknown parameters associated with the jump function are determined by requiring that the sum of the squares of the errors at all control points is a minimum. The accuracy obtained by this approach is about an order of magnitude higher than obtained by the conventional technique of matching boundary conditions at collocation points that are equal in number to the degrees of freedom of the jump function. The degree of over-specification chosen for this problem was 2, i.e. there are 2 times as many control points as the degrees of freedom, i.e. 40 per side. The piezometric contours shown in Fig. 4 correspond
134
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
Fig. 4. Piezometric contours for the case of Fig. 3 under perched conditions; the maximum and minimum values are 58.7 and 33.1 m, respectively, and the contour interval is 0.5 m.
to the same problem, except that the head at the reference point is 5 m rather than 50 m. The values of heads on the contours vary between 33 m outside the inhomogeneity and a maximum of 59 m inside. The contour interval is 0.5 m. The water table inside the inhomogeneity becomes perched; the infiltration causes water to flow from the domain inside the inhomogeneity with its base at 50 m to the domain outside where the base is at 0 m and the head is well below the base inside the inhomogeneity. The head now jumps across the boundary; the boundary condition is that the head inside the inhomogeneity equals the base elevation wherever the head outside is below 50 m and that the head is continuous across the boundary wherever the head outside the inhomogeneity is larger than 50 m. This is illustrated in Fig. 5, which is a detail of Fig. 4. The values of the head along the contours near the boundary of the inhomogeneity
are indicated in the figure. The validity of the use of the Dupuit–Forchheimer solution for perched conditions was demonstrated by Fitts and Strack (1996) by presenting exact solutions for two-dimensional cases of unconfined flow in an aquifer with a stepped base. The flow net is shown in Fig. 6 for a string of leaky curvilinear elements. The curvilinear elements were developed by Strack (1986); the derivation is based on the principles presented in Strack (1989). The elements represent a leaky wall of a width of 1 m (the width of the window displayed in the figure is 400 m) and a hydraulic conductivity of 10 23 m/day; the hydraulic conductivity of the aquifer is 10 m/day. There is a far field of uniform flow from left to right in the figure of discharge 0.5 m 2/day in the aquifer and the flow is unconfined. The hydraulic conductivity is 10 m/day. Note that the flow net consists of equipotentials and streamlines; the potential is given for this
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
135
Fig. 5. Detail of Fig. 4.
case by Eq. (4). The values of the potential vary between 1:283 × 104 and 1:315 × 104 m3 =day: The stream function varies between 2120 and 80 m 3/ day. The contour interval for both functions is 10 m 3/day. Note that the interval in piezometric head between contours is not constant, because the relation between the head and the potential is non-linear. The values in head vary between 30.65 and 51.3 m over the region. The approach outlined above may, in theory, be extended to problems of three-dimensional flow. However, an efficient mathematical formulation for the three-dimensional equivalent of the area-sink, a volume-sink, has not yet been determined and threedimensional applications have so far been limited to specialized three-dimensional models of limited scope (e.g. Fitts, 1990), and to models that combine a two-dimensional formulation with three-dimensional analytic elements (Haitjema, 1982, 1985).
Other noteworthy applications of the approach are the combination of an analytic element model capable of simulating the interchange between surface water and groundwater (Haitjema and Mitchell-Bruker, 1996), and the implementation of variable density flow (Strack, 1995). It may be worth noting that the latter application does not cause the curl of the discharge vector to become rotational, provided that the Dupuit–Forchheimer assumption is adopted. 2.2. Rotational flow The analytic element method has not yet been applied to problems of rotational flow. It is possible to cover rotational flow, such as occurs in aquifers with a sloping base and in aquifers with continuously varying properties (rather than properties that are uniform in sub-domains). The approach is similar to that explained above for irrotational flow, and rests on Helmholtz’s theorem. This theorem states that any
136
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
Fig. 6. Flow net for a curvilinear leaky wall; the contour interval is 10 m 3/day.
vector field may be decomposed in a rotational but divergence-free field and an irrotational field with non-zero divergence. We therefore may add to an irrotational field, as described above, a rotational field that is used to generate the required rotation as a function of position. As for the divergence, we assume that the curl is given; the curl is, in fact, determined by some iterative process, and is fixed during each step of the process. The curl is modeled by means of an area-vortex with a vorticity distribution given by an interpolator function. We represent the curl as 2b , rewrite Eq. (5)
1ij 2i Qj 2b
11
and repeat Eq. (6) 2 i Qi g
12
Note that both equations are linear in terms of the discharge vector Qi as b and g are treated as given functions of position. Thus, superposition of
individual functions Qi remains possible. Non-linearities that may exist in the problem have not been eliminated, but have been hidden in the iterative process for determining the functions b and g . The discharge vector field is now separated into three b
g
h
parts, Q i ; Q i ; and Q i ; where the first vector field has rotation b and is divergence-free, i.e. b
1ij 2i Q j b
13
and b
2i Q i 0
14
The second vector field has a divergence of g and zero curl, and the third vector field has both zero divergence and zero curl. The latter two vector fields may be written in terms of the gradient of a potential, but it is important to note that this potential is no longer explicitly related to the piezometric head;
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
this relation is implicit and heads must be computed by some approximate method from the vector field. The b
vector field Q i is divergence-free and may therefore be written in terms of the curl of a scalar stream function, b
b
Q i 1ij 2j C
15
The curl may be applied to the solution by the use of area-vortices, which may be presented in terms of the sum of a stream function, equal to an interpolator inside a polygon and equal to zero outside, plus an harmonic potential that represents the boundary integrals required to create a vector field that is continuous across the boundary of the area-vortex. The area-vortex is analogous to the area-sink, and the functions that describe the flow generated by the area-vortex may be obtained directly from those of the area-sink with the same density distribution. The main advantage of the approach is that the majority of functions are harmonic, and thus complex variable techniques may be used to determine most of the analytic elements, whereas the functions that contribute to the divergence and curl of the vector field are simple and require minimal computational effort. There is a difficulty associated with the approach for rotational flow that has as yet not been resolved fully. The representation of the vector field as the sum of the gradient of a potential and the curl of a stream function is the cause of this complication. The discharge vector field and the piezometric head are related via Darcy’s law. The application of boundary values in terms of head will require that heads be computed from the discharge vector field, at least at isolated points, which requires integration, either numerical or analytical. The use of an iterative process will be necessary to arrive at a final solution, determining the unknown parameters in the expressions for all functions, including the area-vortices, until some criterion for accuracy has been met. It is not the objective of the current paper to present detailed possible approaches for obtaining solutions for cases of rotational flow, but rather to present the principle of the approach.
3. Transient flow A few comments will be made about transient flow.
137
An analytic element model for transient flow was introduced by Zaadnoordijk (1988) and Zaadnoordijk and Strack (1993), based on the integration of transient wells along boundaries. The disadvantage of that approach is that many parameters need to be remembered as the solution progresses through time, and that the computational effort is significant. An alternative to that approach is to use area-sinks to simulate the release of water from storage using area-sinks of linear density distribution, as applied by Haitjema and Strack (1985). The area-sinks with variable strength are an attractive alternative to the earlier elements that were capable only of constant or linearly varying storage over the element. 4. The superblock approach The superblock approach relies on the notion that the majority of analytic elements in a model are harmonic functions and thus may be represented as the real parts of analytic functions of a complex variable. We subdivide the aquifer into a number of square blocks and determine for each block which elements lie far enough away from the boundary of the block that they may be represented inside the block by a Taylor series. All elements outside the block may thus be combined into a single Taylor series expanded about the center of the block. Computation of the potential in each block thus consists of two parts. The first involves evaluation of all functions that are considered to be inside the block. The second part consists of evaluating the effect of all the elements outside the block, which are combined into a single series, thereby greatly reducing computation time. Superblocks may also be used in determining the unknowns in the solution. 5. Conclusion The analytic element method in its present form is applicable to steady flow in aquifers with properties that are piecewise constant. It has proven to be capable of dealing with large aquifer system, and in obtaining highly accurate solutions. Future developments are promising. One of these is the superblock approach, which is expected to make it possible to deal effectively with models on a
138
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
nation-wide scale (e.g. the Dutch National Groundwater Model, NAGROM, De Lange, 1991) and a state-wide scale (e.g. the Twin Cities Metropolitan Model, Seaberg et al., 1997). Other possibilities for future development are the inclusion of rotational flow, applicable to cases of continuously variable aquifer properties, and of transient flow in the analytic element method. References De Lange, W.J., 1991. A groundwater model of The Netherlands. Note No. 90.066, National Institute for Inland Water Management and Waste Water Treatment, Lelystad, The Netherlands. Detournay, C., 1990. On a Cauchy integral element method for potential flow with corner singularities. BETECH 90, Computational Engineering with Boundary Elements. Fluid and Potential Problems, vol. 1. University of Delaware, pp. 119–130. Fitts, C.R., 1990. Modeling three-dimensional groundwater flow about ellipsoids of revolution using analytic functions. PhD Thesis, Department of Civil Engineering, University of Minnesota, Minnesota. Fitts, C.R., Strack, O.D.L., 1996. Analytic solutions for unconfined groundwater flow over a stepped base. Journal of Hydrology 177, 65–76. Haitjema, H.M., 1982. Modeling three-dimensional flow in confined aquifers using distributed singularities. PhD Thesis, Department of Civil Engineering, University of Minnesota, Minnesota. Haitjema, H.M., 1985. Modeling three-dimensional flow in confined aquifers by superposition of both two- and threedimensional analytic functions. Water Resources Research 21 (10), 1557–1566. Haitjema, H.M., Kreamer, S.R., 1988. A new analytic function for modeling partially penetrating wells. Water Resources Research 24 (5), 683–690. Haitjema, H.M., Mitchell-Bruker, S., 1996. Modeling steady-state conjunctive groundwater and surface water flow with analytic elements. Water Resources Research 32 (9), 2725–2732. Haitjema, H.M., Strack, O.D.L., 1985. An initial study on thermal energy storage in unconfined aquifers. Final Report to Battelle Pacific Northwest Laboratory.
Hardy, R.L., 1971. Multiquadric equations of topography and other irregular surfaces. Journal of Geophysical Reseach 62 (2), 733– 737. Jankovic´, I., Barnes, R., 1999a. Three-dimensional flow through large numbers of spheroidal inhomogeneities. Journal of Hydrology 226, 224–233. Jankovic´, I., Barnes, R., 1999b. High order line elements for twodimensional groundwater flow. Journal of Hydrology 226, 211– 223. Moorman, J.H.N., 1997. Analytic Element Model analysis of the influence of different scenarios for the water level in a future retentions basin, on local and regional hydrology. Proceedings of the International Conference Analytic-based Modeling of Groundwater Flow, Nunspeet, The Netherlands, pp. 105–111. Seaberg, J.K., Hansen, D.D., Block, B.W., Streitz, A.R., de Lange, W.J., Bakker, M., 1997. Development of a regional groundwater model for the Twin Cities Metropolitan Area, USA (paper 1 of 2). Proceedings of the International Conference Analytic-based Modeling of Groundwater Flow, Nunspeet, The Netherlands, pp. 59–75. Strack, O.D.L., 1986. Curvilinear analytic elements. Report on research done for the US Bureau of Mines, Strack Consulting, North Oaks, MN. Strack, O.D.L., 1989. Groundwater Mechanics, Prentice-Hall, Englewood Cliffs, NJ. Strack, O.D.L., 1995. A Dupuit–Forchheimer model for threedimensional flow with variable density. Water Resources Research 31 (12), 1309–1315. Strack, O.D.L., Jankovic´, I., 1999. A multi-quadratic area-sink for analytic element modeling of groundwater flow. Journal of Hydrology 226, 188–196. Strack, O.D.L., Jankovic´, I., Barnes, R., 1999. The superblock approach for the analytic element method. Journal of Hydrology 226, 179–187. Strack, O.D.L., Haitjema, H.M., 1981. Modeling double aquifer flow using a comprehensive potential and distributed singularities. 2. Solutions for inhomogeneous permeabilities. Water Resources Research 17 (5), 1551–1560. Zaadnoordijk, W.J., 1988. Analytic elements for transient groundwater flow. PhD Thesis, Department of Civil Engineering, University of Minnesota, Minnesota. Zaadnoordijk, W.J., Strack, O.D.L., 1993. Area-sinks in the analytic element method for transient groundwater flow. Water Resources Research 29 (12), 4121–4129.
Journal of Hydrology 226 (1999) 128–138 www.elsevier.com/locate/jhydrol
Principles of the analytic element method O.D.L. Strack* University of Minnesota, Department of Civil Engineering, Minneapolis, MN 55455-0220, USA Received 5 September 1997; accepted 7 October 1998
Abstract An overview of the development and current state of analytic element method, as well as a brief presentation of the principles on which the technique is based is presented in this article. Illustrations of a few representative analytic elements are also presented. A brief discussion is added of both rotational flow and transient flow. The latter two aspects of the method require further development before they are ready for implementation in computer programs and can be applied to the modeling of groundwater flow. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Analytic element method; Groundwater modeling; Mathematical modeling
1. Introduction The analytic element method was developed for the mathematical modeling of the flow of groundwater and was originally intended for problems of regional flow. The analytic element method is a numerical method that gives an approximate solution to a problem. Analytic elements are mathematical functions that differ from the many existing classical solutions in that they are not restricted to a single boundary value problem, but rather possess degrees of freedom that makes it possible to combine them. The method is based upon the superposition of analytic functions and might be considered as an extension of classical models of regional flow based on the superposition of elementary functions, for example the Thiem and Theiss solutions for wells, and the potential for uniform flow. Characteristic for such classical models is that each function is unaffected by the other functions in the solution. The discharge * Fax: 11-612-626-7750. E-mail address: [email protected] (O.D.L. Strack).
of each well in a field of wells, for example, is known a priori and is not affected by the discharges of the other ones. The analytic element method presented as a modeling approach by Strack (1989), differs from the classical analytical models in that the elementary solutions possess degrees of freedom that may be chosen in such a way that boundary conditions are met along certain internal boundaries. The technique resembles classical boundary element methods, and differs from these in the nature of the analytic solutions chosen. In boundary element methods, these solutions are boundary integrals of the Cauchy type (double or single layers), the boundaries are usually closed, the integrals are often evaluated numerically, and the formulation used to compute the unknown coefficients is usually different from that used to compute values in the field once the coefficients have been determined. The solutions used in the analytic element method are not necessarily Cauchy integrals; elements may be constructed, for example, by generalizing closed-form analytic solutions obtained by the use of conformal mapping techniques. Examples of such elements are the so-called double-root elements
0022-1694/99/$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S0022-1694(99 )00 144-4
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
and the elements for circular and elliptical inhomogeneities presented by Strack (1989), the analytic elements useful for modeling flow near corners formed by internal boundaries presented by Detournay (1990), the three-dimensional analytic elements for partially penetrating wells presented by Haitjema and Kraemer (1988), and the three-dimensional solutions for inhomogeneities bounded by ellipsoids presented by Fitts (1990). Haitjema (1982) introduced an approach for the regional modeling of the threedimensional flow in confined aquifers based on a combination of two-dimensional and three-dimensional solutions. A major difference between the analytic element method and classical boundary element methods is the use of area-sinks. Area-sinks are functions that represent extraction (or infiltration) over areas bounded by polygons, and are useful for modeling quasi three-dimensional flow. Such flows occur in shallow confined and unconfined aquifers where, at least in part of the domain, the Dupuit–Forchheimer assumption has been adopted. The area-sink, as presented by Strack (1989), is an analytic element that simulates uniform infiltration or extraction over an area that is bounded by a polygon. The use of area-sinks makes it possible to simulate leakage between aquifers, for example, by subdividing the domain into polygons, each with its own uniform rate of leakage. Area-sinks of this kind form the mainstay of the Dutch National Groundwater Model as presented by De Lange (1991) that was constructed by the use of the Multi-Layer Analytic Element Model (MLAEM). The individual leakages of the area-sinks are computed from the difference in head between the aquifers above and below the element. A slightly more sophisticated area-sink, bounded by a triangle and with a leakage that varies linearly, was also presented by Strack (1989) and was used by Haitjema and Strack (1985) to produce a crude model of transient flow. The triangular areasinks were used in the latter model to simulate the release of water from storage. Strack and Jankovic´ (1999) developed the potential for an area-sink that simulates extraction or infiltration that varies over the domain enclosed by a polygon as a radial-basis interpolator (Hardy, 1971). Only available resources, such as computer memory, limit the number of degrees of freedom of this interpolator. The development of this
129
new area-sink makes possible a broader definition of the analytic element method, presented here for the first time. This enables us, in theory, to determine approximate closed-form expressions for the components of general vector fields, i.e. vector fields that have both a non-zero divergence and a non-zero curl. We may use these expressions to apply the analytic element method to aquifers with continuously variable properties, such as hydraulic conductivity, base elevation, and thickness. Another development, based on the combination of the area-sinks with the technique of separation of variables, makes it feasible to create accurate and efficient analytic element models of transient flow. The combination of these advances with the technique of overspecification introduced by Jankovic´ and Barnes (1999a), and the superblock approach discussed briefly at the end of this paper and presented in Strack et al. (1999), will make it possible to develop the analytic element method into an accurate and efficient method for the simulation of general cases of groundwater flow. We describe the analytic element method briefly in what follows. We do this for the two-dimensional steady flow, for the sake of brevity, but the technique is not restricted to two dimensions. A few comments will also be made with regard to transient flow, overspecification, and the superblock approach. It is important to note that some of the recent advances in the analytic element method exist at present only in theory or have been implemented only in part.
2. General approach We formulate the analytic element method in terms of the discharge vector which represents the total flow per unit width that occurs over the saturated thickness of the aquifer. If the flow is irrotational, then the discharge vector may be expressed as the gradient of a scalar potential. The piezometric head may be determined directly from this potential. This is not possible if the flow is rotational; in that case a potential does not exist and the piezometric head must be determined after the discharge vector has been obtained as a function of position. We restrict this discussion to problems of regional
130
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
flow in systems of confined, unconfined, or semiconfined aquifers. The flow in each aquifer is described in terms of its own two-dimensional discharge vector. 2.1. Irrotational flow The discharge vector field is defined by suitable boundary conditions and expressions for the divergence and the curl of the field. We first consider the case of irrotational flow. In this case there exists a potential, F ; the discharge vector may be represented as minus the gradient of this potential, i.e. Qi 22i F
1
where 2i represents the derivatives with respect to the global Cartesian coordinates x1 and x2, or xi, i 1; 2; … i.e. 2i
2 2xi
2
We adopt the Dupuit–Forchheimer assumption, so that the potential can always be expressed in terms of the piezometric head f . The relation between head and potential may be either linear (e.g. for the cases of semi-confined and confined flow) or nonlinear (e.g. for the case of unconfined flow.) We list as examples the expression for the potential for shallow confined flow,
F kH f 2 12 kH 2
3
F
kf
2
1ij 2i Qj 0
5
and 2i Qi g
and the one for shallow unconfined flow, 1 2
layer. This expression for the divergence may be nonlinear in terms of the potential, for example in the case of leakage between an unconfined and a confined aquifer. The discharge vector is expressed in terms of the gradient of the potential rather than head, and therefore an iterative method of solution is required when the expression for the divergence of the discharge vector is non-linear in terms of potential. We approach the problem of determining the head and discharge vector throughout the aquifer in the following manner. We first estimate the divergence throughout the flow domain and represent it as a function of position. This function is taken to be piecewise smooth; the domain is subdivided into areas bounded by polygons and the divergence is represented as either a constant or an interpolator inside each polygon. The divergence jumps across the boundaries of the polygons and is zero in the area not covered by a polygon. The potential function for an area-sink where the leakage is represented by a radial basis interpolator (the multi-quadric interpolator, see Hardy, 1971) is presented by Strack and Jankovic´ (1999). Next we write the basic equations to be solved, adopting the Einstein summation convention. We represent the divergence of the discharge vector as g and the governing equations for Qi become
4
where k represents the hydraulic conductivity, H the thickness of the confined aquifer and f the piezometric head. The divergence of the discharge vector for any given aquifer equals the net extraction rate per unit area for each aquifer. This divergence may be estimated a priori in the cases of rainfall or leakage through the bottoms of rivers or lakes that are above the groundwater table. In many cases of practical interest, however, the divergence is not known a priori but may be expressed in terms of parameters of the problem. For the case of leakage between aquifers, for example, the divergence equals the difference in head across the leaky layer divided by the resistivity of that
6
where 1 ij is the alternating tensor of the second rank, and 1ij 2i Qj represents the curl of the discharge vector. Note that Eq. (5) is satisfied by representing the discharge vector Qi as the gradient of a potential as in Eq. (1). The analytic element method makes use of superposition of solutions; the discharge vector, Qi, is composed of a sum of functions. These functions may be combined in two groups, the first contains g all area-sinks, and may be represented as Q i ; the h
second is written as Q i ; and represents the sum of all remaining functions, which can be expressed as a sum of the gradients of harmonic potentials, whence the superscript h. The extraction rates of the areasinks have been chosen to match the known function
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
131
Fig. 1. Piezometric contours for an area-sink with infiltration rate shown in Fig. 2; maximum and minimum values are 49.970 and 50.085 m, respectively. The contour interval is 0.5 cm. g
g and thus the discharge vector Q i satisfies Eq. (6), i.e. g
2i Q i g
7
The total discharge is represented as the sum of two parts, g
h
Qi Q i 1 Q i
8
and in view of (7) the second part is the gradient of a harmonic function, h
h
Q i 22i F
9
where h
72 F 0
10
Thus, all functions other than the functions that represent the area-sinks are harmonic. The complete solution to the problem thus may be obtained by
superposition of the potentials for area-sinks and additional harmonic potentials that simulate the effect of a variety of internal boundaries upon the flow. Examples of the latter functions are the potentials for linesinks (sources or sinks distributed along lines or curves), line-doublets (double layers, representing a jump in potential along lines or curves) and line dipoles (distributed dipoles with their axes tangent to lines or curves, representing a jump in the stream function along the lines or curves). An efficient manner of representing such functions is by the use of complex variable techniques. The parameters in an analytic element model must be determined via an iterative process when the problem is non-linear. Many of the boundary conditions applied in regional models involve the piezometric head, and the expression for the divergence (often equal to leakage between aquifers) also is a function of the piezometric head in the aquifer. The discharge vector is the gradient of the discharge
132
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
Fig. 2. Contours of constant infiltration rate inside the polygon that bounds an area-sink; the maximum and minimum values plotted are 0:9 × 103 and 20:9 × 1023 m=day; respectively. The contour interval is 0.001 m/day.
potential. If one of the aquifers in the system is unconfined, the relation between the potential and head of that aquifer is non-linear and an iterative procedure must be followed in order to determine all parameters in the solution. The expression for the divergence g is updated at each step in the iteration until sufficient accuracy has been achieved, see Strack (1989). The theory as outlined above has been implemented and is used in computer models of groundwater flow, particularly in large-scale computer models such as those of The Netherlands (De Lange, 1991), of a future retention basin in The Netherlands (Moorman, 1997), and the Twin Cities Metropolitan Area (Seaberg et al., 1997). As an illustration, we present plots of a few representative analytic elements. The effect of an area-sink by itself is illustrated in Fig. 1, which shows the pattern of piezometric contours, produced by an area-sink with a given rate of extraction that varies over the area. The element is bounded by a polygon and resides in an infinite aquifer; there are no other features present to affect the flow. The aquifer is unconfined with a hydraulic conductivity of 10 m/day. The width of the area shown is 400 m and the head at a distance of 1000 m from the center of the area-sink is fixed at
50 m. The values of heads vary between a maximum of 50.085 and a minimum of 49.97 m. The contour interval is 0.5 cm. The maximum value occurs to the left of the center of the figure and the minimum value to the right. The distribution of the extraction rate over the area-sink is shown in Fig. 2. The extraction varies nearly linearly from a value of 210 22 m/day at the extreme left of the area-sink to a value of 10 22 m/day at the extreme right. Note that a rate of infiltration is modeled by means of a negative rate of extraction. The variation in infiltration rate from positive to negative values causes the maximum and the minimum values to be visible in the plot of piezometric contours as shown in Fig. 1. Doublet elements are used to produce a jump in the potential (see Strack and Haitjema,1981; Strack, 1989) that occurs across the boundary of an inhomogeneity in the aquifer properties. An example of the use of doublets to simulate a jump in aquifer properties is shown in Fig. 3. There is uniform infiltration of 0.1 m/day over an area that includes the inhomogeneity, which is bounded by a polygon. The width of the area shown is 400 m. The flow is unconfined everywhere; the base of the aquifer is at elevation 0 m, the hydraulic conductivity is 10 m/day, and the head at a
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
133
Fig. 3. Piezometric contours for flow in an aquifer with an inhomogeneity in aquifer properties; maximum and minimum values are 65.3 and 59.6 m, respectively, and the contour interval is 0.1 m.
reference point 10,000 m away from the inhomogeneity is 50 m. The aquifer is unconfined and the base inside the inhomogeneity is at an elevation of 50 m, i.e. well above the base outside the inhomogeneity. The hydraulic conductivity inside the inhomogeneity is 0.1 m/day, i.e. hundred times lower than that outside. The pattern of piezometric contours is shown in Fig. 3; the heads vary between a minimum of 59.5 m occurring outside the inhomogeneity and a maximum of 65.3 m, visible just to the left of the center of the polygon that binds the inhomogeneity. The contour interval is 0.1 m. Note that the contours are much closer together inside the inhomogeneity than outside as a result of the contrast in the hydraulic conductivity. The jump across the doublet elements is approximated by a polynomial of the order of 20; the number of degrees of freedom associated with each side is 20. The solution has been obtained using the
method of over-specification proposed by Jankovic´ and Barnes (1999b). The technique of over-specification relies on allocating a larger number of control points along the boundary segment than there are degrees of freedom in the doublet function used to model the jump. The unknown parameters associated with the jump function are determined by requiring that the sum of the squares of the errors at all control points is a minimum. The accuracy obtained by this approach is about an order of magnitude higher than obtained by the conventional technique of matching boundary conditions at collocation points that are equal in number to the degrees of freedom of the jump function. The degree of over-specification chosen for this problem was 2, i.e. there are 2 times as many control points as the degrees of freedom, i.e. 40 per side. The piezometric contours shown in Fig. 4 correspond
134
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
Fig. 4. Piezometric contours for the case of Fig. 3 under perched conditions; the maximum and minimum values are 58.7 and 33.1 m, respectively, and the contour interval is 0.5 m.
to the same problem, except that the head at the reference point is 5 m rather than 50 m. The values of heads on the contours vary between 33 m outside the inhomogeneity and a maximum of 59 m inside. The contour interval is 0.5 m. The water table inside the inhomogeneity becomes perched; the infiltration causes water to flow from the domain inside the inhomogeneity with its base at 50 m to the domain outside where the base is at 0 m and the head is well below the base inside the inhomogeneity. The head now jumps across the boundary; the boundary condition is that the head inside the inhomogeneity equals the base elevation wherever the head outside is below 50 m and that the head is continuous across the boundary wherever the head outside the inhomogeneity is larger than 50 m. This is illustrated in Fig. 5, which is a detail of Fig. 4. The values of the head along the contours near the boundary of the inhomogeneity
are indicated in the figure. The validity of the use of the Dupuit–Forchheimer solution for perched conditions was demonstrated by Fitts and Strack (1996) by presenting exact solutions for two-dimensional cases of unconfined flow in an aquifer with a stepped base. The flow net is shown in Fig. 6 for a string of leaky curvilinear elements. The curvilinear elements were developed by Strack (1986); the derivation is based on the principles presented in Strack (1989). The elements represent a leaky wall of a width of 1 m (the width of the window displayed in the figure is 400 m) and a hydraulic conductivity of 10 23 m/day; the hydraulic conductivity of the aquifer is 10 m/day. There is a far field of uniform flow from left to right in the figure of discharge 0.5 m 2/day in the aquifer and the flow is unconfined. The hydraulic conductivity is 10 m/day. Note that the flow net consists of equipotentials and streamlines; the potential is given for this
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
135
Fig. 5. Detail of Fig. 4.
case by Eq. (4). The values of the potential vary between 1:283 × 104 and 1:315 × 104 m3 =day: The stream function varies between 2120 and 80 m 3/ day. The contour interval for both functions is 10 m 3/day. Note that the interval in piezometric head between contours is not constant, because the relation between the head and the potential is non-linear. The values in head vary between 30.65 and 51.3 m over the region. The approach outlined above may, in theory, be extended to problems of three-dimensional flow. However, an efficient mathematical formulation for the three-dimensional equivalent of the area-sink, a volume-sink, has not yet been determined and threedimensional applications have so far been limited to specialized three-dimensional models of limited scope (e.g. Fitts, 1990), and to models that combine a two-dimensional formulation with three-dimensional analytic elements (Haitjema, 1982, 1985).
Other noteworthy applications of the approach are the combination of an analytic element model capable of simulating the interchange between surface water and groundwater (Haitjema and Mitchell-Bruker, 1996), and the implementation of variable density flow (Strack, 1995). It may be worth noting that the latter application does not cause the curl of the discharge vector to become rotational, provided that the Dupuit–Forchheimer assumption is adopted. 2.2. Rotational flow The analytic element method has not yet been applied to problems of rotational flow. It is possible to cover rotational flow, such as occurs in aquifers with a sloping base and in aquifers with continuously varying properties (rather than properties that are uniform in sub-domains). The approach is similar to that explained above for irrotational flow, and rests on Helmholtz’s theorem. This theorem states that any
136
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
Fig. 6. Flow net for a curvilinear leaky wall; the contour interval is 10 m 3/day.
vector field may be decomposed in a rotational but divergence-free field and an irrotational field with non-zero divergence. We therefore may add to an irrotational field, as described above, a rotational field that is used to generate the required rotation as a function of position. As for the divergence, we assume that the curl is given; the curl is, in fact, determined by some iterative process, and is fixed during each step of the process. The curl is modeled by means of an area-vortex with a vorticity distribution given by an interpolator function. We represent the curl as 2b , rewrite Eq. (5)
1ij 2i Qj 2b
11
and repeat Eq. (6) 2 i Qi g
12
Note that both equations are linear in terms of the discharge vector Qi as b and g are treated as given functions of position. Thus, superposition of
individual functions Qi remains possible. Non-linearities that may exist in the problem have not been eliminated, but have been hidden in the iterative process for determining the functions b and g . The discharge vector field is now separated into three b
g
h
parts, Q i ; Q i ; and Q i ; where the first vector field has rotation b and is divergence-free, i.e. b
1ij 2i Q j b
13
and b
2i Q i 0
14
The second vector field has a divergence of g and zero curl, and the third vector field has both zero divergence and zero curl. The latter two vector fields may be written in terms of the gradient of a potential, but it is important to note that this potential is no longer explicitly related to the piezometric head;
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
this relation is implicit and heads must be computed by some approximate method from the vector field. The b
vector field Q i is divergence-free and may therefore be written in terms of the curl of a scalar stream function, b
b
Q i 1ij 2j C
15
The curl may be applied to the solution by the use of area-vortices, which may be presented in terms of the sum of a stream function, equal to an interpolator inside a polygon and equal to zero outside, plus an harmonic potential that represents the boundary integrals required to create a vector field that is continuous across the boundary of the area-vortex. The area-vortex is analogous to the area-sink, and the functions that describe the flow generated by the area-vortex may be obtained directly from those of the area-sink with the same density distribution. The main advantage of the approach is that the majority of functions are harmonic, and thus complex variable techniques may be used to determine most of the analytic elements, whereas the functions that contribute to the divergence and curl of the vector field are simple and require minimal computational effort. There is a difficulty associated with the approach for rotational flow that has as yet not been resolved fully. The representation of the vector field as the sum of the gradient of a potential and the curl of a stream function is the cause of this complication. The discharge vector field and the piezometric head are related via Darcy’s law. The application of boundary values in terms of head will require that heads be computed from the discharge vector field, at least at isolated points, which requires integration, either numerical or analytical. The use of an iterative process will be necessary to arrive at a final solution, determining the unknown parameters in the expressions for all functions, including the area-vortices, until some criterion for accuracy has been met. It is not the objective of the current paper to present detailed possible approaches for obtaining solutions for cases of rotational flow, but rather to present the principle of the approach.
3. Transient flow A few comments will be made about transient flow.
137
An analytic element model for transient flow was introduced by Zaadnoordijk (1988) and Zaadnoordijk and Strack (1993), based on the integration of transient wells along boundaries. The disadvantage of that approach is that many parameters need to be remembered as the solution progresses through time, and that the computational effort is significant. An alternative to that approach is to use area-sinks to simulate the release of water from storage using area-sinks of linear density distribution, as applied by Haitjema and Strack (1985). The area-sinks with variable strength are an attractive alternative to the earlier elements that were capable only of constant or linearly varying storage over the element. 4. The superblock approach The superblock approach relies on the notion that the majority of analytic elements in a model are harmonic functions and thus may be represented as the real parts of analytic functions of a complex variable. We subdivide the aquifer into a number of square blocks and determine for each block which elements lie far enough away from the boundary of the block that they may be represented inside the block by a Taylor series. All elements outside the block may thus be combined into a single Taylor series expanded about the center of the block. Computation of the potential in each block thus consists of two parts. The first involves evaluation of all functions that are considered to be inside the block. The second part consists of evaluating the effect of all the elements outside the block, which are combined into a single series, thereby greatly reducing computation time. Superblocks may also be used in determining the unknowns in the solution. 5. Conclusion The analytic element method in its present form is applicable to steady flow in aquifers with properties that are piecewise constant. It has proven to be capable of dealing with large aquifer system, and in obtaining highly accurate solutions. Future developments are promising. One of these is the superblock approach, which is expected to make it possible to deal effectively with models on a
138
O.D.L. Strack / Journal of Hydrology 226 (1999) 128–138
nation-wide scale (e.g. the Dutch National Groundwater Model, NAGROM, De Lange, 1991) and a state-wide scale (e.g. the Twin Cities Metropolitan Model, Seaberg et al., 1997). Other possibilities for future development are the inclusion of rotational flow, applicable to cases of continuously variable aquifer properties, and of transient flow in the analytic element method. References De Lange, W.J., 1991. A groundwater model of The Netherlands. Note No. 90.066, National Institute for Inland Water Management and Waste Water Treatment, Lelystad, The Netherlands. Detournay, C., 1990. On a Cauchy integral element method for potential flow with corner singularities. BETECH 90, Computational Engineering with Boundary Elements. Fluid and Potential Problems, vol. 1. University of Delaware, pp. 119–130. Fitts, C.R., 1990. Modeling three-dimensional groundwater flow about ellipsoids of revolution using analytic functions. PhD Thesis, Department of Civil Engineering, University of Minnesota, Minnesota. Fitts, C.R., Strack, O.D.L., 1996. Analytic solutions for unconfined groundwater flow over a stepped base. Journal of Hydrology 177, 65–76. Haitjema, H.M., 1982. Modeling three-dimensional flow in confined aquifers using distributed singularities. PhD Thesis, Department of Civil Engineering, University of Minnesota, Minnesota. Haitjema, H.M., 1985. Modeling three-dimensional flow in confined aquifers by superposition of both two- and threedimensional analytic functions. Water Resources Research 21 (10), 1557–1566. Haitjema, H.M., Kreamer, S.R., 1988. A new analytic function for modeling partially penetrating wells. Water Resources Research 24 (5), 683–690. Haitjema, H.M., Mitchell-Bruker, S., 1996. Modeling steady-state conjunctive groundwater and surface water flow with analytic elements. Water Resources Research 32 (9), 2725–2732. Haitjema, H.M., Strack, O.D.L., 1985. An initial study on thermal energy storage in unconfined aquifers. Final Report to Battelle Pacific Northwest Laboratory.
Hardy, R.L., 1971. Multiquadric equations of topography and other irregular surfaces. Journal of Geophysical Reseach 62 (2), 733– 737. Jankovic´, I., Barnes, R., 1999a. Three-dimensional flow through large numbers of spheroidal inhomogeneities. Journal of Hydrology 226, 224–233. Jankovic´, I., Barnes, R., 1999b. High order line elements for twodimensional groundwater flow. Journal of Hydrology 226, 211– 223. Moorman, J.H.N., 1997. Analytic Element Model analysis of the influence of different scenarios for the water level in a future retentions basin, on local and regional hydrology. Proceedings of the International Conference Analytic-based Modeling of Groundwater Flow, Nunspeet, The Netherlands, pp. 105–111. Seaberg, J.K., Hansen, D.D., Block, B.W., Streitz, A.R., de Lange, W.J., Bakker, M., 1997. Development of a regional groundwater model for the Twin Cities Metropolitan Area, USA (paper 1 of 2). Proceedings of the International Conference Analytic-based Modeling of Groundwater Flow, Nunspeet, The Netherlands, pp. 59–75. Strack, O.D.L., 1986. Curvilinear analytic elements. Report on research done for the US Bureau of Mines, Strack Consulting, North Oaks, MN. Strack, O.D.L., 1989. Groundwater Mechanics, Prentice-Hall, Englewood Cliffs, NJ. Strack, O.D.L., 1995. A Dupuit–Forchheimer model for threedimensional flow with variable density. Water Resources Research 31 (12), 1309–1315. Strack, O.D.L., Jankovic´, I., 1999. A multi-quadratic area-sink for analytic element modeling of groundwater flow. Journal of Hydrology 226, 188–196. Strack, O.D.L., Jankovic´, I., Barnes, R., 1999. The superblock approach for the analytic element method. Journal of Hydrology 226, 179–187. Strack, O.D.L., Haitjema, H.M., 1981. Modeling double aquifer flow using a comprehensive potential and distributed singularities. 2. Solutions for inhomogeneous permeabilities. Water Resources Research 17 (5), 1551–1560. Zaadnoordijk, W.J., 1988. Analytic elements for transient groundwater flow. PhD Thesis, Department of Civil Engineering, University of Minnesota, Minnesota. Zaadnoordijk, W.J., Strack, O.D.L., 1993. Area-sinks in the analytic element method for transient groundwater flow. Water Resources Research 29 (12), 4121–4129.