RESEARCH ARTICLE
10.1029/2018WR024025
Key Points:• For a given total thickness of a
two‐layered soil the thickness of theupper layer minimizing or
maximizing evaporation ﬂux is
found• Gardner's, Philip's, and van
Genuchten's soils are studied in 1‐D
ﬂows from a static water table to dry
topsoil• Analytical solutions of Richards
ODEs are juxtaposed with HYDRUS
ﬁnite element computations

Minimizing Evaporation by Optimal Layeringof Topsoil: Revisiting Ovsinsky's Smart
Mulching‐Tillage Technology Via
Gardner‐Warrick's Unsaturated
Analytical Model and HYDRUS
A. R. Kacimov1, Yu.V. Obnosov2, and J. Šimůnek3

1

Department of Soils, Water and Agricultural Engineering, Sultan Qaboos University, Muscat, Oman, 2Institute ofMathematics and Mechanics, Kazan Federal University, Kazan, Russia, 3Department of Environmental Sciences,
University of California, Riverside, CA, USA

Abstract
Correspondence to:A. R. Kacimov,anvar@squ.edu.om;
akacimov@gmail.com

Citation:
Kacimov, A. R., Obnosov, Y. V., &Šimůnek, J. (2019). Minimizing
evaporation by optimal layering of
topsoil: Revisiting Ovsinsky's smart
mulching‐tillage technology via
Gardner‐Warrick's unsaturated
analytical model and HYDRUS. Water
Resources Research, 55. https://doi.org/
10.1029/2018WR024025

Ovsinsky (1899, https://www.rulit.me/books/novaya‐sistema‐zemledeliya‐read‐193251‐1.
html) suggested and tested a water conserving soil no‐till technology for rain‐snow‐fed ﬁeld crops in asemiarid environment in southern Russia. We model Ovsynsky's unsaturated ﬂow fragment, in which 1‐Dsteady evaporation and evapotranspiration through a two‐layered soil from a horizontal static water table toa dry soil surface takes place. Gardner's exponential and algebraic functions are used for the unsaturated
hydraulic conductivity‐suction head relations. The vertical evaporation ﬂux depends on the dyads and triads(correspondingly) of the parameters of these functions, for example, the saturated hydraulic conductivity
and the sorptive number of the two layers. The ﬂux, as a function of the relative thickness of the upper
stratum, is analytically found from the solution of one or two nonlinear equations. This relation can be
nonmonotonic and exhibits either a minimum or maximum depending on whether this stratum is coarser or
ﬁner than the subjacent stratum fed from a horizontal isobar. HYDRUS‐1D simulations conﬁrm these
extrema. This explains the experimental results from the literature on mulching/tillage/soil crusting‐sealing,which can increase, decrease, or have no impact on evaporation from a shallow water table. Alterations of
the soil's homogeneity to reduce evaporation losses can improve the hydrological balance of soil proﬁles.

Received 3 SEP 2018
Accepted 11 MAR 2019
Accepted article online 21 MAR 2019

1. Introduction©2019. American Geophysical Union.
All Rights Reserved.

Crop ﬁelds are sometimes cultivated on soils with a shallow freshwater table (line W1W2 in Figure 1). This
isobar is maintained either naturally by groundwater from an unconﬁned aquifer or by Kornev's (1935) subsurface irrigation. Hydrologically, such agro‐engineering system illustrates the water use efﬁciency‐
constrained dichotomy: soil water should be available to plant roots for transpiration, while evaporationfrom the soil surface (line S2S1 in Figure 1), with ensued secondary salinization, is deleterious to the agro‐
environment. The dilemma, formulated as an optimal control problem, is how to maximize (or maintainat the level of physiological requirements of a particular crop) the former and minimize the latter. In
1899, the Russian agronomist and philosopher Ovsinsky (It is noteworthy that Ovsinsky had to change eventhe title of some chapters of his book because of the animosity from his opponents in the academia. For
example, his original title for section 1 of the book can be translated as “Self‐reﬂection of plants” or even
“Ontological endeavours of plants.” Ovsinsky considered plants as biological creatures having a higher levelof organization than animals (including Aristotelian animals). Ovsinsky drew the following political analogy: plants are republics and animals are monarchies. In the book, Ovsinsky also suggested to upgrade eventhe terminology: agronomy ‐> husbandry of plants as conscious and even intelligent creatures. Ovsinsky perpended on plants as creatures cognizant of farmer's agronomic actions. These Ovsinsky's ideas were vehemently opposed by the contemporaneous crop/plant scientists.) suggested and tested the technology of ano‐till treatment of topsoil in semiarid crop ﬁelds. The main idea was very simple: to convert an initiallyhomogeneous (quasi‐homogeneous) soil body of a given thickness D between W1W2 and S2S1 into a two‐layered composite (Figure 1). The plant roots uptake moisture from this composite, which is agronomicallyfavorable. Evaporation from the soil surface, although relatively less (as compared with surface irrigation), isbad because it reduces the water use efﬁciency and causes secondary salinization, among other negative

KACIMOV ET AL.

1

Water Resources Research

10.1029/2018WR024025

consequences. The Ovsinsky's thickness D of the unsaturated zone in
Figure 1 was high enough to serve as an antievaporation sheath. Thispartially saturated layer also served in Ovsinsky's agronomy as athermal insulator for the plant roots, provided they penetrate deeply
from the superheated soil surface. This dual smart wetting and heat‐stress mitigation of the unsaturated layer has been corroborated bynumerous data on improved harvest/biomass/plant morphologies for avariety of crops.

Figure 1. Two‐layered soil with a vertical ascending evaporative ﬂux fromthe water table at a depth D. Two vertical axes are as follows: Philip's OZ
is oriented downward and OZG is oriented upward. B is the thickness of theupper layer.

Mathematically, Ovsinsky (1899) minimized the evaporative ﬂux, q, by
making the hydrophysical parameters of a layer 2 of thickness B(Figure 1) different from those of the original soil (layer 1). From thisviewpoint, historically, crop cultivation by practices including tillage,viz., plowing, harrowing (post‐plowing dragging of a heavy frame set withteeth or tines to break up soil clods), spudding (rooting‐digging out weeds
by a sharp spadelike tool), mulching, and treatment of the topsoil withpolymers, among others.

After most treatments, the topsoil of thickness B (Figure 1) becomes coarser than a substratum of thickness D‐B (see, e.g., Bodner et al., 2015,
Connolly, 1998, Fuchs & Hadas, 2011, Hillel et al., 1975, Mehari et al.,2011, Minhas et al., 1986, Strudley et al., 2008, Verburg et al., 2012). Rains or intensive sprinkling produce
an opposite effect: a thin (few mm) soil crust (seal) forms due to the impact of water drops on a relatively
ﬁne bare soil surface (see, e.g., Mualem & Assouline, 1989). The natural morphing of soils also results in asimilar vertical stratiﬁcation caused by illuviation, eluviation, lessivage, and cyanomat‐lichen‐moss development, among others (see, e.g., Felde et al., 2014), which evolve over decades to centuries unlike a commonly
annual soil treatment in agronomy. In any type of heterogenization, the engineered or natural porous rectangle M1M2S2S1 has hydraulic and capillary properties that sharply or continuously contrast with thoseof the original subjacent soil in the rectangle W1W2M2M1 (Figure 1).

In this paper, we consider ascending steady unsaturated ﬂow (Figure 1) from an isobaric horizon W1W2
through a two‐layered rectangle W1W2S2S1. We use the unsaturated conductivity functions of Gardner(1958), Willis (1960), Warrick (1988), and van Genuchten (1980) for analyzing the relation q(B). We use
modern computer algebra routines to integrate and solve systems of nonlinear equations and ﬁnd extremaof q(B). We also simulate evaporation for the system in Figure 1 using HYDRUS‐1D (Šimůnek et al.,2016). We assume that ﬂow is Darcian, one‐phase (vapor/gas and solute motions, depositions, and dissolutions are neglected), and isothermal and that the stratiﬁed soils are isotropic (Warrick, 2003).
We answer the following two (main) questions:
1. For a given texture (hydraulic and capillary properties) of each of the two strata in Figure 1 and for a ﬁxed
D, is there an optimal depth B that minimizes the losses q?
2. For ﬁxed B and D, can one optimize the texture of the upper layer with an objective of reducing q?
It is noteworthy that for saturated steady ﬂows for which the governing equation is linear, similar optimization problems have been recently analytically solved for 2‐D ﬂow (Kacimov & Obnosov, 2018). For unsaturated ﬂows studied below, the nonlinearity of the governing (Richards) equations makes 2‐D analyticalsolutions in composite soils prohibitively complex.

2. Analytical Solutions
As in Gardner (1958), Warrick (1988), and Willis (1960) we assume that along W2W1 in the ﬁrst soil layer, thepressure head p1 = 0. Moisture is lifted by the capillarity of two layers from W2W1 to a dry soil surface.
Gravity and the Darcian resistance of both layers oppose the evaporation‐maintained suction.
For most rain‐fed and traditionally irrigated soils in agronomy and land surfaces of arid catchments inhydrology, evaporation is a transient 2–3 stage process (actually, a phase in soil water redistribution). The
boundary conditions for moisture content or suction pressure at S1S2 in real ﬁelds vary with time,
KACIMOV ET AL.

2

Water Resources Research

10.1029/2018WR024025

depending on plants' physiological stages and root evolution, atmospheric conditions (relative humidity of
the air, its temperature, wind speed, and solar radiation), desiccation of the topsoil and ﬂuctuations in thelocus of the water table (see e.g., Boast & Simmons, 2005; Jalota & Prihar, 1990; Stewart & Broadbridge,1999). In models, this requires solving nonlinear transient partial differential equations (PDEs) for each soil
layer, taking into account generally hysteretic constituting relations (phase relative permeabilities and capillary pressure functions) imbedded in those PDEs. Consequently, simpliﬁcations and approximations of thecorresponding boundary value problems (hereafter abbreviated as BVPs) are common (e.g., see Assouline
et al., 2014; De Luca & Cepeda, 2016; Warrick, 1988). Field and laboratory studies of transient evaporation
are performed using lysimeters, tensiometers, theta‐probes, sigma‐probes, thermometers, and other instruments to measure and assess moisture, solute, and temperature dynamics in soils. Mathematical models
combine soil physics, shallow aquifer hydrology, and the near‐surface atmosphere (e.g., see Al‐Shukaili,
2018; Ehlers & Van Der Ploeg, 1976; Geng & Boufadel, 2015; Hillel & Talpaz, 1977; Khan, 1988; Li et al.,2016; Malik et al., 1992; Novák, 2012; Sadeghi et al., 2014; Soylu et al., 2011; Xiao et al., 2010; Zarei et al.,2010). Civil engineers most often ignore the plants and address the transient desiccation of clay liners.
Engineers solve the same BVPs as agronomists but are mainly interested in cracking of layer 2 (inFigure 1); for them, q is of a minor concern. Correspondingly, their models of unsaturated ﬂows with evaporation utilize PDEs that account for soil swelling and shrinkage (see, e.g., Zhou & Rowe, 2005).
As in Gardner (1958), we assume that the soil surface S1S2 has a constant negative pressure head p2 = −ps
and that ps = const > 0. Gardner (1958) suggested several simple equations k(S) for the unsaturated hydraulic conductivity k as a function of the suction head, S (p = ‐S), which were followed by the Brooks‐Corey,
Campbell, Fredlund, Kosugi, van Genuchten‐Mualem, and other empirical k(S) functions. Below, we
answer a question from the Introduction by selecting two Gardner's formulas, viz., exponential and algebraick(S) functions, which we call the G‐E and G‐A equations (in Gardner, 1958, equations (8) and(11), correspondingly).
Gardner and Fireman (1958) stated that “… the steady‐state rate of evaporation from a soil surface mulch
should be inversely proportional to the thickness of the mulch” and experimentally quantiﬁed this monotonic reduction in evaporation (see Figure 7 in their paper). However, in their experiments, D was not ﬁxed,
and the observed decrease in q due to adding more sand on the top of loam seems trivial. Similarly, a decreasein evaporation through homogeneous or layered soils when the water table in Figure 1 drops, i.e., D increases(see also Ripple et al., 1970, their Figure 9), is also trivial. Indeed, with an increase in the thickness of a layerthat conducts water, heat, electricity, and contaminants (according to the Darcy, Fourier, Ohm, Fick, andother resistivity laws), the ﬂux of a corresponding substance always decreases if a layer is thickened, provided
boundary conditions of a ﬂow tube are kept (see Goldshteĭn & Entov, 1994; Polubarinova‐Kochina, 1962).
The theoretical support of a monotonic decrease of q with B at a constant D can be inferred from Gardner(1958), who considered the upper mulch layer to be a purely vapor conducting entity.
Willis (1960) extended Gardner's work and conducted an approximate theoretical analysis of steady evaporation (using the G‐A equation) through a two‐layered soil identical to that in our Figure 1. From Willis' Figure
8, for example, for D = 200 cm, one can already determine the nonmonotonic dependence of evaporation
losses on the mulch thickness (our B), although the ﬁve Willis' curves at this D are cluttered. Figure 9 ofWillis (1960) presents the results of column experiments that also indicate the nonmonotonicity of q(B) ata ﬁxed D. However, Willis did not systematically solve the problem of detecting the value of B that delivers
an extremum of q. That would require a continuous variation of B (or, in Willis' vernacular, the position of
the water table). Willis' experiment also manifests that at ﬁxed B and D the ﬂux q depends on the texturalcomposition of his layered columns in a rather non‐trivial manner.
In the ﬁeld and lab experiments, agronomist and hydrologists found (see, e.g., Assouline et al., 2014; Ehlers& Van Der Ploeg, 1976; Jalota & Prihar, 1990; Klocke et al., 2009; Larson et al., 1983; Price et al., 1998; Prihar
et al., 1996; Schwartz et al., 2010; Unger & Cassel, 1991; Willis & Bond, 1971; Wuest & Schillinger, 2011) thatthe layering of soil (by mulching, harrowing, spudding, tillage, machine‐induced compaction, and by
natural crust‐seal formation) may both decrease and increase evaporation, which is apparently at odds with
Gardner (1958).
We will now prove that the simplest possible Gardner (1958) model based on the 1‐D, steady Richards' equation (ignoring, for example, vapor ﬂow, transiency, intermittency with inﬁltration and redistribution, and
KACIMOV ET AL.

3

Water Resources Research

10.1029/2018WR024025

moisture interception by plant roots), can explain the above‐cited seemingly recalcitrant and incongruous
results. Namely, we illustrate that q(B) can be a single‐extremum function for the ﬂow sketched in
Figure 1. The nonmonotonic function q(B) corroborates the positive, negative, or no‐impact of soil layering
on evaporation.
We utilize modern computer algebra to solve ordinary differential equations (ODEs) and nonlinear equations involving numerical integrations (special functions), which Bakr et al. (1979), Gardner (1958),
Gardner and Fireman (1958), Ripple et al. (1970), Salvucci (1993), Shi et al. (2014), Warrick (1988),
Warrick and Yeh (1990), and Willis (1960), among others, faced when analytically studying evaporation
through homogeneous and layered soils. In a sense, we have realized the program of Hillel et al.
(1975), who, after numerical simulations of two‐layered soils, stated the necessity to optimize heterogeneity of soils, in particular, mulching parameters (see also Al‐Maktoumi et al., 2014; Betti et al., 2016).
Along with the determination of the extrema of q(B) for a two‐layered composite in Figure 1 (with a
G‐E or G‐A constituting relation in both layers), we extend our analysis to continuously
heterogeneous soils.

2.1. G‐E Conductivity Model for Two‐Layered Soils
To be consistent with Philip's (1969, 1991) conventions, in this section, the direction of the vertical axis
OZ is positive downward, with the origin O at the soil surface. We assume that both soil layers in
Figure 1 are characterized by the dyads of the saturated hydraulic conductivity and the sorptive number(both are positive constants within each layer), (K1,α1) and (K2,α2), which are routinely determined in thelab (see, e.g., Wendroth & Wypler, 2008) or ﬁeld (e.g., by tension inﬁltrometry), as well as from pedotransfer functions or stochastic models (e.g., Lu & Zhang, 2004). The unsaturated conductivities obey
the relations
k 1 ðp1 Þ ¼ K 1 exp½α1 p1 ;

k2 ðp2 Þ ¼ K 2 exp½α2 p2 ;(1)where p1(z) and p2(z) are the pressure heads (negative) in the corresponding layers of Figure 1.
The Darcy law states thatv1Z ¼ −k1 ðp1 Þ

dp1
þ k1 ðp1 Þ;dZ

v2Z ¼ −k 2 ðp2 Þ

dp2
þ k 2 ðp2 Þ;dZ(2)where v1Z = v2Z = q are Philip's notations for the vertical Darcian velocities. Ripple et al. (1970, p.34) gave ahydrological balance rather than an agronomist's perspective of evaporation through stratiﬁed soils: “Often,
the only information sought is the dependence of the soil‐limited evaporation, upon the water table depth.”
We will also focus on q.
Along the interface M1M2, the continuities of pressure and ﬂux require that p1 = p2
and v1Z = v2Z = q = const.
Following Philip, we introduce the Kirchhoff potentials in the two layers:p1

ϕ1 ¼ ∫ k1 ðuÞdu ¼−∞

k 1 ðp1 Þ;
α1

p2

ϕ2 ¼ ∫ k 2 ðuÞdu ¼−∞

k2 ðp2 Þ:α2(3)From the conservation of mass, we have two linear ODEs:
d 2 ϕ1dϕ− α1 1 ¼ 0;
dZ
dZ 2
d 2 ϕ2
dϕ2
− α2
¼ 0;
dZ
dZ 2

B < Z < D;
(4)
0 < Z < B:

The Kirchhoff potential is not continuous across M1M2; from the continuities of pressure and ﬂux along thisline, we obtain two equations,

KACIMOV ET AL.

4

Water Resources Research

10.1029/2018WR024025

1
α1 ϕ1 ðBÞ1
α2 ϕ2 ðBÞdϕ ðBÞ
ln¼
ln; α2 ϕ2 ðBÞ− 2
α1
K1
α2
K2
dZ
dϕ1 ðBÞ
¼ α1 ϕ1 ðBÞ−
;
dZ

(5)

which serve for the determination of ϕ1(B) and ϕ2(B). We introduce
dimensionless values: Z* = Z/D, B* = B/D, a*1 ¼ a1 D, a*2 ¼ a2 D, K *2
¼ K 2 =K 1 , p*1 ¼ p1 =D, p*2 ¼ p2 =D, ϕ*1 ¼ ϕ1 =ðK 1 DÞ, ϕ*2 ¼ ϕ2 =ðK 1 DÞ, ϕ*2s ¼
ϕ2s =ðK 1 DÞ, and q* = q/K1, which are slightly different from commonlyused by Philip, Raats, Warrick, and others (see e.g., Blunt, 2017;
Warrick, 2003) because we vary the capillary properties of the layers
keeping D constant. We will drop “*” for the dimensionless quantities.
Figure 2. Gardner's exponential soils; an absolute value of the dimensionless evaporative ﬂux |vz| as a function of B for a coarser upper layer:
K2 = 100, α2 = 4, α1 = 0.4, and ϕ2s = 0.

First, we solve two BVPs for equation (4) with the corresponding boundary conditions ϕ1(1) = 1/α1, ϕ2(0) = ϕ2s, and p1(B) = p2(B):expðα1 ZÞ− expðα1 BÞ þ α1 ϕ1 ðBÞ½ expðα1 Þ− expðα1 Z Þ
; B < Z < 1;
α1 ½ expðα1 Þ− expðα1 BÞ½ expðα2 Z Þ−1ϕ2 ðBÞ þ ϕ2s ½ expðα2 BÞ− expðα2 Z Þ
; 0 < Z < B:
ϕ2 ¼expðα2 BÞ−1

ϕ1 ¼

(6)

From equation (6) we determine the ﬂux:
dϕ1 α1 ϕ1 ðBÞ expðα1 Þ− expðα1 BÞ¼; B < Z < 1;expðα1 Þ− expðα1 BÞdZdϕ
α2 ½ϕ2s expðα2 BÞ−ϕ2 ðBÞq ¼ α2 ϕ2 − 2 ¼; 0 < Z < B:expðα2 BÞ−1dZ

q ¼ α1 ϕ1 −(7)

From equations (5) and (7) we obtainϕ2 ðBÞ ¼ ϕ2s expðα2 BÞ þexpðα2 BÞ−1 α1 ϕ1 ðBÞ− exp½α1 ðB−1Þ:
exp½α1 ðB−1Þ−1α2(8)Substituting (8) into the right‐hand side of the ﬁrst equation (5) we obtain a nonlinear equation with respect
to ϕ1(B). This equation is solved by the FindRoot routine in Wolfram's Mathematica. We placed the foundroot back into the ﬁrst equation (7) and obtained the ﬂux v1Z as a function of (K,α), that is, of hydraulic‐capillary properties of the two soils, thickness B, and the potential at thesoil surface ϕ2s (a constant determined by the atmospheric conditions).
Figure 2 shows the modulus |vz| as a function of B for a mulch material
with K2 = 100, α2 = 4, α1 = 0.4 and an absolutely dry soil surface S1S2 that
corresponds to ϕ2s = 0. In prescribing this boundary condition at the soil
surface, we followed Philip (1991). He also assumed that the soil surface isabsolutely dry that gave the upper bound of evaporation losses. Physically,
evaporation through a very dry soil is prevalently in a vapor form, with aminor contribution of liquid motion through corner zones of pore channels, which Blunt (2017) modeled as triangles. Philip (1991) still used
Richards' equation in domains with a dry‐soil boundary condition.
The minimum of |vz|, qm = 0.47, is attained at Bm = 0.21. The minimum in
Figure 2 is quite pronounced. Similarly, optimal Bs were obtained forother K2, α2, α1 of the G‐E model.
Figure 3. Gardner's exponential soils; an absolute value of the dimensionless evaporative ﬂuxp|v
ﬃﬃﬃ z| as a function of B for a ﬁner upper layer: K 2 ¼ 0:5;
α2 ¼ 4; α1 ¼ 4 2; and ϕ2s = 0.

KACIMOV ET AL.

Figure 3 shows |vz| as a function of B of the topsoil, which is less permeable and has higher capillarity than the substratum: K 2 ¼ 0:5; α2 ¼ 4;
5

Water Resources Research

10.1029/2018WR024025

pﬃﬃﬃ
α1 ¼ 4 2; and ϕ2s ¼ 0 . For this composite, the magnitude of theﬂux attains a maximum, qM = 0.0095, at BM = 0.6. The minimum,qm = 0.0035, is attained at Bm = 0. Consequently, in this case, the ﬁne
top layer of any thickness (e.g., created by the above‐mentioned soil crusting) would evaporationally spoil a coarse homogeneous background soil.
This is contrary to what Ripple et al. (1970) found for a two‐layered(crusted) soil: “… a relatively thin less permeable layer may markedly
decrease evaporation rates” (see their Figure 9). In Ripple's Figure 9, thecurves for q at small B (in our notation) for two‐layered soils are indistinguishable from those of a homogeneous soil. That may be caused by an
insufﬁcient accuracy of Ripple's iterative solutions of the nonlinear equations. Indeed, half a century ago, before the advent of computer algebra, it
would not have been easy to detect the extrema in our Figures 2 and 3.

Figure 4. Gardner's algebraic soils; dimensionless q(B) for D = 20 cm. Layer32
1 is Yolo clay: n1 = 2, a1 = 400 cm /day, b1 = 400 cm , and K1 = 1 cm/day,
5
4and layer 2 is Pachapa sandy loam: n2 = 3, a2 = 3.2 · 10 cm /day,43
b2 = 2.6 · 10 cm , and K2 = 12.3 cm/day (curve 1). Layer 1 is Yolo clay, and8463
layer 2 is Willis' sand: n2 = 4, a2 = 1.7 · 10 cm /day, b2 = 2.6 · 10 cm , andK2 = 68 cm/day (curve 2). The soil surface is extremely dry.

We note that Kumar (1999) also found minima of the function |vz(B)| for acoarse upper stratum, as in our Figure 1. However, Kumar's computations
were cumbersomely indirect; he numerically (using a ﬁnite difference
method) solved a transient Richards' equation and considered a steady‐state limit with the same boundary conditions as in our analytic solution.
He obtained that limit for a few values of B. In our computer algebra code, the q(B) curves and extrema onthem are found in a trice.
2.2. G‐A Conductivity Model for a Two‐Layered SoilIn this section, we select algebraic functions rather than G‐E equation (1):k 1 ðS1 Þ ¼

a1
K1
¼
;
Sn1 1 þ b1 Sn1 1 =b1 þ 1

k2 ðS2 Þ ¼

a2
K2
¼
;
Sn2 2 þ b2 Sn2 2 =b2 þ 1(9)where, for consistency with Gardner‐Warrick's notations, we use the suction pressure heads S1,2 = −p1,2 anda vertical coordinate OZG, with an origin at the water table (Figure 1). This axis is oriented upward such thatq = vz is now co‐oriented with OZG. In equation (9), the triads (a1,2, b1,2, and n1,2) are composed of mathematically arbitrary positive constants that characterize the two soil layers in Figure 1. Equation (2), with (9)inserted, separate and are integrated as
S1

du
;
q=k
ð
1 uÞ þ 1
0

ZG ¼ ∫

0 < Z G < D−B; 0 < S1 < S2B ;
(10)

S2

du
ZG ¼ ∫;
S2B q=k 2 ðuÞ þ 1

D−B < Z G < D; S2B < S2 < SS :Obviously, S2b = S1b (the suction at the interface is continuous). This constant, as well as q, must be found.
We introduce dimensionless values Z *G ¼ Z G =D; B* ¼ B=D; S*1 ¼ S1 =D; S*2 ¼ S2 =D; b*1 ¼ b1 =D2 ; b*2 ¼ b2 =D2 ;q* ¼ q=K 1 ; and K *2 ¼ K 2 =K 1 ; and drop the “*”. Equation (10) is then integrated by the Integrate routine
in Wolfram's (1991) Mathematica as
S1
qSn1 1
F 1; 1=n1 ; 1 þ 1=n1 ; −; 0 < Z G < 1 − B; 0 < S1 < S1B ;1þq
b1 þ b1 q





K2
1
1
−qSn2 2
1
1
− qSn1B2
Z G ðS2 Þ ¼ 1−B þ
S2 F 1; ; 1 þ ;− S1B F 1; ; 1 þ ;;n2
n2 b2 ðK 2 þ qÞn2
n2 b2 ðK 2 þ qÞ
K2 þ q1−B < Z G < 1; S2B < S2 < SS ;Z G ðS1 Þ ¼(11)

where F is the hypergeometric function 2F1.
Since Gardner (1958), (We recall that Gardner himself generalized the earlier Wind and Remson‐Fox
empirical k(S) functions.) soil scientists have studied approximate solutions of equation (10) and exact solutions for various special cases of the G‐A function, for example, n = 1, 3/2, 2, 3, and 4, for which the integrals
KACIMOV ET AL.

6

Water Resources Research

10.1029/2018WR024025

in (10) are evaluated in arctan‐log functions. These quadratures were subject to several inversions ZG(S) → S(ZG) in expression (11); e.g., see Bakr
et al., 1979; Ripple et al., 1970; Salvucci, 1993; Shi et al., 2014; Warrick,1988). Using Mathematica, a systematic solution of nonlinear equations
involving special functions and a continuous variation of all involved
parameters (n, a, and b in (9) and Ss, B, and D in the BVP for an ODE)
is a standard procedure for any n. Consequently, we set S1 = S1B and
ZG = 1‐B and S2 = Ss and ZG = 1 in the corresponding ﬁrst and second
equations (11). We solved the obtained system of two nonlinear equations
with respect to the pair (S1B, q) using the Mathematica FindRoot routine(We recall that at n1 = n2 = 2 for the clay‐like layers in Figure 1, the suction head S1B along the interface is explicitly expressed through q, and thesystem (11) is simpliﬁed to one nonlinear equation with respect to q, anaFigure 5. Gardner's algebraic soils; layer 1 is Diablo loam: n1 = 2,logous to the case of the G‐E functions from section 2.1.).
3
2
a1 = 700 cm /day, b1 = 1450 cm , and K1 = 0.48 cm/day, and layer 2 is
2.2.1.
ExamplesWillis' sand. Dimensionless q(B) for D = 20 cm (curve 1) and D = 15 cm
We ﬁrst consider D = 20 cm and use the Gardner and Fireman (1958) data
(curve 2). The soil surface is extremely dry.
for soil properties; layers 1 and 2 in Figure 1 are Yolo clay and Pachapa
sandy loam, for which the corresponding dimensional parameters in equation (9) are n1 = 2 anda1 = 400 cm3/day, b1 = 400 cm2, and K1 = 1 cm/day, and n2 = 3, a2 = 3.2 · 105 cm4/day,
b2 = 2.6 · 104 cm3, and K2 = 12.3 cm/day.
In Figure 4, curve 1 shows the corresponding dimensionless evaporation losses q(B). Next, we used the same
D and Yolo clay of layer 1 for layer 2. Instead of the Pachapa soil, we used the river bed sand from Willis(1960), for which n2 = 4, a2 = 1.7 · 108 cm4/day, b2 = 2.6 · 106 cm3, and K2 = 68 cm/day. The losses for thistwo‐layered composite are shown as curve 2 in Figure 4.
It is evident from Figure 4 that evaporation is indeed reduced by a thin (approximately 1 cm) layer of a coarser material placed on the top of the Yolo clay. In addition, we have a mild minimum of q(B) for both two‐layered soils. This minimum can easily go unnoticed in experiments. At a given B and D, a coarser mulch(Willis' sand in Figure 4) may—depending on B—evaporate less (at B < Bi) or more (at B > Bi) than a ﬁner
mulch (Pachapa's sandy loam) with the same thickness B. The two curves in Figure 4 intersect at
B = Bi = 0.05. We recall that Yuan et al. (2009) also varied the mulch properties and found in their drying
experiments that the ﬁnest mulch evaporated less than coarser ones.
In the second example, layer 1 is the Diablo loam from Willis (1960), for which n1 = 2, a1 = 700 cm3/day,
b1 = 1450 cm2, and K1 = 0.48 cm/day, and layer 2 is Willis' riverbed sand. We varied D from 20 to 15 cm.
The corresponding curves of q(B) are shown as curves 1 and 2 in Figure 5. Again, as in Figure 4, both curves
in Figure 5 have minima (we computed the minima using the FindMinimum Mathematica routine) at
approximately the same dimensionless Bm = 0.02. The minimum qm
increased signiﬁcantly (from 2.4 to 3.4) with an decrease in D. We recall
Hadas' (1997) conclusion: “Tillage‐produced ‘soil mulch,’ a structured soil
layer, reduces soil water evaporation in a similar manner to plant residues
mulch, provided its aggregate sizes range between 8.3 and 5.0 mm diameter and its depth is between 5 and 10 cm.”
In the third example, we selected two claylike layers, the Diablo and Yolo
soils, as layers 1 and 2. In Figure 6, curves 1 and 2 show q(B) at D = 20 cm
and 15 cm, respectively. In comparison with Figure 5, much more pronounced minima were attained at Bm = 0.54 and 0.52.

Figure 6. Gardners' algebraic soils; layer 1 is Diablo loam; layer 2 is Yolo
clay. Dimensionless q(B) at D = 20 cm (curve 1) and D = 15 cm (curve 2).
The soil surface is extremely dry.

KACIMOV ET AL.

In the fourth example, a homogeneous Diablo clay with D0 = 48 cm had
q0 = 0.84. The subscript here indicates the soil which is either thickened
or heterogenized in the following manner. We compared two types of
melioration of D0: thin mulching with the addition of Willis' sand of thickness B = 2 cm on the top of the clay (making D = 50 cm) and hallowing,
which converted a layer of thickness B = 10 cm into a pseudo‐Yolo clay

7

Water Resources Research

10.1029/2018WR024025(retaining D = 48 cm). From equation ((11)) we get for the former q1 = 0.7 and for the latter q2 = 0.74, thesecomposites reduce evaporative losses by approximately 17 and 12%, respectively, compared with ahomogeneous soil.
2.3. Average Suction in the Composite
The distributions S(Z) in Figure 1 can be immediately obtained by inverting equation (11) using
Mathematica routines. Here we introduce the vertically averaged values of the capillary pressure head in
the strata of Figure 1:
1−B

CP1 ¼

1

1
1
∫ S1 ðZ G ÞdZ G ; CP2 ¼ ∫ S2 ðZ G ÞdZ G ; CP ¼ CP1 þ CP2
1−B 0
B 1−B

(12)

We computed the integrals in (12). For the example of the Diablo‐Yolo clay (D = 20 cm, curve 1 in Figure 6),
CPm = 34.8 for the minimal‐q composite.
2.4. Continuous Variation in Soil Properties With Depth
Do the minima for q occur in soils with continuous variations in hydraulic properties with Z? In this section,
we return to the G‐E model and Philip's vertical coordinate convention in Figure 1. Now we assume thatk(Z) = K0(Z)exp[α(Ζ) p(Z)], where K0(Z) and α(Z) are two given functions. For simplicity, we assume thatthese two functions are correlated according to an empirical function from Communar and Friedman(2015):αðZ G Þ ¼ Ca

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
K 0 ðZ Þ;(13)

where K0(Z) is in cm/hr, α(Z) is in 1/cm, and Ca is in h1/2/cm3/2. There are numerous pedotransfer functions
and other empirical relations between K0(Z) and α(Z), which can be used instead of equation (13). For arbitrary K0(Z) and α(Z) in equation (13), the ﬁrst‐order ODEs (2) cannot be integrated in the form of Z(p) as in
Gardner (1958), Willis (1960), and Warrick (1988). Instead, we must numerically solve a BVP for an ODE(e.g., see Salvucci, 1993):
hi
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
d pðZ ÞvZ
þexp −Ca K 0 ðZ ÞpðZ Þ −1 ¼ 0;
dZ
K 0 ðZ Þ

0 < Z < D;pð0Þ ¼ ps ;pðDÞ ¼ 0:(14)This is easily done using the Mathematica NDSolve routine, but integration in quadratures is not possible.
There is still some room for explicit integration of (14) if we expand the exponent there into a Taylor series
and retain the ﬁrst two terms (also see Barontini et al., 2007, for a similar linearization of Richards' equation
with the G‐E k(S) in transient evaporation models). This linearization reduces (14) to a linear ODE:dpðZ Þ
C a vZvZ− pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pðZÞ ¼ 1−;
dZ
K 0 ðZ Þ
K 0 ðZ Þ

0<Z<D;

pð0Þ ¼ ps ;pðDÞ ¼ 0:(15)

A general solution of BVP (15) is
!
! !

Z
t
dt
vZ
dτ
c þ ∫ 1−
exp − Ca vZ ∫ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ dt :pðZÞ ¼ exp Ca vZ ∫ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
K 0 ðt Þ
K 0 ðt Þ
K 0 ðτ Þ
0
0
0
Z

(16)

The constant of integration c = ps here is due to the boundary condition at Z = 0. Using the boundary condition at Z = D in (15), we obtain the following nonlinear equation with respect to vZ:
!

t
vZ
dτ
exp −Ca vZ ∫ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ dt ¼ −pS :∫ 1−
K 0 ðt Þ
K 0 ðτ Þ
0
0
D

(17)

For a given particular K0(z), we can reduce a nonlinear integral equation (17) to an algebraic one. For
example, if
KACIMOV ET AL.

8

Water Resources Research

10.1029/2018WR024025

Figure 7. Dimensional evaporation losses as a function of β1 for Gardner's exponential unsaturated conductivity function.
The saturated hydraulic conductivity K0 varies with Z according to equation (18); the sorptive number is correlated to K0
according to (13). Curves 1–3 correspond to three values of suction, ps = −1,000, −10,000, and −100,000 cm (curves 1–3,
respectively).

K 0 ðZÞ ¼ ðβ1 Z þ β2 Þ2

ðβ1 Z þ β2 ≠0;∀Z∈½0; DÞ;(18)

we getβ2vz =β2þ− pS
β1 −Ca vz β1 þ Ca vz



β1 D þ β2β2

Ca vz =β1¼

β1 D þ β2 vz =ðβ1 D þ β2 Þþ:
β1 þ Ca vz
β1 −Ca vz(19)

We used the FindRoot routine in Mathematica to solve equation (19) with respect to vz.
Figure 7 shows the magnitude of the dimensional evaporative ﬂux in cm/hr as a function of β1 at D = 20 cm,
β 2 = 1 (cm/hr)1/2, and Ca = 0.04 h1/2/cm3/2 for three values of suction, ps = −1,000, −10,000, and−100,000 cm (curves 1–3, respectively).
If we retain the third (quadratic) term in the Taylor series expansion of the exponential function in (14), thenthe BVP for the nonlinear ODE is!
dpðZ Þ
0:04
1
2− vZ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pðZ Þ þ 0:0008pðZ Þ þ−1 ¼ 0; 0 < Z < D; pð0Þ ¼ ps ; pðDÞ ¼ 0:
dZ
K 0 ðZ Þ
K 0 ðZ Þ(20)

We solved (20) using the NDSolve Mathematica routine in the same manner as the nonlinearized BVP (14).
It is noteworthy that NDSolve and FindRoot use numerical evaluations (see Wolfram, 1991) for integration
and root ﬁnding, which fail if the initial guesses are poor (we skip over examples of such types of crashes forthe routines in our computations).

3. HYDRUS‐1D Simulations
In this section, we use HYDRUS‐1D (Šimůnek et al., 2016), a public domain software that solves the governing water ﬂow (Richards) equation using the ﬁnite element method. In all transient ﬂow simulations, time, t,
is in days and the vertical size of two‐layered columns is in centimeter. We use several default HYDRUS
options: no vapor ﬂow, the VGM capillary pressure and phase permeability functions, and the HYDRUS
Soil Catalogue.
The column in Figure 1 is composed of clay (stratum 1) and sand (stratum 2).
As in the analytical solution of section 2 we ﬁxed the total thickness D = 20 cm (HYDRUS‐1D works withdimensional quantities). An ﬁnite element mesh discretizes it into 500 elements in the vertical, ZH, direction.
The bottom and top boundary conditions are constant pressure heads p = 0 and ps = −10,000 cm,
KACIMOV ET AL.

9

Water Resources Research

10.1029/2018WR024025

respectively. Initially, the pressure head is distributed hydrostaticallybetween these two values. Simulations are run for t = 20 days during
which steady‐state evaporation from the water table is established.
Because of the very steep retention curve of the default sand, when smallchanges in the water content correspond with large changes in the pressure head (for sands the pressure head can ﬂuctuate from −100 to−10,000 while still being within the default HYDRUS‐1D water content
tolerance limit of 0.001), the water content tolerance was reduced to 1.e−5. We also disabled the internal interpolation tables to get more precise
results. The steady‐state and transient surface ﬂuxes at the top of the soilcolumn computed using HYDRUS‐1D are plotted in Figure 8 as a functionof the thickness B (top) and for three values of B = 1, 6, and 18 cm (bottom), respectively. In notations of Figure 2, the results are |vz| = 0.114,
0.1003, and 0.320 cm/day. Therefore, at the mulch thickness B = 6 cm aminimum (actually, a pseudo‐minimum for these three discrete valuesof B) is attained. We surmise that the minimum is unique, global, andclose to the pseudo‐minimum. Thus, the analytically found extremum in
Figure 2 for Gardner's bare soils and ﬁnite difference method extrema
found by Kumar (1999) are qualitatively corroborated by FEM simulations for VGM two‐component composites.

Figure 8. (top) The dimensional steady‐state evaporation rate (ﬂux) as afunction of the thickness B and (bottom) the transient evaporation rates
for thicknesses B = 1, 6, and 18 cm. The evaporation rate is calculated using
HYDRUS‐1D from an isobar p = −10,000 cm, S1S2 (Figure 1) for a two‐layered soil of the total thickness D = 20 cm, and varying thickness B of the
top mulch made of a coarse layer (sand) with a subjacent clay layer andwater supply from an isobar p = 0 cm, W1 W2.

The effect of uptake by plants' roots on water ﬂow and the water content
distribution along the unsaturated soil proﬁle often prevails over a purely
capillary ﬂow. In order to quantify this juxtaposition of evaporation, transpiration, and variation of the thickness of our composite we used the root
water uptake (RWU) model in HYDRUS‐1D, assuming a grass with a uniform root distribution throughout the proﬁle (0 < ZG < 20 cm in Figure 1),
Feddes' stress response function, and potential transpiration of 0.5 cm/
day. The results are shown in Figure 9 (for discretely varied B = 2, 4, …
16, 18 cm) and evidence the following. The total evapotranspirative ﬂux
still has a global minimum with B as for pure evaporation.
For small thicknesses of sand, there is no evaporation and only RWU
declining with the sand thickness. Once the proﬁle is mostly sand (14–18 cm), there is a further reduction in RWU, while the proﬁle starts evaporating. It is noteworthy that the RWU curve in Figure 9 exhibits a startling, hardly expected behavior: a global end‐maximum at B = 2 cm, a
global end‐minimum at B = 18 cm, and two internal extrema, viz., local
maximum and minimum.

4. Concluding Remarks

Figure 9. Dimensional steady‐state evaporation ﬂuxes as functions of B for
the same composite as in Figure 8 but root water uptake by Feddes' grass
added in HYDRUS simulations: total evapotranspiration, evaporation, andtranspiration.

KACIMOV ET AL.

We examined a 1‐D ascending ﬂow through an unsaturated two‐layeredsoil in Figure 1. Analytically, we did not study soils' transient multistage
drying but rather analyzed a steady‐state evaporation from a stationary
isobar (water table), as Gardner, Willis, and Warrick did. For each stratumin Figure 1, we used the Gardner (1958) exponential and algebraic functions for the unsaturated conductivity‐suction pressure head k(S), which
are characterized by dyads and triads of empirical soils' constants. We
employed modern computer algebra routines (solutions of systems ofnonlinear equations, integration, ﬁnding minima, and evaluation of
hypergeometric functions) to determine the constant evaporation ﬂux q
as a function of the thickness B of the upper stratum of the unsaturated
zone. We interpreted this texturally contrasted stratum asmulching/tillage and illustrated that q(B) might have nontrivial minima
10

Water Resources Research

10.1029/2018WR024025(qm, Bm) if the total composite thickness D is ﬁxed, a fact found numerically by Kumar (1999). Numericalsimulations in HYDRUS‐1D for the van Genuchten soil properties qualitatively corroborate the analyticalresults, viz., for a ﬁxed total thickness of a two‐layered soil a certain thickness of the upper coarse layer minimizes an evaporative ﬂux from bottom isobar to a top one. A HYDRUS‐simulated example with Feddes'
grass uptaking moisture from the soil composite illustrates nontrivial extrema of the global ﬂux from a shallow water table and of the transpiration due to water uptake by plants' roots.
The reduction in q by shallow (1–2 inches of B in Figure 1) spudding/harrowing (without traditional tillage)
was tested by Ovsinsky (1899) and later replicated and developed for various climatic and agronomic conditions (Allen, 1981; Faulkner, 1943; Friedrich et al., 2014; Kassam et al., 2018; Maltsev, 1954; Novakovska
et al., 2018; Verhulst et al., 2010; Yadav et al., 2018). Along with many other advantages (including a better
aeration of the root zone, circadian condensation of the dew in the B‐layer of Figure 1, better nitriﬁcationand ensued reduction or the elimination of the application of agrochemicals, and more proliﬁc activity of soil
microbes), this no‐till technology in semiarid climates preserves more moisture in the root zone. No‐till orshallow soil treatment technology is, of course, not a panacea in agronomy (constraints and drawbacks are
discussed, for example, by Garmashov, 2018; Pykhtin, 2017; and Turusov et al., 2018). Indeed, even if solely
evaporation losses are a mathematical criterion, the minima in Figures 2 and 8 are rather mild, that is, even
optimal mulching only moderately reduces the ascending ﬂux.

Abbreviations
BVP
G‐E
G‐A
ODE
PDE
VGM

Acknowledgments
This work was funded by a grant fromthe Sultan Qaboos Higher Center for
Culture and Science – Diwan of Royal
Court and the Research Council of
Oman (TRC) [RC/AGR/SWAE/17/01]
and by the subsidy allocated to Kazan
Federal University for the state assignment in the sphere of scientiﬁc activities, project1.12878.2018/12.1. Helpfulcomments by three anonymous referees
are appreciated. The paper is theoretical, and no data are used.

KACIMOV ET AL.

= boundary value problem= Gardner's exponential unsaturated conductivity relation k(S)= Gardner's algebraic function for unsaturated conductivity k(S)= ordinary differential equation= partial differential equation= van Genuchten‐Mualem

References
Allen, H. P. (1981). Direct drilling and reduced cultivation. Ipswich: Farming Press.
Al‐Maktoumi, A., Al‐Ismaily, S., Kacimov, A., Al‐Busaidi, H., Al‐Saqri, S., & Al‐Haddabi, M. (2014). Soil substrate as a cascade of capillarybarriers for conserving water in a desert environment: Lessons learned from arid nature. Journal of Arid Land, 6(6), 690–703. https://doi.
org/10.1007/s40333‐014‐0068‐7
Al‐Shukaili, A., (2018). The effect of stratiﬁcation on sand dunes hydrology in arid zone environments: Case of Oman. MSc dissertation.
Sultan Qaboos University, Oman.
Assouline, S., Narkis, K., Gherabli, R., Lefort, P., & Prat, M. (2014). Analysis of the impact of surface layer properties on evaporation fromporous systems using column experiments and modiﬁed deﬁnition of characteristic length. Water Resources Research, 50, 3933–3955.
https://doi.org/10.1002/2013WR014489
Bakr, H. M. A., Fathi, A. M., & Gheibeh, A. R. S. (1979). Evaluation of the parameters involved in Gardner's relation between unsaturated
hydraulic conductivity and soil water matric suction. Agricultural Water Management, 2(1), 25–35.
Barontini, S., Ranzi, R., & Bacchi, B. (2007). Water dynamics in a gradually nonhomogeneous soil described by the linearized Richards
equation. Water Resources Research, 43, W08411. https://doi.org/10.1029/2006WR005126
Betti, G., Grant, C. D., Murray, R. S., & Churchman, G. J. (2016). Size of subsoil clods affects soil‐water availability in sand–clay mixtures.
Soil Research, 54(3), 276–290. https://doi.org/10.1071/SR15115
Blunt, M. J. (2017). Multiphase ﬂow in permeable media: A pore‐scale perspective. Cambridge: Cambridge University Press. https://doi.org/
10.1017/9781316145098
Boast, C. W., & Simmons, F. W. (2005). Evaporation of water from bare soil. In D. Hillel (Ed.), Encyclopedia of Soils in the Environment, (pp.
494–502). Oxford: Elsevier. https://doi.org/10.1016/B0‐12‐348530‐4/00358‐1
Bodner, G., Nakhforoosh, A., & Kaul, H. P. (2015). Management of crop water under drought: A review. Agronomy for Sustainable
Development, 35(2), 401–442.
Communar, G., & Friedman, S. P. (2015). Cylindrically conﬁned ﬂows of water from a point source and to a sink, above a water table.
Vadose Zone Journal, 14(8). https://doi.org/10.2136/vzj2014.12.0183
Connolly, R. D. (1998). Modelling effects of soil structure on the water balance of soil–crop systems: A review. Soil and Tillage Research,
48(1‐2), 1–19. https://doi.org/10.1016/S0167‐1987(98)00128‐7
De Luca, D. L., & Cepeda, J. M. (2016). Procedure to obtain analytical solutions of one‐dimensional Richards' equation for inﬁltration in
two‐layered soils. Journal of Hydrologic Engineering, 21(7), 04016018.
Ehlers, W., & Van Der Ploeg, R. R. (1976). Evaporation, drainage and unsaturated hydraulic conductivity of tilled and untilled fallow soil.
Zeitschrift für Pﬂanzenernährung und Bodenkunde, 139(3), 373–386. https://doi.org/10.1002/jpln.19761390312
Faulkner, E. H. (1943). Plowman's Folly. Norman: Univ. of Oklahoma Press.

11

Water Resources Research

10.1029/2018WR024025

Felde, V. J. M. N. L., Peth, S., Uteau‐Puschmann, D., Drahorad, S., & Felix‐Henningsen, P. (2014). Soil microstructure as an under‐exploredfeature of biological soil crust hydrological properties: Case study from the NW Negev Desert. Biodiversity and Conservation, 23(7),1687–1708. https://doi.org/10.1007/s10531‐014‐0693‐7Friedrich, T., Kassam, A., Corsi, S., Jat, R. A., Sahrawat, K. L., & Kassam, A. H. (2014). Conservation agriculture in Europe. In R. A. Jat, K.
L. Sahrawat, & A. H. Kassam (Eds.), Conservation agriculture: Global prospects and challenges (pp. 127–170). Boston: CABI Press.
Fuchs, M., & Hadas, A. (2011). Mulch resistance to water vapor transport. Agricultural Water Management, 98(6), 990–998. https://doi.org/
10.1016/j.agwat.2011.01.008
Gardner, W. R. (1958). Some steady solutions of the unsaturated moisture ﬂow equation with application to evaporation from a water table.
Soil Science, 85(4), 228–232. https://doi.org/10.1097/00010694‐195804000‐00006
Gardner, W. R., & Fireman, M. (1958). Laboratory studies of evaporation from soil columns in the presence of a water table. Soil Science,85(5), 244–249. https://doi.org/10.1097/00010694‐195805000‐00002
Garmashov, V. M. (2018). Principles and methods of optimization of basic soil treatment and rehabilitation of fertility of common chernozem by crop rotation in the Central Chernozem Zone. DSc. Dissertation, Kamennaya Step, Russia (in Russian).
Geng, X., & Boufadel, M. C. (2015). Numerical modeling of water ﬂow and salt transport in bare saline soil subjected to evaporation.
Journal of Hydrology, 524, 427–438. https://doi.org/10.1016/j.jhydrol.2015.02.046
Goldshteĭn, R. V., & Entov, V. M. (1994). Qualitative methods in continuum mechanics. Pitman monographs and surveys in pure and applied ,athematics, (Vol. 72). Harlow: Longman.
Hadas, A. (1997). Soil tilth—The desired soil structural state obtained through proper soil fragmentation and reorientation processes. Soiland Tillage Research, 43(1‐2), 7–40. https://doi.org/10.1016/S0167‐1987(97)00033‐0
Hillel, D., & Talpaz, H. (1977). Simulation of soil water dynamics in layered soils. Soil Science, 123(1), 54–62. https://doi.org/10.1097/
00010694‐197701000‐00007
Hillel, D. I., Van Bavel, C. H. M., & Talpaz, H. (1975). Dynamic simulation of water storage in fallow soil as affected by mulch ofhydrophobic aggregates. Soil Science Society of America Journal, 39(5), 826–833. https://doi.org/10.2136/
sssaj1975.03615995003900050015x
Jalota, S. K., & Prihar, S. S. (1990). Bare‐soil evaporation in relation to tillage. In B. A. Stewart (Ed.), Advances in Soil Science, (Vol. 12, pp.
187–216). New York: Springer.
Kacimov, A., & Obnosov, Y. V. (2018). Hydraulically optimal porous liner around a porous lens: Strack's problem revisited. Journal of
Hydraulic Engineering ISH, 1–11. https://doi.org/10.1080/09715010.2018.1516576
Kassam, A., Friedrich, T. and Derpsch, R., 2018. Global spread of conservation agriculture. International J. of Environmental Studies,
https://doi.org/10.1080/00207233.2018.1494927
Khan, A. R. (1988). Assessing the evaporation zone in the bare soil from the soil water ﬂux and soil heat ﬂux measurements. Journal ofAgronomy and Crop Science, 161(4), 234–237. https://doi.org/10.1111/j.1439‐037X.1988.tb00662.x
Klocke, N. L., Currie, R. S., & Aiken, R. M. (2009). Soil water evaporation and crop residues. Transactions of the ASABE, 52(1), 103–110.
https://doi.org/10.13031/2013.25951
Kornev, V. G. (1935). Subsurface irrigation. Moscow‐Leningrad: Selhozgiz. (in Russian)
Kumar, C. P. (1999). Evaporation from shallow water table through layered soil proﬁles. ISH Journal of Hydraulic Engineering, 5(2), 65–75.
https://doi.org/10.1080/09715010.1999.10514654
Larson, W. E., Swan, J. B., & Shaffer, M. J. (1983). Soil management for semiarid regions. Agricultural Water Management, 7(1–3), 89–114.
https://doi.org/10.1016/0378‐3774(83)90077‐X
Li, H., Wu, F., Zhan, H., Qiu, F., & Wang, W. (2016). The effect of precipitation pulses on evaporation of deeply buried phreatic water in
extra‐arid areas. Vadose Zone Journal, 15(5). https://doi.org/10.2136/vzj2015.09.0127
Lu, Z., & Zhang, D. (2004). Analytical solutions to steady state unsaturated ﬂow in layered, randomly heterogeneous soils via Kirchhofftransformation. Advances in Water Resources, 27(8), 775–784. https://doi.org/10.1016/j.advwatres.2004.05.007
Malik, R. S., Anlauf, R., & Richter, J. (1992). A simple model for predicting evaporation from bare soils. Zeitschrift für Pﬂanzenernährung
und Bodenkunde, 155(4), 293–299. https://doi.org/10.1002/jpln.19921550409
Maltsev, T. (1954). New system of soil cultivation and sowing. Moscow: Znanie. (in Russian)Mehari, A., Van Steenbergen, F., & Schultz, B. (2011). Modernization of spate irrigated agriculture: A new approach. Irrigation andDrainage, 60(2), 163–173. https://doi.org/10.1002/ird.565
Minhas, P. S., Khosla, B. K., & Prihar, S. S. (1986). Evaporation and redistribution of salts in a silt loam soil as affected by tillage induced soil
mulch. Soil and Tillage Research, 7(4), 301–313. https://doi.org/10.1016/0167‐1987(86)90018‐8
Mualem, Y., & Assouline, S. (1989). Modeling soil seal as a nonuniform layer. Water Resources Research, 25(10), 2101–2108. https://doi.org/
10.1029/WR025i010p02101
Novák, V. (2012). Evapotranspiration in the Soil‐Plant‐Atmosphere System. Dordrecht: Springer. https://doi.org/10.1007/978‐94‐007‐3840‐9
Novakovska, I., Bulgakov, V., Rucins, A. & Dukulis, I., (2018). Analysis of soil tillage by ploughs and optimization of their aggregation.
Engineering for Rural Development, 335–341 DOI: https://doi.org/10.22616/ERDev2018.17.N401
Ovsinsky, I. E. (1899). A new system of agronomy. Kiev. (in Russian)Philip, J. R. (1969). Theory of inﬁltration. Advances in Hydroscience, 5, 215–296. https://doi.org/10.1016/B978‐1‐4831‐9936‐8.50010‐6
Philip, J. R. (1991). Upper bounds on evaporation losses from buried sources. Soil Science Society of America Journal, 55(6), 1516–1520.
https://doi.org/10.2136/sssaj1991.03615995005500060002x
Polubarinova‐Kochina, P. Y. (1962). Theory of ground water movement, (2nd edition of the book ed.). Princeton: Princeton University Press.
(in Russian) was published in 1977, Nauka, Moscow
Price, J., Rochefort, L., & Quinty, F. (1998). Energy and moisture considerations on cutover peatlands: Surface microtopography, mulch
cover and Sphagnum regeneration. Ecological Engineering, 10(4), 293–312. https://doi.org/10.1016/S0925‐8574(98)00046‐9
Prihar, S. S., Jalota, S. K., & Steiner, J. L. (1996). Residue management for reducing evaporation in relation to soil type and evaporativity.
Soil Use and Management, 12(3), 150–157. https://doi.org/10.1111/j.1475‐2743.1996.tb00536.x
Pykhtin I.G., (2017). Tillage: Reality and myths. Zemledelie. No 1. 33–36 (in Russian) https://cyberleninka.ru/article/n/obrabotka‐pochvy‐
deystvitelnost‐i‐mify
Ripple, C. D., Rubin, J., & Van Hylckama, T. E. A. (1970). Estimating steady‐state evaporation rates from bare soils under conditions of highwater table. Menlo Park, CA: USGS, Water Resources Division.
Sadeghi, M., Tuller, M., Gohardoust, M. R., & Jones, S. B. (2014). Column‐scale unsaturated hydraulic conductivity estimates in coarse‐
textured homogeneous and layered soils derived under steady‐state evaporation from a water table. Journal of Hydrology, 519, 1238–1248.
https://doi.org/10.1016/j.jhydrol.2014.09.004

KACIMOV ET AL.

12

Water Resources Research

10.1029/2018WR024025

Salvucci, G. D. (1993). An approximate solution for steady vertical ﬂux of moisture through an unsaturated homogeneous soil. Water
Resources Research, 29(11), 3749–3753. https://doi.org/10.1029/93WR02068
Schwartz, R. C., Baumhardt, R. L., & Evett, S. R. (2010). Tillage effects on soil water redistribution and bare soil evaporation throughout aseason. Soil and Tillage Research, 110(2), 221–229. https://doi.org/10.1016/j.still.2010.07.015
Shi, W. J., Shen, B., & Wang, Q. J. (2014). Effect of ignoring parameter b on the maximum upward ﬂux calculated from the Gardner rationalmodel. Canadian Journal of Soil Science, 94(1), 97–103. https://doi.org/10.4141/cjss2013‐029
Šimůnek, J., van Genuchten, M. T., & Šejna, M. (2016). Recent developments and applications of the HYDRUS computer software
packages. Vadose Zone Journal, 15(7), 25. https://doi.org/10.2136/vzj2016.04.0033
Soylu, M. E., Istanbulluoglu, E., Lenters, J. D., & Wang, T. (2011). Quantifying the impact of groundwater depth on evapotranspiration in asemi‐arid grassland region. Hydrology and Earth System Sciences, 15(3), 787–806. https://doi.org/10.5194/hess‐15‐787‐2011
Stewart, J. M., & Broadbridge, P. (1999). Calculation of humidity during evaporation from soil. Advances in Water Resources, 22(5), 495–505.
https://doi.org/10.1016/S0309‐1708(98)00024‐4
Strudley, M. W., Green, T. R., & Ascough, J. C. (2008). Tillage effects on soil hydraulic properties in space and time: State of the science. Soiland Tillage Research, 99(1), 4–48. https://doi.org/10.1016/j.still.2008.01.007
Turusov, V.I., Garmashov, V.M., Kornilov, I.M. & Nuzhnaya, N.A., (2018). Agro‐bioenergy efﬁciency of various pea soil basic management
practices under conditions of the south‐east of the Central Chernozem Zone. Proceedings Michurin's State Agrarian University, 6‐10.
http://www.mgau.ru/sciense/journal/PDF_ﬁles/vestnik_2_2018.pdf#page=10
Unger, P. W., & Cassel, D. (1991). Tillage implement disturbance effects on soil properties related to soil and water conservation: A literature review. Soil and Tillage Research, 19(4), 363–382. https://doi.org/10.1016/0167‐1987(91)90113‐C
van Genuchten, M. Th. (1980). A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society ofAmerica Journal, 44, 892–898.
Verburg, K., Bond, W. J., & Hunt, J. R. (2012). Fallow management in dryland agriculture: Explaining soil water accumulation using apulse paradigm. Field Crops Research, 130, 68–79. https://doi.org/10.1016/j.fcr.2012.02.016
Verhulst, N., Govaerts, B., Verachtert, E., Castellanos‐Navarrete, A., Mezzalama, M., Wall, P., et al. (2010). Conservation agriculture,
improving soil quality for sustainable production systems. In Advances in soil science: Food security and soil quality, (pp. 137–208). Boca
Raton: CRC press.
Warrick, A. W. (1988). Additional solutions for steady‐state evaporation from a shallow water table. Soil Science, 146(2), 63–66. https://doi.
org/10.1097/00010694‐198808000‐00001
Warrick, A. W. (2003). Soil water dynamics. New York: Oxford University Press.
Warrick, A. W., & Yeh, T. C. J. (1990). One‐dimensional, steady vertical ﬂow in a layered soil proﬁle. Advances in Water Resources, 13(4),207–210. https://doi.org/10.1016/0309‐1708(90)90042‐3
Wendroth, O., & Wypler, N. (2008). Unsaturated hydraulic properties: Laboratory evaporation. In M. R. Carter & E. G. Gregorich (Eds.),
Soil sampling and methods of analysis (Chap. 81, pp. 1089–1127). Boca Raton, FL: CRC Press.
Willis, W. O. (1960). Evaporation from layered soils in the presence of a water table. Soil Science Society of America Journal, 24(4), 239–242.
https://doi.org/10.2136/sssaj1960.03615995002400040009x
Willis, W. O., & Bond, J. J. (1971). Soil water evaporation: Reduction by simulated tillage. Soil Science Society of America Journal, 35(4),526–529. https://doi.org/10.2136/sssaj1971.03615995003500040016x
Wolfram, S. (1991). Mathematica. A system for doing mathematics by computer. Redwood City: Addison‐Wesley.
Wuest, S. B., & Schillinger, W. F. (2011). Evaporation from high residue no‐till versus tilled fallow in a dry summer climate. Soil Science
Society of America Journal, 75(4), 1513–1519. https://doi.org/10.2136/sssaj2010.0368
Xiao, B., Zhao, Y. G., & Shao, M. A. (2010). Characteristics and numeric simulation of soil evaporation in biological soil crusts. Journal of
Arid Environments, 74(1), 121–130. https://doi.org/10.1016/j.jaridenv.2009.06.013
Yadav, G. S., Saha, P., Babu, S., Das, A., Layek, J., & Debnath, C. (2018). Effect of no‐till and raised‐bed planting on soil moisture conservation and productivity of summer maize (zea mays) in Eastern Himalayas. Agricultural Research, 7(3), 300–310. https://doi.org/
10.1007/s40003‐018‐0308‐8
Yuan, C., Lei, T., Mao, L., Liu, H., & Wu, Y. (2009). Soil surface evaporation processes under mulches of different sized gravel. Catena,78(2), 117–121. https://doi.org/10.1016/j.catena.2009.03.002
Zarei, G., Homaee, M., Liaghat, A. M., & Hoorfar, A. H. (2010). A model for soil surface evaporation based on Campbell's retention curve.
Journal of Hydrology, 380(3‐4), 356–361. https://doi.org/10.1016/j.jhydrol.2009.11.010
Zhou, Y., & Rowe, R. K. (2005). Modeling of clay liner desiccation. International Journal of Geomechanics, 5(1), 1–9. https://doi.org/
10.1061/(ASCE)1532‐3641(2005)5:1(1)KACIMOV ET AL.

13