## Abstract

Steady two-dimensional groundwater flow in a porous low-permeable trough is studied by the method of boundary-value problems for holomorphic functions. In the overlying highly permeable aquifer the hydraulic head varies linearly, i.e. flow is unidirectional. The exchange of groundwater between aquifer and trough does not affect the flow in the aquifer. It is assumed that through a horizontal aquifer–trough interface the head is transmitted into the trough, where the bounding effect of the bed causes circulatory seepage. Triangular troughs are studied and an isosceles form with a base angle of 38° is proved to have the highest circulation rate at a given cross-sectional area. In the class of arbitrary forms, solution to this optimal shape design problem is obtained by tackling the Schwartz and Signorini singular integrals and a unique and global maximum of the rate is found. The extreme curve coincides with an optimal soil channel of minimal seepage losses having depth to width ratio of 0.371. The global maximum differs from that one for the triangular class in 3% only that corroborates stability and robustness of the optimization criterion and makes possible isoperimetric estimates of the seepage intensity through arbitrary troughs.

## 1. Introduction

In hydrogeology, a standard conceptualization of aquifers and aquitards as hydrostratigraphic units assumes groundwater flow along the layers in the former and across the layers (leakage) in the latter (Freeze & Cherry 1979). This quasi-one-dimensional model is called hydraulic or Dupuit–Forchheimer (Polubarinova-Kochina 1962). The model works well in laminated formations, i.e. when the boundaries between layers of contrasting conductivity are nearly parallel (not necessary horizontal). Quite often, however, the confining layers of aquifers have disconformities or breaches appearing either as natural imperfections or caused by mining/petroleum activities (Brikowski 1993; Anderson 2003). Groundwater flow becomes two-, three-dimensional in the vicinity of these flaws (Strack 1989, pp. 393–396). We consider one possible type of this complication, a trough BMC in the underlying bedrock (figure 1*a*), which Freeze & Cherry (1979, fig. 4.3) called a ‘bedrock valley aquifer’. Such bedrock depressions are typical in coastal aquifers of the Batinah region in Oman, where the trough depths reach 500–600 m as compared with points B and C (Eastern Batinah Resource Assessment 1995). In the trough, the premises of the Dupuit–Forchhemer theory do not hold because seepage there is two-dimensional although commensurate with one-dimensional flow in the overlying homogeneous isotropic aquifer ABCD.

With the increase of the depth of burial, porosity of the rock decreases due to consolidation and cementing. This reduction of the interstitial space makes the trough massif in figure 1*a* less permeable to groundwater than the aquifer.

The basic assumption in this paper is that the flow in the trough does not affect the flow in the aquifer, and that this amounts to the assumption that the flow through the trough, *Q*, is small relative to the flow rate in the aquifer. Figure 1*b* depicts the streamlines of steady groundwater movement, which is induced in the trough by the ambient flow in the aquifer.

Although this movement is relatively minor in rate and slow as compared with the main flow in the aquifer, the exchange of water between the trough and aquifer can cause contamination of the aquifer. Indeed, deep groundwater in the trough often has much higher salinity (EC of 80 000–90 000 μS cm^{−1} is common in Omani trough marls as compared with 500–1000 μS cm^{−1} in the aquifer). Groundwater samples in deep trough-screened wells (Eastern Batinah Resource Assessment 1995) showed distinctly different chemical constitution with prevailing Cl^{−} anions as compared with ones in shallow (near-wadi and interfluve) wells tapping the aquifer. This indicates higher residence time of groundwater in the trough as compared with the aquifer. Topologically, the ‘eddy’ shown in figure 1*b* departs from a horizontal route of bulk groundwater flow in the aquifer, deliquesces minerals in the trough, advects them back to the aquifer and, hence, intensifies solute migration. It causes dispersion that, if no circulation in a delitescent trough, would be controlled by slow diffusion only.

Mathematical modelling undertaken in this paper becomes a cheap supplement to field geophysical methods (Eastern Batinah Resource Assessment 1995). As we shall show, modelling gives a simple hydrodynamic clue to the pattern of water circulation through the aquifer–trough interface.

We study analytically groundwater movement in the trough considering it as a Toth-type element (Freeze & Cherry 1979), where the hydraulic head is a harmonic function. In congruity with the Toth model, we assume that the boundary condition on the interface BOC between the aquifer and trough is imposed by the ambient (aquifer) flow. This is not an ideal simplification even for the Toth elements recharged from an isobaric boundary (Kacimov 1996) and the refraction boundary condition (Obnosov 1998) would be more adequate. However, the advantage of the Toth model stems from its simplicity and versatility of shaping the contour BMC. From the model, the complex potential in the element, in particular, the leakage intensity along the interface is reconstructed for an arbitrary trough geometry. We further use the trough potential as the feedback to the aquifer Dupuit–Forcheheimer model with a distributed losses at BOC.

We solve explicitly and rigorously the problem of trough shaping by considering the trough bed as a control variable in an optimal shape design problem (Pironneau 1984). The cross-sectional area of the trough and ambient flow in the aquifer are given and the trough of maximal *Q* is searched. Two boundary-value problems for the complex position and complex potential are solved and concatenated by the help of the Schwartz and Signorini formula. A unique and global maximum of *Q* is attained at a fascinatingly simple curve of an optimal irrigation channel found earlier by Ilyinsky & Kacimov (1984). This closed-form maximum in the class of arbitrary trough boundaries is compared with one obtained in the class of triangular troughs and is shown to differ in only 3% in terms of the criterion of optimization and less than 5% in the trough depth–width ratio.

## 2. Toth triangular element in a confined aquifer

We consider steady, fully saturated, Darcian flow in a confined aquifer of a constant thickness *b* and conductivity *K* (figure 1*a*). A trough BMC is filled with a porous medium of hydraulic conductivity *K*_{T}, *K*_{T}≪*K*. We assume that the trough is an isosceles ‘triangle’ of a width 2*L* and base angle *απ*, 0≤*α*≤1. If 0≤*α*≤1/2 (figure 1*a*,*b*) the ‘triangle’ is finite and common. If 1/2≤*α*≤1 point *M* is at infinity (figure 1*c*). The case *α*=1 (figure 1*d*) represents a breach in an impermeable layer separating two aquifers (Strack 1989, fig. 5.21).

The interface BOC in figure 1*a* is parallel to the aquifer top and bottom and we set the origin of Cartesian coordinates *x*O*y* at point O.

The incident groundwater flow in the aquifer is uniform and we assume that the aquifer specific discharge far from the trough (at points *A*_{1} and *D*_{1}) has only a horizontal component *u*=*Kϵ* where *ϵ*=const.>0 is the hydraulic gradient magnitude. The hydraulic head *h*(*x*) in the aquifer and *h*_{t}(*x*, *y*) in the trough are measured from the line MON, which is a constant head boundary for both media due to the symmetry of the flow pattern. If *K*_{T}→0 then *h*=*h*_{0}=−*ϵx*. As we have already assumed, the trough does not affect aquifer even for *K*_{T}≠0. As BMC is impervious the trough acts as a half-eddy as depicted in figure 1*b*, i.e. water seeps into the trough through BO and exudes back into the aquifer through OC. The specific discharge in the trough has both components. Point O in figure 1*b* is a hinge point at which the vertical component *v*_{T} of is zero. A longitudinal pressure drop in the aquifer induces this circulation, which increases with *K*_{T} and *ϵ*.

The questions are: what is the seepage rate *Q*? Which trough generates maximal *Q* at a given cross-sectional area of BMC?

In order to answer these questions we consider a two-dimensional flow in BMCOB induced by a linear variation of *h*_{T}=*h* along BOC. We introduce a complex physical coordinate *z*=*x*+i*y* and the complex potential *w*_{T}=*ϕ*_{T}+i*ψ*_{T} where *ϕ*_{T}=−*K*_{T}*h*_{T} is the seepage velocity potential and *ψ*_{T} is the stream function. The specific discharge is . Both *ϕ*_{T}(*x*, *y*) and *ψ*_{T}(*x*, *y*) are harmonic in the ‘triangle’. Along BMC, we set *ψ*_{T}=0 and the complex potential domain is sketched in figure 2*a* with the shape of BOC unknown. The stream function reaches its maximum, *Q*, at point O.

Along BOC the horizontal component of velocity *u*_{T}=*ϵ* because we assumed *ϕ*_{T}=*K*_{T}*ϵx*_{T} due to pressure continuity along the leaky interface BOC. Consequently, the hodograph domain *V*_{T}=*u*_{T}+i*v*_{T} is also a common triangle shown in figure 2*b* for *α*<1/2, corresponding to figure 1*b*, and an infinite ‘triangle’ in figure 2*c* for *α*≥1/2, corresponding to figure 1*c*,*d*. We can solve the flow problem by a conformal mapping of a mirrored hodograph (domain in the *u*_{T}−i*v*_{T} plain), which is a holomorphic function, onto the physical domain ‘triangle’. However, we implement an alternative, boundary-value problem technique (Kacimov 1996; Kacimov & Youngs 2005), bearing in mind our ultimate goal to study arbitrary troughs for which the hodograph domain is generally not known *a priori*.

We map conformally BMC onto the half-plane Im *ζ*≥0 of an auxiliary plane *ζ*=*ξ*+i*η* (figure 3*a*) by the function(2.1)i.e. with the correspondence of points B, C→±1, M→∞.

From equation (2.1) at *ζ*=*ξ*, |*ξ*|≤1 (section BC),(2.2)where *Γ* is the Gamma-function and _{2}*F*_{1} is the hypergeometric function.

Now in the half-plane of figure 3*a*, we determine a holomorphic function *w*_{T} by the following boundary conditions:(2.3)By the Signorini formula (see Kacimov & Youngs 2005) the solution to the mixed boundary-value problem (2.3) is(2.4)where *x*[*τ*] in the singular integral is taken from equation (2.2) and the branch of positive at *ζ*≥1 is selected.

In the limit *ζ*→*ξ*, |*ξ*|≤1 from equation (2.4) we obtain the dimensionless stream function *Ψ*_{T}=*ψ*_{T}/(*K*_{T}*ϵL*) along BOC,(2.5)From equations (2.2) and (2.5) using Mathematica (Wolfram 1991) routines *ParametricPlot* and *CauchyPrincipalValue* we evaluated *ψ*_{T}(*X*) functions, which are shown in figure 4 for *α*=0.2, 0.4, 0.99 (curves 1–3, correspondingly).

As we have mentioned, *Q*=*ψ*_{T}(0) and figure 5*a* presents a dimensionless flow rate *Q*/(*K*_{T}*ϵL*) as a function of *α*, 0≤*α*≤1. If we use the data for Wadi Al-Maawel (Eastern Batinah Resource Assessment 1995) with *L*≈15 000 m, *α*≈0.01, *K*_{T}=1 m day^{−1} assessed from geophysical surveys and *ϵ*=0.075 from a network of monitoring wells then *Q*≈1.8 m^{2} day^{−1}. From the piezometric surface in this catchment, we calculated that the flow rate in the aquifer is an order of magnitude higher, i.e. our main assumption on relatively low *Q* is valid. We note also that *Q*→0 as *α*→0, i.e. the trough vanishes.

Another dimensionless form of *Q* is shown in figure 5*b* as , where *S*=*L*^{2} tan *απ* is the trough cross-sectional area. Clearly, this form is meaningful at *α*≤0.5 only because at *α*≥1/2 the area *S*=∞ but *Q* is finite. The curve *μ*(*α*) reveals a unique and global peak *μ*_{m}≈0.5057 at *α*_{m}≈0.21, i.e. *απ*≈38°. It implies that in the class of isosceles triangles of a given area the trough of an angle *α*_{m} provides a maximum of *Q*(*α*).

Arbitrary ‘triangles’ of base length 2*L* and base angles *βπ* (point B) and *γπ* (point C) were studied by the same technique. For non-isosceles ‘triangles’, we shift the origin of *xy* coordinates to point C. The only difference with the isosceles subclass is that, from the Schwartz–Christoffel formula, the kernel in (2.5) isUsing *ContourPlot* of Mathematica we plotted *μ* as a function of *β* and *γ* (figure 5*c*). A unique and global maximum of *μ*(*β*, *γ*) is attained at *β*=*γ*=*α*_{m}=0.21, i.e. in the class of arbitrary triangles an isosceles one is the best in the sense of our isoperimetric pair *Q*−*S*. If *β* is fixed, then a local maximum of *μ* as a function of 0<*γ*<1−*β* is also unique and can be readily found.

We can readily solve the Signorini problem for trapezoidal, semi-elliptical or any other troughs, for which an explicit conformal mapping analogous to equation (2.1) can be used in the kernels of the integral representation (2.4). Then *Q* can be maximized by varying the geometrical parameters (side slope and depth to width ratio for a trapezium or axes ratio for an ellipse).

The trough in figure 1*a* is located at the aquifer bottom as is common in real situations. However, a leaky bulge can be placed at the top of the aquifer (figure 1*e*). Then, prior to implementing of our results, we have to check that flow remains confined in BMC. If pressure head *p*_{T}=*h*_{T}−*y* is not high enough then a phreatic surface *F*_{1}*F*_{2} can develop in the bulge and one arrives at seepage problems studied in geotechnics (Pavlovskii 1922) for the so-called Senkov dams.

## 3. General optimum

In this section, we solve the following optimization problem in the class of arbitrary troughs shown in figure 6*a*.

*Problem*. At a given trough area *S* determine the shape BMC providing maximum *Q*.

We assume that *K*_{T} and *ϵ* are fixed. We implement the same technique as in Kacimov (2001*a*) with an important difference that now none of canonical (e.g. complex potential or Zhukovskii function) domains are known. So, we shall vary BMC directly in the auxiliary half-plane Im *ζ*≥0 shown in figure 3*b* and solve a Dirichlet and mixed boundary-value problems there. We note that the correspondence of points in figure 3*b* is C, B→±1, O→∞, M→0 that is different from figure 3*a* for a ‘triangular’ trough.

First, the control of the shape of BMC is introduced by(3.1)where *U*_{n}=sin [*n* arccos (*ξ*)] are the Chebyshev polynomials of the second kind, *f*(*ξ*)is an arbitrary symmetric (*y*(−*ξ*)=*y*(*ξ*)) function belonging to the Holder class with the only constraint that *y*(±1)=0 that ensures that the trough tips lay on the bedrock. The coefficients *b*_{n} are to be found from optimization. Below we drop the limits of summation. We note that neither the trough width 2*L* nor depth *H* are fixed.

Since *y*=0 at BOC, i.e. at |*ξ*|≥1 in figure 3*b*, using equation (3.1) as a boundary condition for a holomorphic function *z*(*ζ*) we can write an integral representation,(3.2)where *x*_{O}=0 because our trough is symmetrical.

From equation (3.2),(3.3)where *T*_{n}=cos [*n* arccos *ξ*] and(3.4)Similarly to Kacimov (2001*a*), we use equations (3.1) and (3.3) and orthogonality of the Chebyshev polynomials. Then the trough area is(3.5)i.e. the constraint in the *problem* is a quadratic form of the Fourier coefficients of the control function.

In order to express the criterion through the same coefficients we search in the half-plane of figure 3*b* for a holomorphic function , which satisfies the following boundary conditions:(3.6)where *x*(*ξ*) is taken from equation (3.4). Similarly to (2.4) the Signorini formula yields(3.7)In the limit *ζ*→∞ we have *w*^{*}→−*Q* and from equation (3.7),(3.8)Now we assemble a Lagrange function *Φ*=*Q*+*λS*:(3.9)where *λ* is the Lagrange multiplier to be found. The necessary extremum condition is ∂*Φ*/∂*b*_{n}=0 for all *b*_{n} which by differentiating equation (3.9) gives(3.10)Substituting coefficients (3.10) into the constraint (3.5) we determine the multiplier(3.11)where *ζ*[ ] is the Riemann zeta-function. We note that in the quadratic equation for *λ* we retained the positive root (3.11) only because the negative *λ* gives positive *b*_{2k−1} and negative *Q* that does not make sense physically.

Now we put *λ* from equation (3.11) into (3.10) and (3.8) and arrive at(3.12)Due to for all at *m*≠*K* and , we infer that the second variation of the Lagrange functional is strictly positive, i.e. the sufficient extremum condition holds, that means we have found a unique and global maximum of *Q*. Using equation (3.12) in (3.1) and (3.3) and converting the series using eqns (5.4.6.9) and (5.4.6.10) from Prudnikov *et al*. (1986), we obtain parametric equations of the optimal BMC as(3.13)where *θ*=arccos *ξ*, 0≤*θ*≤*π*. Curve (3.13) is shown in figure 6*b*. Amazingly, equation (3.13) is exactly what Ilyinsky & Kacimov (1984) found (see their eqn (1.8)) as an optimal soil channel possessing minimal seepage losses at a given cross-sectional area.

By inserting the coefficients from equation (3.12) into the Plemelj–Sokhotsky limit of equation (3.7), we obtain after some algebra,(3.14)Unfortunately, we could not convert the last equation to an explicit form as in Ilyinsky & Kacimov (1984). Using equation (3.14) *ψ*_{T}/*Q* was plotted in figure 6*c* (the right half of the distribution) as a function of *x*/*L* by Mathematica routines.

From equation (3.12) *μ*=0.5209, i.e. the global maximum is only 3% better than one we found in the class of triangular shapes. From equation (3.13) the ratio *H*/(2*L*)≈0.371 as compared with 0.388 for the optimal triangle. The same robustness of optima was found by Kacimov (1985, 1992) for soil channels. In a hydrogeological context this is even more important than in irrigation engineering. Indeed, the trough contour is hardly known precisely although *H*/(2*L*) and *S* can be assessed from lithological profiles. Then *Q* for an arbitrary trough can be bounded from above by the global optimum and can be approximated by a geometrically close triangle.

The admissible class of curves BMC in figure 6*a* is not limited to monotonic *x*(*y*) functions or simple ones (although hydrogeologically bedrock shapes can hardly make fancy loops). If we consider an adjoint problem, i.e. search for a trough of minimal *S* at a given *Q*, then disconnected contours such as the ones in figure 1*c* are included directly into the admissible class. As Kacimov (1991) showed in a mathematically similar optimal shape design problem for trenches, the best curve can even be self-intersecting which signifies practically (but not mathematically!) an ill-posed optimization problem.

## 4. Unconfined aquifer

In this section, we consider a Dupuit–Forchheimer phreatic surface flow in the aquifer where the hydraulic head coincides with the thickness *h*(*x*) of the saturated layer above O*x* and, correspondingly, *h*=*H*_{0} at AB and *h*=*H*_{1} at CD as shown in figure 7. Again, we assume here that the trough does not influence the aquifer. A phreatic surface AD in the aquifer is the Dupuit parabola and the Dirichlet boundary condition along BOC of our ‘triangular’ trough is(4.1)where (*H*_{0L}, *H*_{1L})=(*H*_{0}, *H*_{1})/*L*.

In the auxiliary plane of figure 3*a* we set *ϕ*_{T}(*ξ*, 0)=−*Kh*_{T}[*x*(*ξ*)], |*ξ*|≤1, where *h*[*x*] is taken from the parabolic function of equation (4.1). Dropping the final Signorini integral solution we present the results of calculations as graphs *Ψ*_{T}(*X*) along BOC where *Ψ*_{T}=*ψ*_{T}/(*K*_{T}*L*). Figure 8 shows these functions for *α*=0.3, *H*_{0L}=1 and *H*_{1L}=0.9, 0.5, 0.1 (curves 1–3, correspondingly).

Clearly, the nonlinear variation of the potential on BOC generates asymmetric leakage as compared with the situation in figure 1*b*, i.e. the hinge point *G* is now displaced with respect to the geometrical centre O of the trough. The flow rate *Q* and the abscissa *x*_{G} of G were found by the *FindMinimum* routine from the maximum of *Ψ*_{T}(*x*) along BOC. Figure 9*a* shows *Q*/(*K*_{T}*L*) as a function of *α* for *H*_{0L}=1 and *H*_{1L}=0.1, 0.5 (curves 1–2) and figure 9*b* illustrates *μ*=*Q*/(*K*_{T}*S*) for the same *H*_{1L} revealing similar extrema as for confined flows in the aquifer. One can intuit that the optimal shape design problem has a unique global extremum for an arbitrary monotonic potential distribution along BOC. In the studied case of triangular contours, the maximum of *Q* is attained at *β*≠*γ*. Figure 9*c* illustrates *X*_{G}(*α*). It is clear from figure 9*c* that flow in flatter troughs is more skewed than in deep ones.

We attempted to determine a global maximum of *Q* at a given *S* similar to the confined flow scenario but failed to reveal a closed form solution, although numerical optimization for truncated *b*_{K} series is tenable. Remembering the class-robustness of the extrema in confined flows we expect that it holds for phreatic flows as well and, therefore, an optimal triangle here is close to the unknown global one.

## 5. Trough feedback

In this section, we assay the trough effect on the aquifer flow (neglected in previous sections). We consider a confined flow in a rectangular element ABCD (figure 1*a*) where the piezometric surface curves into convex and concave segments in response to non-uniformly distributed losses/gains into/from a ‘triangular’ trough. We shall check how the head *h*_{0}=−*ϵx* is altered by *h*_{T}(*x*, *y*).

Space-distributed leakage *E*(*x*) through the segment BOC is modelled in terms of the Dupuit–Forchheimer approximation (head is constant along any vertical line in ABCD), i.e.(5.1)where *E*>0 generates head losses along BO and *E*<0 reflects head recovery along OC. We assume that in equation (5.1) *E*(*x*)=*ψ*(*x*)/*L*, where *ψ*(*x*) is taken from the Toth model according to equation (2.4). Since *ψ* and *x* along BOC are expressed through an auxiliary variable in figure 3*a*, we transform equation (5.1) to *ξ* as(5.2)The nonlinear ordinary differential equation (5.2) was solved with the boundary conditions *h*(±1)=±*ϵ* using *NDSolve* of Mathematica and a shooting algorithm. Computations showed negligible alteration of *h*(*x*) even for high *ϵ*. We can take the solution of equation (5.2) and impose it as a new Dirichlet boundary condition in the Toth trough element and carry on the next iteration by equation (2.4) with a modified Cauchy-integral kernel *x*(*τ*).

## 6. Conclusion and discussion

By conformal mappings and the Signorini formula, we solved the seepage flow problem in a triangular trough, where the Darcian saturated continuum is agitated by a unidirectional flow in an overlying highly permeable aquifer. Along a straight line of contact between the aquifer and trough the hydraulic head drops linearly that generates a hinging topology of streamlines in the trough.

We solved an optimal shape design problem in the class of triangular isosceles (non-isosceles) troughs with the flow rate as a criterion, the trough cross-sectional area as a constraint and the base angle (angles) as control variable. The same problem was solved in the class of arbitrary trough shapes but the control function, used as a kernel of singular integrals, involved an infinite number of series coefficients. These coefficients were found explicitly because the Lagrange objective functional is a quadratic form of the coefficients.

The explicitly determined shape coincides with the optimum found by Ilyinsky & Kacimov (1984) for a free boundary problem in irrigation engineering. This discovery corroborates the Maxwell statement: ‘…Geometry itself is a part of the science of motion, and that it treats, not of the relation between figures already existing in space, but of the process by which these figures are generated by the motion of a point or line.’ Optimal forms in general and the studied hydrogeological unit in particular do fit this category of figures inherently—although often subtly—related to the creeping movement of water in porous media, which in the case of optimal regimes is generated by simple hydrodynamic singularities (Kacimov 2001*a*,*b*, submitted).

In geological systems, the trough shape cannot be changed by our own will but our optimal shapes lead to interesting engineering, philosophical and even tralatitious analogies. For example, we found that an optimal triangular trough is symmetrical under uniform aquifer flow conditions. If flow is non-uniform, i.e. water accelerates in the aquifer—the simplest example studied in this paper was a phreatic surface Dupuit flow—then the triangle of maximal seepage rate at a given area is not symmetrical any more. It illustrates explicitly the suppleness of the trough to ambient agitation in the framework of the constructal law, which states (Bejan 2005): ‘For a finite-size flow system to persist in time (to live), it must evolve in such a way that it provides easier access to the currents that flow through it.’

## Acknowledgments

This study was supported by Sultan Qaboos University, project IG/AGR/SOIL/02/04. Comments and criticism of two anonymous referees and Yu. V. Obnosov are highly appreciated.

## Footnotes

- Received April 26, 2005.
- Accepted November 17, 2005.

- © 2006 The Royal Society