## Abstract

The pioneering solution of Zhukovskii for a steady two-dimensional flow of an ideal heavy fluid with a nonlinear free boundary condition is extended to a Darcian flow of groundwater encumbered by an impermeable barrier. The stoss or/and lee sides of the barrier are covered by a macrovolume of a liquid contaminant. Explicit parametric equations of the sharp interface are obtained by inversion of the hodograph domain. Zhukovskii's gas-finger shape is shown to be a particular case of our new class of free surfaces. For a cap of a light liquid, partially covering the roof, from the given cross-sectional area of the cap, the affixes of the conformal mapping are found as a solution of a system of two nonlinear equations. The horizontal width and vertical height of the cap are determined. If the dimensionless incident velocity is higher than the density contrast, then the interface (cap boundary) cusps at its apex. For a relatively small velocity, the interface spreads to the vertices of the barrier, the apex zone remaining blunt shaped. We depict all the relevant domains and plot the flow nets using computer algebra routines.

## 1. Historical motivation from Zhukovskii's solution

There are very few analytical solutions to the problem of steady two-dimensional motion of an ideal irrotational heavy fluid with a free surface. This is caused by the necessity to satisfy the Bernoulli equation
1.1aalong an unknown isobaric streamline (free surface). In equation (1.1a), *V* is the magnitude of the velocity vector ** V**=

*V*e

^{iθ},

*θ*is the angle that

**makes with the horizontal axis**

*V**x*,

*g*=9.8 m s

^{−2}and

*y*is the vertical coordinate oriented upwards, against gravity.

Equation (1.1a), written in a system of coordinates moving with a constant horizontal velocity of the wave (with respect to a laboratory frame), is the main mathematical hindrance of the theory of nonlinear gravity-controlled waves, which stems from Stokes (1880) and which attracted towering figures in mechanics: Longuett-Higgins, Michell, Nekrasov, Sretenskii and Stoker, among others (for recent reviews, see Craik 2005; Gandzha & Lukomsky 2007; Dallaston & McCue 2010). Another branch of classical fluid mechanics dealing with equation (1.1a) is the theory of jets (Gurevich 1961), which sprang from Zhukovskii (1891), who himself was motivated by Michell. In both fields born to Stokes and Zhukovskii, the nonlinearity of equation (1.1a) makes a standard algorithm of problem solving impossible in linear theories, which is based on conformal mapping of a complex potential domain on the hodograph (Kirchhoff's function) domain (Gurevich 1961), because the latter is not known owing to the second term in equation (1.1a). Zhukovskii (1891) invented an ‘inverse’ trick: he specified a distribution of sin*θ* as a function of an auxiliary variable along the free surface such that equation (1.1a) is automatically satisfied. Zhukovskii (1891) found explicit expressions of the free surface, BOC (figure 1), which is sustained by two vertical walls CA_{dr} and BA_{dl}.

These walls and BOC confine an air chamber where the gas pressure is constant and high enough to withstand pressure exerted by the heavy fluid, which flows past this solid-gas ‘finger’ protruding as a static cap above the walls. The Zhukovskii (1891) idea was extended by Bervi (1894) and Richardson (1920) to other types of flows (e.g. an open-channel flow with a non-planar bed). Gurevich (1961) classified these solutions as ‘artificial’ because the impermeable walls (e.g. the channel bed or weir contour over which a free surface is formed) are ‘inversely’ reconstructed as part of the solution. Unfortunately, Zhukovskii (1891) presents a purely analytical solution but does not specify the incident velocity at infinity, which emerges as a part of his solution. Tremendous progress has been made in tackling directly stated and directly solved, gravity-controlled, free boundary problems by a combination of analytical and numerical (boundary-element method) techniques (Daboussy *et al*. 1998; Dias & Vanden-Broeck 2011). However, to the best of our knowledge, since Zhukovskii (1891), no non-trivial but explicit analytical expressions for the free surface, flow net and velocity field have been derived for boundary-value problems involving equation (1.1a) in its full form.

One fact, which seems to be not well exploited (overlooked?) in the literature on fluid dynamics, is that the Zhukovskii (1891) free surface (see his eqn (11)) satisfies the condition
1.1bwhere *V* _{0}/2 is the magnitude of the ambient velocity far away from the ‘finger’ in figure 1. In groundwater hydrology, equation (1.1b), with a replacement *V* _{0}→*k*, where *k*=const. is a hydraulic conductivity of soil, is a standard boundary condition on the so-called phreatic surface. This condition follows from physics, rather than a free ‘inverse’ choice as in Zhukovskii (1891), although the phreatic surface is—exactly as in Zhukovskii (1891)—an isobar and streamline (Polubarinova-Kochina 1962–1977; Strack 1989). It is often called the ‘water table’. For seepage flows, *V* in equation (1.1b) is the magnitude of the Darcian velocity. The characteristic functions in Darcian flows, the velocity potential and stream function obey Laplace's equation as in Stokes (1880) and Zhukovskii (1891). Consequently, the Zhukovskii (1891) solution of the flow problem depicted in figure 1 gives an exact solution of an interesting free-surface seepage problem to be studied here. Injection of air through a pressurized chamber as in Zhukovskii (1891) is hardly practical in applications to subsurface mechanics. So, we slightly modify the Zhukovskii (1891) flow scheme to adjust it to a situation often encountered in contaminant hydrology and environmental engineering (remediation of oil-contaminated aquifers), with BOC in figure 1 separating a lens of a light non-aqueous phase liquid (LNAPL) from groundwater descending from above (for a review of non-aqueous phase liquids (NAPLs), see Fetter 1999).

It is noteworthy that Zhukosvkii—simultaneously with his analytical studies of ideal fluids—worked on phreatic surface Darcian flows and—by modifying the Kirchhoff technique—introduced a function that now bears his name (e.g. Polubarinova-Kochina (1962–1977)). Notice that the model by Zhukovskii (1889) for two-dimensional seepage flow past a vertical barrier is mathematically similar to the flow in figure 1*a*. However, according to Vedernikov (1939) (see also Polubarinova-Kochina (1962–1977)), the Zhukovskii model gives rise to some inconsistencies in conjugation of the free surface and barrier edge.

In both ideal fluid and seepage flows, Zhukovskii reconstructed the characteristic functions (complex potential and complex velocity) by what was later called Chaplygin's method of singularities (Gurevich 1961), and did not always portray the complex potential and hodograph domains and the flow nets. In this study, we use the results of Zhukovskii (1891), and obtain his solution as a special case of a broad class of directly stated and solved problems in subsurface mechanics. We exploit the so-called inversion of the hodograph, which has been widely used in studies of Darcian free-surface flows by (Polubarinova-Kochina (1962–1977)) and her students (Polubarinova-Kochina 1962–1977; Emikh 1993). We depict all the relevant domains and plot the flow nets using Wolfram's (1991) Mathematica routines, which were, of course, not available in Zhukvoskii's days.

## 2. Engineering motivation in contaminant hydrology

NAPLs are categorized as LNAPLs that are lighter than fresh water and dense NAPLs (DNAPLs, denser than water). Crude oil, diesel, trichloroethylene, bitumen and other NAPLs (Fetter 1999; Oostrom *et al*. 1999) contaminate the groundwater/aquifer matrix/vadose zone soil and require costly and protracted environmental engineering techniques in remediation, monitoring and control. Often, NAPLs in a porous/fractured subsurface make macrovolumes (lenses, hydrocarbon traps, ganglia, blobs, streaks, etc.; Kacimov & Obnosov 2001*a*,*b*, 2008; Kacimov 2009, in press) separated from groundwater by sharp interfaces. The locus and shape of the interface between two liquid-saturated domains of contrasting density depends on the ambient groundwater flow and the geological boundaries (e.g. bed-, cap-rock unconformities, troughs, bulges, anticline flanks, etc.; Fetter 1999). In relatively coarse aquifers, commonly considered as homogeneous hydrostratigraphic units, the so-called clayey lenticular bodies (Schieber 2011; Van der Meer & Menzies 2011) or man-made structures (e.g. tunnels, subsurface waste repositories or missile silos; Kacimov & Nicolaev 1992; Kacimov 2007; Kacimov & Youngs 2009; Kacimov *et al*. 2011*a*,*b*) occur. They act as hydrodynamic encumbrances with respect to a natural or induced groundwater flow. Lehr & Wright (1963), Strack (1989), Kung (1990) and Oostrom *et al*. (1999) studied the winding seepage past such barriers, in a sense similar to cascades of aerofoils or other barriers in classical fluid mechanics (Milne-Thomson 1958). A vertical cross section of such a coarse-material aquifer with imbedded obstacles is shown in figure 2*a*. In the case of aquifer contamination by NAPLs, a common remediation technique is to either induce a flushing flow from surface ponds, injection–abstraction wells, well galleries, etc. (Fetter 1999; Al-Maamari *et al*. 2011) or to rely on ‘natural attenuation’ (i.e. to leave the contaminated soil/rock to periodic leaching by rainfall-induced infiltration).

Detection of ‘alien entities’ (baffles in figure 2*a*), which obstruct flow due to their contrast in texture with the main porous medium, is in most remediation projects problematic as a result of imprecision of geophysical tools and limited core samples from the boreholes. Therefore, analytical and numerical modelling of multi-phase flow and transport in the subsurface (Pinder 2008) is a vital step in many surveys, containment projects and decontamination by pump-and-treat, air-sparging, permeable reactive barriers, LNAPL skimming, DNAPL scraping, etc. (Fetter 1999).

In managed aquifer recharge projects, an intensive vertical groundwater gradient is artificially generated by injection of fresh water from infiltration ponds. In these hydrodynamically agitated clean-up schemes, NAPLs can be flushed away (ideal scenario) or the NAPL–groundwater interface can take a non-trivial (even counterintuitive) static shape (deleterious scenario). In this paper, we study analytically the latter when an LNAPL is attached as a macrovolume to a solid barrier, with a strong influence of the conjugated dynamic groundwater on the interface configuration.

We consider a homogeneous porous matrix of hydraulic conductivity *k* in a descending, fully saturated, Darcian seepage with an incident velocity *V* _{0} far above the obstacles (figure 2*a*). Figure 2*a* presents a ‘late flushing’ regime of a clean-up operation. At , the medium in figure 2*a* was filled by LNAPL and intensively leached afterwards. The major part of the LNAPL was removed (e.g. entrained by a strong groundwater motion as in Al-Maamari *et al*. 2011). What has remained is shown in a grey tone as LNAPL caps (figure 2*a*) sitting on the rooves of the baffles. These caps, as static lenses, are hydrodynamically stable because the seepage flow just pushes them to the stoss side of the barriers. If there is no ambient flow of groundwater in figure 2*a* (i.e. *V* _{0}=0), the ‘alien’ caps would rise to the regional water table (or bottom of an infiltration pond) owing to the Archimedian force. On the other hand, if *V* _{0} is large enough, then the caps would be flushed downwards, similar to the so-called hydrodynamic hydrocarbon traps in anticlines (Kacimov & Obnosov 2008).

Although the size of the light-phase caps in figure 2*a* may be relatively small, their insidious environmental effect is caused by (i) dispersive stripping of contaminants from the lens–groundwater interface to the main flow (Chrysikopoulos *et al*. 1994), and hence a protracted groundwater contamination, and (ii) retarded resurfacing of LNAPLs on/in the regional water table/infiltration pond (wherefrom LNAPLs are often skimmed; Fetter 1999) if *V* _{0} in figure 2*a* is sufficiently high. In our remediation of LNAPL contaminated sites (Al-Maamari *et al*. 2011), we faced a recurring spiking of contaminant concentration in monitoring wells: after several months of a continuous air-sparging and pump-and-treat operation, the contaminant level was reduced to a target value, but after a following period of hydrodynamic dormancy (no artificial agitation of the aquifer), LNAPLs re-emerged as a free product. We already attributed this to ‘hydrodynamic’ accumulation of ‘floating’ caps of LNAPLs on the phreatic surface (Kacimov *et al*. 2011*a*,*b*). In this paper, we warn the environmental engineers, treating heterogeneous aquifers, of potential LNAPL caps and pendants like in figure 2.

We consider a tall two-dimensional barrier I in figure 2*a*. Our objective was to find the shape of the sharp interface BOC as a function of the incident velocity, quantity of LNAPL, densities of groundwater and LNAPL, barrier sizes and soil permeability, and to determine the flow field near the cap barrier.

## 3. Zhukovskii's flow regime and description of method

We zoom in on barrier I, a rectangle EFGH, in figure 2*a* and depict its upper part in figure 2*b*. We assume that the barrier is not tilted, i.e. aligned with the vector of gravity acceleration and that the middle points of the barrier (A_{dr}, A_{lr}) are so far from E, F, G and H that flow there uniformizes and acquires the incident velocity *V*_{0}↓↓** g**, i.e. this vector is collinear with the vector of gravity acceleration.

We also assume that other barriers are far apart from the selected one such that the flow domain, *G*_{z}, in figure 2*b* is unbounded from all sides, i.e. the bafflers in figure 2*a* are relatively small when compared with the distance between the neighbours. In other words, we study the near field of a solitary barrier. We neglect any dispersion or diffusion of LNAPL from a perfectly confined lens. The capillary fringe between the groundwater and LNAPL is also neglected.

In this section, we consider the ‘Zhukovskii’ (or ‘critical’) case when the cap covers the whole roof EF of the barrier, i.e. points B and C of the interface coincide with E and F, respectively (figure 2*b*). Following Zhukovskii (1891), we select the origin O of a Cartesian coordinate system *x*O*y* coinciding with the apex O of the LNAPL cap, as illustrated in figure 2*b*. The horizontal dimension of the barrier is L and the cap is confined from below by a horizontal segment EF. From above, LNAPL is bounded by an abrupt interface BOC whose height is *b* above the middle of the barrier roof.

Darcy's law states ** V**(

*x*,

*y*)=−

*k*∇

*h*, where the hydraulic head

*h*(

*x*,

*y*) is counted from point O, i.e.

*h*(0,0)=0. The head is related to groundwater pore pressure

*p*

_{w}via the relation

*p*

_{w}=

*ρ*

_{w}

*g*(

*h*−

*y*+

*C*

_{p}), where

*ρ*

_{w}is the water density and

*C*

_{p}is a pressure head at point O (this not essential constant vanishes from the analysis). A complex potential is defined as

*w*=

*φ*+i

*ψ*, where

*φ*=−

*kh*is the velocity potential and

*ψ*is a complex-conjugated stream function. As in Zhukovskii (1891), both

*φ*and

*ψ*are harmonic functions. LNAPL density is

*ρ*

_{o}=const,

*ρ*

_{o}<

*ρ*

_{w}and the pressure of LNAPL inside the cap is hydrostatic. There is no pressure jump across BOC and from elementary analysis (Polubarinova-Kochina 1962–1977 or Kacimov & Obnosov 2001

*a*,

*b*, 2008), the water potential on the interface satisfies the condition 3.1Obviously, BOC, CA

_{dr}and BA

_{lr}are no-flow lines. Point O is a stagnation point and BOC is symmetric with respect to O

*y*. Therefore, the streamline A

_{u}O splits into two branches OCA

_{dr}and OBA

_{lr}, and we select

*ψ*=0 along this streamline. The complex potential domain

*G*

_{w}is shown in figure 2

*c*, where the potential at points B and C, located on two opposite banks of the cut, is

*φ*

_{0}=

*C*

_{d}

*b*, in congruency with equation (3.1).

A complex Darcian velocity *V* =*u*+i*v* is defined, where *u*(*x*,*y*) and *v*(*x*,*y*) are its horizontal and vertical components. It is well known (Polubarinova-Kochina 1962–1977; Kacimov & Obnosov 2001*a*,*b*) that the image of BOC in the hodograph domain, *G*_{V}, shown in figure 2*d*, is a circle of diameter *C*_{d} and centre at the point (0,−i*C*_{d}/2). The equation of the circle is *u*^{2}+*v*^{2}+*C*_{d}*v*=0. In *G*_{V}, the image of the barrier walls is a cut such that points A_{u}, A_{dr} and A_{lr} (different in *G*_{z}) merge into a single point (tip of a finite straight cut). Obviously, from figure 2*d*, 0<*V* _{0}<*C*_{d}, which determines the hydrodynamic constraint on the regime of figure 2*b*.

Thus, we have two domains, shown in figure 2*c*,*d* for two characteristic flow functions. The Zhukvoskii linear relation between the magnitude of velocity and sin of the angle, which the velocity vector makes with a horizontal axis, results in a circular arc as the image of the free surface in the hodograph plains. Like in classical fluid mechanics of a weightless fluid, the Kirchhoff function makes the circular hodograph domain boundaries straight (Gurevich 1961), in what follows, we will invert figure 2*d* onto standard polygons conformally and map them onto the polygons in figure 2*c*.

It is well known that *V* (*x*,*y*) is an antiholomoprhic function, and *u*−i*v*=d*w*/d*z* is holomorphic. We invert the disc with a cut in figure 2*d* into a half-plane *G*_{ω} with a cut in the plane d*z*/d*w* shown in figure 2*e* (see Polubarinova-Kochina (1962–1977) on the theory of this method). Obviously, *G*_{ω} is a tetragon. Next, we map *G*_{ω} onto *G*_{w} via an auxiliary plane *ζ*=*ξ*+i*η* (figure 2*f*).

The domain *G*_{w} is mapped on the half-plane *G*_{ς} (*η*>0) shown in figure 2*f* by the function
3.2with the correspondence of points B→−1, A→1, O→0, .

By the Schwarz–Christoffel formula, *G*_{w} is mapped on *G*_{ς} by the function
3.3where the branch of the radical is fixed in the upper half-plane and is selected to be positive at −1<*τ*<1. We emphasize that at 1<*τ*, the radical is and at *τ*<−1, it is . For the sake of brevity, we omit the calculation of the index of this problem. The mapping constant *A*_{1}>0 in equation (3.3) is evaluated from the locus of point A (=−i*V* _{0}) in *G*_{V}. Then, upon integration in (3.3), we obtain
3.4From equations (3.2) and (3.4),
3.5The function (3.5) satisfies the symmetry condition because arcsin *ζ* and satisfy the same condition.

Owing to *z*=*L*/2−i*b* at *ζ*=1, from equation (3.5),
3.6We introduce dimensionless variables as , and the superscript asterisks for them is dropped.

We separate the real and imaginary parts in equation (3.5) and, by taking into account equation (3.6), obtain the parametric equations of the right half of the sharp interface, OC, 3.7The cross-sectional area of the cap is 3.8

Obviously, the Zhukovskii (1891) flow (e.g. see parametric eqns (13) of CO) is a particular case of the class defined by equations (3.2) and (3.5). Namely, if *C*_{d}/*V* _{0}=3/2 in equations (3.2)–(3.7), then we arrive at exactly the Zhukovskii (1891) solution.

In figure 3*a*, we plot the interfaces (3.8) for *C*_{d}=3/4, *V* _{0}=2/3, 1/2 and 1/9 (curves 1–3, correspondingly). Curve 2 is, obviously, the Zhukovskii (1891) contour. The flow net is shown in figure 3*b* for *C*_{d}=0.2, *V* _{0}=0.15, where the solid curves are streamlines (*ψ*=0,0.01,0.02,…,0.07 for the right half of figure 3*b*, counted from the axis of symmetry rightwards), and dashed lines are equipotentials (*ϕ*=−0.05,0,0.05,…,0.25 counted downwards from the uppermost one in figure 3*b*). As is evident from figure 3*a*, the increase in the incident velocity causes a drastic increase in the size of the lens of a fixed LNAPL density.

By contrast to the ‘inverse’ solution of Zhukovskii (1891), ours is obtained by a physically and mathematically direct method, all boundary conditions and characteristic domains are determined by the conceptual model rather than by manipulations with boundary conditions of analytic functions.

As we have already mentioned, if *V* _{0}>*C*_{d}, then the hodograph domain in figure 2*d* cannot be realized such that the interface is sustained by exactly the vertices of the rectangle in figure 2*a*. Obviously, the ambient velocity can be higher than *C*_{d}. Then, what? The answer to this question is given in §4.

## 4. Partially covered barrier

Now, we study the situation when the roof is only partially covered by the LNAPL; i.e. points B and C (the cap tips) do not coincide with E and F. The cross-sectional area *S* of the cap is given, and our objective is to predict how this given ‘icing’ of the LNAPL is spread over a horizontal barrier EF. The width of the whole oil lens is *l*, *l*<1 (we keep the same dimensionless quantities as in §3). Both *l* and *b* are to be found.

The *G*_{w} domain corresponding to *G*_{z} is an upper half-plane *ψ*>0 with a cut A_{dr}FCMONBEA_{dl} shown in figure 4*b*, where M and N are two symmetric inflexion points on the interface. The potentials at points (M, N, E and F) are not known in advance (and practically not needed). The *G*_{V} domain in figure 4*c* is a lower half-plane with a circular cut MC(B)NOM and a straight cut FAE. There are two stagnation points (O and C) and two points of infinite velocity (F and E). The horizontal segments CF and BE in figure 4*a* are not covered by the cap. Owing to symmetry, we consider only the right half of *G*_{z} in figure 4*a* and corresponding halves of *G*_{V} and *G*_{w}.

The complex velocity at point M, *V* _{M}=*u*_{M}+i*v*_{M}, is not known in advance (not needed in applications). We emphasize again that at the sharp edge, F in figure 4*a*, the velocity blows up and therefore the LNAPL–groundwater interface cannot be sustained there (see Vedernikov 1939 who commented on Zhukovskii (1889)). In §4 and in Zhukovskii (1891), the infinite velocity was avoided by forcing the cap to cover the whole roof (figure 2*a*) or—in the vernacular of Zhukovskii (1891)—the gas bubble to be attached to the walls of the nozzle (figure 1).

As in §3, we invert the circular tetragon (right half of figure 4*c*) into a regular tetragon *G*_{ω} shown in figure 4*d*. *G*_{w} of the right half of *G*_{z} is an upper half-plane, and we map it onto the half-plane *G*_{ς} (*η*>0) by the function
4.1*G*_{ω} is mapped on *G*_{ς} by the function
4.2where the radical's branch in the integrand is fixed in the upper half-plane by the condition of its positiveness at *τ*>−1 (then, at *τ*<−1, this radical becomes i|*τ*+1|^{3/2}). The Schwartz–Christoffel constant *A*_{2}>0 is to be found. The correspondence of points in mappings (4.1) and (4.2) is C→1, M→*m*, O→−1, *A*_{dr}, *A*_{u}→−*a*, , where 1<*m*<−1, *a*>1 are the affixes of points M and A, which are also to be found.

Integration of equation (4.2) gives
4.3where the vanishing at infinity branch of the function *f*(*ζ*) is fixed in the upper half-plane.

We use equation (4.3) and the relation for point M in *G*_{ω} *f*(*m*)=−(*u*_{M}/*v*_{M}+i)/*C*_{d} (figure 4*d*), where and obtain the following equations (−1<*m*<1,*A*_{2}>0):
4.4and
4.5Next, by using equation (4.3) and the condition *f*(−*a*)=−i/*V* _{0}, we obtain the equation
4.6By using (4.1), (4.3) and (4.5), we integrate,
4.7The last geometrical condition in equation (4.7) is , i.e.
4.8At point C in figure 4*a*, *z*(1)=*l*/2−i*b*. Then, from equations (4.7) and (4.8),
4.9As in §4, *y*(*ξ*) and *x*(*ξ*) are the imaginary and real parts of the equation of BOC obtained from *z*(*ζ*) in equation (4.7). The dimensionless cross-sectional area of the cap (LNAPL volume) is integrated explicitly as
where
and
From the last three relations, we obtain
4.10We use the routine *FindRoot* of Mathematica to solve the system of nonlinear equations (4.6) and (4.10) with respect to *m* and *a* at given *C*_{d}, *V* _{0} and *S*. The found values (*a*,*m*) are inserted into equations (4.8) and (4.9) that gives the cap sizes. Table 1 shows the results of computations for *C*_{d}/*V* _{0}=2 and *S* varying from 0.001 with the increment of 0.005*k*, *k*=0,1,…,7.

Now we analyse what happens if points M and N in figure 4*c* approach each other and merge into the point −i*C*_{d}. In this limit, there are two different situations, depending on *V* _{0} and *C*_{d}. In congruity with terminology of Polubarinova-Kochina (1962–1977) for similar problems in horizontal drainage (see also Kacimov *et al*. 2011*a*,), these limits are called the ‘critical’ regimes.

If *V* _{0}>*C*_{d} (as in figure 4*c*), then merging of points N and M requires point O to jump from the origin of coordinates to −i*C*_{d}. Consequently, in *G*_{V}, the two-tip circular cut (figure 4*d*) degenerates into a full circle (figure 4*f*). In *G*_{z}, the apex O of the cap becomes a cusp where the magnitude of the Darcian velocity equals *C*_{d}. The half of *G*_{V} in this ‘first critical’ regime is a circular triangle. Riesenkampf (see Polubarinova-Kochina (1962–1977), pp. 178–179) and Kacimov (in press) mapped conformally this triangle onto a half-plane.

If in figure 4*c* *V* _{0}<*C*_{d}, then merging of the points M and N causes a different criticality. Point O remains at the origin of coordinates, but points B and C in figure 4*c* jump from the origin to the point −i*C*_{d}. Similarly, points E and F jump from infinity to the same point −i*C*_{d}. Consequently, the hodograph in figure 4*c* transforms into one presented in figure 2*d*, i.e. into the interior of a circle with a straight cut CAB. Therefore, the solution in §4 is the ‘second critical’ regime. It is noteworthy that the sudden jumps of characteristic points and corresponding changes of *G*_{V} (resulting in a reduced number of vertices of the corresponding circular polygons or sudden vanishing of a whole subdomain, such as the circle in transition from figure 4*d*–*f*) are theoretically tackled by the Carathéodory theory.

The analysis of criticality can be also conducted for fixed *C*_{d} and *V* _{0}, but for varying *S*. For the ‘second critical’ regime in the limit *S*→*S*_{cr}, where *S*_{cr} is given by equation (3.8), the solution of this section degenerates into one of §4. The maximal static area, which can be sustained by the barrier for this particular LNAPL and relatively small incident velocity, is bounded by the contour (3.7). In the ‘first critical’ regime with a cusping apex of the cap, *S*_{cr} can be easily obtained as in Riesenkampf (Polubarinova-Kochina (1962–1977)) and Kacimov (in press) (we drop the corresponding derivations for the sake of brevity).

Upon determining the couple (*a*,*m*), we expressed *ς*(*w*) from equation (4.1), put it into equation (4.7) and using the routines of Mathematica, plotted the cap contour, walls and flow net. Figure 5*a* shows the results for the transition to the ‘first critical’ regime at *C*_{d}=0.1, *V* _{0}=0.2 and *S*=0.001, 0.041, 0.081, panels from left to right, correspondingly. As is evident from this figure, with the increase in the area of the LNAPL cap, its middle zone tapers and eventually makes a cusp at point O. The tips B and C of the cap flanks do not, however, spread to the edges E and F of the barrier.

Figure 5*b* illustrates the flow nets and transformation of the interface at *C*_{d}=0.2, *V* _{0}=0.1 and *S*=0.006, 0.016, 0.040, panels from left to right, correspondingly. As we can see from figure 5*b*, with an increase in *S*, the lens shape bulges from a ‘pancake’ at small *S* to a pronounced protuberance approaching contour (3.7), for which the ‘second critical’ regime is attained. In the limit, the whole roof of the barrier is covered by the LNAPL, the barrier edges coincide with the two ends of the free surface, and the velocity at these two points becomes finite, as in the Zhukvosky scheme of figure 1.

We also used an alternative approach to find the affixes. Instead of *S*, we fixed the cap bottom width l, found the pair (*a*,*m*) and, eventually, (*b*,*S*). In this way, only one nonlinear equation rather than a system of two was solved by the *FindRoot* routine of Mathematica. To save space, we drop the corresponding details.

## 5. Light non-aqueous phase liquid pendants

The barriers in figure 2*a* can also hydrodynamically detain LNAPL macrovolumes on their lee sides (i.e. from below). The interface then confines a light ‘alien’ volume, which is attached to the bottom GH of the barrier in figure 2*a*. Figure 6*a* shows such a Zhukovskii (1891) type regime with a fully covered floor of a barrier. *G*_{z} for this regime is geometrically symmetric with one in figure 2*b*. *G*_{w} and *G*_{V} are shown in figure 6*b*,*c*. Obviously, the solution obtained in §4 holds—with simple changes of notations—for the pendant GOH in figure 6.

Obviously, LNAPL volumes can accumulate simultaneously on both the barrier roof (figure 2) and floor (figure 6). Moreover, an LANPL bubble without any barrier can be suspended in a descending flow, which hydrodynamically balances the upward-oriented Archimedian force acting on the bubble. This problem, analogous to another famous Zhukovskii (1891) bubble problem (see also Gurevich 1961), will be studied in our next paper.

The very presence of attached LNAPL lenses as in figures 2 and 6 is counterintuitive to many hydrogeologists (e.g. Fetter 1999 showed only DNAPL lenses in aquifuge troughs). This makes our solutions practically deleterious. Indeed, a large infiltration pond (a common situation in managed aquifer recharge schemes in Oman) seeping for a long time into an allegedly clean subjacent porous medium may hydrodynamically—by maintaining *V* _{0} in figure 2—entrap LNAPLs originally present in the vadose zone. Then, a long-term degradation of quality of water stored/recovered in or from the aquifer due to dispersion–diffusion (Fetter 1999) will occur as a result of insidious lenses of LNAPL.

## 6. Conclusions and future directions

In this study, we found a striking analogy between a seminal Zhukovskii solution for a heavy fluid flow with a gas–liquid interface and a Darcian flow past a barrier with a groundwater–LNAPL interface. A static NAPL cap or pendant clings to the roof/bed of the barrier. Interfaces emerge in many engineering applications of classical fluid mechanics (Stone 2010) and the old solutions—such as Zhukovskii (1891)—can be revisited and used in hydrogeology, contaminant hydrology and reservoir engineering. For practitioners dealing with aquifers, pays and vadose zones, the very existence of non-trivial multi-phase flow regimes with entrapped macrovolumes of a static ‘alien’ phase is a surprise because they are used to simple models of one-phase, non-hysteretic Darcian flows. In ideal flows of free water (Verhoff 2010), observations of attached gas cavities is a field–laboratory routine, and the corresponding two-phase macro-scale complexity is evident. In subsurface mechanics, a direct monitoring of LNAPLs such as in figures 2 and 6 is, obviously, not possible, and mathematical modelling is a valuable tool complementing laboratory tomography of soil columns and field resistivity surveys. One of the main results of this paper is the effect of an incident velocity on static LNAPLs within a saturated soil/rock. Namely, a more intensive leaching (higher *V* _{0}) from an infiltration pond located above the baffle lens in figure 2*a* causes an alteration of the interface as in figure 3*a*. This is totally opposite to what a remediation engineer may expect. Indeed, instead of an improved clean-up, a larger cap and, consequently, more hydrodynamic entrapment of LNAPL occur.

Fluid dynamics of Darcian flows in porous media is fascinatingly rich with beautiful analytical curves, which have been revealed and analysed in various potential fields: the Taylor–Saffman bubble as an optimal waste repository and concrete dam contour, the Preissman channels as an optimal low-permeability sedimentation trough, ellipses, parabolas and parabolloids as heterogeneities generating collimated ‘internal’ velocity fields, the Nicomedes conchoids as isobars in a seepage flow to a semicircular drainage trench, catenary as a baffling contour in a tension-saturated flow, among others (Kacimov & Nicolaev 1992; Kacimov & Obnosov 2000; Kacimov 2001, 2006*a*,*b*, 2008; Kacimov & Youngs 2009). A beautiful property of the family of contours described by equation (3.7), which imbeds the Zhukovskii (1891) curve, is its trochoidal form. We recall Kozeny's trochoids (Polubarinova-Kochina (1962–1977)) as contours of soil channels seeping with phreatic surfaces. Similarly, the streamlines of an ideal, circulation-free fluid flow past a cylinder are the conchoids of Nicomedes—a fact ignored in classical textbooks on fluid mechanics.

We can easily extend our analysis to other types of barriers (similar to Verhoff 2010, who studied flows without porous media) or to tension-saturated flows where the potential is also a harmonic function. Instead of a tall barrier I in figure 2*a*, examined in this study, we can tackle a thin but broad barrier II. The recent progress in finding analytical solutions to flow problems with arbitrary driving singularities in piece-wise homogeneous porous media (Obnosov 2006, 2009, 2011) makes the extension to seepage regimes with a macro-volume of LNAPL attached to low-permeable (rather than impermeable) baffles and to ambient groundwater moving from injection to abstraction wells promising, rather than to vertical infiltration as in this study.

## Acknowledgements

This work has been supported by the ‘Estimating natural groundwater recharge and discharge in North Oman’ (grant no. SR/SCI/ETHS/11/01), His Majesty Research Trust Fund (Oman) and the Russian Foundation of Basic Research (grant no. 12-01-97015-r_povolgh'e_a). The helpful comments of three anonymous referees are appreciated.

- Received May 24, 2012.
- Accepted July 26, 2012.

- This journal is © 2012 The Royal Society