Published online 27 February 2007
Published in Vadose Zone J 6:168-174 (2007)
DOI: 10.2136/vzj2005.0141
© 2007 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
Dipole-Generated Unsaturated Flow in Gardner Soils
Anvar Kacimov*
Dep. of Soils, Water and Agricultural Engineering, P.O. Box 34, Al-Khod 123, Sultan Qaboos Univ., Sultanate of Oman
* Corresponding author (anvar{at}squ.edu.om)
Received 4 December 2005.
 |
ABSTRACT
|
|---|
The Bystrov explicit analytical solution for viscous, low-Reynolds number flow in layers of variable thickness is interpreted as infiltration in a Gardner homogeneous soil obstructed by a subterranean cavity, stone, or other obstacle. Mathematically, the two-dimensional advectiondispersion equation for the Kirchhoff potential is solved by combination of a term responsible for incident unidirectional infiltration, and a term describing a dipole. Superposition results in a separatrix (a cavity or stone contour), outside of which streamlines are deflected from vertical lines, and constant potential (pressure, moisture content) lines demarcating lobe-shaped domains. The physical impedance of the obstacle causes a buildup of moisture near the leading edge and a dry zone near the trailing edge of the obstacle. The criticality conditions of the model were also tested (i.e., that the moisture content in the flow domain is less than porosity but greater than zero).
Abbreviations: AEM, analytical element method MHE, modified Helmholtz equation SSADE, steady-state advectiondispersion equation
 |
INTRODUCTION
|
|---|
ON A LARGE SCALE, flow of moisture in the vadose zone is often modeled as purely descending (infiltration) or ascending (capillary rise and evaporation) through a soil that is homogeneous in its hydraulic properties. These assumptions are useful in constructing the hydrological balance of a catchment, in particular for assessing total groundwater recharge from the vadose zone. However, real soils are heterogeneous and interspersed by constructed (e.g., the Yucca mountain system of tunnels and shafts) and natural (e.g., stones, mole or worm holes) inclusions. For example, traditional water supply systems in Oman use falaj (qanat) tunnels that are partially submersed under the water table and partially span through the vadose zone. Infiltration either bypasses the tunnel or breaks through, thus adding to the overall water supply of the tunnel. In such cases, uniform moisture inflow through a flat soil surface leads to a meandering of flow paths through more permeable zones. This funneling causes a more rapid arrival of solute particles at some locations. The tunnel or other obstacles will also produce certain areas in the flow system where saturation can be substantially higher (e.g., perched "water table") or lower (reduced root development) than that of the surrounding locations.
Philip et al. (1989a, 1989b) enunciated a conceptual model of unsaturated flow obstructed by a subterranean cavity with an impermeable boundary. The lack of any water entering a cavity in the unsaturated zone can be quite puzzling to those who are used to saturated flow systems where water will leak into any subsurface opening, such as an isobaric (tunnel, cavity) or constant head (drain) boundary. Although the physical phenomenon of flow impedance by empty (or, generally, coarser-material) inclusions embedded in a fine-textured soil matrix has long been known to petroleum engineers as the "end effect" (Dullien, 1992), and to geotechnical engineers as the "capillary barrier effect" (Kachinskii, 1970), the analytical description of Philip elucidated fascinating spatial distributions of moisture (e.g., lobe-shaped isobars "hanging" from the cavity, and a "dry shadow" behind the obstacle).
Philip's analysis was based on the Gardner exponential soil conductivity function (Gardner, 1958), which leads to the steady-state advectiondispersion equation (SSADE) for the Kirchhoff potential, or the modified Helmholtz equation (MHE). Since the resulting equations are linear, the model is often called quasilinear. Linearity allows one to superpose particular solutions (e.g., point singularities) into new solutions. For example, Philip et al. (1989a) combined multipoles of the SSADE and matched the no-flow condition at the boundary of an impermeable circular cylinder. For more complicated shapes (e.g., an ellipse or polygon) rigorous solutions become awkward and a boundary-layer approximation (Philip, 1990a) or numerical method must be employed. The method of distributed singularities of the SSADE has recently become popular (e.g., Bakker and Nieber, 2004; Warrick and Knight, 2002) in studies of piece-wise homogeneous porous media in which the pressure head and the normal Darcian velocity are matched along lines or surfaces separating the different soil materials. In subsurface mechanics this method is known as the analytical element method (AEM), which has been widely used also in groundwater hydrology (Strack, 1989). An advantage of the quasilinear model within the AEM stems from the absence of free boundaries (Philip, 1990b), which otherwise would have required tedious delineation of the exact flow domain in the GreenAmpt model (Youngs et al., 2004) that describes the impeding disturbance of infiltration by cavities.
Warrick and Fennemore (1995) suggested a different method to study two- or three-dimensional flow past impermeable obstacles based on original mid-19th century ideas by Rankine, who designed ship hulls by placing hydrodynamic singularities inside the ship body. Later in the 20th century, the entire theory of airfoil lift commenced with a Rankine-type combination of a simple vortex inside a wing and ambient unidirectional flow (e.g., Elizarov et al., 1997). Warrick and Fennemore (1995) also determined the shape as a stream line of flow generated with a Gardner sinksource pair placed in a uniform infiltration flow field. This idea was expanded further by Kacimov (2006), who found an isobaric depression generated by two singularities. In this paper we pursue the same WarrickFennemore method of semi-inversely determined shapes.
After the quasilinear model for unsaturated flow was introduced into the literature, several soil physicists recognized the beauty of the model and tried to employ methods involving complex analysis. Such methods work well for two-dimensional saturated flow where the characteristic functions (hydraulic head and stream function) satisfy the Laplace equation. Therefore the complex potential is a holomorphic function (Strack, 1989). Philip et al. (1989a)trying to embark on the machinery of these functionswrote that "Knight (unpublished manuscript, 1988) has shown that Cauchy-Riemann-like equations...yield simple relations between the constituent multipoles of...two quasi-conjugate functions." However, few if any approaches have been reported for the quasilinear model that were as efficient as the methods of conformal mapping, hodograph, or Signorini formulabased techniques, which have proven to be very versatile in analyses of saturated flow (Polubarinova-Kochina, 1977).
A group of Soviet scientists led by K.N. Bystrov (1956, 1959) (scanned copies of papers in Russian referenced in the following two paragraphs are available from the author) developed a systematic approach to describe viscous laminar flow (without any porous media) in what they called "layers of variable thickness." Bystrov implemented the mathematical technique of Bers and Gelbart (1944) for generalized analytic functions, which he formally assembled from metaharmonic functions, in particular those which satisfy the MHE. He defined the operations of differentiation and integration for quasianalytic functions and Cauchy-type formula. Using this approach, explicit rigorous solutions were obtained in terms of the Kirchhoff potential and/or stream function for a source, sink, vortex, dipole, multipole, and their superpositions (Bystrov, 1959; Yourisov, 1965). Later some of these solutions were rediscovered in soil physics (Philip, 1969) and heat transfer theory (Concer, 1959). The work of Bystrov and his associates preceded similar mathematical developments by Duffin (1971).
Unfortunately, none of Bystrov's definitive papers have been translated or referenced in the Western literature. Moreover, even major Russian books on subsurface mechanics (e.g., Polubarinova-Kochina, 1977) did not report Bystrov's elegant solutions. This is extremely surprising since Soviet scientists studied two- and three-dimensional unsaturated flow using quasilinear models before Philip, Raats, Zachman, or anyone else in the West. For example, Polubarinova-Kochina and Kulabukhova (1959) and Kulabukhova (1959) analytically investigated transient quasilinear seepage from point and line sources (subsurface emitters). We note here that Frenkel (1949) was first to use multipoles of the SSADE, while Borzykh and Cherepanov (1978) pioneered the use of Mathieu function expansions for SSADE problems involving convective heat transfer. Despite the conspicuous lacuna in the use (and even referencing) of Bystrov's clear results in soil physics, they were widely implemented in the Soviet spacecraft engineering, contemporaneously with the completely independent analytical endeavors by the abovementioned Western soil physicists and hydrologists dealing with the SSADE and MHE equations. For example, the cover of the book by Novikov et al. (1991) has Bystrov's (1959) flow net for the case of a source (sink) (a less general solution was obtained by Philip, 1969), while the content of the book itself reflects the development and wide applicability of many of Bystrov's ideas.
The objective of this study was to address this lacuna and to resuscitate Bystrov's legacy by fine-tuning it to water flow in the vadose zone. Bystrov's dipole solution is exploited by adjusting it to the quasilinear model. From Bystrov's explicit and simple expressions for the Kirchhoff potential and the stream function, and with the help of Mathematica (Wolfram, 1991), one can construct the flow nets and inspect systematically the shape of the corresponding subterranean cavities.
 |
ANALYTICAL DESCRIPTION OF FLOW SYSTEM
|
|---|
Consider an impermeable cavity BMCB (Fig. 1A
shows a vertical cross-section) placed in an ambient unsaturated flow field. The moisture content, pressure head and Darcian velocity vector at infinity are ma, pa*,
a*, respectively. A system of Cartesian coordinates (x*,y*),y** = y* originates at point O where the singularity occurs.
The quasilinear model assumes the hydraulic conductivity k to be of the form
 | [1] |
where ks is the saturated hydraulic conductivity,
is the sorptivity number, and p* is the pressure head. The Kirchhoff potential
* is defined as
 | [2] |
The potential should satisfy the inequalities 0 <
* <
s*, where
s* = ks/
corresponds to the complete saturation limit (Philip, 1992). If mathematically we obtain
* >
s* on the contour of the cavity, then this means physically that either the cavity leaks and flow becomes supercritical (Kacimov, 2000) or, provided that a no-flow condition is kept at the boundary of the obstacle (e.g., a stone; see Warrick and Fennemore, 1995), a "perched" saturated zone is formed. The saturated zone should be matched with the ambient unsaturated flow field as, for instance, in Philip (1992) and Kacimov (2003) by conjugating the pressure head and streamlines along the boundary between the two zones. If
* < 0 (as is the case in the vicinity of a sink or in the sink portion of a doublet singularity), then the corresponding domain should be excluded as physically meaningless. We recall that in the common theory of horizontal drains or vertical wells in groundwater hydrology a sink is excluded from the flow domain by contouring a certain area around it by an effective drain or well quasi-circumference where the positive pressure head is fixed (see Polubarinova-Kochina, 1977, for details).
Define the dimensionless variables (x, y, y*,
,
a) =
x*,
y*,
y**,
*/
s*,
a*/ks). The potential
(x,y) satisfies the SSADE:
 | [3] |
We now introduce a stream function
(x,y) according to the conditions (Warrick, 2003, p. 251253)
 | [4] |
where u(x,y) and v(x,y) are the horizontal and vertical components of
a. As is well-known (Warrick, 2003),
satisfies the same SSADE given by Eq. [3].
 |
CAVITY-HAMPERED INFILTRATION
|
|---|
Bystrov (1959) obtained the following solution to Eq. [3]
 | [5] |
in terms of the stream function and
 | [6] |
in terms of the potential, where q is the dipole strength (constant), r =
, while Ko, K1, K2 are the MacDonald (or modified Bessel) functions of zero, first, and second order, respectively. The first term on the right-hand sides of Eq. [5] and [6] corresponds to ambient unidirectional infiltration with 
= v
=
a, whereas the second term represents the dipole. The dipole center coincides with O, while the dipole axis lies on Oy. For completeness, and since many readers may not have easy access to Bystrov's work (in Russian), we summarize in an Appendix of this paper the derivation of Eq. [5] and [6] following the original development by Bystrov (1959).
The velocity components calculated from Eq. [4] and [5] are:
 | [7] |
Streamlines (two of which,
=
1 and
=
2,
1 >
2 > 0, are depicted in Fig. 1 as dotted and dashed lines, respectively) typically consist of two disconnected branches. One branch stretches from +
to
in the y direction, while the other branch loops and touches the Oy axis inside the cavity. The streamlines
= c < 0 are symmetric with those shown in Fig. 1 about the Oy axis. The streamline
= 0 is a mathematical separatrix. It consists of an upper ray CA1, which bifurcates at the leading edge B of the obstacle, then coincides with the cavity contour between B and C, and converges back to a vertical line BA2 at the trailing edge B. The separatrix divides the whole flow domain in Fig. 1 into two zones. The external zone is the soil surrounding the cavity, while the internal zone is analogous to the Rankine ship-interior (i.e., this zone is not relevant to our infiltration problem).
The ordinates yB and yC are found from the conditions
B =
C = 0 and vB = vC = 0. These conditions result in the following nonlinear equation obtained from Eq. [5]:
 | [8] |
We solved Eq. [8] using the FindRoot routine of Mathematica (Wolfram, 1991). The positive root of Eq. [8] corresponds to point B and the negative root to point C. Obviously, the total height of the obstacle is h = yB yC. Figure 2A
shows h,yB, and yC as functions of
a/q (curves 13). We note that all graphs here and below are plotted in (x,y*) coordinates.
Next we substituted yC resulting from Eq. [8] into Eq. [6] and tested whether or not
C < 1. This inequality is necessary to ensure that flow is entirely unsaturated. Similarly, we calculated
B and checked if
B > 0. The last limitation follows from the definition of the Kirchoff potential (i.e.,
< 0 would imply a negative hydraulic conductivity). Figure 2b shows
B/q and
C/q (curves 1 and 2) as functions of
a/q. As can be seen from the first curve, if we assume q = 1, then flow becomes supercritical near the leading edge C (if we interpret BMC as a cavity whose leading edge wets up and starts dripping), or a saturation zone will develop there (if we interpret BMC as a physical stone) at
a
0.28.
We now determine the second characteristic size of our obstacle, the coordinates of the most protruding point M in Fig. 1. Two conditions at this point should hold: uM = 0 and
M = 0. We write them as a system of two nonlinear equations in polar coordinates (O, r,
) as follows:
 | [9] |
By eliminating
Eq. [9] reduces to one equation in terms of r:
 | [10] |
Equation [10] was solved also with the FindRoot routine of Mathematica to obtain the root rM. From this we obtained
M, xM, yM, and the cavity prolateness e = 2xM/h. Figure 2c shows plots of xM, yM, and e (curves 13) as functions of
a/q. The curve for e shows that the obstacle becomes less ovalic as the ambient moisture content or infiltration rate increases. We recall that for saturated conditions a dipole placed against a unidirectional groundwater flow field generates a circular barrier as a separatrix (Kacimov and Nicolaev, 1992).
We used the ContourMap option of Mathematica and plotted in Fig. 3
the flow net for
a = 0.1, q = 1. Curve 1 indicates the cavity contour
= 0. As mentioned earlier, the interior of this contour is physically irrelevant to our problem. However, we have to verify that the Kirchhoff potential along the exterior of the separatrix satisfies the inequality 0 <
< 1. Curve 2 is the equipotential
= 1. All equipotentials
> 1 lie inside this critical loop and, as we can inspect from comparing curves 1 and 2, well inside the cavity contour. Curve 3 is the equipotential for
= 0. All equipotentials
< 0 are inside that critical loop and inside contour 1. Hence, all nonphysical equipotentials are well confined inside the separatrix. Curve 4 indicates
=
a. Between two branches of this equipotential and the trailing edge we have Philip's "dry shadow," that is, the zone where the moisture content is less than that of the incident infiltration flow field far from the obstacle. Curve 5 is the equipotential for
= 0.08.

View larger version (25K):
[in this window]
[in a new window]
|
Fig. 3. Flow net for a = 0.1, q = 1. Curve 1 is the cavity contour ( = 0). Curves 2 through 8 are the equipotential lines for = 1, 0, 0.1, 0.08, 0.12, 0.16, and 0.2, respectively. Curves 9 through 11 represent the stream lines = 0.1, 0.2, and 0.3, respectively.
|
|
The equipotential lines for
= 0.12, 0.16, and 0.2, shown in Fig. 3 by curves 6, 7, and 8, respectively, demarcate progressively wetter zones as compared with undisturbed infiltration. As was discovered by Philip et al. (1989a) and Warrick and Fennemore (1995), and as is obvious from our Fig. 3, these curves have two lobes beetling around and beneath the cavity. Two-branched streamlines
= 0.1, 0.2 and 0.3 are shown as curves 9, 10, and 11, respectively. Corresponding negative values of the stream function would produce streamlines (not shown in Fig. 3) symmetrical with curves 9 through 11. The branches of curves 9 through 11 contained within the separatrix have again only mathematical interest. We note that the branches 9 through 11 inside and outside the obstacle lie on the left and right side, respectively, with respect to the y axis in Fig. 3.
 |
CONCLUSIONS
|
|---|
We tried to reinvoke the Bystrov dipole solution by applying it to the quasilinear model (a Gardner soil). Infiltration in an infinite homogeneous soil is shown to be hampered by an impermeable ovalic obstacle, whose shape is not given a priori but reconstructed using Rankine's method. The more recent methods for solution of quasilinear flow with barriers apply not only to two spatial dimensions but also to more realistic three-dimensional barriers, which were studied by Fennemore and Warrick (1997) using Rankine's method.
The beauty of the dipole solution is that both the stream function and Kirchhoff potential are expressed explicitly and in closed form. We recall that in the case of a circular cylinder (Philip et al., 1989a) the solution was obtained as an infinite series. Warrick and Fennemore (1995) considered a more general case than here (a finite distance between the source and sink), but did not obtain an elegant stream function expression. To the best of our knowledge, the Bystrov solution brings about the first essentially two-dimensional infiltration regime for which the entire flow net is expressed in a rigorous, yet simple mathematical form. We recall also that even for a seemingly simpler singularity (e.g., a line sink or source) the expression for streamlines was not reported by Bystrov (1959) or Philip (1969), but obtained recently in series form with terms involving both MacDonald and trigonometric functions (Knight, 2006).
We plotted the vertical and horizontal sizes of the oval-shaped separatrix as functions of the incident infiltration intensity and the dipole strength. From the streamlines deflected by the obstacle, and isobaric contours delineating zones of constant moisture content, further geotechnical and agricultural implications can be inferred. Examples are the potential dripping of tunnels designed to store radioactive waste, or the moistening of certain areas of irrigated agricultural plots where soils are interspersed with horizontal worm holes.
An array of dipoles placed horizontally, as in Yourisov (1965), can be used to model identical obstacles set as a row. In that case the strength of all dipoles should be equal. Similarly to Warrick and Fennemore (1995), dipoles can be arranged in a double-periodic pattern, in which case the dipole strengths in each row are not necessarily equal (an elementary cell of this pattern can be a quincunx).
Because of the full mathematical analogy between the Gardner soil model and solute transport in confined aquifers (Batu, 1987), we can import all dipole solutions to equivalent contaminant transport problems. In that case one simply turns Fig. 1 90° and considers only the upper part. The
= 0 line will then represent an aquifer bed with an area BMC through which solutes cannot pass. Instead of the Kirchhoff potential we now have a certain solute concentration in the liquid phase of a fully saturated rock, while solutes will build up near the upstream part of the obstacle and concentrations will decrease near the back part (as compared with the background level in a unidirectional flow field far from the obstacle).
 |
APPENDIX
|
|---|
In this Appendix we derive the basic Eq. [5] and [6] and show how they follow from Bystrov (1959). We keep all of the notations of Bystrov (1959), but add subscript "B" to all variables from Bystrov (1959). We reference equations from Bystrov (1959) as [B1], [B2], etc., and his figures as Fig. B1, Fig. B2, etc. Equations in this Appendix are labeled as [R1], [R2], etc. We assume that the readers will have access to Bystrov (1959) (scanned images and further details can be obtained from the author at anvar@squ.edu.om or downloaded from http://www.ksu.ru/eng/anvar/ [verified 14 Feb. 2007]).
Bystrov introduced his velocity potential
B and the stream function
B using his relation [B4]. These two functions satisfy the partial differential equations [B5] and [B6], respectively. We shall immediately assume that the parameter µ, which determines the exponential growth of the Bystrov layer [B1'], has a value of 1/2. This is done solely to relate Philip's quasilinear model and notations in soil physics (i.e., our dimensionless coordinates) with Bystrov's hydrodynamics. We notice also that the orientation of Bystrov's Cartesian coordinates (xB, yB) coincides with ours (i.e., x, y* in our Fig. 1).
Bystrov introduced next other functions
B* and
B*, given by Eq. [B7] and [B8], respectively. We will need only Eq. [B8], in which we assume without any loss of generality that y0 = 0 (i.e., the origin of the coordinate system coincides with the dipole to be introduced later). As Eq. [B6'] shows,
B* satisfies the MHE
 | [R1] |
We now transform our stream function
, which, as stated in the main text, satisfies the SSADE:
 | [R2] |
wherewe stressthe orientation of the y axis is opposite to the one used by Bystrov, while the x axis is the same as that of Bystrov. While it may not make much sense to have different axes for the same problem, this was done here only to keep our notations in congruity with those by Bystrov and those common in soil physics.
As Philip commonly did, we introduce now the function H(x,y) =
exp[y/2]. Upon substitution into Eq. [R2], this leads to the MHE:
 | [R3] |
This shows that H satisfies the same MHE as
B*. Considering Bystrov's solution for
B as given by Eq. [B18], one may notice that the Bystrov dipole strength M is related to our Q as M/4
= q, while cos
= xB/r according to Bystrov's Eq. [B9] (we also recall that µ = 1/2). Note that the radial coordinate r in our and Bystrov's notation is the same. Thus, Bystrov's dipole stream function is
 | [R4] |
Next, from Eq. [R4] we obtain
 | [R5] |
We previously noted already that H coincides with
B*; therefore, from Eq. [R5] we conclude that a single dipole in terms of the Philip H function is described by
 | [R6] |
We now return to our own stream function and recall that xB = x to obtain
 | [R7] |
Bystrov did not study infiltration and consequently limited his interest to pure singularities. In our case we have a background descending infiltration flow field whose stream function in the absence of a cavity or tunnel is given by
 | [R8] |
In the quasilinear model we can assemble stream functions (and Kirchhoff potentials) using a linear combination of known functions. Hence we add Eq. [R7] and [R8] and obtain Eq. [3] in the main text.
To derive Eq. [6] in the paper, consider the first part of Eq. [4] as follows:
 | [R9] |
This equation may be integrated to yield
 | [R10] |
where
a is the Kirchoff potential at infinity (far from the dipole). Attempts to integrate the result of the differentiation in Eq. [R10] directly using Mathematica failed since Mathematica does not evaluate Eq. [R10] in an explicit fashion. We instead changed the order of integration and differentiation of Eq. [R10] to obtain with [R7] the expression
 | [R11] |
The integration in Eq. [R11] can be done using Mathematica (any other package can be used as well) to yield 2exp[y/2]K0(r/2). Substituting the calculated integral back into Eq. [R11] and differentiating with respect to y finally leads to Eq. [6] in the main text.
 |
ACKNOWLEDGMENTS
|
|---|
This study was supported by project CL/SQU-UAEU/0/3/02, Sultan Qaboos University. Helpful comments by three anonymous referees, J. Knight, J. Nieber, and R. van Genuchten, are appreciated.
 |
REFERENCES
|
|---|
- Bakker, M., and J.L. Nieber. 2004. Analytic element modeling of cylindrical drains and cylindrical inhomogeneities in steady two-dimensional unsaturated flow. Available at www.vadosezonejournal.org. Vadose Zone J. 3:10381049.[Abstract/Free Full Text]
- Batu, V. 1987. Introduction of the stream function concept to the analysis of hydrodynamic dispersion in porous media. Water Resour. Res. 23:11751184.
- Bers, L., and A. Gelbart. 1944. On a class of functions defined by partial differential equations. Trans. Am. Math. Soc. 56:6793.[CrossRef]
- Borzykh, A.A., and G.P. Cherepanov. 1978. Plain problem in the theory of convective heat and mass transfer. (In Russian.) Prikladnaya Mathematika I Mekhanika 42:848855.
- Bystrov, K.N. 1956. On determination of sources, vortices and multipoles in curved fluid layers of varying thickness. (In Russian.) Trans. Moscow Distr. Pedagog. Inst. 75(43):N3.
- Bystrov, K.N. 1959. On two-dimensional steady flows in a layer with exponentially varying thickness. (In Russian.) Trans. Moscow Distr. Pedagog. Inst. 75:3159.
- Concer, D.B. 1959. Heat flow towards a moving cavity. Q. J. Mech. Appl. Math. 12:222232.[Abstract/Free Full Text]
- Duffin, R.J. 1971. Yukawan potential theory. J. Math. Anal. Appl. 35:105130.[CrossRef]
- Dullien, F.A.L. 1992. Porous media. Fluid transport and pore structure. Academic Press, San Diego, CA.
- Elizarov, A.M., N.B. Ilinskiy, and A.V. Potashev. 1997. Mathematical methods of airfoil design. Akademie Verlag, Berlin.
- Fennemore, G.G., and A.W. Warrick. 1997. Simulation of unsaturated water flow around obstructions: Three-dimensional Rankine bodies. Adv. Water Resour. 20:1522.
- Frenkel, Ya.I. 1949. Theory of phenomena of atmospheric electricity. (In Russian.) GITTL, Leningrad.
- Gardner, W.R. 1958. Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Sci. 85:228232.
- Kachinskii, N.A. 1970. Soil physics. (In Russian.) Vysshaya Shkola, Moscow.
- Kacimov, A.R. 2000. Circular isobaric cavity in a descending unsaturated flow. J. Irrig. Drainage ASCE. 126:172178.[CrossRef]
- Kacimov, A.R. 2003. Unsaturated quasi-linear flow analysis in V-shaped domains. J. Hydrol. 279:7082.[CrossRef]
- Kacimov, A.R. 2006. Analytic element solutions for seepage towards topographic depressions. J. Hydrol. 318:262275.[CrossRef]
- Kacimov, A.R., and A.N. Nicolaev. 1992. Steady seepage near an impermeable obstacle. J. Hydrol. 138:1740.[CrossRef]
- Knight, J.H. 2006. Analytical solutions for linearised unsaturated soil water flow. p. 183194. In Proc. 5th Int. Conf. on the Analytic Element Method, Kansas State University, Manhattan. 1417 May 2006.
- Kulabukhova, I.I. 1959. Two problems of transient seepage at incomplete soil saturation. Izvestiya Akad. Nauk SSSR. (In Russian.) Mekhanika I Mashinostroenie. N3:196201.
- Novikov, P.A., L.Y. Lyubin, and V.I. Novikova. 1991. Flow and heat and mass transfer in slot systems. (In Russian.) Nauka i Tekhnika, Minsk.
- Philip, J.R. 1969. Theory of infiltration. Adv. Hydrosci. 5:215296.
- Philip, J.R. 1990a. Conjectures on certain boundary-layer equations and natural coordinates. Proc. R. Soc. Lond. A 428:307324.
- Philip, J.R. 1990b. How to avoid free boundary problems. p. 193207. In K.H. Hoffman and J. Sprekels (ed.) Free boundary problems: Theory and applications. Longman, New York.
- Philip, J.R. 1992. What happens near a quasi-linear point source? Water Resour. Res. 28:4752.[CrossRef]
- Philip, J.R., J.H. Knight, and R.T. Waechter. 1989a. Unsaturated seepage and subterranean holes: Conspectus, and exclusion problem for circular cylindrical cavities. Water Resour. Res. 25:1628.
- Philip, J.R., J.H. Knight, and R.T. Waechter. 1989b. The seepage exclusion problem for parabolic and paraboloidal cavities. Water Resour. Res. 25:605618.
- Polubarinova-Kochina, P.Ya. 1977. Theory of ground-water movement. (In Russian.) Nauka, Moscow.
- Polubarinova-Kochina, P.Ya., and I.I. Kulabukhova. 1959. On transient seepage at incomplete saturation of soil. (In Russian.) Izvestiya Akad. Nauk SSSR, Mekhanika i Mashinostroenie N2:5765.
- Strack, O.D.L. 1989. Groundwater mechanics. Prentice Hall, Englewood Cliffs.
- Warrick, A.W. 2003. Soil water dynamics. Oxford Univ. Press, Oxford, UK.
- Warrick, A.W., and G.G. Fennemore. 1995. Unsaturated water flow around obstructions simulated by two-dimensional Rankine bodies. Adv. Water Resour. 18:375382.[CrossRef]
- Warrick, A.W., and J.H. Knight. 2002. Two-dimensional unsaturated flow through a circular inclusion. Water Resour. Res. 38:181186.
- Wolfram, S. 1991. Mathematica. A system for doing mathematics by computer. Addison-Wesley, Redwood City, CA.
- Youngs, E.G., A.R. Kacimov, and Yu.V. Obnosov. 2004. Water exclusion from tunnel cavities located in the saturated capillary fringe with uniform precipitation flowing to a water bearing substratum. Adv. Water Resour. 27:237243.[CrossRef]
- Yourisov, V.A. 1965. On one problem of flow in a plane-variable layer. (In Russian.) Trans. Moscow Distr. Pedagog. Inst. 142(5):6991.