J. Hydrol. Hydromech., 66, 2, 170–180
DOI: 10.2478/johh-2018-0002

Inverse estimation of soil hydraulic properties and water repellency following
artificially induced drought stress
Vilim Filipović1*, Thomas Weninger2, Lana Filipović1, Andreas Schwen3, Keith L. Bristow4,
Sophie Zechmeister-Boltenstern5, Sonja Leitner5
1

University of Zagreb, Faculty of Agriculture, Department of Soil Amelioration, Svetošimunska 25, 10000 Zagreb, Croatia.
University of Natural Resources and Life Sciences Vienna (BOKU), Institute of Hydraulics and Rural Water Management, Muthgasse 18,1190 Vienna, Austria.
3
Austrian Agency for Health and Food Safety (AGES), Institute for Plant Protection Products, Spargelfeldstraße 191, 1220 Vienna, Austria.
4
CSIRO Agriculture & Food, PMB Aitkenvale, Townsville, QLD 4814, Australia.
5
University of Natural Resources and Life Sciences Vienna (BOKU), Institute of Soil Research, Peter-Jordan-Straße 82, 1190 Vienna,
Austria.
*
Corresponding author. Tel.: 00385 1 2393711. E-mail: vfilipovic@agr.hr
2

Abstract: Global climate change is projected to continue and result in prolonged and more intense droughts, which canincrease soil water repellency (SWR). To be able to estimate the consequences of SWR on vadose zone hydrology, it isimportant to determine soil hydraulic properties (SHP). Sequential modeling using HYDRUS (2D/3D) was performed on
an experimental field site with artificially imposed drought scenarios (moderately M and severely S stressed) and a control plot. First, inverse modeling was performed for SHP estimation based on water and ethanol infiltration experimental
data, followed by model validation on one selected irrigation event. Finally, hillslope modeling was performed to assesswater balance for 2014. Results suggest that prolonged dry periods can increase soil water repellency. Inverse modelingwas successfully performed for infiltrating liquids, water and ethanol, with R2 and model efficiency (E) values both >0.9. SHP derived from the ethanol measurements showed large differences in van Genuchten-Mualem (VGM) parameters for the M and S plots compared to water infiltration experiments. SWR resulted in large saturated hydraulic conductivity (Ks) decrease on the M and S scenarios. After validation of SHP on water content measurements during a selectedirrigation event, one year simulations (2014) showed that water repellency increases surface runoff in non-structured
soils at hillslopes.
Keywords: Inverse modeling; Water and ethanol infiltration; SHP estimation; Water dynamics; HYDRUS (2D/3D).
INTRODUCTION
Soil water repellency (SWR) is a reduction in the rate ofwetting and retention of water in soil caused by drying and thepresence of various hydrophobic coatings on soil particles. The
physical background is not yet fully understood, but it is widely
accepted that this phenomenon is most likely caused by waterrepellent compounds that can coat soil mineral particles or bepresent as interstitial matter in soil pores (Doerr et al., 2000).
SWR can increase substantially due to seasonal events, e.g.
drought periods and/or wildfires (Jordan et al., 2013; Schwen etal., 2015). With the projected continuation of climate change,
frequency and severity of drought events will intensify in mostregions of the globe and the relevance of SWR effect on soil
water dynamics is expected to increase in the future (Fischer
and Knutt, 2014; Stocker et al., 2013).
SWR is not a stationary soil property but highly variable
over time, depending on soil water content. Generally, SWR
increases with decreasing water content (Jordan et al., 2013).
Furthermore, although dry soils may be very water repellent
initially, the repellency effect can disappear after prolonged
contact with water. The duration of this process is described as
repellency persistence. Besides persistence, SWR is also defined by severity of repellency, both of which can be expressed
quantitatively (Chau et al., 2014). The intensity of SWR ischaracterized by the contact angle between soil surface andinfiltrating water (Subedi et al., 2013). Hydrophobic conditions
are present if contact angle above 90° is present, however eventhe contact angle in between 0° and 90° can affect water infil-

170

tration (Hallett et al., 2001). The SWR cannot be measured
directly as physical value in the field, but has to be obtained
indirectly, e.g. by observing the difference in flow behavior
between water and fully-wetting fluid. Due to specific physicochemical properties of ethanol (i.e. lower surface tension),ethanol is considered to be complete wetting fluid which iscommonly used in SWR estimation (i.e. zero repellency; Lamparter et al., 2010; Watson and Letey, 1970). Further established methods (laboratory) are based on the influence which
SWR has on other soil physical parameters, e.g. contact angle
between water and soil surface, water drop penetration time, or
capillary effect (Letey et al., 2000; Shang et al., 2008).
In comparison to non-repellent soils, water repellent soils
have different infiltration patterns (initially postponed infiltration which increases after contact angle between water and soil
particles decreases e.g. Bughici and Wallach, 2016; Debano,1975) and increased fractions of preferential flow (e.g. Ritsema
et al., 1993; 2000) or surface runoff (Lemmnitz et al., 2008).
Soil infiltration experiments with a KBr tracer performed by
Clothier et al. (2000) demonstrated (i) transient behavior of
fingered preferential flow during the breakdown of hydrophobicity as a result of increasing soil water content, and (ii) solutepenetration of the whole soil pore space after complete wetting.
Hence, SWR appears to be reversible and very dynamic in timeand space, thus making it difficult to predict. In a review by
Jordan et al. (2013) the authors show that the majority of previous SWR studies are focused on the relationship between SWR
and different soil properties (texture, organic matter content,
soil chemical characteristics) or microbiological activity as a

Inverse estimation of soil hydraulic properties and water repellency following artificially induced drought stressresponse to soil management. However, SWR dynamics under
different climatic scenarios and how it can affect soil moisturein the long term remains unclear.
To assess the water dynamics under various field conditions
as well as for predictions of future scenarios, modeling has
proven to be an appropriate tool (Šimůnek et al., 2016; Vereecken et al., 2016). However, in available modeling software
applications for plot or profile scale SWR is not accounted for
as a separate parameter, but may be preferably expressed during
the procedure of soil hydraulic properties (SHP) estimation.
Several studies highlight that the influence of SWR on SHP isevident in the hysteresis effect (Bauters et al., 1998; Czachor etal., 2010). Generally, the hysteresis is highly related to SWR,
and the SWR effect is primarily detectable on the wetting curve(Hardie et al., 2013; Stoffregen and Wessolek, 2014). In alaboratory study focusing on SWR-influenced SHP, Diamantopoulos et al. (2013) performed multistep inflow/outflow
experiments with water and ethanol on four substrates, where
they gradually induced water repellency by adding water repellent material (hydrophobic sand) in different ratios to soil. The
experiments were performed with initially dry or initially saturated conditions to account for hysteresis, and inverse parameter estimation was performed to obtain SHP. Their results
showed that SWR affects SHP on the wetting curve, contributing to the hysteresis effect, and that the artificial mixtureswith a higher fraction of water repellent substances had a larger
effect of SWR on SHP compared to naturally repellent soils.
Therefore, during SHP estimation in water repellent soils, it isimportant to account for the most severe SWR effects expectedduring the initial soil wetting process (infiltration).
Hysteretic forms of SHP are sometimes implemented inmodeling applications to account for SWR, e.g. Nieber et al.
(2000) simulated infiltration in wettable and water repellent
sand with a 2D finite-element model where SWR was taken
into account by including hysteresis in the water retention curve
and two slightly different equations for the unsaturated hydraulic conductivity function. Ganz et al. (2014) modeled water
infiltration patterns in water repellent soils as well, but using a
3D-simulation in HYDRUS (2D/3D). They emphasized astrong need for the inclusion of hysteresis (model implemented
in HYDRUS by Lenhard and Parker, 1992) and a scaling procedure based on independently measured contact angle data(method by Bachmann et al., 2007).
To investigate the impact of different rainfall distribution
patterns on soil water dynamics, modeling using data from an
artificially induced drought stress field experiment was performed in order to investigate the full effect of SWR on thehysteretic wetting curve. The objective of this study was (i) toestimate SHP from disc infiltrometer experiments with waterand ethanol using inverse modeling approach and (ii) to further
assess the impact of different artificially induced drought scenarios on SWR and consequently on local vadose zone hydrology using HYDRUS (2D/3D).
MATERIALS AND METHODS
Field site description
The experimental field was set up at the iLTER-site (International Long Term Ecological Research) in the Rosalian Mountains, Austria (47°42´26.33˝ N, 16°17´54.5˝ E, 600 m a.s.l.;
Leitner et al., 2017). The mean annual temperature is 6.5°C andthe mean annual precipitation is 796 mm at this location. The
experimental site was situated in a forested hillslope with mature beech trees (Fagus silvatica L.) and no understory on aplateau with a sloping angle of 16°. The soil type was classified

as Podsolic Cambisol according to the WRB (World reference
base for soil resources, IUSS, 2014) covering impermeable
granitic bedrock at 75–80 cm below the soil surface following
the hillslope curvature. The soil profile was covered with anorganic matter O horizon (0–7 cm), followed by an eluvial
humus Aeh-horizon (7–25 cm), a cambic, slightly humusosesquioxidic Bhs-horizon (25–50 cm) over weathered granitic rockdebris (C-horizon 50–75 cm) (Schwen et al., 2014). To assess
the impact of changed rainfall distribution patterns on various
soil properties, this experimental trial was established in 2013
(Leitner et al., 2017). Briefly, two artificially induced drought
stress scenarios were applied during the vegetation period(May-October): a moderately (M) stressed scenario which had
six consecutive cycles of four weeks drying followed by anintensive 75 mm irrigation, and a severely (S) stressed scenario
which had three cycles of eight weeks drying followed by alarger irrigation event with 150 mm of irrigation. Stressed plots
were protected from natural rainfall during vegetation periods
by a plastic roof 1.20 m above the soil surface (each treatment
having 4 plots, 2 m × 2 m, Fig. 1). Plots were irrigated withsprinkler irrigation systems with axial-flow full cone nozzles(Series 460, Lechler GmbH) installed under the roofs using
descaled tap water from a nearby field station. The duration of
the irrigation events were 2h each. To be able to compare theresults with natural conditions, four control plots (C) received
only natural rainfall. Drought plots had additional trenches (20cm deep) at the upper end of the plot to avoid any lateral flow
and/or surface runoff from elevated ground to enter the plots. In
each plot, soil volumetric water content was measured at 10 cm
depth in the Aeh horizon (VWC, TDR theta ML2x probes,
UMS, Germany), with measurement intervals of 30 minutes.
Climatic data were collected from a meteorological station
located 500 m from the field site and used to calculate evapotranspiration according to Penman-Monteith (Monteith, 1981).

Fig. 1. iLTER experimental site scheme with bold lines as bridge
pathways, dashed lines as plot roofs, colored squares as experimental units (2 × 2 m), black bold dashes as TDR-probes, X-ed
circles as locations of performed infiltration experiments andbrown points as trees. C stands for control, M for moderate, and S
for severe stress. Infiltration experiments for control areas in
Schwen et al. (2015) were conducted outside of the equippedsquares to minimize influences to the soil system, especially to soil
biology by ethanol.

Measurements of soil hydraulic parameters
Infiltration experiments were performed during September
2014 in all three scenarios (M, S, and C) in four repetitions as
described elsewhere (Schwen et al., 2014, 2015) using a self-

171

Vilim Filipović, Thomas Weninger, Lana Filipović, Andreas Schwen, Keith L. Bristow, Sophie Zechmeister-Boltenstern, Sonja LeitnerTable 1. Basic soil physical properties at the iLTER experimental site (Austria) for Podsolic Cambisol soil profile: clay, silt and sand fraction, total porosity φ, bulk density ρb, and soil organic carbon (OC) content, Values were derived from a soil profile where samples weretaken from incremental 5-cm-layers (uppermost sample from depth = 0–5cm), averaged for soil horizons (standard deviation in brackets).
More detailed in Schwen et al., 2014, Fig. 1(b).
Horizon

Depth
cm

O

0–7

Aeh

7–25

Bhs

25–50

C

50–75

Clay

Silt

Sand

OC

–1

–1

–1

–1

gg
0.120
(–)
0.090
(0.0125)0.088(0.0092)
0.083(0.0382)

gg
0.520
(–)
0.252
(0.0343)0.240(0.0524)0.191(0.0204)

constructed tension disc glass infiltrometer which allowed
using both water and ethanol as infiltrating liquids (Schwen etal., 2015). These data served as an input for SHP estimation
using inverse modeling (the procedure is explained in the next
section). Before the infiltration measurements were made, theorganic litter was removed to ensure that the measurements
were performed at the top of the mineral Aeh horizon.
Additionally, a thin layer of uniform glass beads was used at
the soil surface to ensure good contact between soil and theporous disc (Dragonite, Jaygo Inc.; diameter: 0.45 mm).
Infiltration measurements (water and ethanol) were conductedusing different pressure heads (e.g., –10, –5, –3, –1 cm), in four
replicates (one per plot) with each liquid. At the end of thewater infiltration, the soil was covered to prevent any physicaldisturbance. Two days later, infiltration experiments were
performed at the same specific spots using ethanol as theinfiltration liquid with the glass beads layer replaced where
necessary. Differences in the dynamic viscosity (η) of water (η= 1.0 mPa) and ethanol (η = 1.2 mPa) results in different liquid
infiltration rates, even at identical liquid contents. To be able tocompare the two infiltration experiments, the infiltration ratesof ethanol were corrected for the difference in viscosity
between water and ethanol using a factor of 1.2 (Jarvis et al.,2008). Considering different physicochemical properties ofwater and ethanol, the ethanol pressure head values were scaled
based on the capillary rise equation which takes into accountthe difference between surface tension and density of particularliquids:hi =2σ i cosγ
r ρi g(1)where σ is the surface tension (mN m−1), γ is the contact angle(°), r is the equivalent capillary radius (m), ρ is the density of
the liquid (g cm−3), and g is the acceleration due to gravity (ms−2). The subscript i refers to water (w) or ethanol (e). With thewater and ethanol surface tension at 20°C of 72.7 mN m−1 and22.4 mN m−1, and a density of 0.998 g cm−3 and 0.789 g cm−3,
respectively, a correction factor between he and hw of 2.5 wasassumed (Diamantopoulos et al., 2013; Lamparter et al., 2010).
Multiplying he with 2.5 results in the effective supply pressure(he, eff) giving the applied ethanol pressure heads of –25, –12.5,–7.5, and –2.5 cm. The correction factor assumes that the contact angles of water/ethanol and soil surface are identical (γe =γw), despite the fact that water contact angle in the field was
oscillating in time. However, the differences in the contact
angle reflect SWR through the different infiltration volumes foreach liquid and consequently inversely estimated SHP. Initialwater content was determined using the gravimetric method on

172

gg
0.360
(–)
0.659
(0.0438)
0.673(0.0611)0.726(0.0565)

gg
0.051
(–)
0.016
(0.0095)
0.002(0.0027)
0.000(–)

Porosity,φcm3 cm–3
0.630
(–)
0.510
(0.0326)
0.428(0.0268)
0.378(0.0164)

Bulk density, ρbg cm–3
0.980
(–)
1.300
(0.0941)1.516(0.0712)1.644(0.0391)

undisturbed soil cores of 250 cm3 volume (n = 4 per treatment).
Average initial water content values were 0.17, 0.19, and 0.31
cm3 cm–3 for M, S and C treatments prior to the water infiltration experiments. The same water contents were assumed to bepresent prior to ethanol infiltration experiments (no additional
sampling/measurements were done in order to prevent disruption of the infiltration spots). Using the same soil core samples,
SHP were estimated using the evaporation method (Schindler etal., 2010; device: HYPROP, UMS GmbH, Munich, Germany).
Although SHP was measured by the evaporation method, only
part of the measured data was used to obtain certain hydraulicparameters (porosity, θs) while the rest of the curve fitting was
performed using inverse modeling based on the data from theinfiltration experiments. SWR tend to be more expressed indrier soils and decreases with increased soil moisture (Dekker
and Ritsema, 1994; Diamantopoloulos et al., 2013; Liu et al.,2012), thus starting with initially dried soil and performing
infiltration experiments to investigate SWR in different treatments seemed appropriate. Particle size distribution was determined by a combination of sieving and sedimentation experiments according to Gee and Or (2002). Basic soil physicalproperties for the iLTER experimental site are given in Table 1.
Numerical modeling
Numerical modeling was performed with the HYDRUS(2D/3D) model (Šimůnek et al., 2016) using a three-step simulation process:1. Inverse modeling based on water and ethanol field
infiltration data – to obtain SHP;
2. Simulation of a particular irrigation event using obtained
SHP – to validate the model using field TDR measurements;
3. Seasonal simulation (2014) – to assess the effect of SWR onwater dynamics in hillslope areas.
This approach is explained in more detail in the following
sections.
Inverse modeling to estimate soil hydraulic properties
Tension disc infiltration measurement data (average of four
repetitions) with water and ethanol for the S, M, and C plots(Schwen et al., 2015) were used to obtain SHP using inverse
modeling (Hopmans et al., 2002). A numerical solution of the
Richards’ equation coupled with the Levenberg-Marquardt
nonlinear minimization method implemented in the HYDRUS(2D/3D) model was used. The program solves the equation
numerically using a quasi-three-dimensional axisymmetricfinite element code. The Richards’ equation, which describes
isothermal Darcian flow in a variably saturated rigid porous

Inverse estimation of soil hydraulic properties and water repellency following artificially induced drought stressmedium, is used in the model in its modified form (Šimůnek etal., 1998):∂θ 1 ∂  ∂h  ∂  ∂h  ∂K=
 rK  +  K  +∂t r ∂r  ∂r  ∂z  ∂z  ∂z(2)where θ is the volumetric water content [L3 L–3], h is the pressure head [L], K is the unsaturated hydraulic conductivity[L T–1], r is a radial coordinate [L], z is vertical coordinate [L],positive upwards, and t is time [L]. Equation (2) was solved
numerically for the following initial and boundary conditions
which reflect the initial and boundary conditions of the tension
disc infiltrometer experiment:θ ( r , z,t ) = θi t = 0(3)

h ( r , z , t ) = h0

(4)0 < r < r0 , z = 0∂h ( r , z , t )= −1 r > r0 , z = 0∂z(5)h ( r , z , t ) = hi(6)

r2 + z2 → ∞3θ s −θ r(1 + α h )n mfor h < 0(7)

θ ( h ) = θ s for h ≥ 0 
1 m lm
K ( h ) = K s Se 1 −  1 − Se  2(8)

θ − θrθs − θr(9)1
m =1− ; n > 1
n

(10)

Se =

Modeling water dynamics in the M, S and C scenarios–3

where θi is the initial soil water content [L L ], h0 is the timevariable supply pressure head imposed by the tension disc
infiltrometer for water (–10, –5, –3, –1 cm) and ethanol (–25,–12.5, –7.5, –2.5 cm) [L], and r0 is the disc radius (porous discradius of 2.9 cm) [L]. The SHP, estimated from ethanol
infiltration volumes (scaled to match water physicochemical
properties), were assumed to reflect the water infiltration inhydrophilic soil.
Soil hydraulic functions θ(h) and K(h) used in the inverseand direct simulations (next section) were described using thevan Genuchten-Mualem model (VGM, van Genuchten, 1980)
defined as follows:θ ( h ) =θ r +head, providing initial condition in this form ensures a morestable and unique solution of the inverse problem. Soil surface
boundary conditions below the disc infiltrometer and the remaining soil surface are represented by Eqs. 4 and 5, respectively. Eq. 6 assumes that all subsurface boundaries are distantfrom the supply source and do not influence the results in any
way. The inverse solution was obtained using a combination of
cumulative infiltration data and observed initial/final watercontent after minimization of the objective function. The simulated axisymmetrical domain was 15 cm wide and 20 cm longsoil block with 2501 nodes and increased density along theupper boundary due to the tension disc infiltrometer placement.
The soil hydraulic parameters (θr, α, n, and Ks) were initiallyderived from particle size distribution and bulk density datausing the ROSETTA pedotransfer functions (Schaap et al.,2001) (See Table 2). The θr parameter was not modified, asŠimůnek et al. (1998) and González et al. (2015) found that thisparameter had little effect on the simulated θ and h time series.
The inverse modeling approach proposed by Šimůnek and van
Genuchten (1996) was then used to calibrate α, n, and Ks in the
top soil layer of each treatment starting with the initial Aeh
horizon properties (Table 2).

where θr and θs denote residual and saturated volumetric watercontent [L3 L–3], respectively, Ks is the saturated hydraulicconductivity [L T–1], Se is the effective saturation [–], α [L–1]
and n [–] are shape parameters, and l [–] is a pore connectivity
parameter. Pore connectivity parameter (l) was fixed to 0.5 asrecommended by Mualem et al., (1976) to avoid optimizationof large number of parameters. Please note that the initial condition was given in terms of the soil water content (values presented in field site description chapter). Šimůnek and van
Genuchten (1997) showed that, compared to the use of pressure

After performing inverse VGM parameters estimation, directmodeling was performed. Simulations included selected irrigation event (on June 24th 2014, starting 60 hours before and after
irrigation was performed) in the M and S treatments during one
year period (2014, on a daily time frame) in all three scenarios(M, S, and C). Simulations were performed for twodimensional variably saturated porous media using Richards’
equation:∂θ∂   A ∂h
 K  K ij=
+ K izA   − S∂t ∂xi  
∂x j



(11)

where θ represents the soil volumetric water content [L3 L–3], hpressure head [L], xi (i = 1, 2) the spatial coordinates [L], t time[T], KijA is the components of the dimensionless hydraulicconductivity anisotropy tensor (KA) in the two main spatial
directions xi, K is the unsaturated hydraulic conductivity [L T–1],
and S accounts for root water uptake [L3 L–3 T–1]. The root
water uptake was calculated using the Feddes et al. (1978)
equation. The potential root water uptake was calculated taking
into account seasonal dynamics of LAI and interception of the
canopy at the forest sites (beech trees). A constant rooting
depth (75 cm) with a root distribution adapted from Huang etal. (2011) was assumed. Effective precipitation was further
calculated by subtracting losses due to interception from grossprecipitation. Seasonal variations in beech canopy were estimated using dynamic LAI with a maximum value of 5.8 andinterception capacity of 2.0 mm (Armbruster et al., 2004; Breuer et al., 2003). Based on the assumptions stated above, potential root water uptake was calculated in the study of Schwen etal. (2014), performed on the same site.
The simulated domain in 2D vertical space (for specific irrigation event simulation and one year modeling) was 0.75 m
deep and 2 m long (corresponding to one plot; see Figure 6 in
the results section). Atmospheric boundary conditions wereselected at the top and seepage conditions at the right side
(down slope) to mimic the possible lateral subsurface movement as the soil profile was located on an impermeable sloped173

Vilim Filipović, Thomas Weninger, Lana Filipović, Andreas Schwen, Keith L. Bristow, Sophie Zechmeister-Boltenstern, Sonja LeitnerTable 2. Van Genuchten-Mualem (VGM) soil hydraulic parameters derived from pedotranfer soil functions (PTFs, Rosetta) basedon soil texture and bulk density (Table 1) with the measured saturated water content value θs based on evaporation experiments(Schwen et al., 2014).


R 2 = Horizon  n ( S − O )2 ii
RMSE =  i =1
n



Aeh
Bhs
C

θr
cm3 cm–3
0.04330.04040.0391θs
cm3 cm–3
0.47
0.35
0.32

α
cm–1
0.026
0.03350.0427n–1.45541.45631.4735Ks
cm day–1
88.37
44.22
36.8

l
–
0.5
0.5
0.5

bedrock. The simulation domain had 14882 nodes with theincreased density at the top boundary (with 29350 2D elements). The soil layering and hydraulic properties were selectedaccording to the Table 2, with the Aeh horizon extended to thesoil surface.
Although Richards’ equation is not considered applicable forhydrophobic medium (e.g., Diamantopoulos and Durner, 2013),
because SWR is a reversible process and not a constant state,some of the classical physical approaches are still suitable whencritical water content is exceeded.
Numerical simulations (i.e., infiltration volumes and TDR
measurements) were evaluated using the coefficient of determination (R2), root mean square error (RMSE) and model efficiency coefficient (E) (Nash and Sutcliffe, 1970):




n
n
 i =1(Oi − O )  i =1( Si − S ) 

 i =1(Oi − O ) ( Si − S )
n

(12)

0.5

2

E =1−

n
 i =1( Oi − Si )2n
 i =1( Oi − O )(13)where Oi and Si are observed and simulated values, respectively, O and S represent the averages of observed and simulated
values, respectively, and n is the number of observed/simulated
points.
RESULTS AND DISCUSSION
Inverse simulations of tension disc infiltrometer data

Tension disc infiltration data were used to estimate SHP by
performing inverse modeling with HYDRUS (2D/3D). Modeled
data is compared with field data in Fig. 2. Water infiltration for
the M and S scenarios had a very steady slow inflow rate

Fig. 2. Measured (circle) and simulated (solid line) infiltration cumulative fluxes for the M (moderately stressed), S (severe stressed) and C(control) scenarios with fitting performed using inverse optimization for water and ethanol liquids in HYDRUS (2D/3D) with an indicationof the imposed pressure head and its duration (dotted line). Ethanol curves (_ethanol) are scaled in order to take account the different physicochemical properties and are directly comparable to water curves (_water).

174

Inverse estimation of soil hydraulic properties and water repellency following artificially induced drought stress
Table 3. VGM parameters derived with inversion procedure using HYDRUS(2D/3D) from field infiltrometer data performed with water(_w) and ethanol (_e) liquids and statistical parameters (R2, E, RMSE) describing goodness of model fitting for Aeh horizon.
Scenario
M_w
M_e
S_w
S_e
C_w
C_e

α
(cm–1)
0.036
0.100
0.025
0.055
0.101
0.056

n
(–)
1.717
1.713
1.389
2.178
1.318
1.393

Ks
cm day–1
23.45
1454.2627
440
257.02
460.8

R2

E

0.998
0.981
0.995
0.999
0.989
0.982

0.998
0.985
0.995
0.999
0.9910
0.981

RMSE
(cm3)
2.77
21.6
4.14
4.25
5.11
11.03

dynamics (e.g. van Genuchten, 1980). Table 3 shows a large
difference between soil hydraulic properties with water andethanol for all optimized parameters. A large increase in Ks
values in the M and S scenarios was found when ethanol as acomplete wetting liquid was used (Table 3) compared to water
infiltration where Ks was lower. Decreases of 98.3% and 93.8%
in Ks values were noticed when comparing ethanol and water
infiltration for M and S scenarios. In contrast, the control scenario shows a lower increase of 44.2% compared to the water
infiltration measurements. This difference also indicates that

Fig. 3. Plotted retention curves for different VGM soil hydraulicparameter sets for various precipitation manipulation scenarios anddifferent liquids (w, water vs. e, ethanol). Ethanol curves (_e) are
scaled in order to take account the different physicochemical properties and are directly comparable to water curves (_w). M, moderately stressed scenario; S, severely stressed scenario; C, control
scenario.

during the entire infiltration experiment and did not respond tothe change in supply pressure. By contrast, a typical water
infiltration curve found for non-repellent soils was observed in
scenario C (although with reduced effect e.g. Šimůnek and van
Genuchten, 1997), showing increased infiltration volumes atlower pressure heads (close to saturation). On the contrary, for
M and S scenarios data showed reduced infiltration at lower
pressure heads (h = –3 cm and –1 cm) and indicated that thelarger pores were more hydrophobic than the smaller pores(Leue et al., 2015; Schwen et al., 2015). The final infiltration
volumes were similar with 200.4, 185.7 and 175 ml for M, S
and C scenarios, respectively. However, if these results arecompared to the ethanol experimental data, it can be seen thatthe ethanol infiltration volumes are larger, with the final volumes of 519, 531 and 303 ml for M, S and C scenarios, respectively, indicating increased water repellency in the M and S
scenarios. These data also show that the C scenario showed thesmallest difference in total volume and infiltration curve behavior between water and ethanol. Similar results were observed by
Jarvis et al. (2008) where they compared grassland to arable
land with water and ethanol measurements and found that water
repellency is smaller when water and ethanol infiltration volumes are similar (and vice versa).
The inverse optimization modeling worked well in both cases, showing a good fit for both infiltrating liquids, which can beobserved visually and through statistical indicators (Fig. 2;
Table 3). This shows that the VGM model describes soil hydraulic parameters well at this site. Observations through theprofile indicate the presence of few macropores with a rather
uniform structure (Schwen et al., 2014), suggesting that a single
porosity model would be sufficient to describe the soil water

Fig. 4. Pressure head distribution for M (moderately stressed) S(severely stressed) and C (control) scenarios at the end of the tension infiltration field experiment using water (_w) and ethanol (_e).

175

Vilim Filipović, Thomas Weninger, Lana Filipović, Andreas Schwen, Keith L. Bristow, Sophie Zechmeister-Boltenstern, Sonja LeitnerFig. 5. Measured and simulated soil volumetric water content at 10 cm depth (Aeh horizon, TDR) using optimized parameters from water
infiltration experiments from M (moderately stressed) S (severely stressed) and C (control) scenarios for the irrigation event at 24th of June2014 (60 h before and after the event) with a 75 mm after 4 weeks of drought (M plots) and 150 mm irrigation after 8 weeks of drought (Splots) of irrigation for the event duration of 120 h.

Fig. 6. Pressure head distribution for M (moderately stressed) S(severely stressed) and C (control) scenarios on the August 20th
2014 after 150 mm of the irrigation event on the previous day.

176

the control scenario had a certain degree of soil water repellency, which was most likely linked to the natural rainfall distribution e.g. critically low soil water content during longer periods
with no rain (Schwen et al., 2015).
Fig. 3 shows fitted water retention curves based on the Table
2 parameter set for measurements derived from water (fulllines) and ethanol infiltration (dotted lines). Again greater
differences in the shape of the retention curve at the larger
tension (dryer soil condition) were recorded between the artificially induced drought stress scenarios, while the control plot
had a similar shape of the retention curve for both, water andethanol obtained SHP. It should be noted that the small differences in water retention curves of control plot might be also
connected to the fitting procedure in HYDRUS since morethan one set of SHP can be realistic. In addition to changes in
Ks values, α and n values were increased significantly during
the optimization process for ethanol data (M and S plots), resulting in a different curve shape.
Pressure head snapshots were taken at the final time of thesimulated disc tension infiltration field experiment performedin HYDRUS (2D/3D) (Fig. 4). These data show that the water
infiltration plume in M and S scenarios is significantly reduced
compared with the ethanol infiltration plume, resulting in a very
dry bottom part of the soil block (20 cm depth) even after 1.5
hours of infiltration. This part of the soil block had a pressure
head of approximately –150 and –650 cm for M and S scenarios,
respectively, while the ethanol treated soil had a pressure headaround –75 cm. The control plot also showed differences in final
pressure head values, but the lowest pressure heads were –41.58
cm for water and –54.21 cm for ethanol infiltration, respectively.
The differences in duration and infiltration volumes between
the two applied liquids should also be taken into consideration.
Occurrence of such extensive infiltration reduction increases
the risk of preferential flow in structured soils (Jarvis et al.,2008) and the potential of surface runoff as well (Cerdà and
Doerr, 2007), which is further evaluated through simulations ofartificially induced drought stress scenarios in the next section.

Inverse estimation of soil hydraulic properties and water repellency following artificially induced drought stressFig. 7. Simulated cumulative seepage face flux (at the right side boundary, a), cumulative infiltration (b), cumulative runoff (c) and volumetric water content (d) at the M (moderately stressed) S (severely stressed) and C (control) scenarios during 2014 using natural rainfall
and irrigation events.

Simulations of artificially induced drought stress scenarios

After the SHP were obtained by inverse modeling, a direct
simulation of a particular irrigation event (on June 24th 2014)
was performed for the M and S scenarios and also for the control scenario without any irrigation (Fig. 5). The selected event
was simulated with a temporal resolution of 1h to provide detailed information regarding the model performance using theobtained SHP. The simulation was set to start and finish 60
hours before and after the irrigation event. Because the plots
were covered during the rest of the season, simulating such
shorter periods before and after the irrigation event in order tovalidate the model was assumed sufficient. The simulations
were performed using the properties of the Aeh horizon derived
from the water infiltration inverse modeling and compared towater contents measured by TDR probes at 10 cm depth (Aehhorizon) while the rest of the horizons had the same values foreach layer as presented in the Table 2. Simulation of watercontent reproduced the field observation well with R2 values
above 0.9 and confident model efficiency values, while themodel efficiency was low (–0.50) in the control scenario due tothe absence of any variation during this limited period. Still, theRMSE of 0.005 indicated a good fit of the data of the control
scenario. These model performance results for water content
simulations are in line with previous studies using the HYDRUS model (e.g., Ajdary et al., 2007; Kandelous et al., 2011;
Nakhaei and Šimůnek, 2014).
Fig. 5 clearly shows a dependence between soil water content and the amount of applied irrigation water for the M and S
scenarios, e.g. an instant increase of soil water content in the M
scenario from 0.12 to 0.29 cm3 cm–3 at the time 60.5 hours isobserved, and in S scenario from 0.15 to 0.34 cm3 cm–3, bothcorresponding to the 75 and 150 mm of irrigation, respectively.
The control scenario did not show any increase in soil watercontent due to the absence of natural rainfall and irrigation.

After obtaining hydraulic properties from the infiltration experiment and validating the model by comparing it to fieldwater content measurements, one year simulations were performed for 2014 using a daily time step. The simulation included plant uptake and natural rainfall and irrigation events, inorder to maximize the difference between the scenarios (interms of water balance) and to reveal potential downsides of thewater repellency (e.g., surface runoff, low water content). Fig. 6
shows the pressure head distribution on August 20th, 2014 (oneday after irrigation events on both irrigated plots) in the three
scenarios. The differences between the scenarios are seemingly
negligible for the M and S scenarios. The low conductivity of
the Aeh layer (23.45 cm day-1 and 27 cm day–1 for the M and S
scenario, respectively) delayed the infiltration of the irrigation
plume and induced surface runoff. On the contrary, due to its
non-repellent state which was expressed in larger Ks values(257.02 cm day–1), the irrigation plume saturated simulated
profile till approximately 50 cm depth.
Our simulations suggest that in this particular case, the effectof hydrophobicity on large scale hydrology can have a substantial impact on water balance and its distribution between theinfiltration and surface runoff. Fig. 7 shows differences in waterbalance (e.g., seepage face flux, infiltration rate and surface
runoff) between artificially induced drought stress scenarios
and the control plot. Due to similar SHP for the two imposed
scenarios, the water balance and its distribution is almost identical. However, large seepage face flux (which in this casemimics subsurface lateral flow because of the impermeable
bedrock at 75 cm) can be observed in the control plot (34.7 cm
vs. 15 cm in the drought scenarios). Small Ks values in the M
and S scenarios resulted in decreased infiltration and subsurface
flow but increased surface runoff (77 cm in drought stressscenarios vs. 59 cm in control). The simulated water content(Fig. 7d) followed the rainfall and irrigation inputs very welland reflected the increased water levels at the top as well. The

177

Vilim Filipović, Thomas Weninger, Lana Filipović, Andreas Schwen, Keith L. Bristow, Sophie Zechmeister-Boltenstern, Sonja Leitnerextent of the oscillations was linked directly to the estimated
SHP, e.g. a quicker response of volumetric water content
changes can be seen in a control plot with the applied irrigation(simulated). Our results are in accordance with previous studies(e.g. Lemmnitz et al., 2008) which indicate that the occurrenceof soil water repellency on hillslopes is important when addressing larger scale soil hydrology on a seasonal basis.
CONCLUSIONSThe HYDRUS (2D/3D) model was used to estimate soil hydraulic properties from field tension disc infiltration experiments and to quantify soil water repellency effects using datafrom water and ethanol infiltration measurements. Additionally,
simulations to fit the TDR field moisture measurements and to
perform 2D water balance modeling on an artificially induced
drought stress field experiment were conducted. The water andethanol infiltration experiments showed a large variation among
the treatments and control scenario, revealing the importance of
prolonged soil drying on soil water repellency. The inverse
modeling was performed successfully with R2 and model efficiency (E) values above 0.9, indicating good fit with the fieldmeasured infiltration data for both liquids (water and ethanol).
Soil hydraulic properties derived from the ethanol measurements showed significantly greater Ks values for the M and S
scenarios, thus suggesting linkage between water repellency
and reduced infiltration. Direct simulation of irrigation events
showed good reliability of the model to fit water contents
measured at 10 cm depth using TDR probes. One year simulations (2014) showed that the non-structured water repellent
soils have a potential to produce increased surface runoff, as
well as reduced subsurface lateral flow (if impermeable or low
conductivity layer is present) or vertical drainage. Climaticchange scenarios are predicting more intense and prolonged
droughts, as well as more extreme rainfall events, which canlead to increased soil water repellency and result in changed
water flow patterns at the plot scale. Further studies are needed
to clarify the occurrence, non-linear nature and impact of SWR.
Acknowledgements. This project was partly funded by the Austrian Climate and Energy Fund (ACRP 6–DRAIN–
KR13AC6K11008). Sonja Leitner received a PhD fellowship of
the AXA Research Fund. We also acknowledge travel funding
from the bilateral agreement project HR20/2016 funded by theAustrian Agency for International Cooperation in Education(OeAD) and the Ministry of the Science and Education of the
Republic of Croatia.
REFERENCES

Ajdary, K., Singh, D.K., Singh, A.K., Khanna, M., 2007. Modelling of nitrogen leaching from experimental onion fieldunder drip irrigation. Agr. Water Manage., 89, 15–28.
Armbruster, M., Seegert, J., Feger, K.H., 2004. Effects ofchanges in tree species composition on water flow dynamics– model applications and their limitations. Plant Soil, 264,
13–24.
Bachmann, J., Deurer, M., Arye, G., 2007. Modeling water
movement in heterogeneous water-repellent soil: 1. Development of a contact angle–dependent water-retention model.
Vadose Zone J., 6, 436–445.
Bauters, T.W.J., DiCarlo, D.A., Steenhuis, T.S., Parlange,
J.-Y., 2000. Soil water content dependent wetting front characteristics in sands. J. Hydrology, 231–232, 244–254.

178

Breuer, L., Eckhardt, K., Frede, H.G., 2003. Plant parametervalues for models in temperate climates. Ecol. Model., 169,
237–293.
Bughici, T., Wallach, R., 2016. Formation of soil-water repellency in olive orchards and its influence on infiltration pattern. Geoderma, 262, 1–11.
Cerdà, A., Doerr, S.H., 2007. Soil wettability, runoff and erodibility of major dry-Mediterranean land use types on calcareous soils. Hydrol. Process., 21, 17, 2325–2336.
Chau, H.W., Biswas, A., Vujanovic, V., Si, B.C., 2014. Relationship between the severity, persistence of soil water repellency and the critical soil water content in water repellent
soils. Geoderma, 221–222, 113–120.
Clothier, B.E., Vogeler, I., Magesan, G.N., 2000. The breakdown of water repellency and solute transport through a hydrophobic soil. J. Hydrol., 231–232, 255–264.
Czachor, H., Doerr, S.H., Lichner, L., 2010.Water retention of
repellent and subcritical repellent soils: new insights frommodel and experimental investigations. J. Hydrol., 380,
104–111.
Debano, L.F., 1975. Infiltration, evaporation, and water
movement as related to water repellency 1. In: Gardner,
W.R., Moldenhauer, W.C. (Eds.): Soil Conditioners. SSSA
Spec. Publ. 7. SSSA, Madison, WI., pp. 155–164.
Dekker, L.W., Ritsema, C.J., 1994. How water moves in awater repellent sandy soil, 1. Potential and actual water repellency. Water Resour. Res., 30, 2507–2517.
Diamantopoulos, E., Durner, W., 2013. Physically-based modelof soil hydraulic properties accounting for variable contact
angle and its effect on hysteresis. Adv. Water Resour., 59,
169–180.
Diamantopoulos, E., Durner, W., Reszkowska, A., Bachmann,J., 2013. Effect of soil water repellency on soil hydraulicproperties estimated under dynamic conditions J. Hydrol.,486, 175–186.
Doerr, S.H., Shakesby, R.A., Walsh, R.P.D., 2000. Soil water
repellency: its causes, characteristics and hydrogeomorphological significance. Earth Sci. Rev., 51, 33–65.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of
Field Water Use and Crop Yield. John Wiley & Sons, NewYork.
Fischer, E.M., Knutti, R., 2014. Detection of spatially aggregated changes in temperature and precipitation extremes. Geophys. Res. Lett., 41, 2, 547–554.
Ganz, C., Bachmann, J., Noell, U., Diujnisveld, W.H.M., Lamparter, A., 2014. Hydraulic modeling and in situ electricalresistivity tomography to analyze ponded infiltration into awater repellent sand. Vadose Zone J., 13, 1, 1–14.
Gee, G.W., Or, D., 2002. Particle-size analysis. In: Dane, J.H.,
Topp, G.C. (Eds.): Methods of Soil Analysis. Part 4 – Physical Methods. SSSA Book Series, No. 5. SSSA, Madison,
WI, USA, pp. 1381–1402.
González, M.G., Ramos, T.B., Carlesso, R., Paredes, P., Petry,
M.T., Martins, J.D., Aires, N.P., Pereira, L.S., 2015. Modelling soil water dynamics of full and deficit drip irrigated maize cultivated under a rain shelter. Biosyst. Eng., 132, 1–18.
DOI: 10.1016/j.biosystemseng.2015.02.001.
Hallett, P.D., Baumgartl, T., Young, I.M., 2001. Subcriticalwater repellency of aggregates from a range of soil management practices. Soil Sci. Soc. Am. J., 65, 184–190.
Hardie, A.H., Lisson, S., Doyle, R.B., Cothing, W.E., 2013.
Evaluation of rapid approaches for determining the soil water retention function and saturated hydraulic conductivity ina hydrologically complex soil. Soil Till. Res., 130, 99–108.

Inverse estimation of soil hydraulic properties and water repellency following artificially induced drought stress

Hopmans, J.W., Šimůnek, J., Romano, N., Durner, W., 2002.
Inverse methods. In: Dane, J.H., Topp, G.C. (Eds.): Methodsof Soil Analysis. Part 4. SSSA Book Series, No. 5. SSSA,
Madison, WI, pp. 963–1008.
Huang, M., Barbour, S.L., Elshorbagy, A., Zettl, J.D., Si, B.C.,
2011. Water availability and forest growth in coarse-texturessoils. Canad. J. Soil Sci., 91, 199–210.
IUSS, 2014. World reference base for soil resources. FAO,
Rome.
Jarvis, N., Etana, A., Stagnitti, F., 2008. Water repellency, nearsaturated infiltration and preferential solute transport in amacroporous clay soil. Geoderma, 143, 223–230.
Jordán, A., Zavala, L.M., Mataix-Solera, J., Doerr, S.H., 2013.
Soil water repellency: origin, assessment and geomorphological consequences. Catena, 108, 1–5.
Kandelous, M.M., Šimůnek, J., van Genuchten, M.Th, Malek,K., 2011. Soil water content distributions between two emitters of a subsurface drip irrigation system. Soil Soil Sci. Soc.
Am. J., 75, 488–497.
Lamparter, A., Bachmann, J., Deurer, M., Woche, S.K., 2010.
Applicability of ethanol for measuring intrinsic hydraulicproperties of sand with various water repellency levels. Vadose Zone J., 9, 445–450.
Leitner, S., Minixhofer, P., Inselsbacher, E., Keiblinger, K.M.,
Zimmermann, M., Zechmeister-Boltenstern, S., 2017. Shortterm soil mineral and organic nitrogen fluxes during moderate and severe drying-rewetting events. Appl. Soil Ecol.,114, 28–33.
Lemmnitz, C., Kuhnert, M., Bens, O., Güntner, A., Merz, B.,
Hüttl, R.F., 2008. Spatial and temporal variations of actual
soil water repellency and their influence on surface runoff.
Hydrol. Process., 22, 1976–1984.
Lenhard, R.J., Parker, J.C. 1992. Modeling multiphase fluid
hysteresis and comparing results to laboratory investigations.
In: Genuchten, M.Th., Leij, F.J., Lund, L.J. (Eds.); Proc.
Intl. Workshop on Indirect Methods for Estimating the Hydraulic Properties of Unsaturated Soils. University of California, Riverside, CA.
Letey, J., Carrillo, M.L.K., Pang, X.P., 2000. Approaches tocharacterize the degree of water repellency. J. Hydrol., 231–
232, 61–65.
Leue, M., H.H. Gerke, Godow, S.C., 2015. Droplet infiltrationand organic matter composition of intact crack and biopore
surfaces from clay-illuvial horizons. J. Plant Nutr. Soil Sci.,
178, 250–260.
Liu, H., Ju, Z., Bachmann, J., Horton, R., Ren, T., 2012. Moisture-dependent wettability of artificial hydrophobic soils andits relevance for soil water desorption curves. Soil Sci. Soc.
Am. J., 76, 342–349.
Monteith, J.L., 1981. Evaporation and surface temperature. Q.
J. R. Meteorol. Soc., 107, 1–27.
Mualem, Y., 1976. A new model for predicting the hydraulicconductivity of unsaturated porous media. Water Resour.
Res., 12, 3, 513–521.
Nakhaei, M., Šimůnek, J., 2014. Parameter estimation of soil
hydraulic and thermal property functions for unsaturated porous media using the HYDRUS-2D code. J. Hydrol. Hydromech., 62, 7–15.
Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting throughconceptual models. Part I. A discussion of principles. J. Hydrol., 10, 282–90.
Nieber, J., Bauters, T.W.J., Steenhuis, T.S., Parlange, J.Y.,
2000. Numerical simulation of experimental gravity-driven
unstable flow in water repellent sand. J. Hydrol., 231–232,
295–307.

Ritsema, C.J., Dekker, L.W., 2000. Preferential flow in water
repellent sandy soils: principles and modeling implications.
J. Hydrol., 231–232, 308–319.
Ritsema, C., Dekker, L.W., Hendrickx, J.M.H., Hamminga, W.,
1993. Preferential flow mechanism in a water repellent
sandy soil. Water Resour. Res., 29, 2183–2193.
Schaap, M.G., Leij, F.J., van Genuchten, M.T., 2001. ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol.,
251, 163–176.
Schindler, U., Durner, W., von Unold, Georg., Müller, L., 2010.
Evaporation method for measuring unsaturated hydraulicproperties of soils: extending the measurement range. Soil
Sci. Soc. Am. J., 74, 1071–1083.
Schwen, A., Zimmermann, M., Bodner, G., 2014. Verticalvariations of soil hydraulic properties within two soil profiles and its relevance for soil water simulations. J. Hydrol.,
516, 169–181.
Schwen, A., Zimmermann, M., Leitner, S., Woche, S.K., 2015.
Soil Water Repellency and its Impact on Hydraulic Characteristics in a Beech Forest under Simulated Climate Change.
Vadose Zone J., 14, 12, 1–11.
Shang, J., Flury, M., Harsh, J.B., Zollars, R.L., 2008. Comparison of different methods to measure contact angles of soil
colloids. J. Colloid Interface Sci., 328, 299–307.
Šimůnek, J., van Genuchten, M.Th., 1996. Estimating unsaturated soil hydraulic properties from tension disc infiltrometerdata by numerical inversion. Water Resour. Res., 32, 2683–
2696.
Šimůnek, J., van Genuchten, M.Th., 1997. Parameter estimation of soil hydraulic properties from multiple tension disc
infiltrometer data. Soil Sci., 162, 383–398.
Šimůnek, J., Angulo-Jaramillo, R., Schaap, M.G., Vandervaere,
J.P., van Genuchten, M.T., 1998. Using an inverse method
to estimate the hydraulic properties of crusted soils from tension disc infiltrometer data. Geoderma, 86, 61–81.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2016. Recent
developments and applications of the HYDRUS computer
software packages. Vadose Zone J., 15, 7. DOI:
10.2136/vzj2016.04.0033
Stocker, T.F., Qin, D., Plattner, G.-K. Tignor, M.M.B., Allen,
S.K., Boschung, J., Nauels, A., Xia, Y., Bex, V., Midgley,
P.M. (Eds.), 2013. Climate Change. The Physical Science
Basis. Summary for Policymakers. Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge Univ. Press,
Cambridge and New York, pp. 1–30.
Stoffregen, H., Wessolek, G., 2014. Scaling the hydraulic functions of a water repellent sandy soil. Int. Agrophys., 28,
349–358.
Subedi, S., Kawamoto, K., Komatsu, T., Moldrup, P., Wollesende Jonge, L., Müller, K., Clothier, B., 2013. Contact angles
of water-repellent porous media inferred by tensiometerTDR probe measurement under controlled wetting and drying cycles. Soil Sci. Soc. Am. J., 77, 1944–1954.
van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil
Sci. Soc. Am. J., 44, 892–898.
Vereecken, H., Schnepf, A., Hopmans, J.W., Javaux, M., Or,
D., Roose, T., Vanderborght, J., Young, M.H., Amelung,W., Aitkenhead, M., Allison, S.D., Assouline, S., Baveye,
P., Berli, M., Brüggemann, N., Finke, P., Flury, M., Gaiser,T., Govers, G., Ghezzehei, T., Hallett, P., Hendricks
Franssen, H.J., Heppell, J., Horn, R., Huisman, J.A., Jacques,D., Jonard, F., Kollet, S., Lafolie, F., Lamorski, K., Leitner,

179

Vilim Filipović, Thomas Weninger, Lana Filipović, Andreas Schwen, Keith L. Bristow, Sophie Zechmeister-Boltenstern, Sonja LeitnerD., McBratney, A., Minasny, B., Montzka, C., Nowak, W.,
Pachepsky, Y., Padarian, J., Romano, N., Roth, K., Rothfuss,Y., Rowe, E.C., Schwen, A., Šimůnek, J., Tiktak, A., Van
Dam, J., van der Zee, S.E.A.T.M., Vogel, H.J., Vrugt, J.A.,
Wöhling, T., Young, I.M., 2016. Modeling soil processes:

Review, key challenges, and new perspectives. Vadose Zone
J., 15, 5. DOI: 10.2136/vzj2015.09.0131.
Watson, C.L., Letey, J., 1970. Indices for characterizing soilwater repellency based upon contact angle-surface tension
relationships. Soil Sci. Soc. Am. Proc., 34, 841–844.
Received 17 July 2017
Accepted 29 December 2017

Note: Colour version of Figures can be found in the web version of this article.

180