Transp Porous Med (2018) 121:643–664
https://doi.org/10.1007/s11242-017-0978-x

Steady Flow from an Array of Subsurface Emitters:
Kornev’s Irrigation Technology and Kidder’s Free
Boundary Problems Revisited
A. R. Kacimov1· Yu. V. Obnosov2· J. Šimůnek3

Received: 18 July 2017 / Accepted: 23 November 2017 / Published online: 4 December 2017
© Springer Science+Business Media B.V., part of Springer Nature 2017

Abstract Kornev’s (Subsurface irrigation, Selhozgiz, Moscow-Leningrad, 1935) subsurface irrigation with a periodic array of emitting porous pipes is analytically modeled as asteady potential Darcian flow from a line source generating a phreatic surface. The hodograph method is used. The complex potential strip is mapped onto the triangle of theinverted hodograph. An analogy with the Deemter (Theoretische en numerieke behandelingvan ontwaterings-en infiltratie stromings problemen (in Dutch). Theoretical and numericaltreatment of flow problems connected to drainage and irrigation. Ph.D. dissertation, Delft
University of Technology, 1950) drainage problem and Kidder (J Appl Phys 27(8):867–869,1956) free-surface flow toward an array of oil wells underlain by a “wavy” oil–water interface
is drawn. For a half-period of Kornev’s flow, the “wavy” phreatic surface has an inflection
point. The “waviness” of the phreatic surface is controlled by the spacing between emitters,
the strength of line sources, and the pipe pressure and radius. Numerical modeling with
HYDRUS involved two factors which constrained the saturated–unsaturated flow: the positive pressure head at the outlet of the modeled domain and lateral no-flow boundaries, witha qualitative corroboration of analytical solutions for potential (fully saturated) and purely
unsaturated flows. HYDRUS is also applied to a generalized Philip’s regime of an unsaturated
flow past a subterranean hole, which is impermeable at its top and leaks at the bottom.

B A. R. Kacimov

anvar@squ.edu.om; akacimov@gmail.com
https://www.squ.edu.om/agr/Departments/Soils-Water-and-Agricultural-Engineering/Faculties/ProfAnvar-Kacimov
Yu. V. Obnosov
yobnosov@kpfu.ru
J. Šimůnek
jsimunek@ucr.edu

1

Department of Soils, Water and Agricultural Engineering, Sultan Qaboos University, Al-Khod 123,
PO Box 34, Muscat, Sultanate of Oman

2

Institute of Mathematics and Mechanics, Kazan Federal University, Kazan, Russia

3

Department of Environmental Sciences, University of California Riverside, Riverside, CA, USA

123

644

A. R. Kacimov et al.

Keywords An array of subsurface irrigation emitters · Line sources versus oil production
wells · “Wavy” phreatic surface · Conformal mappings of hodograph-complex potential
domains · HYDRUS-2D simulations of “backwater” bottom conditions“Singularity is almost invariably a clue.”
Sherlock Holmes

Abbreviations
BVP
K-35
PP
SI

Boundary value problem
Reference to the book Kornev, V.G., 1935. Subsurface irrigation. Selhozgiz,
Moscow-Leningrad (in Russian)Porous pipe
Subsurface irrigation

1 Introduction
Kornev (1935) (hereafter abbreviated K-35) proposed and tested a new technology of subsurface irrigation (SI) that uses an array of porous pipes (PP) that are buried at the same
depth d p (practically ranging from 0.2 to 1 m), with a period of 2L (0.5–2 m), and a constantdischarge 2Q (up to several tens of l/h/m) per a unit length in the direction of the pipe axis.
Fig. 1 (see also the original Fig. 38 in K-35) shows a vertical cross section of one period of
this SI system. A homogenous soil has the saturated hydraulic conductivity K 1 . The emitting
PP of a radius rC (2–5 cm) is modeled as a line source at point A. In this paper, we study
steady flows at a constant, positive, and relatively high pressure in PP (the pressure heads instandard SI range between a few tens of centimeters to a few meters).
Several decades after the publication of K-35, numerous SI projects used PP in France,
Iran, Israel, Japan, Oman, Russia, Spain, Saudi Arabia, the USA and other countries witharid/semi-arid agriculture (see, e.g., Al-Rawahy et al. 2004; Ashrafi et al. 2002; Communarand Friedman 2015; El-Nesr et al. 2014; Honari et al. 2017; Iwama et al. 1991; Kato andTejima 1982; Lazarovitch et al. 2005; Moniruzzaman et al. 2011; Shani et al. 1996; Siyal
and Skaggs 2009; Tanigawa et al. 1988; Warrick and Shani 1996; Yabe and Tanigawa 1992).
Kacimov and Obnosov (2016, 2017) and Obnosov and Kacimov (2017) revisited K-35
and modeled an isolated subsurface emitter in homogeneous and two-layered soils, at both
positive and negative pipe pressures. If the pressure in PP is positive and permeability of walls
is not too small, in particular, for mole emitters which have no walls (Bobchenko 1957), thenseepage from an array of these emitters into the soil forms a “wavy” water table F2 FF 1 witha capillary fringe, which is a “subdued replica” of F2 FF 1 (see Kacimov 2006). Similarly,
undulating phreatic surfaces emerge above periodic tile drains, although the troughs andcrests appear above and between the drains (Deemter 1950; Ilyinsky and Kacimov 1992).
The soil surface S2 S1 and curve F2 FF 1 in Fig. 1 sandwich a tension-saturated and unsaturated zone. The plant roots uptake moisture from this zone. Evaporation from the soil surface,although relatively small (as compared with surface irrigation), is undesirable because it
reduces water use efficiency and causes secondary salinization, among other negative consequences (see K-35 for details). The effective thickness D of the unsaturated zone (20–25 cm
in Kornev’s SI projects) in Fig. 1 is high enough to serve as an anti-evaporation sheath. This

123

Steady Flow from an Array of Subsurface Emitters: Kornev’s…

645

Fig. 1 One period of subsurface irrigation from an array of porous pipes maintained at a positive pressure
and Kidder’s (1956) free boundary oil flow. The mathematical source is placed at a depth, d, under the summit
F of a “wavy” phreatic surface F2 FF 1 ; m/2 is an amplitude of the “wave” as compared with an average
horizon W2 W1 ; the flow rate from a constant total head emitter is 2Q; the emitter’s radius is rC . The soil
above W2 W1 is of a thickness D, the top soil layer of thickness B contrasts with the main soil in texture,
hydraulic conductivity, and capillarity

partially saturated layer also serves as a thermal insulator for the plant roots, provided they
penetrate deeply from the superheated soil surface. This dual “smart” wetting and heat-stressmitigation of the unsaturated layer has been confirmed by numerous data on improved harvest/biomass/plant morphologies for a variety of crops cultivated by Kornev in the USSR
and France (K-35).
In Sect. 2 of this paper, we consider a saturated 2-D steady flow from a source A cappedby a free boundary F1 FF 2 . Correspondingly, we solve a free boundary value problem (BVP)
for a potential flow that is laterally confined by the rays F2 D2 and F1 D1 in Fig. 1. For thispurpose, we use complex variables (the hodograph method), which gives the positions ofpoints F and F1 and therefore the value of D (the locus of horizon W1 W2 ) for a given d p
and emitter characteristics. This free BVP is mathematically equivalent to Deemter’s (1950)
tile drainage problem with a phreatic surface and Kidder’s (1956) coning problem of oil flow
to an array of producing wells (sinks) in an inclined formation with a subjacent stationary“wavy” oil–water contact (a free boundary as a sharp interface between a moving oil andstatic groundwater in the oil formation). Our Fig. 1 also sketches Kidder’s formation, forwhich the angle of a dip is a scaling constant. Kidder only studied a limiting case of a criticalproduction rate, for which the sink strengths assume the largest value for which a stable front
free boundary can exist without water being drawn into the well, other factors remain thesame. Within one period, Kidder’s free surface has a cusp–crest at point F beneath the welland zero-slope troughs at points F1 and F2 . Kornev’s flow from the sources in Fig. 1 is amathematical inversion of Deemter’s and Kidder’s flow to sinks.

123

646

A. R. Kacimov et al.

For a general (subcritical) flow to sinks and SI flow from Kornev’s sources, the freeboundary within one flow period (Fig. 1) has two inflection points, N1 and N2 , where theslope of the free boundary is maximal.
In Sect. 3, we use the HYDRUS-2D finite element (FE) code (Šimůnek et al. 2016) andsimulate a saturated–unsaturated system shown in Fig. 1. We illustrate that numerical solutions of Richards’ equation are in agreement with the analytical results in Sect. 2. We also
simulate unsaturated flows past subterranean holes with a small water level at the hole bottom, combining unsaturated and saturated regimes from Philip et al. (1989) and Warrick andZhang (1987). These HYDRUS applications are motivated by the design and operation ofmole drains-emitters and tunnels for which the boundaries of contact between the ambient
soil (rock) and porous-fractured medium can be either impermeable to the descending infiltration due to the capillary barrier phenomenon or permeable (e.g., Wang and Bodvarsson
2003).
We assume that all flows are Darcian, one-phase (vapor/gas and solute motions, depositions/dissolutions on/of the soil matrix are neglected), and isothermal and that the bare(rootless) stratified soils are isotropic.
We answer the following questions:(a) How do values of L , K 1 , the pressure in PP, and the flow rate, 2Q, in the saturated
fragment of the free boundary flow affect the position of F2 FF 1 in Fig. 1? In particular,
how large is the trough-crest difference, m, in the ordinates of points F and F1 ?
(b) What are the results of HYDRUS simulations for flow in a positive-pressure zone consisting of a “mound” created by seepage from a cascade of positive-pressure PP on the
top of an already existing unconfined aquifer, which is hydraulically drainable to a thickhighly permeable substratum?
(c) Is there agreement between the numerical results for these “nonstandard” drainage
conditions of the irrigated topsoil and the results obtained by simplified analytical models
and solutions, which ignore capillarity and consider flow domains not confined frombelow?

2 Steady 2-D Saturated Flow from Array of Sources
Strack (1989) modeled an arbitrary number of interfering emitters and drains (line sources
and sinks), which generate a positive pore pressure zone around the emitters (line sources),with the flow domain bounded by a phreatic surface. Kacimov and Obnosov (2016) revisited
the Riesenkampf solution for a single emitter (see Polubarinova-Kochina 1962), used theVedernikov (1939) model of a capillary fringe and tension-saturated flow, and modeled asingle PP. This model just rescales the pressure head at F2 FF 1 in Fig. 1 by a summand − pc ,in comparison with the water table pressure p = 0. Mathematical models in the Strack(1989) and Riesenkamf solutions are the same. Therefore, without any loss of mathematicalgenerality, we ignore the fringe in Fig. 1.
To connect our analytical solution with Kidder’s (1956), we assume that F2 FF 1 is a
streamline, i.e., evapotranspiration losses from the water table are neglected. The “saturated
bulbs” around the line sources of the whole array (laterals) of SI in Fig. 1 coalesce, and—sufficiently deep under the sources—a descending, 1-D, positive-pressure flow with a Darcian
velocity vD = Q/L occurs, where 2Q is the strength of the line source in Fig. 1. Similarflows are shown in Fig. 7.38 of Strack (1989).

123

Steady Flow from an Array of Subsurface Emitters: Kornev’s…

647

Fig. 2 A vertical cross section for a half-period of the flow domain in Fig. 1. A “wavy” phreatic surface has
an inflection point N1 (a), the complex potential plane (b), and the hodograph plane (c), corresponding to thephysical flow domain in a

2.1 Conformal Mappings
Due to symmetry, we follow Deemter (1950) and Kidder (1956) and consider only the righthalf of the period, a physical domain G z , as shown in Fig. 2a. The origin of the Cartesian
coordinates (x, y) coincides with the apex F of the free surface. At this point, the Darcian total
hydraulic head h(x, y) = 0, h = p + y, where p is the pressure head. The total (hydraulic)head in the pipe is H p = const > 0. The source is located at z = − id, where d = const > 0
is to be found. Obviously, a posteriori, we must check that d < d p (Fig. 1). Otherwise, theassumed flow scheme is not valid and, practically, the water table seeps out as a “wet spot”
on S1 S2 in Fig. 1 (the situation to be avoided in SI, see K-35).
We introduce the complex potential w = φ + iψ, where φ = − K 1 h is the velocity
potential and ψ is the stream function. At point F, the potential φ = 0, and ψ = 0 along
AD1 . Therefore, ψ = Q along AFN 1 F1 D1 . We also introduce a complex Darcian velocity
V = u + iv. The complex potential domain G w , which corresponds to G z , is then the stripof a width Q (Fig. 2b), and the hodograph domain G V is the right half-plane with a circular
cut FN 1 F1 (Fig. 2c). At the inflection point N1 (the tip of the cut in Fig. 2c), the magnitudeof the Darcian velocity vector attains its maximum along FN 1 F1 . At this point, the angle

123

648

A. R. Kacimov et al.

πγ of this vector with the x-axis is also maximal. The complex velocity at point D (infinity)
is given and is V D = − iv D = − i Q/L. At a certain depth, dh , flow becomes almost 1-D
vertical, i.e., a horizontal line Dl Dr in Fig. 2a (a dashed line) is almost an equipotential line,
which is imaged by the corresponding almost vertical segment in Fig. 2b (dashed line). Along
Dl Dr in Fig. 2a, the pressure head is also almost constant. We note that in Kidder’s (1956)
solution, the incident Darcian velocity (far upstream of his production wells) is also constant(actually, a scaled parameter a in Kidder’s parametrization). McCarthy (1994) studied flows toa solitary line sink with free boundaries, but with an alternative to Deemter (1950) and Kidder(1956) condition at infinity, viz. McCarthy’s flow domains are laterally unbounded, whereas
Deemter–Kidder’s and ours are vertically unbounded. First, the circular quadrilateral G V ismirror-imaged with respect to the real axis. The corresponding domain G ω (Fig. 3a) is invertedwith respect to the coordinate origin to obtain the domain G ω1 (Fig. 3b), which has a straightcut FN 1 F1 . The function ω1 (z) = 1/(u− iv) = dz/dw is holomorphic. PolubarinovaKochina (1962) gives full details about the hodograph method, which irrigation engineers(e.g., Swamee and Chahar 2015) use for designing of irrigation channels. Hodographs similar
to one in Fig. 2c and their inversions are used for analytical studies of flow from rectangular
irrigation-MAR channels (e.g., Choudhary and Chahar 2007).
We map G w onto an auxiliary half-plane G ζ (Fig. 3c) by the functionw(ς) = i Q −

Q
ς +μ
log
.
π
1+μ

(1)

Similarly, as for Strack’s and Riesenkampf’s emitters-drains, flow determined by Eq. (1) isradial 1-D in the vicinity of point A in Fig. 2a. The affixes −μ, ν of points D1 , N1 (Fig. 3c)will be found later. At point C, i.e., the apex (the contact of the exterior of the PP wallwith the ambient soil), the hydraulic head HC is easily calculated from H p , the given wall
thickness (approximately 1 cm in K-35) and the permeability (see Kacimov and Obnosov
2016, 2017). In the simplest situation of a mole drain (an emitter with no seepage-impeding
wall, Bobchenko 1957), HC = H p . Then, at point C, w(c) = −K 1 HC + i Q, and affix c is
c = (1 + μ) exp(π K 1 HC /Q) − μ.

(2)

We map G ω1 onto G ζ using the function
ω1 (ς) =i 2ς + cosecπγ − 1
i
.
−
√
2K 1
K1
ς(ς − 1)

(3)

The angle πγ (Fig. 3a) will be found later. The mapping function (3) is the inversion of theelementary function:
1 + ς1 (ω1 )(K 1 ω1 + i)2 − cot2 πγ
1 − sinπγς1 (ω1 ) =,ς(ω1 ) =. (4)K 1 ω1 + i
2
1 − sinπγ ς1 (ω1 )In Eq. (4), ς1 (ω1 ) maps G ω1 onto the upper half of the ς1 -plane, with the correspondenceof points (F1 , F, A) → (− 1,1,1/ sin πγ ), and ς(ς1 ) is the automorphism of the upper
half-plane, which gives the correspondence of points (F1 , F, A) → (0, 1, ∞).
Note that ω1 = (cotanπγ − i)/K 1 at point N1 in the ω1 -plane. From Eq. (4), it thenfollows that the affix ν = (1 − sin πγ )/2.
Using Eqs. (1) and (3), we get
⎛ ς
⎞
ς

dw
ς + μ⎠i Q ⎝ 2u + cosecπγ − 1z(ς) = ω1 (u)du − 2 log
du = −. (5)√du2K 1 π1+μ(u + μ) u(u − 1)
1

123

1

Steady Flow from an Array of Subsurface Emitters: Kornev’s…

649

Fig. 3 Mirrored hodograph domain (a) and its inversion with respect to the origin (b); an auxiliary halfplane (c); Deemter–Kornev’s flow with evaporation from the phreatic surface (d); the hodograph plane for anevaporating free surface (e)The last integral is evaluated in elementary functions:
ς
√
1

√

du
ς + ς − 1 = log(ς + μ) + 2 log= 2 logu(u − 1)

ς
+
ς +μ

ς −1
,
ς +μ

and
ς
1

du
2
log
= √
√
μ(μ + 1)
(u + μ) u(u − 1)

That transforms (5) intoz(ς) =√
√(μ + 1)ς + μ(ς − 1)
.
√
ς +μ

√
√
ς(μ + 1) + μ(ς − 1)(cosecπγ − 1 − 2μ)
log
√
√
μ+ς
μ(μ + 1)√√
(μ + 1)ς + (μ + 1)(ς − 1).
+ 2 log
√
ς +μ
−iQ
K1π



(6)

To find γ and μ, we use the coordinates of points A and C:z(∞) = − id, and z(c) = − i(d − rC ).

(7)

123

650

A. R. Kacimov et al.

In accordance with Eqs. (6) and (7), we obtain√(cosecπγ − 1 − 2μ)Qlogd=μ + 1 + μ + log[4(1 + μ)] .
√K1πμ(1 + μ)From the second condition (7) and (2), (6), we get√√c(μ + 1) + μ(c − 1)Q
(cosecπγ − 1 − 2μ)log
d = rC +
√
√
μ+c
K1π
μ(μ + 1)√√c − 1 + c − HC ,+ 2 log(8)(9)

where μ = −ς(−i/v D ) = −ς(−i L/Q).
We introduce dimensionless quantities z∗ = z/L , d ∗ = d/L , m ∗ = m/L , HC∗ =HC /L , rC∗ = rC /L , w∗ = w/(K 1 L), v ∗D = v D /K 1 and Q ∗ = Q/(K 1 L) and drop the“*” for the sake of brevity.
From (3), we deriveQ − 1 + (Q − 1)2 + Q 2 cot2 πγ
1 − sin πγ

.
(10)
μ=
2
1 − Q + sin πγ (Q − 1)2 + Q 2 cot2 πγEliminating d from Eqs. (8) and (9) and taking into account (2) and (10), we obtain thefollowing nonlinear equation to determine γ :√√ 
μ + 1 + μ exp(π HC /2Q)
Q (cosecπγ − 1 − 2μ)
HC − rC +log √√πμ(μ + 1)c + μ(exp(π HC /Q) − 1)
√


2 μ+1
+ 2 log √
= 0.
(11)
√
c−1+ cWe used the FindRoot Mathematica (Wolfram 1991) routine to solve this and subsequent
nonlinear equations.
The separation of the real and imaginary parts in Eq. (6) gives a parametric equation of
the free boundary:
cosecπγ − 1 − 2μarccos2 arccos ξ +√μ(1 + μ)x(ξ ) =Qπy(ξ ) =

ξ +μQlog
, 0 < ξ < 1.
π1+μ(1 + μ)ξ
μ+ξ

,
(12)

The depth of point F1 is m = −Imz(0) = −y(0), and the second Eq. (12) gives
m=

Q
1+μ
log
.
π
μ

(13)

We note that along with the hodograph method, one can use the complex potential andZhukovsky function and the Zhukovsky-Chaplygin method (see e.g. Anderson 2013; Strack1989).
We used the ParametricPlot Mathematica routine to plot the free surfaces FF 1 . For thesake of further comparisons with HYDRUS, in Fig. 4a, we show dimensionless free surfaces
and loci of sources for the following dimensional input parameters: rC = 0.025 m, HC =0.85 m, K 1 = 0.25 m/day, fixed L = 1.0 m with varying Q = 0.265, 0.275, 0.285, 0.295,
0.305 0.315 m2 /day (curves 1–6) with the corresponding position of the source indicated bya solid dot. In Fig. 4b we fixed Q = 0.3 m2 /day and plotted three free surfaces for varying L:

123

Steady Flow from an Array of Subsurface Emitters: Kornev’s…

651

Fig. 4 Phreatic surfaces and loci of sources in dimensionless coordinates for Kornev’s subsurface irrigation.
Dimensional input parameters for a are rC = 0.025 m, HC = 0.85 m, K 1 = 0.25 m/day, L = 1.0 m, and
Q = 0.265, 0.275, 0.285, 0.295, 0.305, and 0.315. For b rC = 0.025 m, HC = 0.85 m, K 1 = 0.3 m/day,
Q = 0.3 m2 /day, and 1 L = 0.8 m, 2 L = 1.2 m, and 3 L = 2.4 m(1) L = 0.8 m, (2) L = 1.2 m, (3) L = 2.4 m (curves 1–3). The corresponding dimensionlessd and m are: (1) d = 0.87, m = 0.16, (2) d = 0.50, m = 0.58, (3) d = 0.25, m = 0.81.
As is evident from Fig. 4a, an increase in Q makes the phreatic surface steeper. Also, therange of variation of Q is quite narrow (we recall that HC and L are fixed for the examples
computed in Fig. 4a). Curve 1 in Fig. 4a illustrates the case of an “almost confined” flow, forwhich FF1 is an almost horizontal straight line. A mathematically equivalent problem for a
confined flow to horizontal sinks placed above a horizontal bedrock and fed from a ponded
soil surface has been analyzed by Vedernikov (1939) (see also Polubarinova-Kochina 1962).
Curve 6 in Fig. 4a corresponds to an “almost critical” flow regime. In case of L > Q/K 1 ,the neighboring emitters do not interfere with each other and one arrives at Riesenkampf’s
critical regime for a solitary emitter (Polubarinova-Kochina 1962).
For comparisons, we used Kidder’s Eqs. (18) and (19) as the degree of waviness of his freeboundary. We also computed a dimensionless area Ar between FF 1 and a horizontal oil–watercontact (a dashed line in Kidder’s flow, our Fig. 1). For Q → ∞ in Kidder’s solution, when
his free surface degenerates into a cycloid, the maximum waviness parameters are m →2/π, Ar → 0.15 (obviously, lim Q→0 (m, Ar) = (0,0), which corresponds to a horizontal
oil–water interface). Therefore, even in the limit of the “most perturbed equilibrium”, theKidder cycloid is only moderately wavy. Kidder’s “critical” flow regime gives the mostundulating free boundary because this regime is a precursor to water breakthrough into theoil wells i.e. at subcritical oil flows the free boundary (with inflection points) differs evenless from a straight oil–water interface.
In Fig. 5 we use the ParametricPlot, Re, and Im routines of Mathematica applied toEq. (6) and show the flow net for rC = 0.025 m, HC = 0.85 m, K 1 = 0.25 m/day, L =1 m, and Q = 0.295 m2 /day, for which m = 0.32 m, γ = 0.155, μ = 0.74. In Fig. 5,
the streamlines are plotted by solid lines for ψ = 0.11 j Q, where j = 1,2,…9; as thephreatic surface is shown by a dashed-dotted line; the equipotential lines are plotted for
φ = 0.05 j, j = − 2, −1, 0, 1, . . .8 by dashed lines. The lowermost equipotential almost
coincides with a horizon Dl Dr depicted in Fig. 2a, b.

123

652

A. R. Kacimov et al.

Fig. 5 A flow net for
rC = 0.025 m, HC = 0.85 m,
K 1 = 0.25 m/day, L = 1.0 m,
and Q = 0.295 m2 /day. The
streamlines are for ψ = 0.11 j Q
where j = 1, 2, . . . , 9 (solidlines) and the equipotential lines(dashed lines) are for φ = 0.1 j,
j = − 2, − 1, 0, 1, . . ., 5

For the situation of a fixed L and varying Q in Fig. 4a, the mathematical problem is solvable in the range Q min < Q < Q max (curves 1 and 6 illustrate the regimes close to thelower and upper bounds of Q). Physically, for a given circular PP operated at a given pipepressure in a given soil, i.e., for a given triad K 1 , rC , and HC , the value of Q is a part ofthe solution (Kacimov and Obnosov 2016, 2017). Modeling of a real PP by a mathematicalsource bypasses the involvement of the real geometry of the emitter contour into the analysis.
Namely, the equipotential lines close to the source are believed to be “almost circular” andone of these small-radius (rC ) “circles” is considered as a PP-soil contact boundary (see PK77, Fujii and Kacimov 1998). If the “deep drainage” conditions at infinity (point D1 –D2 ) in
Fig. 1 are known by specification of the pore pressure for saturated flows (as for example, in
Kidder 1956; Kacimov and Obnosov 2016) or moisture content for unsaturated flows (as forexample, in Kacimov and Obnosov 2017), then Q for a “circular” (or any other shape) emitteris determined by the above-mentioned triad. Vedernikov (1939) analyzed the shape of realequipotential lines around mathematical sinks-sources and found when they deviate from prescribed circles, although this affects flow in the very vicinity of the mathematical singularities.
In the conceptual model of Fig. 2a, the pore pressure at infinity (point D1 ) is not known,
which is often common in practical SI when farmers do not have details about the subsurface

123

Steady Flow from an Array of Subsurface Emitters: Kornev’s…

653

deeper than a couple of meters. Therefore, in our model with the “infinity” of z in Fig. 2a, Q
is “conditioned” by the triad. This “conditioning” will be rectified in Section C, where we
consider finite modeling domains with prescribed pressure boundary conditions along theinlet and outlet of a “stream tube” (flow domain) of the Richards’ equation.
With the emphasized mathematical commonality between our and Deemter–Kidder’s solutions, there are some physical differences. Namely, the solvability of Kidder’s (1956) problem
is determined by his Eq. (16); for a given difference in the ordinates of the sink and points
F1 and F2 (Kidder’s parameter b), the oil production rate Q cannot exceed a critical value.
The Darcian velocity in the oil far upstream of Kidder’s wells (Kidder’s parameter a) can
be arbitrarily high at a given distance between neighboring wells. The infinite increase in
Kidder’s a requires an infinite increase in his b.
In Kornev’s flow (Fig. 2a), Q can be mathematically arbitrarily high if HC is allowed to riseto infinity for a given emitter radius (of course, in this case, the source in Fig. 1 should be deep
enough such that the free surface does not “bounce” the soil surface from below). Kornev’s
Q and v D can be, however, bounded from below. As in Strack’s (1989) “free drainage” flows
from a cluster of sinks and sources, the Darcian velocity in the saturated plume at infinity(our v D ) is K 1 . If the pore pressure at infinity (point D1 in Fig. 2a) is atmospheric (we recall
that capillarity in this Section is neglected), then the criticality of Kornev’s flow commencesat Q = Q m such that at Q < Q m the saturated plumes of neighboring emitters in Fig. 1
are disconnected (the dotted line in Fig. 1 shows one isolated plume for this case). In otherwords, at this pressure, for any Q < Q m , Kornev’s flow degenerates into Riesenkapmf’s
flow (Kacimov and Obnosov 2016).

3 HYDRUS-2D Simulations
In this Section, we use HYDRUS (2D/3D) (Šimůnek et al. 2016) and compare numerical(FE) and analytical solutions. In all transient flow simulations, time, t, is in days. We use
several default HYDRUS options: iteration criteria, time step controls, an initial (at t = 0)
pressure head of − 100 cm throughout the modeled 2D domains, the VGM capillary pressure
and phase permeability functions, and the soil catalog.

3.1 Flow From Emitters
Ashrafi et al. (2002), El-Nesr et al. (2014), Siyal and Skaggs (2009), Wang et al. (2017), among
many others (see also Šimůnek et al. 2016), used HYDRUS for modeling SI from a PP either
as a solitary emitter, i.e., without overlapping of flows generated by laterally neighboring
emitters, or as a system of parallel emitters. In their conceptual HYDRUS models, they
assumed a “free drainage” or specified negative pressure head boundary condition at thebottom of the flow domain. Hanson et al. (2008) modeled in HYDRUS SI with a shallow
water table, the case analytically investigated by Philip (1989). Both Hanson et al. (2008)
and Philip (1989) assumed a zero-pressure (water table) boundary condition at the outlet of
their flow domain that takes place if drainage water flows out of the domain without raising
the water table, assuming that sufficient natural groundwater drainage occurs. Hanson et al.
(2008) solved Richards’ equation for an infinite number of emitters 2L apart as in Fig. 1
with no-flow boundary conditions on both the left and right sides (at a distance of L) of the
HYDRUS domains. Here, we also consider a half-period of an array of emitters, but our
flow domain is bounded by a horizontal positive-pressure interface I1 I2 between a soil layer
where the emitter is placed and a highly permeable “water-bearing” substratum, into which

123

654

A. R. Kacimov et al.

the emitted water drains. Vedernikov (1939) (see also PK-77, Bouwer 1978) distinguishedpre-irrigation hydrogeological conditions with this “groundwater bearing” hydrostratigraphic
unit, which is also called the “main”, confined aquifer. The aquifer creates a “backwater”
condition above I1 X H (Fig. 6), such that in the soil layer with the emitter or another irrigation
source (Ilyinsky and Kacimov 1991) the water table “bulges” even if the flow rate from theemitters is small. In other words, even for Kornev’s “negative pressure” emitters the originallyflat water table “waves” due to the accretion from the superjacent soil.
We start with a homogeneous rectangle 100 cm × 200 cm, shown in Fig. 6a and use the“2D - Simple” (rectangular) type of geometry of HYDRUS. Finite element discretization is
as follows: 40 × 40 nodes in both horizontal and vertical directions, with 160 1-D and 3200
2-D finite elements. The soil is loam with K 1 = 24.96 cm/day. We model the right half ofthe subsurface emitter by a rectangle ACC 1 C2 C3 A of a height and width of 5 and 2.5 cm,
respectively, point A having X H = 0, Z H = L v = 112.5 cm in the HYDRUS Cartesian
coordinate system. The interior of this rectangle consists of sand of the hydraulic conductivity
K s = 713 cm/day. This large conductivity ensures that the total head in the emitter rectangle
is not lost from line CAC 3 to the sand–loam interface, i.e., the “sand” filling of the rectangle isalmost equivalent to no-filling in an irrigation pipe with a zero Darcian resistance of the pipe
walls (mole hole emitters, Bobchenko 1957). The boundary conditions are shown in Fig. 6a:
the upper and left faces of the rectangle are no-flow boundaries; at point C3 the pressure
head p = 90 cm; along CAC 3 the pressure head is hydrostatic, i.e., decreases from 90 cm topo = 85 cm at point C (i.e., the total head, h, is constant along C3 C2 ).
Along the bottom I1 I2 of the flow domain (Fig. 6a), the choice of the boundary condition
depends on the knowledge of the “deep subsurface” into which water is drained. In common
applications of HYDRUS (see, e.g., El-Nesr et al. 2014; Siyal and Skaggs 2009), the “freedrainage” boundary condition (in a steady-state flow this implies a unit gradient), i.e., freegravitational flow at the bottom, I1 I2 , of the modeled flow domain. In other words, the watertable is postulated to be “somewhere deep”.
We assume that prior to irrigation from emitters, groundwater is static, the water table
is horizontal and positioned at an elevation pi > 0 above I1 I2 in the “minor” aquifer, as
indicated in Fig. 6b. We repeat that the boundary condition p = pi along I1 I2 holds also at
any t > 0 during recharge from the emitters. Vedernikov (1939) explained that this unalteredisobaricity of the line I1 I2 is realized if the substratum is “mighty” enough. Then, the irrigationwater seeping from a furrow, channel, emitter, or other sources is well drained either laterally(this is indicated by the horizontal arrow in the substratum, Fig. 6c) or vertically into an evendeeper aquifer. A steady saturated–unsaturated flow (at t > 0) in the soil layer (above I1 I2 )makes a 2-D groundwater mound (“unconfined aquifer”) as compared with an undisturbed
situation of Fig. 6b. We repeat that this 2-D flow is occluded by pi and L, viz. a decrease of
L or an increase of pi (provided all other parameters are fixed) reduces Q and increases thelocus of the phreatic surface. This follows from both our HYDRUS simulations and from thevariational principles (Gol’dsheteı̆n and Entov 1994). Moreover, if one increases pi above
a certain threshold value, the constant total head line CAC 3 works as a drain (mathematicalsink), i.e., groundwater seeps from the minor aquifer into this line (the situation modeledby Deemter 1950). This draining action of the subsurface pipes and mole drains after heavy
rains was detected and described in K-35 and Bobchenko (1957).
What if a hypothetical emitter is not a small square (as in Fig. 6a) but a “stretched” thin stripof a length L shown in Fig. 6c? In this case, if we ignore capillarity—recalling the prescribedno head loss along the sand strip from C3 C to C2 C1 —we arrive at another 1-D flow regime.
Above CC 1 , groundwater stands still at an elevation pi . From C3 C2 to I1 I2 , a fully saturated
1-D vertical flow takes place, with a Darcian velocity of v Z H = K 1 ( po + L v − pi )/L v .

123

Steady Flow from an Array of Subsurface Emitters: Kornev’s…

ZH

HYDRUS2D
mesh 40*40
layers655“Minor” unconfined
aquifer, relatively
less permeable than
the substratum

ZH

100 cm200 cmNo-flowboundary

HydrostaticpressureheadCboundary
C3AC
A

No-flowboundaries

Van Genuchten’sloam, K1=25 cm/day

C1
5cm

A

XH
Highly permeable
substratum,“main”, confined
aquifer, static
groundwater

ZH

Lv=112.5cm

No-flowboundary

Constant pressure headboundary, p=pi
XH

Z

I2

po

C
C3

X

I1

I1

Pre-irrigation
static watertable horizon

I2

(b)

C3
C2
2.5cm
Van Genuchten’s“sand lens”,
Kd =713 cm/day, models
an emitter centered at
XH=0, ZH=112.5 cm andof a radius rC=2.5 cm

I1

pi

C1
C2
Lv

Pressurehead p=po
Thin-strip“aquifer”

I2
XH

Project kidder-kornev March 10 - saturated-unsaturated flow from periodic emitters
Domain Properties, Material Distribution

(a)

(c)

“aquitard”

Pressurehead p=pi

Inflexion point N1

F
m=40 cmd=53.5 cm

F1
A

Free drainage
outlet I1I2

Z

I1
-10.000-8.182-6.364-4.545

I2

X-2.727-0.909

0.909

2.727

4.545

Pressure Head - h[cm], Min=-112.146, Max=91.307

6.364

8.182 10.000(d)Project kidder-kornev March 10 - saturated-unsaturated flow from periodic emitters
Results, Pressure Head, Time 5 - 10 days

Fig. 6 HYDRUS modeling results: a the finite element mesh and a rectangular “sand lens” as a model ofa constant total head subsurface emitter. b Undisturbed hydrogeological conditions with a static horizontal
water table in an unconfined aquifer and a subjacent highly permeable confined aquifer. c A “degenerate”
case of the constant total head sand lens stretched into a strip. d Isobars close to the phreatic surface for the“free drainage” boundary condition along the bottom I1 I2 of the modeled rectangle, the maximum pressure
head at the “lens” apex is po = 85 cm (t = 10 days). e A colored map of isobars in the whole flow domain( po = 85 cm; t = 30 days). f Isobars close to the phreatic surface for the “backwater” boundary conditionpi = 110 cm along I1 I2 and po = 85 cm in the rectangular “lens” (t = 50 days). g Isobars close to
the phreatic surface for the “backwater” boundary condition pi = 45 cm along I1 I2 and po = 85 cm inthe rectangular “lens” (t = 50 days). h Isobars close to the phreatic surface for the “backwater” boundary
condition pi = 0 cm along I1 I2 and po = 85 cm in the rectangular “lens” (t = 50 days). i Isobars close to
the phreatic surface for the “backwater” boundary condition pi = 45 cm along I1 I2 and po = 85 cm along
the semi-circular emitter (t = 50 days). j Isobars around a semi-circular emitter: rC = 2.5 cm, po = 85 cm,
pi = 138 cm, d p = 85 cm, an emitter axis is 188 cm above a horizon Dl Dr . k Isobars for flow past an empty
hole: rC = 15 cm, pi = 10 cm, d p = 75 cm, a hole axis is 147.5 cm above a “backwater” horizon (the outlet
side of the rectangle), the inlet of the rectangle receives infiltration at a rate of qi = 3 cm/day. l Isobars forflow past a hole with a small quantity of water at its bottom: rC = 15 cm, pi = 10 cm, d p = 75 cm, the holeaxis is 147.5 cm above a “backwater” horizon, qi =6 cm/day

Flow in Fig. 6c reflects a standard conceptualization and a hydraulic connection between
two highly permeable units, the strip CC 1 C2 C3 (an upper “aquifer”) and the “main aquifer”.
The rectangle C3 I1 I2 C2 (Fig. 6c) acts as a “thick aquitard“. We recall that in agreement

123

656

A. R. Kacimov et al.

Free drainage
outlet I1I2

Z

X-100.000-82.608-65.216-47.824-30.431-13.0394.353

21.74539.13756.52973.92191.314Pressure Head - h[cm], Min=-100.000, Max=91.314(e)Project kidder-kornev March 10 - saturated-unsaturated flow from periodic emitters
Results, Pressure Head, Time 3 - 30 days

Isobaricoutlet I1I2
at pi=110 cm
Z

X-10.000-8.182-6.364-4.545-2.727-0.909

0.909

2.727

4.545

6.364

8.182

10.000(f)

Pressure Head - h[cm], Min=-200.000, Max=110.000

Project kidder-kornev March 10 - saturated-unsaturated flow from periodic emitters
Results, Pressure Head, Time 5 - 50 days

Isobaricoutlet I1I2
at pi=45 cm
Z

X-10.000-8.182-6.364-4.545-2.727-0.909

0.909

2.727

4.545

6.364

8.182

10.000Pressure Head - h[cm], Min=-200.000, Max=91.516

Project kidder-kornev March 10 - saturated-unsaturated flow from periodic emitters
Results, Pressure Head, Time 5 - 50 days

Fig. 6 continued

123

(g)

Steady Flow from an Array of Subsurface Emitters: Kornev’s…

657

F

A
F1

Isobaricoutlet I1I2
at pi=0Z

I1-10.000-8.182-6.364-4.545

I2

X-2.727-0.909

0.909

2.727

4.545

6.364

8.182

10.000Pressure Head - h[cm], Min=-200.000, Max=90.529(h)Project kidder-kornev March 10 - saturated-unsaturated flow from periodic emitters
Results, Pressure Head, Time 5 - 50 days

po=85cm
Semicircularemitter atpressure head
po=85 cm
(emitter’sapex) andisobaricoutlet I1I2
at“backwater”
pressure head
pi=45 cmC
2rC= 5cm

C1 C1
C3Constant
total head
semicircle,
emitter fullof water at
positive pressureZX

0.242 0.260 0.277 0.294 0.311 0.328 0.345 0.362 0.379 0.396 0.413 0.430(i)Water Content - th[-], Min=0.242, Max=0.430

Project circular emitter Kornev's irrigation
Results, Water Content, Time 5 - 50 days

Isobars around asemi-circularemier: rC=2.5 cm,
po =85 cm, pi=138
cm, dp=85 cm, anemier axis is 188
cm above a horizon
DlDr.

Z

X

-12.000 -10.000-8.000-6.000-4.000-2.000

0.000

2.000

4.000

6.000

8.000

10.000Pressure Head - h[cm], Min=-47.308, Max=138.000

Project Oct 20 circular emitter-comparisons with analytics
Results, Pressure Head, Time 5 - 50 days(j)Fig. 6 continued

123

658

A. R. Kacimov et al.

Isobars for ﬂow past an
empty hole in Philip’s
regime: rC=15 cm, pi=10cm, dp=75 cm, the holeaxis is 147.5 cm above a“backwater” horizon, theinlet of the rectangle
receives inﬁltraon at a
rate of qi =3 cm/day
Z

X-20.000 -17.273 -14.545 -11.818-9.091-6.364-3.636-0.909

1.818

4.545

7.273

10.000Pressure Head - h[cm], Min=-25.241, Max=10.000(k)Project circular emitter-4
Results, Pressure Head, Time 5 - 50 days

Isobars for flow past a hole with asmall quantity of water at its
bottom: rC=15 cm, pi=10 cm,
dp=75 cm, the hole axis is 147.5
cm above a “backwater” horizon, qi=6 cm/day.
ZX-9.050-7.318-5.586-3.855-2.123-0.391

1.341

3.073

4.805

6.536

8.268 10.000

Pressure Head - h[cm], Min=-9.050, Max=10.000

Project circular emitter-5
Results, Pressure Head, Time 5 - 50 days(l)Fig. 6 continued

with common assumptions of groundwater hydraulics (PK-77), in aquifers the water flow isprevalently horizontal while in aquitards it is mostly vertical.
We used HYDRUS to simulate the two 1-D regimes sketched in Fig. 6b, c as asymptotic limits. HYDRUS confirmed these hydrogeologically “degenerate” scenarios. Next, we
simulated 2-D flows for realistic rC (Fig. 6a). A steady-state regime is established in about
10 days.
In the scenario shown in Fig. 6d, e, we assumed po = 85 cm and a HYDRUS “freedrainage” condition along I1 I2 . The colored images of the pressure heads in the vicinity ofthe phreatic surface p = 0 and in the whole rectangle are shown in Fig. 5d, e, respectively. The
computed values of d and m are 53.5 and 40 cm, respectively. The distribution of velocity
along I1 I2 showed a variation from 25.1 cm/day at point I1 (X H = 0, Z H = 0) to 24.8
cm/day at point I2 (X H = 100 cm, Z H = 0).
In scenarios shown in Fig. 6f–h, we assumed po = 85 cm and pi = 110, 45, and 0 cm,
respectively, along I1 I2 . The obtained pairs of (d, m) and the range of variations of the verticalcomponent of the Darcian velocity vector along I1 I2 are: (60, 10 cm), (13–14.6 cm/day);(50, 30 cm), (22.2–25.2 cm/day); and (43, 90 cm), (27.3–32.2 cm/day), respectively. If pi isdecreased to negative values (as in Ashrafi et al. 2002), then the water table drops further,123

Steady Flow from an Array of Subsurface Emitters: Kornev’s…

659

Q increases and the saturated zone around the emitter shrinks (and may even convolute into
Philip’s “saturated bulb”).
Next, we used a “2D - General” type of geometry of HYDRUS and modeled a semicircular arc C3 C1 C of a diameter of 5 cm (Fig. 6i) as the emitter boundary. This constanttotal head contour was placed at the same L v as the rectangular “sand lens” (in Fig. 6d–h).
The emitter boundary condition corresponds to po = 85 cm at the semicircle’s apex C. The
colored map of the moisture content (Fig. 6i) illustrates the water table position, in particular,
d = 46 cm, m = 40 cm. Similarly to Ashrafi et al. (2002) Siyal and Skaggs (2009), andWang et al. (2017), we can include a thin pipe wall made of a low permeable clay.
The results of HYDRUS simulations qualitatively agree with the analytical solutions in
Fig. 4. To illustrate this, we compared the case shown in Fig. 5 with numerical results for aloamy soil, a circular emitter of a radius of 2.5 cm subject to a total head of 85 cm, and placed
188 cm above a horizon Dl Dr (Fig. 2a) where a positive pressure head (“backwater”) is 138
cm. The upper no-flow boundary of a rectangular flow domain is placed 85 cm (d p in Fig. 1)
above the emitter axis (this “cap” segment can be selected at any sufficiently high elevation
above the emitter). Fig. 6j shows the computed isobars in the range of −10 cm< p <10 cm.
The middle isobar is a phreatic surface. Its point F (see Fig. 2a) is about 48 cm above thedrain axis as compared to 64 cm in the analytical solution (Fig. 5). Point F1 is about 34 cm
above the drain axis (32 cm in the analytical solution, Fig. 5). The presence of the inflexion
point on all isobars in Fig. 6j (including the phreatic line) is also apparent, similarly to Khan
and Rushton (1996a, Figs. 1a, 2, 4, 8a) and Rushton (2003, Figs. 4.23, 4.25a).
We note that even at pi = 0 (Fig. 6f), i.e., in the case of no “backwater” in the substratum(Fig. 6b), the pore pressure is positive along I2 F1 . If we keep pi = 0 and increase L above
100 cm, then at a certain value L = L cr , a purely unsaturated flow will take place close to
X H = L cr . In the analytical solution, this corresponds to the regime discussed in Sect. 2:
the free surfaces of the neighboring sources do not intersect (Kacimov and Obnosov 2016),
the corresponding water tables are “ lighthouse-shaped” (dotted line in Fig. 1a). If, however,pi > 0 (“backwater” regime), then the phreatic surfaces in both analytical and numericalsolutions always intersect and within one half-period (Fig. 2a) an inflection point exists.

3.2 Infiltration Over Subterranean Holes in the Philip Regime
HYDRUS has been successfully used in studies of 2D and 3D flows toward horizontal drains,
including mole drains (e.g., Boivin et al. 2006; Brunetti et al. 2017; Filipović et al. 2014;
Karandish et al. 2017). In these studies, drains were represented with a standard seepage
face boundary condition along the drain contour (i.e., the drain was assumed to be empty).
Therefore, the drains acted as water sinks, as in Khan and Rushton (1996a, b) and Rushton (2003) who used saturated flow models in their analysis. However, Philip et al. (1989)
showed that an empty cavity (hole) may act as an obstruction to an unsaturated descending
flow if its intensity is small enough (see also Stormont and Zhou 2005). In other words,
the hole hydrodynamically behaves as neither a sink or a source but as a dipole (Kacimov2007). Therefore, if one models an empty hole in this regime, then the HYDRUS seepage
face boundary condition along the contour can be expected to be equivalent to the no-flow
condition. This capillary barrier effect has been used in design and hydrological analysisof the Yucca Mountain tunnels as elements of subsurface nuclear waste repositories (e.g.,
Wang and Bodvarsson 2003). Birkholzer et al. (1999) investigated a descending unsaturated flow in the regime of Philip et al. (1989), without considering an opportunity of water
dripping from the roof of the tunnel and accumulation at the bottom of the drift. In thissubsection, we consider a flow regime past a hole with a shallow water level accumulated

123

660

A. R. Kacimov et al.

at the hole bottom that is a common situation in mole drains and emitters (Bobchenko
1957).
In Fig. 6k, we selected the same soil and the rectangular domain as in Fig. 6j, except
the hole size (rC = 15 cm) much larger than that of the emitter in Fig. 6j. We set d p =75 cm and the bottom draining boundary Dl Dr at the depth of 147.5 cm under the axis. The
boundary condition along the outlet is pi = 10 cm (“mild “backwater”). The infiltration rate
qi = 3 cm/day was assumed along the top of the rectangle. Fig. 6k shows for a steady-stateflow the pressure head contours in the range of − 20 cm< p < 10 cm. They, as well as
the moisture and velocity fields computed by HYDRUS, are qualitatively similar to what is
depicted in Fig. 1 of Philip et al. (1989). The pressure head profile shows a roof-drip lobe, a
retarded zone on the upflow side of the hole, and a dry shadow on the downflow side. Alongthe vertical axis x = 0 (Fig. 6k), the minimum of the pressure head (the driest point) ofp ≈ − 25 cm is attained at the bottom of the hole, while the maximum of p ≈ − 5.8 cm isattained at the hole’s apex, i.e., the upper stagnation point (Philip et al. 1989), similarly towhat Wang and Bodvarsson (2003) and Birkholzer et al. (1999) reported for large tunnels.
Due to the encumbrance of our hole, groundwater mounds in the right bottom corner of theflow domain (Fig. 6k). An increase of the infiltration rate qi to approximately 5 cm/day would
result in a breakthrough of moisture into the hole through its apex. For higher values of qi ,water would drip through the hole’s upper segment (see Kacimov 2000).
Figure 7 shows the case of a hole partially filled with a static water and placed into
an ambient unsaturated flow of a small qi . The segment Ch2 Ch1 is a constant total headboundary, the segments Ch1 Sf and SfAp are impermeable. If the water level in the hole issmall and the draining horizon Dl Dr (Fig. 2a) is deep enough with a small “backwater”,
then a “perched” water table (as in Warrick and Zhang 1987), rather than a “dry shadow”
of Philip et al. (1989) is formed just under the hole. The hole in Fig. 7 again does not drain
the incident infiltration but replenishes the “regional” water table by seepage from Ch2 Ch1 .
Obviously, the regional water table (not shown in Fig. 7) mounds in response to this seepage
just close to the hole axis, unlike mounding at the corner in Fig. 6k. Figure 6l shows HYDRUS
simulations for a steady-stated flow and the same soil and geometry as in Fig. 6k, with thefollowing difference in boundary conditions along the hole contour: the topsoil is subject toan infiltration flux qi = 6 cm/day, the water level in the hole is 2.9 cm. A constant total headboundary condition was specified along Ch2 Ch1 , and a seepage face boundary condition was
imposed along Ch1 Ap. The computations in Fig. 6l are qualitatively in congruency with theflow topology and the soil wetness field from Philip et al. (1989) and Warrick and Zhang(1987).

4 Conclusions and Perspectives
We decoupled a saturated–unsaturated steady-state Darcian flow from an infinite line ofequally spaced PPs that emit water under a positive pressure, each at the same rate. For
the first flow fragment of a potential flow, we followed Holmes (see the epigraph) and used
Deemter’s (1950) and Kidder’s (1956) “singularity clues” when applying the hodographmethod for line sources. From an explicit analytical solution to a free BVP in a homogenous
and isotropic soil stratum, the computed phreatic surface is wavy with inflection points. Thissurface bulges above the emitters and is depressed between two neighboring linear sources.
In the limiting case, the circular cut in the hodograph domain (Deemter 1950) degenerates
into a full half-circle (see Kidders’ Fig. 2) that corresponds to cusps of the phreatic surface.

123

Steady Flow from an Array of Subsurface Emitters: Kornev’s…

661

qi
Isobar p=-p2,
0<p1<p2

Isobar p=-p1,

Ap
Philip’sroof-drip
lobe

No-flow segment of the hole
contour (streamline)

Air inside the hole atatmospheric pressure

Sf
Ch1

Free static water
filling the hole
Positive soilpressure zone

Ch2

Streamlines

Segment of the hole contour
at a constant total head, water
seeps from the hole
Perched phreatic surface,isobar p=0

Fig. 7 Descending infiltration flow past a hole with an impermeable roof and constant-head bottom section
seeping into the soil

In our Kidder-type phreatic flow, we ignored evaporation from F1 F (Fig. 2a), i.e., weassumed that unsaturated flow does not affect potential flow. The Deemter (1950) andPolubarinova-Kochina (1962) evaporation model takes into account coupling between two
flows by assuming that along the phreatic surface ψ = ex + Q, where e = const > 0 is
a given intensity of evaporation. Flow from the source (Fig. 3d) then splits into ascending
and descending parts, with a separatrix AE and a stagnation point E (see Fig. 3 of Deemter1950). Along this curve, as well as along EF 1 and ED1 , the stream function ψ = Q E , where
Q E < Q is the part of the source-emitted water that descends as deep percolation. The
hodograph domain (Fig. 3e) remains the same as in Fig. 2a, modulo translated upward to a
distance μ.
Our HYDRUS simulations substantiate the analytical findings. For saturated–unsaturated
flows, we followed Hanson et al. (2008), i.e., we considered a HYDRUS rectangular domain
with an emitter modeled as a constant total head contour (a semicircle or rectangle). However,in order to compare the numerical and analytical solutions, we bounded the bottom of therectangle from below by a positive-pressure horizon (“backwater”) similarly to Choudhary
and Chahar (2007), which represents a common hydrogeological “outlet drainage” (see,also Vedernikov 1939, PK-77). This “backwater” is routinely encountered in regular surface
irrigation projects with seepage from earth channels and furrows (Bouwer 1978; Ilyinsky

123

662

A. R. Kacimov et al.

and Kacimov 1991; Kacimov and Youngs 2005), but was not previously modeled in SI. We
also modeled unsaturated–saturated zonation more general than in Philip et al. (1989) fora subterranean hole placed in a descending infiltration-induced ambient flow, with a mixed
boundary condition along the hole contour, which is impermeable near hole’s roof and leaks
through a constant total head segment at hole’s bottom. The future HYDRUS analysis might
be warranted to show how flow fields from emitters and past holes are affected by soil layeringand fluxes.
Acknowledgements This work was supported by TRC (Oman) Grant ORG/EBR/15/002 “Artificial capillarybarriers: A smart agro-engineering technique for saving irrigation water in Oman”, and through a special
program of the Russian Government supporting research at Kazan Federal University. Helpful commentsfrom three referees are highly appreciated.

References
Al-Rawahy, S., Rahman, H.A., Al-Kalbani, M.S.: Cabbage (Brassica oleracea L.) response to soil moistureregime under surface and subsurface point and line applications. Int. J. Agric Biol. 6, 1093–1096 (2004)
Anderson, E.I.: Stable pumping rates for horizontal wells in bank filtration systems. Adv. Water Resour. 54,
57–66 (2013)
Ashrafi, S., Gupta, A.D., Babel, M.S., Izumi, N., Loof, R.: Simulation of infiltration from porous clay pipe insubsurface irrigation. Hydrol. Sci. J. 47(2), 253–268 (2002)
Birkholzer, J., Li, G., Tsang, C.F., Tsang, Y.: Modeling studies and analysis of seepage into drifts at Yucca
Mountain. J. Contam. Hydrol. 38(1), 349–384 (1999)
Bobchenko, V.I.: Subsurface Irrigation. Selkhozgiz, Moscow (1957). (in Russian)
Boivin, A., Šimůnek, J., Schiavon, M., van Genuchten, M.T.: Comparison of pesticide transport processes inthree tile-drained field soils using Hydrus-2D. Vadose Zone J. 5(3), 838–849 (2006)
Bouwer, H.: Groundwater Hydrology. McGraw, New York (1978)
Brunetti, G., Šimůnek, J., Turco, M., Piro, P.: On the use of surrogate-based modeling for the numericalanalysis of low impact development techniques. J. Hydrol. 548, 263–277 (2017)
Choudhary, M., Chahar, B.R.: Recharge/seepage from an array of rectangular channels. J. Hydrol. 343(1),71–79 (2007)
Communar, G., Friedman, S.P.: Cylindrically confined flows of water from a point source and to a sink, above
a water table. Vadose Zone J. (2015). https://doi.org/10.2136/vzj2014.12.0183
Deemter, J.V.: Theoretische en numerieke behandeling van ontwaterings—en infiltratie-stromings problemen(in Dutch). Theoretical and numerical treatment of flow problems connected to drainage and irrigation.
Ph.D. Dissertation, Delft University of Technology (1950)
El-Nesr, M.N., Alazba, A.A., Šimůnek, J.: HYDRUS simulations of the effects of dual-drip subsurface irrigation and a physical barrier on water movement and solute transport in soils. Irrig. Sci. 32, 111–125(2014)
Filipović, V., Mallmann, F.J.K., Coquet, Y., Šimůnek, J.: Numerical simulation of water flow in tile and moledrainage systems. Agric. Water Manag. 146, 105–114 (2014). https://doi.org/10.1016/j.agwat.2014.07.
020
Fujii, N., Kacimov, A.R.: Analytically computed rates of seepage flow into drains and cavities. Int. J. Numer.
Anal. Methods Geomech. 22, 277–301 (1998)
Gol’dsheteı̆n, R.V., Entov, V.M.: Qualitative Methods in Continuum Mechanics. Longman, London (1994)
Hanson, B.R., Šimůnek, J., Hopmans, J.W.: Leaching with subsurface drip irrigation under saline, shallow
groundwater conditions. Vadose Zone J. 7(2), 810–818 (2008)
Honari, M., Ashrafzadeh, A., Khaledian, M., Vazifedoust, M., Mailhol, J.C.: Comparison of HYDRUS-3D
soil moisture simulations of subsurface drip irrigation with experimental observations in the South ofFrance. J. Irrig. Drain. Eng. 143(7), 04017014 (2017)
Ilyinsky, N.B., Kacimov, A.R.: Estimates of backwater and drying levels in hydrodynamic model. Fluid Dyn.
26(2), 224–231 (1991)
Iwama, H., Kubota, T., Ushiroda, T., Osozawa, S., Katou, H.: Control of soil water potential using negative
pressure water circulation technique. Soil Sci. Plant Nut. 37(1), 7–14 (1991)
Ilyinsky, N.B., Kacimov, A.R.: Problems of seepage to empty ditch and drain. Water Resour. Res. 28(3),871–877 (1992)

123

Steady Flow from an Array of Subsurface Emitters: Kornev’s…

663

Kacimov, A.R.: Circular isobaric cavity in descending unsaturated flow. J. Irrig. Drain. Eng. 126(3), 172–178(2000)
Kacimov, A.R.: Capillarity and evaporation exacerbated seepage losses from unlined channels. J. Irrig. Drain.
Eng. 132(6), 623–626 (2006)
Kacimov, A.R.: Dipole-generated unsaturated flow in Gardner soils. Vadose Zone J. 6, 168–174 (2007)
Kacimov, A.R., Obnosov, Y.V.: Tension-saturated and unsaturated flows from line sources in subsurface
irrigation: Riesenkampf’s and Philip’s solutions revisited. Water Resour. Res. (2016). https://doi.org/10.
1002/2015WR018221
Kacimov, A., Obnosov, Y.: Analytical solution for tension-saturated and unsaturated flow from wicking porous
pipes in subsurface irrigation: the Kornev-Philip legacies revisited. Water Resour. Res. 53(3), 2542–2552(2017). https://doi.org/10.1002/2016WR019919
Kacimov, A.R., Youngs, E.G.: Steady-state water-table depressions caused by evaporation in lands overlying
a water-bearing substratum. J. Hydrol. Eng. 10(4), 295–301 (2005)
Karandish, F., Darzi-Naftchali, A., Šimůnek, J.: Application of HYDRUS (2D/3D) for predicting the influenceof subsurface drainage on soil water dynamics in a rainfed-canola cropping system. Irrig. Drain. J. https://doi.org/10.1002/ird.2194
Kato, Z., Tejima, S.: Theory and fundamental studies on subsurface irrigation method by use of negative
pressure. Trans. Jpn. Soc. Irrig. Drain. Reclam. Eng. 101, 45–54 (1982). (in Japanese)
Khan, S., Rushton, K.R.: Reappraisal of flow to tile drain I. Steady state response. J. Hydrol. 183(3–4), 251–366(1996a)Khan, S., Rushton, K.R.: Reappraisal of flow to tile drains III. Drains with limited flow capacity. J. Hydrol.
183(3–4), 383–395 (1996b)
Kidder, R.E.: Flow of immiscible fluids in porous media: exact solution of a free boundary problem. J. Appl.
Phys. 27(8), 867–869 (1956)
Kornev, V.G.: Subsurface Irrigation. Selhozgiz, Moscow-Leningrad (1935). (in Russian)Lazarovitch, N., Šimůnek, J., Shani, U.: System-dependent boundary condition for water flow from subsurface
source. Soil Sci. Soc. Am. J. 69(1), 46–50 (2005)
McCarthy, J.F.: Analytical solutions of 2D cresting models using the hodograph method. Transp. Porous Media
15(3), 251–269 (1994)
Moniruzzaman, S.M., Fukuhara, T., Ito, M., Ishii, Y.: Seepage flow dynamics in a negative pressure difference
irrigation system. J. Jpn. Soc. Civ. Eng. Ser. B 67(4), 97–102 (2011)
Obnosov, Y.V., Kacimov, A.R.: Steady Darcian flow in subsurface irrigation of topsoil impeded by substratum:
Kornev–Riesenkampf–Philip legacies revisited. Irrig. Drain. (2017, in press)
Philip, J.R.: Multidimensional steady infiltration to a water table. Water Resour. Res. 25, 109–116 (1989)
Philip, J.R., Knight, J.H., Waechter, R.T.: Unsaturated seepage and subterranean holes: conspectus, and exclusion problem for circular cylindrical cavities. Water Resour. Res. 25(1), 16–28 (1989)
Polubarinova-Kochina, P.Y.: Theory of Ground Water Movement. Princeton University Press, Princeton. The
second edition of the book (in Russian) was published in 1977, Nauka, Moscow (1962)
Rushton, K.R.: Groundwater Hydrology: Conceptual and Computational Models. Wiley, Chichester (2003)
Shani, U., Xue, S., Gordin-Katz, R., Warrick, A.W.: Soil-limiting flow from subsurface emitters. I: pressure
measurements. J. Irrig. Drain. Eng. 122(5), 291–295 (1996)
Siyal, A.A., Skaggs, T.H.: Measured and simulated soil wetting patterns under porous clay pipe sub-surface
irrigation. Agric. Water Manag. 96, 893–904 (2009)Šimůnek, J., van Genuchten, MTh, Šejna, M.: Recent developments and applications of the HYDRUS computer
software packages. Vadose Zone J. 15(7), 25 (2016). https://doi.org/10.2136/vzj2016.04.0033
Stormont, J.C., Zhou, S.: Impact of unsaturated flow on pavement edgedrain performance. J. Transp. Eng.
131(1), 46–53 (2005)
Strack, O.D.L.: Groundwater Mechanics. Prentice Hall, Englewood Cliffs (1989)
Swamee, P.K., Chahar, B.R.: Design of Canals. Springer, New Delhi (2015)
Tanigawa, T., Yabe, K., Tejima, S.: Comparison of prediction of actual measurement about dynamic distributionof soil moisture Tension. Trans. Jpn. Soc. Irrig., Drain. Reclam. Eng. 137, 9–16 (1988). (in Japanese)
Vedernikov, V.V.: Theory of Seepage and Its Applications to Problems of Irrigation and Drainage. Gosstroiizdat,
Moscow-Leningrad (1939). (in Russian)Wang, J.S.Y., Bodvarsson, G.S.: Evolution of the unsaturated zone testing at Yucca Mountain. J. Contam.
Hydrol. 62, 337–360 (2003)
Wang, J.J., Huang, Y.F., Long, H.Y., Hou, S., Xing, A., Sun, Z.X.: Simulations of water movement and solutetransport through different soil texture configurations under negative-pressure irrigation. Hydrol. Process.
31, 2599–2612 (2017)
Warrick, A.W., Shani, U.: Soil-limiting flow from subsurface emitters. II: Effect on uniformity. J. Irrig. Drain.
Eng. 122(5), 296–300 (1996)

123

664

A. R. Kacimov et al.

Warrick, A.W., Zhang, R.: Steady two-and three-dimensional flow from saturated to unsaturated soil. Adv.
Water Resour. 10(2), 64–68 (1987)
Wolfram, S.: Mathematica. A System for Doing Mathematics by Computer. Addison-Wesley, Redwood City
(1991)
Yabe, K., Tanigawa, T.: Studies on the shape of porous pipe under cultivated condition. Experimental studies
on sub-irrigation method (VIII). Trans. Jpn. Soc. Irrig. Drain. Reclam. Eng. 162, 43–48 (1992). (inJapanese)123