Profiling ponded soil surface in saturated seepage into drain-line sink: Kalashnikov’s method of lateral leaching revisited

Abstract

Two boundary value problems are solved for potential steady-state 2D Darcian seepage flows towards a line sink in a homogeneous isotropic soil from a ponded land surface, which is not flat but profiled. The aim of this shaping is ‘uniformisation’ of the velocity and travel time between this surface and a horizontal drain modelled by a line sink. The complex potential domain is a half-strip, which is mapped onto a reference plane. Either the velocity magnitude or a vertical coordinate along the land surface are control variables. Either a complexified velocity or complex physical coordinate is reconstructed by solving mixed boundary-value problems with the help of the Keldysh-Sedov formula via singular integrals, the kernel of which are the control functions. The flow nets, isotachs and breakthrough curves are found by computer algebra routines. A designed soil hump above the drain ameliorates an unwanted ‘preferential flow’ (shortcut) and improves leaching of salinised soil of a cropfield during a pre-cultivation season.

1. Introduction

Systematic tile drains (mathematically modelled by line sinks S of intensity 2q ^* [m²/s]) are commonly installed at the depth d ^* [m] (commonly 1–3 m) with a period 2L ^* [m] (commonly few meters-tens of meters) under a flat soil surface BMC (Figure 1(a)). The leachate (water collected by the drains) is conveyed to the so-called collector tubes laid perpendicular to the drains, commonly 0.5–1 m below them. In humid-semihumid climates (e.g. Canada, Holland, UK, Minnesota, Northern Russia), after occasional heavy Summer rains/regularly in Springs, the drains remove infiltrated rainwater/thawed snow from the soil of the root zone of cropfields (see e.g. [39, 44, 50, 52]. Drain pipes flow gravitationally (without pumping) perpendicular to the vertical planes, as shown in Figure 1(a)–(c).

In Central Asia and other arid and semi-arid countries [4], the periodic drains in Figure 1(a) have another hydromeliorating purpose: in Winters and early Springs, the soil surface is ponded by a fresh river water, which leaches the salts (commonly targeted anions are Cl- and SO₄ ²-). They are accumulated in the root zone due to evapotranspiration during the hot cultivation seasons (late Springs – Autumns). For leaching, the cropland is divided into rectangular bays demarcated by bunds. The pioneer of this technology, the Russian agricultural engineer Bushuev, experimented with such bays as seasonal infiltration basins and found that desalination of the root zone is increased up to 10 times if the infiltration rate of ponded water is high, as compared with the scenarios with mild infiltration when soil salts are only diffused into the bay water, rather than advected downwards [35–37].

On another hand, the salt-amended leachate, generated by ponding-seepage, should not percolate to shallow groundwater (itself often saline, with salinities sometimes exceeding 5000 mg/l) subjacent to the drains [3]; [5]; [15]; [40]; [45]). Otherwise, the water table rises to the depth of few tens of cm (from 2 m and less the depth of groundwater is considered dangerously low) that is detrimental even to relatively salt-tolerant cotton plants. In other words, leaching by ponding without drainage causes both waterlogging (with an ensued anoxia), intensive evaporation and adverse secondary salinisation. For example, in the Hungry Steppe of Uzbekistan (the area of the Steppe is almost 1 mln. ha) of the 1910s, the croplands, which had no drain trenches or tile drains, were lost to badlands due to waterlogging after 1–2 years of border irrigation. Recuperation (reducing the soil salt content by 2–3 times) by smart drainage leaching schemes in the 1930s took decades [34, 35].

After the last drought in California, a heavy winter ponding of orchards was also introduced as the so-called Ag-MAR technology (see e.g. [14]). The hydrodynamic goal is opposite to leaching a thin vadose zone in Uzbekistan, viz. the surface water is forced as much as possible to a deep water table through a thick vadose zone, hoping to replenish an unconfined aquifer subjacent to a cropfield.

Figure 1. (a) Systematic Vedernikov’s drains seeping under a flat ponded soil surface; (b) one period of Kalashnikov’s system of furrows with seepage into ‘lateral’ line sinks; (c) proposed ponded ‘humped’ soil with seepage to a line sink; (d) sketch of an advective travel time along streamlines from BMC towards the sink.

The disadvantage of the leaching system in Figure 1(a) is caused by the topology of seepage: above the drain (near point M in Figure 1(a)) the velocity of pore water motion is relatively high as compared with the peripheral zones (near points B and C), i.e. a hydrodynamic shortcut takes place between the vicinity of point M on the ponded surface and the sink below. Heavy leaching norms are required to wash out soil in Figure 1(a), unless drains are buried deep that is also costly to install and maintain. Deep drainage also requires longer leaching times, which, as we already mentioned, are limited by the Winter – early Spring season.

Analytical solutions for potential steady 2D seepage towards the array of sinks in Figure 1(a) were obtained by Vedernikov [48] (see also [42], hereafter abbreviated as PK-62) and later by Dutch, British, Danish and American irrigation-drainage engineers and applied mathematicians (see e.g. [26] for references to Western sources).

In the 1960s, the Soviet irrigation engineer Kalashnikov ([30], [31], [32], after decades of designing and constructing agricultural drainage in Uzbekistan, patented ¹ and published the so-called ‘lateral leaching’ technology, which was applied to recuperation of salinised soils in Central Asia. The crux of Kalashnikov’s method is depicted in Figure 1(b), which shows a vertical cross-section of one period of a system of furrows (channels) into which fresh water is supplied to seep towards horizontal drains (trenches or buried tile drains). In this scheme of seepage, less fresh water is needed than in Vedernikov’s flat-surface drainage (Figure 1(a)), and only a part of the soil surface is ponded viz. the segment F _B F _C remains dry. The peripheral soil zones near points B and C in Figure 1(c) are better leached, albeit a phreatic line F _B MF _C is formed. Soil above this curve remains unsaturated, saline and poorly leached [36]. Analytical solutions to several free boundary problems related to Kalashnikov’s lateral leaching are obtained in a set of papers and a book by Emikh (e.g. [7–11]). Emikh ignored the shape depth of furrows in Figure 1(b) and found that for this simplest, zero-depth geometry of the leaching channel the phreatic line in Figure 1(b) ‘overhangs’ the sink (see [26] for terminology) and has (one branch) up to 2 inflexion points.

In this paper, we propose the following:

• A construction method of a drained cropfield, totally opposite to a standard one (Legostaev, 1953, [33]) when the land surface is levelled (and bunded) prior to ponding. In our case, this surface is undulated that minimises the disadvantages of the seepage topology in the Vedernikov and Kalashnikov (Figure 1(a) and (b)) drainage patterns.

• An analytical ‘constructal’ (in terms of Adrian Bejan) design technique, which utilises the method of inverse boundary value problems (hereafter, BVPs) [19, 20, 47] and reconstructs a non-flat soil boundary BMC, with a porous medium bulging above the drain.

Specifically, in order to make a hummocky ground surface, one should bulldoze soil from the left (point C) and right (point B) to make a porous hump BMC (Figure 1(c)). The flow path along MS (the shortest streamline in Figure 1(a)) is then increased. Moreover, along BMC, we specify a desired distribution of velocity or stream function depending on the vertical coordinate i.e. – in a sense – we control the flow net by shaping a boundary of the flow domain.

Historically, Zhukovsky [53] championed in analytical aerodynamics by superposing a dipole at infinity (a unidirectional ambient flow of an ideal fluid) with a line vortex that gave a special streamline, viz. the Zhukovsky airfoil (flow inside the airfoil zone is obviously discarded). Vedernikov’s potential flow (seepage) in Figure 1(a) is a superposition of an infinite number of dipoles, for each of which the sink and source are distance 2d apart, that generates a flat equipotential horizontal soil surface BMC (and seepage from the sources to this line is, obviously, discarded). In the inverse method [6], one starts with the Zhukovsky or Vedernikov combinations of hydrodynamic singularities and then ‘perturbs’ the result by an ‘enforced’ control of a part of the boundary of the flow domain. The Zhukovsky-Chaplygin method of superposed singularities is also used in modelling of potential seepage flows with standard free boundaries [2, 25]; [21].

2. Inverse BVPs and their solutions

Similarly to Kacimov et al. [27], Kacimov and Obnosov [23, 25], and Obnosov and Kacimov, 2018 [41] we apply the theory of holomorphic functions to potential 2D Darcian flows. We generalise the Vedernikov [48], pp.174–178) solution for seepage flow shown in his Figure 72 and our Figure 1(a).

We select a system of Cartesian coordinates Ox ^* y ^* and introduce a complex physical variable z ^* = x ^* +i y ^* such that points B and C have complex coordinates -L ^* and L ^* . A line sink is placed at z ^* = -id ^* . We fix this depth, albeit Kalashnikov, Legostaev and other irrigation engineers debated d as an important control variable. A homogeneous, isotropic and capillarity-free soil of a hydraulic conductivity k [m/s] and effective porosity m extends deep in the -y ^* direction. A constant ponding depth of free water above points B and C in Figure 1(c) is $h^{^{\ast}}_p$ [m]. Seepage obeys the Darcy law:

(1) \begin{align}{\rm{ }}{\overrightarrow V}^{^{\ast}}({x^*},{y^*}) = - k\nabla {h^*}({x^*},{y^*})\end{align}

where ${\rm{ }}{\overrightarrow V ^*} = {u^*} + {\rm{i }}{v^{^{\ast}}}$ [m/s] is a complexified vector of specific discharge (u and v are the velocity components), and h ^* [m] is the piezometric head. The complex potential is ${w^*} = {\phi ^*} + {\rm{i}}{\psi ^*}$ , ${\phi ^*} = - k{\kern 1pt} {h^*}$ , where ${\phi ^*}$ [m²/s] is the velocity potential and ${\psi ^*}({x^*},{y^*})$ [m²/s] is a stream function, conjugated with ${\phi ^*}$ via the Cauchy–Riemann conditions. The functions ${w^*}({z^*})$ and ${V^*}({z^*})$ are holomorphic and antiholomorphic everywhere in the flow domain G _z (symmetric with respect to the ordinate axis), except at point S. In the vicinity of this point, we follow Emikh [10], PK-62 and Vedernikov [48]: the singularity is isolated by selecting an equipotential line NS _B S _C, which is a quasi-circle of a small radius r. This contour, comprising S, models a perfect horizontal drain, i.e. a perforated tube, fully occupied by water, with a gravel pack wrapping the tube such that the pressure head, p ^* (x,y) [m], at point N is $p_N^*$ [m]. Vedernikov [48], p. 161 calls this regime of drainage ‘the worse’, i.e. seepage is impeded by the piezometric head inside the drain and collector tubes. Mityushev and Adler [38] studied the case of an empty circular drain placed between two semi-infinite Kalashnikov’s channel rays (Figure 1(b)). This regime is the ‘best’ in terms of Vedernikov. The empty drains of Mityushev and Adler are relevant rather to mole drains, which intercept infiltration generated by drizzles, rather than tube drains operated in hydraulically ‘bad’, surface-ponded circumstances (see [13, 18] for layered soils and various boundary conditions along the drain contours, both circular and modelled as Zhukovsky’s slots). Generally speaking, $p_N^*$ varies in the direction perpendicular to the vertical plain in Figure 1 [17], but in our model below we follow PK-62 and Vedernikov [48], i.e. ignore the Darcy–Weisbach drop of the pressure head in the free water moving inside the drain pipe.

A nonflat soil surface BMC is a part of solution.

We count ${\phi ^*}$ and ${h^*}$ from BMC and assume that ${\psi ^*} = 0$ along the segment MONS (Figure 1c), i.e. point M is fiducial (w _M = 0). We make a cut BCS in G _z . Then along the rays CC _i S _B S and BB _i S _B S, the stream function is ${\psi ^*} = \pm {\kern 1pt} q$ . The complex potential domain is a half-strip shown in Figure 2(a). With such a choice of the fiducial point, the pressure head in G _z is ${p^*} = {h^*} - {y^*} + h^{^{\ast}}_p$ . The functions ${p^*}$ , ${h^*}$ , ${\phi ^*}$ , ${\psi ^*}$ ${u^*}$ , ${v^*}$ are harmonic.

Figure 2. (a) Complex potential domain G _w for a solitary drain under a profiled soil surface; (b) reference plain onto which G _w is mapped to; (c) the hodograph plane for a class of bulged soil surfaces obeying the phreatic line boundary condition along the soil surface.

The drop of piezometric head between BMC and the drain contour is $h_d^* = {d^*} + h^{^{\ast}}_p - {r^*} - p_N^*$ . This constant has to be positive to ensure seepage into the drain. Otherwise, a line source (subsurface emitter) will be modelled. At infinity (points B _i and C _i ), the velocity potential $\phi _i^* = k\,h_i^*$ is another positive constant, a part of solution. The complex potential domain, a half-strip G _w , is shown in Figure 2(a).

We introduce dimensionless quantities:

$\begin{array}{l}\qquad\qquad\qquad(x,y,L,r,h,{h_i},{h_p},p,{p_N}) = ({x^*},{y^*},{L^*},{r^*},{h^*},h_i^*,h^{^{\ast}}_p,{p^*},p_N^*)/{d^*},{\rm{ }}\\[5pt] \qquad\qquad\qquad (w,q) = ({w^*},{q^*})/(k{\rm{ }}{d^*}),\quad \left( {V,{\rm{ }}u,{\rm{ }}v} \right) = ({V^*},{u^*},{v^*})/k,\quad T = T{}^*k/({d^*}{\rm{ }}m).\end{array}$

We map conformally a reference half-plane $\eta \ge 0$ of an auxiliary variable $\zeta = \xi + {\rm{i}}\eta $ onto G _w by the function:

(2) \begin{align}{\rm{ }}w = - {\rm{ }}\frac{{2{\rm{i}}q}}{\pi }\arcsin \zeta \end{align}

We note that the seepage scenario in Figure1(c) and the corresponding complex potential domain in Figure 2(a) follow [48] assumption that soil is bounded from below by a deep bedrock, i.e. all leached surface water is intercepted by the drains. Alternative hydrogeological settings with a horizontal highly permeable substratum are investigated by Emikh [10]. In this case, the drains intercept only a portion of infiltrated water.

In the inverse method, a part of the boundary of G _z is not given but found from a boundary condition, which is different from a standard free boundary condition. Specifically, while a free boundary condition is determined by physics and cannot be changed by a designer (e.g. the conditions along a phreatic line or sharp interface between fresh and saline groundwater, PK-62, [46], in inverse problems an intelligent designer of an hydro-object (channel or pit bed, slope of an earth dam slope, etc.) selects a real or imaginary part of a holomorphic function as a control variable [19].

2.1 BVPs with control of velocity along the soil surface

In this subsection, we start with the simplest inverse problem. Namely, we consider the case of the velocity magnitude along BMC imaged by a circle (a solid curve in Figure 2c) of a diameter $({V_M} - {V_B})/2$ centred at the point $ - {\rm{i}}({V_M} + {V_B})/2$ in the plane u+iv. Here ${V_M}$ and ${V_B} = {V_C}$ are assumed to be the absolute values of the velocities at points M and B (C), respectively, such that ${V_M} \gt {V_B}$ . This domain G _V is exactly the same as in the Riesenkampf problem of seepage from a line source (subsurface emitter) with a phreatic line without evaporation (see PK-62, [24].

G _V in Figure 2(c) is mirrored with respect to the u-axis to get a domain in the plane of u-iv. The reference half-plane $\eta \ge 0$ of the plane $\zeta = \xi + i\eta $ is mapped onto the mirrored G _V, first by squaring the inverse Zhukovsky function and then by the corresponding linear function:

(3) \begin{align}{V_0}(\zeta ) = {\rm{i}}\left( {\frac{{{V_M} + {V_B}}}{2} - \frac{{{V_M} - {V_B}}}{2}{{\left( {\zeta + \sqrt {{\zeta ^2} - 1} } \right)}^2}} \right)\end{align}

where the branch of the square root is fixed in the half-plane $\eta \ge 0$ by the condition of its positivity for $\zeta = \xi \gt 1$ . In solution (3), velocity is zero at points $ \pm 1/\lambda $ , where $\lambda = \sqrt {1 - V_B^2/V_M^2} $ .

Obviously,

\begin{align*}{\rm{d}}z{\rm{ = }}\frac{{{\rm{d}}w}}{{{\rm{d}}\zeta }}\frac{{{\rm{d}}\zeta }}{{{V_0}(\zeta )}},\end{align*}

which upon integration yields:

(4) \begin{align}{\rm{ }}z(\zeta ) = - \frac{{2q}}{{{\rm{i}}\pi }}\int\limits_0^\zeta {\frac{1}{{{V_0}(\tau )}}} \frac{{{\rm{d}}\tau }}{{\sqrt {1 - {\tau ^2}} }}{\rm{ + i}}{y_M} = \frac{{2{\rm{i}}q}}{{\pi ({V_M} + {V_B})}}\ln \frac{{\sqrt {{\zeta ^2} - 1} - {V_B}/{V_M}\zeta }}{{\sqrt {{\zeta ^2} - 1} + \zeta }}{\rm{ + i}}{y_M}\end{align}

The logarithmic function in equation (4) has four branch points $ \pm 1$ and $ \pm {V_M}/\sqrt {V_M^2 - V_B^2} $ . The branch of function (4) is fixed in the $\zeta$ -plane with a cut along the segment $\left( { - {V_M}/\sqrt {V_M^2 - V_B^2} ,{V_M}/\sqrt {V_M^2 - V_B^2} } \right)$ . The fixed branch is real for $\varsigma = \xi \gt {V_M}/\sqrt {V_M^2 - V_B^2} \gt 1$ and equals ${\rm{i}}\pi + \log \left( {{V_B}/{V_M}} \right)$ at the point $\varsigma = 1$ on the top side of the cut.

From the condition z(1) = -L, we get two real equations (the locus of point C in the physical plane). From the equality $z(\infty ) = - {\rm{i}}$ , we get the third real equation (the locus of point S). Thus, we got three solvability conditions:

(5) \begin{align}\frac{{2q\ln ({V_M}/{V_B})}}{{\pi ({V_M} + {V_B})}} = {y_M},\;\quad \frac{{2q}}{{{V_M} + {V_B}}} = L,\quad \frac{{2q\ln (2{V_M}/({V_M} - {V_B}))}}{{\pi ({V_M} + {V_B})}} = {y_M} + 1\end{align}

Equation (5) is a system of three equations with respect to three unknowns y _m , V _M and V _B . This system has the following exact solution:

\begin{align*}{V_B} = \frac{{q{{\rm{e}}^{\pi /L}}}}{{L(1 + {{\rm{e}}^{\pi /L}})}},\quad {V_M} = \frac{q}{{L(1 + {{\rm{e}}^{\pi /L}})}} + \frac{q}{L},\quad {y_M} = \frac{L}{\pi } \ln(2 + {{\rm{e}}^{\pi /L}}) - 1.\end{align*}

which we put back into equations (3) and (4).

In the latter, we use the Re and Im Mathematica routines and the routine ParametricPlot in the interval -1< $\xi$ < 1 that gives BMC and the flow net. Figure 3(a) shows the found hump shapes for q = 1 and L = 2.5, 3 and 3.5. Figure 3(b) shows the flow net for q = 1, L = 3. The dashed lines in Figure 3(b) are equipotentials. Clearly, close to the sink the equipotentials become almost circular. The computed triad $({y_M},{V_M},{V_B})$ for the case in Figure 3(b) is (0.508, 0.42, 0.247).

Figure 3. (a) Hump BMC for q = 1 and L = 2.5, 3 and 3.5; (b) the flow net for the case q = 1, L = 3.

Example. Legostaev [36], p.108) summarised numerous soil leaching experiments in Central Asia and wrote: ‘The most effective soil leaching is attained at the infiltration rates 20-30 cm/day<V<130–150 cm/day’. Legostaev [35] estimated that for light, initially fully saturated soils the minimum leaching norm is 2200 m³/hectar. Based on these recommendations, we specify the following dimensional quantities: a sandy loam having k = 1 m/day, the drain depth is d = 0.5 m and its flow rate 2q =1 m²/day, the distance between two neighbouring drains is 2L = 10 m. Solution of equation (5), converted to dimensional quantities gives V _M = 0.49 and V _B = 0.31 m/day that is well within the Legostaev recommended bounds. The hump height is y _M = 0.35 m.

An arbitrary (not circular) case of the image of BMC in the hodograph plane is depicted in Figure 2c by a dashed curve. This curve is assumed to be symmetric with respect to the v-axis. The curve passes through the points B, C and M. The corresponding function $V(z(\zeta ))$ stands on the shoulders of ${V_0}(\zeta )$ in the following sense. Both the general $V(\zeta )$ (to be found) and already found ${V_0}(\zeta )$ (see equation (3)) have the same ${V_M}$ and ${V_B}$ . Moreover, both functions equal zero at the points $ \pm 1/\lambda $ . Next, we introduce a function

(6) \begin{align}{\rm{ }}\omega {\kern 1pt} {\kern 1pt} (z)\,{\rm{ = }}\,{\rm{ln }}\left( {V(z)/{V_0}{\rm{(}}\zeta {\rm{(}}z{\rm{)}}} \right)\,\end{align}

This auxiliary function eliminates singularities common for both $V(\zeta )$ and ${V_0}(\zeta )$ .

The functions $V(z) = u(z) - {\rm{i}}v(z) = {\rm d}w/{\rm d}z$ and ${V_0}{\rm{(}}\zeta {\rm{(z))}}$ are holomorphic and do not vanish in the domain G _z and are complex conjugate to the corresponding complexified velocity vectors. The single-valued holomorphic branch of function (6), fixed in the domain G _z, satisfies the symmetry condition $\overline {\omega ( - \bar z)} \equiv \omega (z)$ . This is a direct consequence of the identity $\overline {f( - \bar z)} \equiv - f(z)$ , which is, obviously, true for function ${V_0}{\rm{(}}\zeta {\rm{(}}z{\rm{))}}$ and both for the velocity and $V(z)$ .

We have ${\rm{ Re }}\,\omega {\kern 1pt} \,{\kern 1pt} {\rm{ = }}\,{\rm{ln |}}V/{V_0}{\rm{|}}$ and ${\rm{Im}}{\kern 1pt} \,\omega {\kern 1pt} \,{\rm{ = }}\,\,{\theta _0} - \theta $ where $\theta $ and ${\theta _0}$ are the angles that the velocity vectors $\overrightarrow V $ and $\overrightarrow {{V_0}} $ make with the x-axis.

The vectors $\overrightarrow V $ and $\overrightarrow {{V_0}} $ are directed downwards along the sides $B{B_i}$ , $C{C_i}$ and upward along ${B_i}S$ , ${C_i}S$ (see Figure 1(c)). Hence, we can fix the branch of the function in (6) for which ${\rm{Im}}{\kern 1pt} \,\omega \,\,{\rm{ = }}\,\,0$ along ${B_i}S$ and $B{B_i}$ . Then, due to the identity $\overline {\omega ( - \bar z)} \equiv \omega (z)$ , we get ${\rm{Im}}{\kern 1pt} \,\omega \,\,{\rm{ = }}\,\,0$ along ${C_i}S$ and $C{C_i}$ .

We introduce the function $\,g(\xi ) = {\rm{ln |}}V(\xi )/{V_0}(\xi ){\rm{|}}$ for $ - 1 \le \xi \le 1$ . Clearly, $g(\xi )$ should satisfy the conditions

(7) \begin{align}g(0) = 0,\quad {\rm{ }}g( \pm 1) = 0\end{align}

Below, we select shape-controlling functions $g(\xi )$ that satisfy conditions (7). Consequently, we obtain the so-called mixed BVP for $\omega {\kern 1pt} {\kern 1pt} (\zeta )$ :

(8) \begin{align}\begin{array}{l}{\rm{Re }}\,\omega {\kern 1pt} {\kern 1pt} (\xi ) = g(\xi ),\quad {\rm{ }} - 1 \le \xi \le 1,\\{\mathop{\rm Im}\nolimits}\, \omega {\kern 1pt} {\kern 1pt} (\xi ) = 0,\quad{\rm{ }}\xi \le - 1,\quad \xi \ge 1\;.{\rm{ }}\end{array}\end{align}

We have to solve this BVP in the class of functions bounded at the points $\zeta = \pm 1$ and at infinity. A unique solution to this problem is given by the Keldysh-Sedov formula (see e.g. [16]:

(9) \begin{align}\omega (\zeta ) = \frac{{\sqrt {1 - {\zeta ^2}} }}{{\pi {\rm{i}}}}\int\limits_{ - 1}^1 {\frac{{g(\tau ){\rm{d}}\tau }}{{\sqrt {1 - {\tau ^2}} (\tau - \zeta )}}}, \end{align}

where the branch of the function ${\rm{ }}\sqrt {1 - {\zeta ^2}} $ is fixed in the upper half of the $\zeta $ – plane to be positive at ${\rm{ - 1 \lt }}\zeta = \xi \lt 1$ .

The chosen branch of the radical ${\rm{ }}\sqrt {1 - {\zeta ^2}} $ satisfies the symmetry condition

(10) \begin{align}\overline {f( - \bar \zeta )} \equiv f(\zeta )\end{align}

Therefore, function (9) meets condition (10) if the function $g(\tau )$ satisfies the condition $g(\tau ) \equiv g( - \tau )$ , which we enforce on our class of shape controls.

The integral of the Cauchy type in equation (9) is singular at ${\rm{ }}\zeta = \xi $ , for $ - 1 \lt \xi \lt 1$ . It is calculated using the Plemelj–Sokhotski formulas with the integral term evaluated in the sense of v.p. [16]. For instance, the PrincipalValue routine of Wolfram’s (1991) Mathematica can be used.

From equations (2), (3), (6), (9) we get:

(11) \begin{align}{\rm{ }}z(\zeta ) = \frac{{2q}}{{{\rm{i}}\pi }}\int\limits_0^\zeta {\frac{{\exp ( - \omega (\tau ))}}{{{V_0}(\tau )}}} \frac{{{\rm{d}}\tau }}{{\sqrt {1 - {\tau ^2}} }}{\rm{ + i}}{y_M}\end{align}

Note that the solution (11) satisfies the condition $\overline {z( - \bar \zeta )} \equiv - z(\zeta )$ , as it should be. Besides, function (11) should satisfy the conditions

(12) \begin{align}z( \pm 1) = \mp L,\,\,\,z(\infty ) = - {\rm{i}}\end{align}

and the condition $z( \pm 1/\lambda ) = \infty $ which is obviously met, since function (3) has simple zeros at the points $ \pm 1/\lambda $ , if $\lambda = \sqrt {1 - V_B^2/V_M^2} $ .

For example, as a control function $\,g(\xi ) = {\rm{ln |}}V(\xi )/{V_0}(\xi ){\rm{|}}$ , we choose the following family:

(13) \begin{align}g(\xi ) = a{(1 - {\xi ^2})^\alpha }|\xi {|^{2\beta }}\end{align}

where $\alpha $ , $\beta $ and a are positive parameters. Clearly, all functions of the family (13) obey conditions (7).

Equation (9) with kernel (13) is integrated as follows:

(14) \begin{align}{\rm{ }}\omega (\zeta ) = \frac{{{\rm{i}}a\sqrt {1 - {\zeta ^2}} }}{{\pi \zeta }}{\mathop{\rm B}\nolimits} (1/2 + \alpha ,1/2 + \beta ){\mathop{\rm F}\nolimits} \left( {1,1/2 + \beta ;1 + \alpha + \beta ;{\zeta ^{ - 2}}} \right)\end{align}

for $|\zeta |\gt 1$ , and

(15) \begin{align}{\rm{ }}\omega (\zeta ) = {\rm{ - }}\frac{{{\rm{i}}a}}{\pi }{\mathop{\rm B}\nolimits} ({\textstyle{1 \over 2}} + \alpha ,\beta - {\textstyle{1 \over 2}})\zeta \sqrt {1 - {\zeta ^2}} {\mathop{\rm F}\nolimits} (1,1 - \alpha - \beta ;{\textstyle{3 \over 2}} - \beta ;{\zeta ^2}) + a(1 - {\rm{i}}\tan \pi \beta ){\zeta ^{2\beta }}{(1 - {\zeta ^2})^\alpha }\end{align}

for $|\zeta | \lt 1$ due to equation (15.3.7) from [1]. Here ${\mathop{\rm B}\nolimits} (x,y)$ and F are the Beta function ([1], equations 6.2.1,2) and the hypergeometric function, respectively. In equation (15), the branches of the analytic functions $\sqrt {1 - {\zeta ^2}} $ , ${(1 - {\zeta ^2})^\alpha }$ are fixed in the $\zeta$ -plane cut along the rays $( - \infty , - 1)$ and $(1, + \infty )$ . Both functions, $\sqrt {1 - {\zeta ^2}} $ , ${(1 - {\zeta ^2})^\alpha }$ , are positive for $ - 1 \le \xi \le 1$ . The branch of the function ${\zeta ^{2\beta }}$ is fixed in the $\zeta$ -plane cut along the ray $(0,\infty )$ . This fixation is attained by the condition $0 \lt \arg \varsigma \lt 2\pi $ . Thus, our formula (15) in the MS is correct in the unit circle cut along the interval (0,1).

The conditions $z( \pm 1) = \mp L$ have to be fulfilled. We select $z(1) = - L$ . We note, that the condition $z( - 1) = L$ is also met due to the symmetry property of the solution (11). This complex condition leads to two real ones:

\begin{align*}{\mathop{\rm Re}\nolimits} \left( {\frac{{2q}}{\pi }\int\limits_0^1 {\frac{{\exp ( - \omega (\tau ))}}{{{V_0}(\tau )}}} \frac{{{\rm{d}}\tau }}{{\sqrt {1 - {\tau ^2}} }}} \right) = {y_{M,}}\quad {\mathop{\rm Im}\nolimits} \left( {\frac{{2q}}{\pi }\int\limits_0^1 {\frac{{\exp ( - \omega (\tau ))}}{{{V_0}(\tau )}}} \frac{{{\rm{d}}\tau }}{{\sqrt {1 - {\tau ^2}} }}} \right) = - L.\end{align*}

which, upon substitution $\tau = {\mathop{\rm Cos}\nolimits} \varphi $ , are transformed to

(16) \begin{align}\begin{array}{l}\dfrac{{2q}}{\pi }\int\limits_0^{\pi /2} {\dfrac{{\exp ( - a{{{\mathop{\rm Cos}\nolimits} }^{2\beta }}\varphi {\textrm{Sin}^{2\alpha }}\varphi )}}{{\sqrt {V_M^2{\textrm{Sin}^2}{\kern 1pt} \varphi + V_B^2{\textrm{Cos}^2}{\kern 1pt} \varphi } }}} {\mathop{\rm Sin}\nolimits} \gamma (\varphi )d\varphi = {y_{M,}}\quad \\[20pt] \dfrac{{2q}}{\pi }\int\limits_0^{\pi /2} {\dfrac{{\exp ( - a{{{\mathop{\rm Cos}\nolimits} }^{2\beta }}\varphi {\textrm{Sin}^{2\alpha }}\varphi )}}{{\sqrt {V_M^2{\textrm{Sin}^2}{\kern 1pt} \varphi + V_B^2{\textrm{Cos}^2}{\kern 1pt} \varphi } }}} {\mathop{\rm Cos}\nolimits} \gamma (\varphi )d\varphi = L,\end{array}\end{align}

where

\begin{align*}\begin{array}{l}\gamma (\varphi ) = \dfrac{a}{{2{\kern 1pt} {\kern 1pt} \pi }}\textrm{Sin}2\varphi {\mathop{\rm B}\nolimits} ({\textstyle{1 \over 2}} + \alpha ,\beta - {\textstyle{1 \over 2}}){\mathop{\rm F}\nolimits} (1,1 - \alpha - \beta ;{\textstyle{3 \over 2}} - \beta ;\textrm{Co}{\textrm{s}^2}{\kern 1pt} \varphi ) + \\[8pt] \quad \quad \quad \textrm{aTan}\pi \beta {{\mathop{\rm Cos}\nolimits} ^{2\beta }}\varphi \textrm{Si}{\textrm{n}^{2\alpha }}\varphi + {\mathop{\rm ArcSin}\nolimits} \dfrac{{{V_M}\textrm{Sin}{\kern 1pt} \varphi }}{{\sqrt {V_M^2\textrm{Si}{\textrm{n}^2}{\kern 1pt} \varphi + V_B^2\textrm{Co}{\textrm{s}^2}{\kern 1pt} \varphi } }} - \varphi .\end{array}\end{align*}

The second condition (12) can be written as follows:

(17) \begin{align}\frac{{2q}}{\pi }\int\limits_0^{\log (1 + \sqrt 2 )} {\left( {\frac{{\exp ( - {\gamma _1}(\varphi ))}}{{{V_M}{\mathop{\rm Cosh}\nolimits} \varphi - {V_B}{\mathop{\rm Sinh}\nolimits} \varphi }} + \frac{{\exp ( - {\gamma _2}(\varphi )){\rm Tanh}(\varphi /2)}}{{{V_M}{\mathop{\rm Cosh}\nolimits} \varphi - {V_B}}}} \right)} d\varphi = {y_M} + 1\end{align}

where

\begin{align*}{\gamma _1}(\varphi ) = \frac{a}{{2\pi }}{\rm Sinh}2\varphi {\mathop{\rm B}\nolimits} ({\textstyle{1 \over 2}} + \alpha ,\beta - {\textstyle{1 \over 2}}){\mathop{\rm F}\nolimits} (1,1 - \alpha - \beta ;{\textstyle{3 \over 2}} - \beta ; - {{\mathop{\rm Sinh}\nolimits} ^2}{\kern 1pt} \varphi ) + \frac{{a{{{\mathop{\rm Cosh}\nolimits} }^{2\alpha }}\varphi {{\rm Sinh}^{2\beta }}\varphi }}{{{\mathop{\rm Cos}\nolimits} \pi \beta }} + \varphi ,\end{align*}

\begin{align*}{\gamma _2}(\varphi ) = \frac{a}{\pi }{\rm Cosh}\varphi {\mathop{\rm B}\nolimits} ({\textstyle{1 \over 2}} + \alpha ,{\textstyle{1 \over 2}} + \beta ){\mathop{\rm F}\nolimits} (1,{\textstyle{1 \over 2}} + \beta ;1 + \alpha + \beta ; - {{\mathop{\rm Sinh}\nolimits} ^2}{\kern 1pt} \varphi ).\end{align*}

Similarly to the above described special case of eqns. (4)-(5) we – upon specification of the three parameters in equation (13) – solved a system of three nonlinear equations by the FindRoot and evaluated y _M , V _M and V _B . After that we plotted BMC.

Generally, we can consider the system of equations (16), (17) as a system with respect to two triads of six parameters: (a, $\alpha $ , $\beta $ ) and (y _M, V _M , V _B ). By fixing either triad, we get a system of equations with respect to another triad. If we fix (a, $\alpha $ , $\beta $ ), then the system (16), (17) directly gives (y _M, V _M , V _B ) without any problems, albeit reconstruction of the flow net using the corresponding function (11) requires a significant computation time in Mathematica.

If we choose special values of $\alpha = 1$ , $\beta = 1$ in equation (13), then all functions in eqns. (14), (15) and $\gamma (\varphi )$ , ${\gamma _1}(\varphi )$ , ${\gamma _2}(\varphi )$ in equations (16), (17) are expressed in elementary functions. Namely,

\begin{align*}\omega (\zeta ) = - \frac{a}{2}\zeta \sqrt {{\zeta ^2} - 1} {\left( {\zeta - \sqrt {{\zeta ^2} - 1} } \right)^2},\end{align*}

\begin{align*}\gamma (\varphi ) = - \frac{a}{8}{\rm Sin}\,4\varphi + {\mathop{\rm ArcSin}\nolimits} \frac{{{V_M}{\rm Sin}{\kern 1pt} \varphi }}{{\sqrt {V_M^2{\textrm{Sin}^2}{\kern 1pt} \varphi + V_B^2{\textrm{Cos}^2}{\kern 1pt} \varphi } }} - \varphi .,\end{align*}

\begin{align*}{\gamma _1}(\varphi ) = {a \over 8}{\rm Sinh}\,4\varphi - {a \over 4}{\textrm{Sinh}^2}\,2\varphi + \varphi ,\,\,{\gamma _2}(\varphi ) = {a \over 2}{\mathop{\rm Cosh}\nolimits} \varphi /{\left( {1 + {\mathop{\rm Cosh}\nolimits}\, \varphi } \right)^2}\end{align*}

As an illustration, we solved the system (16), (17) for the parameters L = 3, q = 1 and a = 1, 2, 3, 4. The corresponding surfaces BMC and flow nets for the first three values of a are plotted in Figure 4.

Figure 4. Upper panel: the soil surface profiles for the control function (13), L = 3, q = 1 and a = 1, 2, 3, 4. Lower panel (from left to right): flow nets at L = 3, q = 1 and a = 1, 2, 3, correspondingly.

As Figure 4 illustrates, the family of controls (13) generates both a single-bump BMC (curve 1) and BMC with a maximum and two minima, the soil contours resembling a cross-section of a propelling sting ray (curves 2–4), which would agroengineeringly mean an excessive bulldozing of the original flat soil surface.

2.2 Mixed BVP for $z(\zeta )$

In this subsection, we consider the limit of $L \to \infty $ that corresponds to Vedernikov’s solitary drain shown in his Figure 74 and solution in his pp. 178–181. Similarly to Kacimov [21], we shape the ponded soil surface. We formulate a mixed BVP:

(18) \begin{align}\begin{array}{l}{\rm{ }}y(\xi ) = G(\xi ),\;\quad G( - \xi ) = G(\xi ),{\rm{ }} - 1 \le \xi \le 1,\\[4pt] {\rm{ }}y(0) = {y_M},{\rm{ }}y{\rm{(}} \pm 1) = 0,\quad y\prime ( \pm 1) = 0,\\[4pt]{\rm{ }}x(\xi ) = 0,{\rm{ }}1 \lt |\xi |,\\[4pt]{\rm{ }}x( \pm 1) = m {\rm{ }}\infty ,\end{array}\end{align}

for an analytic function $z(\zeta )$ in the upper half-plane of Figure 2(b). The shape-controlling function $G(\xi )$ belongs to a broad class (e.g. Holder’s one) that suffices for existence of singular integrals below. The dimensionless height of point M is a given constant y _M . Equations (18) show that $z(\zeta )$ should have singularities $z(\zeta ):\; {(1 - {\zeta ^2})^{ - 1/2}}\;{\rm{ at }}\;\zeta \to \pm 1$ . The index of this problem is 1 in the class of functions holomorphic in the upper half-plane and bounded at infinity.

According to the Keldysh-Sedov formula [16], a general solution of BVP (18) is

(19) \begin{align}{\rm{ }}z(\zeta ) = \frac{1}{{\pi \sqrt {1 - {\zeta ^2}} }}\left( {\int\limits_{ - 1}^1 {G(\tau )} \frac{{\sqrt {1 - {\tau ^2}} }}{{(\tau - \zeta )}}{\rm{d}}\tau + {c_1}\zeta + {c_2}} \right)\end{align}

where ${c_1},\;{c_2}$ are arbitrary real constants. These constants should be found from two conditions:

a given locus of the sink $z(\infty ) = - {\rm{i}}$ and the symmetry identity $\overline {z( - \bar \zeta )} \equiv - z(\zeta )$ . The former yields ${c_1} = - \pi $ , and the latter gives ${c_2} = 0$ . Thus, from (19) we get

(20) \begin{align}{\rm{ }}z(\zeta ) = \frac{1}{{\pi \sqrt {1 - {\zeta ^2}} }}\int\limits_{ - 1}^1 {G(\tau )} \frac{{\sqrt {1 - {\tau ^2}} }}{{(\tau - \zeta )}}{\rm{d}}\tau - \frac{\zeta }{{\sqrt {1 - {\zeta ^2}} }}\end{align}

We selected the following classes of shape-controlling functions:

(21) \begin{align}G(\xi ) = {y_M}{(1 - {\xi ^2})^{1 + \gamma }}({\rm{a}}),\,\,G(\xi ) = {y_M}{(1 - {\xi ^2})^{1 + \gamma }}{\rm{Cos [}}\pi {\rm{ }}c{\rm{ }}\xi /2{\rm{]}}\ ({\rm{b}}),\end{align}

where $\gamma $ is a positive parameter. Substituting the kernel function (21) (a) into equation (20) and integrating, we obtain

(22) \begin{align}{\rm{ }}z(\zeta ) = - {y_M}\frac{{\Gamma (5/2 + \gamma ){\mathop{\rm F}\nolimits} (1/2,1;3 + \gamma ;1/{\zeta ^2})}}{{\sqrt \pi \Gamma (3 + \gamma )\zeta \sqrt {1 - {\zeta ^2}} }} - \frac{\zeta }{{\sqrt {(1 - {\zeta ^2})} }}\end{align}

in the vicinity of infinity, and

(23) \begin{align}{\rm{ }}z(\zeta ) = - {y_M}\left( {\frac{{2\Gamma (5/2 + \gamma )}}{{\sqrt \pi \Gamma (2 + \gamma )}}\frac{{\zeta {\mathop{\rm F}\nolimits} (1, - 1 - \gamma ;3/2;{\zeta ^2})}}{{\sqrt {1 - {\zeta ^2}} }} - {\rm{i(1}} - {\zeta ^2}{)^{1 + \gamma }}} \right) - \frac{\zeta }{{\sqrt {(1 - {\zeta ^2})} }}\end{align}

in the vicinity of $\zeta = 0$ . In equations (22), (23), and elsewhere, $\Gamma $ is the gamma function. Obviously, in the limit y _M = 0 our solution (22)–(23) degenerates into the Vedernikov [48] flat soil surface geometry, his equation (384), i.e. a combination of a line sink and line source (e.g. a dipole), located distance 2d apart for which the equipotential lines and streamlines make two families of mutually orthogonal circles.

Figure 5(a) shows stream (solid) and equipotential (dashed) lines found using equations (22)–(23) for $z(\zeta (\phi + {\rm{i}}\psi ))$ with q = 1, ${y_M} = 0.5$ , $\gamma = 1$ and $\psi = 0.3n$ , $n = \overline { - 9,9} $ ; $\phi = 0.3n$ , $n = \overline {1,6} $ . Similarly to Figure 3(b), in Figure 5(a), the near-sink equipotentials are almost circular.

Figure 5. Flow net for the control function (21) with ${y_M} = 0.5$ , $\gamma = 1$ .

For the class of kernels given by equation (21)(a), we involved the NIntegrate routine of Mathematica, with the v.p. (principal value) option in evaluation of singular integrals for plotting the soil surface contour BMC. For reconstruction of the flow net, we used a standard NIntegrate routine. As an example, in equation (21)(b), we selected the parameters q = 3, $\gamma = 1$ , y _M = 1, c = 1. In Figure 5(b), the contour BMC (the line $\phi $ = 0) ‘caps’ the flow domain G _z . For the right half of G _z , the flow net is plotted. The streamlines $\psi $ = 0.5, 1, 1.5 and 2 are solid lines. Obviously, the streamlines $\psi $ = 0 and 3.0 are the straight segment [-1,1] and the ray $[ - 1, - \infty ],$ correspondingly. The equipotential lines $\phi $ = 0.5, 1, 1.5 are plotted as dotted lines. Next, in equation (21)(b) we selected parameters q = 3, $\gamma = 1$ , y _M = 1, c = 2. The corresponding contour BMC is shown in Figure 5(b) as a dashed contour. It has two minimum points (compare with BMC in Figure 4). The flow net for this ‘fancy’ soil surface is not plotted in Figure 5(b) to avoid cluttering.

2.3 Travel time along streamlines and advective breakthrough curves

It is well known (see e.g. [22, 28, 29] that a purely advective travel time (dimensionless, introduced above) of a marked particle moving along any streamline US (Figure 1(c)) is evaluated by the integral along this flow path:

(24) \begin{align}{\rm{ }}{T_{US}}(\psi ) = \int\limits_0^\infty {\frac{{{\rm{d}}\phi }}{{|V(\phi ,\psi ){|^2}}}} ,{\rm{ }} - q \le \psi \le q\end{align}

This function determines the dispersion-free break-through curve, BTC (Figure 1(d)) of the drain. Unfortunately, even bulging of BMC in Figure 1(c) cannot eliminate the equality ${\rm{ }}T( \pm q) = \infty $ , i.e. the long tail of BTC ² . However, an intelligent designer can try to make BTC as flat as possible by selection of the kernel functions, e.g. (13) or (21). We formulate two criteria of tracer advection and corresponding leaching efficiency: one is local, i.e. applied to a single (the only straight) streamline MS, and another is integral, based on ‘uniformity’ of the breakthrough curve:

(25) \begin{align}\begin{array}{l}{T_{MS}} = \int\limits_0^\infty {\dfrac{{{\rm{d}}\phi }}{{|V(\phi ,0){|^2}}}} = \int\limits_{{y_M}}^{ - 1} {\dfrac{{{\rm{d}}y}}{{|V(0,y){|^{}}}}} ,{\rm{ }}\\[16pt] Cr{\rm{ = }}\dfrac{1}{{2q}}\int\limits_{ - {\psi _c}}^{{\psi _c}} {T(\psi )} \,{\rm{d}}\psi ,\end{array}\end{align}

where ${\psi _c}$ is a given constant (less than q).

For comparisons, for a single drain-sink under a flat ponded soil surface, the potential of which is given by formula (6.2) in PK-77, p.353 we have from equations (22)–(23):

(26) \begin{align}{z_0}(w) = - \tan \frac{{{\rm{i}}\pi w}}{{2q}}\end{align}

The Vedernikov [48] equation (384) (in our dimensional notations) reads:

(27) \begin{align}{z_{0V}}\left( w \right) = - \sqrt {|MC|(|MC| + 2r)} \tan \frac{{{\rm{i}}\pi w}}{{2{\rm{q}}}}\end{align}

where r (commonly few cm) is the radius of a circular equipotential contour (drain) and |MC| is the distance between the drain apex and a flat soil surface. Obviously, at practical values of d and r, equations (26) and (27) almost coincide. Then from equation (26) along the only straight and the shortest streamline MS the vertical component of velocity is $v\left( y \right) = V(0,y) = - \frac{{2q}}{{\pi (1 - {y^2})}}$ . Consequently, from equation (25), the advection time for the Vedernikov case in Figure 1(a) is T _MSV = $\pi$ /(3 q). For instance, for q = 3, this time is T _MSV = 0.35.

Using equation (24), in Figure 6(a) we plotted ${T_{US}}(\psi )$ for soils profiled by the control (21) with $\gamma = 1,\,2,\,3$ (curves 1–3, correspondingly), q = 3 and a fixed y _M = 0.5. As is evident from Figure 6(a), the travel time is well uniformised for $|\psi | \lt 1$ . Larger values of the stream function, $1 \lt {\kern 1pt} \,|\psi |\,\lt 3$ , determine the ‘tail’ of the BTC. Figure 6(b) presents Cr( $\gamma $ ) according to equation (25) for the same class of controls, y _M = 0.3, 0.4, 0.5 (curves 1–3, correspondingly) and ${\psi _c} = 2$ .

Figure 6. Contour BMC for q = 3, $\gamma = 1$ , y _M = 1, c = 1 in equation (21) (b) and the corresponding flow net. A dashed line shows BMC for q = 3, $\gamma = 1$ , y _M = 1, c = 2.

Figure 7. (a) Advective travel time distribution along streamlines ${T_{US}}(\psi )$ , water particles seeping from the protruding ponded soil surface to the drain for the control function equation (21), $\gamma = 1,\,2,\,3$ (curves 1–3, correspondingly), q = 3 and y _M = 0.5; (b) the uniformity coefficient Cr( $\gamma $ ) for y _M = 0.3, 0.4 and 0.5 (curves 1–3, correspondingly), ${\psi _c} = 2$ .

We used the NIntergate routine of Mathematica to compute the integrals involved in making Figure 6.

3. Concluding remarks

Legostaev [36], p. 142) emphasised that ‘…only an impeccably functioning drainage facilitates intensive desalination of soils and groundwater’. Therefore, an adequate mathematical modelling of pore water motion towards agricultural drains is a must for keeping high and time-wise stable (sustainable for decades of irrigation) WUE.

Micro-profiling of a flat soil surface in orchards and cropfields with perennials is already used as a technology for increasing WUE. Indeed, a tree-focused soil relief funnels the runoff to the near crown zone (see e.g. [43]. However, in adverse climatic, irrigation and precipitation-runoff conditions, this micro-relief can spur both secondary salinisation (or ‘return of leached salts’ in terminology of [49] and waterlogging [36]. In terms of labour-cost expenses, profiling of the soil surface suggested in this paper is similar to one reported by Postel [43].

Our technique is applicable to ponded-leached and drained cropfields in hot-arid climates and is potentially capable to improve the topology of Darcian seepage towards line sinks by oppressing the vice of a ‘short cut’ (‘preferential flow’) between two equipotential lines, viz. the ponded soil surface and the gravel pack of the drain. During the cultivation stage, crops can be planted in the topographic depressions (Figure 1), i.e. two birds (‘uniform leaching’ through desalinised soil and ‘funneled irrigation’ by fresh water) can be killed by the same profiling bullet. We suggested two mathematical methods of designing the shape of a ponded soil surface and ensued seepage flow net: by variation of the magnitude of the Darcian velocity and of the vertical coordinate of a bulged land surface. The obtained explicit analytical solutions for characteristic anti- and holomorphic functions (complexified velocity, complex potential and complex physical coordinate) allow for reconstruction of the shape of the land surface and evaluation of the advective travel time along streamlines in a piston-type displacement of pore water pushed through the soil by its surface ponding. If the flow net in the flow domain, capped by the designed surface ‘hump’, is found, then the hydrodynamic dispersion along stream lines can be easily added (see e.g. Frenkel, 1978). Other straightforward generalisations of our solutions are

• The water table with groundwater supplied from a highly permeable substratum (PK-62) can be taken into account. In this case, G _W is a half-strip with two vertical cuts. Such polygon (enneagon) can be easily mapped onto a reference half-plane by the Schwartz-Christoffel formula.

• Mathematical shaping of BMC can be done by specification of $\theta (\xi )$ and $x(\xi )$ in Sections 2.1 and 2.2, respectively, that leads to the Dirichlet rather than mixed BVPs for characteristic functions. Integral solutions to these BVPs are easily expressed [16], analogously to equations (11) and (20).

• If the ponding depth h _p > y _M in Figure 1(c), then our solution covers the surface water regimes when the head drop between the soil surface and drain tube decreases with time. Indeed, for a rigid porous skeleton and incompressible pore water, as long as the soil remains fully saturated, the Laplace equation is valid for characteristic functions and the transience in h _d does not alter our analytical formulae.