Geoderma 380 (2020) 114699

Contents lists available at ScienceDirect

Geodermajournal homepage: www.elsevier.com/locate/geoderma

Estimation of vineyard soil structure and preferential flow using dye tracer,
X-ray tomography, and numerical simulations

T

Vilim Filipovića, , Jasmina Defterdarovića, Jiří Šimůnekb, Lana Filipovića, Gabrijel Ondrašeka,Davor Romića, Igor Bogunovićc, Ivan Mustaća, Josip Ćurićd, Radka Kodešováe
⁎a

University of Zagreb Faculty of Agriculture, Department of Soil Amelioration, Zagreb, Croatia
Department of Environmental Sciences, University of California Riverside, Riverside, CA, USAc
University of Zagreb, Faculty of Agriculture, Department of General Agronomy, Zagreb, CroatiadUniversity of Zagreb, School of Medicine Clinical Hospital Dubrava, Department of Diagnostic and Interventional Radiology, CroatiaeUniversity of Life Sciences Prague, Faculty of Agrobiology, Food and Natural Resources, Department of Soil Science and Soil Protection, Prague, Czech Republicb

ARTICLE INFOABSTRACT

Handling Editor: Yvan Capowiez

The appearance and distribution of soil pores have a significant influence on water flow and solute transport in
the soil vadose zone. The pore system is highly variable in arable soils where crop rotation, tillage, trafficking,soil amendments, and various management practices are commonly implemented. The aim of this study was toassess the porous system and preferential flow pathways in a vineyard soil using undisturbed soil columns, andby combining laboratory and numerical methods with dye staining and X-ray imaging. It was hypothesized thatthe integration of various methods could reveal more information about soil structure, and flow and transport
behavior of structured arable soil. Soil water retention and hydraulic conductivity curves were obtained using
the evaporation method, while water flow was assessed using intermittent leaching experiments. Water flow andthe transport of Brilliant Blue were simulated using HYDRUS-1D. A single-porosity model of soil hydraulicproperties provided a good description of data collected during the evaporation experiments. Data collectedduring leaching experiments did not provide enough experimental evidence for the occurrence of nonequilibrium flow patterns and the differentiation between the single- and dual-permeability models of soil hydraulic properties. However, dye staining and X-ray imaging revealed a complex pore-architecture network withlarge vertical and horizontal biopores. The staining patterns (Brilliant Blue FCF) within the vertical column
sections documented the extent of preferential flow. The study showed that the bi-modal character of pore
structure could often be hidden when a limited number or non-adequate methods are applied for its quantification from water flow behavior. The impact of preferential pathways on dye transport can be investigated withobservations and simulations. A combination of various methods enabled us to adequately assess vineyard soil
structure and fine-tune the description and extent of preferential water flow.

Keywords:
Arable soil structure
Dye staining
Nonequilibrium flow
X-ray imaging
Numerical modeling
HYDRUS

1. Introduction
Most pronounced soil heterogeneity in terms of structural changes ispresent in the topsoil layer of arable soils, which are heavily modified
during each crop growing season by tillage and trafficking (Jirků et al.,2013; Strudley et al., 2008). Changes in climate, biological activity, andsoil management practices may have a strong influence on soil structure(Fér et al., 2016, 2019; Kodešová et al., 2011). Pagliai et al. (2004)
demonstrated deterioration of the soil structure in arable soils due toenhanced oxidation of organic matter by aeration, mechanical dispersion because of moist soil compression, or disaggregation of bare soil by
rainfall. A newly developed pore system (e.g., after the tillage) in arable⁎

soils, even with a well-connected network of larger pores, is often unstable and changes rapidly. As reported by Horn (2004), the temporalvariability of soil structure and its evolution is influenced by various
biotic and abiotic factors, and it is often challenging to properly estimate its influence on subsurface flow and transport properties. In the
Mediterranean, vineyards are one of such agricultural systems withenhanced soil structure degradation due to topography, climatic conditions, intensive management, and trafficking practices (Bogunović
et al., 2020; Galati et al., 2015; Prosdocimi et al., 2016a, 2016b). Viticulture areas register some of the highest rates of soil and water losses
in Europe as a result of employing unsustainable practices and management (e.g., large interrow widths, low leaf cover, steep slopes, andCorresponding author.
E-mail address: vfilipovic@agr.hr (V. Filipović).

https://doi.org/10.1016/j.geoderma.2020.114699
Received 27 April 2020; Received in revised form 26 August 2020; Accepted 29 August 2020
0016-7061/ © 2020 Elsevier B.V. All rights reserved.

Geoderma 380 (2020) 114699

V. Filipović, et al.

conventional tillage; Cerdan et al., 2010). Thus, vineyard soils present
an interesting research topic about soil structure heterogeneity andpossible occurrence of preferential flow.
Preferential flow, i.e., unequal or nonequilibrium flow, is a phenomenon in which liquid, gas, and dissolved chemicals primarily movethrough a small fraction of the porous media, such as wormholes, root
holes, and cracks, thus allowing much faster transport (Gerke, 2006;
Ranjit Kumar et al., 2017). The most common types of preferential flow
are macropore flow (through biopores or erosion cracks), finger flow
(unequal water flow due to the hydraulic instability or heterogeneity),
and funnel flow (flow above sloping soil layers). At some sites, thecombination of different flow patterns can occur as well (Guo and Lin,2018). Preferential flow is mostly dependent on pore characteristics
and soil properties, especially soil texture and structure, as well as on
hydraulic properties, and the human impact through soil management
practices in agricultural soils. Macropores have a positive influence on
preferential flow due to their number, size, linkage, and continuity(Jarvis, 2007). The presence of preferential transport reduces the contact time of solutes with the soil matrix, making their adsorption difficult and resulting in significant environmental issues, such as
groundwater pollution. Additionally, preferential flow influences thehydrological cycle, soil erosion, agricultural production, and various
biochemical and ecological processes within the vadose zone. It has apositive influence over flood regulation, erosion prevention, andaquifer recharge, but a negative impact in terms of faster contaminanttransport, thus affecting water quality (Guo and Lin, 2018).
The effects of vegetation are two-sided. On the one hand, the rootsystem increases the abundance of soil pores. On the other hand, thepresence of the vegetation cover has resulted in a decrease of preferential flow (Stumpp and Maloszewski, 2010). Soil fauna has the same
positive influence as the vegetation. Generally, soil fauna influences
many physical and chemical processes in soil (Kooistra and Pulleman,2018). It has an essential role in altering various soil processes, such asorganic matter decomposition and soil structure formation. Soil compaction from trafficking, which frequently occurs in intensively managed arable landscapes, can increase lateral preferential flow andgenerate the destruction of the macropore system (Jarvis, 2007; Sander
and Gerke, 2007).
To decrease nutrient and pesticide losses, and thus to protect surface
and subsurface water resources, extensive understanding of soil water
preferential flow and transport of dissolved substances is fundamental(Salvador et al., 2011). Nonequilibrium flow and transport in structured
soil are highly complex. The characterization of these processes thus
requires a combination of various methods, which can result in their
better quantification.
Dye staining is an efficient method for assessing preferential flow insoils at different scales (e.g., column, profile). Dyes with the following
characteristics are used for the monitoring of water flow within the soil
profile and the estimation of soil pore system appearance, size, anddistribution (Flury and Flühler, 1994, 1995): visibility, non-toxicity,
non-accumulation in soil flora and fauna, and minimal subsequent degradation in the soil. The most commonly used dye tracer in soil hydrology on both field and column scale is Brilliant Blue FCF (e.g.,
Fuentes et al., 2014; Schneider et al., 2018). For example, Bebej et al.
(2017) described the interactions between the Brilliant Blue dye solution with soil particles and their effect on the mobility of substances
around the zones of preferential flow. Sorption of Brilliant Blue FCF in
soils was found to be affected by soil pH, ionic strength, and soil
composition (Germán-Heins and Flury 2000). Kodešová et al. (2015)
confirmed that dye distribution is significantly affected by the presenceof large pores due to the earthworm population, root system spreading,and organic matter in the soil. It was also noted that the influence of
topology, macropores, and soil water content is important.
However, dye staining is a destructive method when used on soil
samples. A dye can sorb to macropore walls, resulting in an overestimation of nonequilibrium flow patterns. Thus, an advanced non-invasive and non-destructive X-ray computed tomography (CT) method
is often used when precise quantification of the soil pore system needs
to be achieved. It is a method superior to dyeing because of a highspatial resolution, which clearly exposes and differentiates soil pores,
roots, and solid-phase fractions. In soil science, the CT method is applied for monitoring of time-dependent soil physical properties such as
porosity, pore characteristics, density, meandering, and volumetricwater content, as well as solute transport parameters, fractal characteristics, soil aggregation, and retention and hydraulic conductivity
curves under unsaturated conditions (Cercioglu, 2018; Périard et al.,2016; Rogasik et al., 2003; Tracy et al., 2015). The method has considerable potential for the quantification of different processes within
the soil profile, such as root development, soil structure, water flow,
and solute transport (Koestel, 2017). However, the implementation of
X-ray computed tomography alone is often not sufficient to study nonlinear processes (i.e., flow and transport) in heterogenous arable soils.
Numerical models have been developed to describe such nonequilibrium flow processes in structured soils. Commonly used, physically-based numerical models for simulating water flow and solutetransport in soils use either single continuum or bi- and multi-continuum approaches (Kodešová et al., 2010). Bi-modal concepts fornonequilibrium flow assume that the porous soil system is divided into
two domains, in which one domain explains matrix (micropore) flow
and the other fracture (macropore) flow, with two sets of transport
properties (Gerke and van Genuchten, 1993). Haghverdi et al. (2020)studied the performance of five unimodal water retention models andeleven combinations of alternative (e.g., Peters–Durner–Iden PDI) andbi-modal expressions on 94 soil samples with high-resolution data setsfrom Turkey and the United States. Overall, Haghverdi et al. (2020)
showed that alternative expressions provided a better fit than unimodal
expressions. To correctly characterize the soil and to parameterize more
complex models, extensive and reliable data sets need to be analyzed byvarious techniques. Guo et al. (2019) combined time-lapse ground-penetrating radar, controlled infiltration, and high-frequency moisturemonitoring with a 2D/3D numerical modeling to investigate preferential flow in saprock. As concluded by Guo et al. (2020), the pairingof a useful combination of geophysical methods is needed to advanceour understanding of the complex interactions between plant (e.g., treeroots) and soil. Guo and Lin (2018) suggested that the magnitude andspatial pattern of preferential flow are regulated by interactions among
soil architecture, soil properties, landscape setting, and land use/landcover, which change with seasonal soil wetness conditions and water
input characteristics.
The overall goal of this study was to assess the soil porous systemand preferential flow pathways in vineyard soil by combining laboratory and numerical methods with dye staining and X-ray imaging. It
was hypothesized that the integration of various methods could reveal
more information regarding soil structure and flow behavior in structured arable soil. The specific objectives were i) to identify soil hydraulic parameters for two conceptual models (i.e., single-porosity anddual-permeability models) by combining laboratory (evaporation andleaching experiments) and numerical approaches, ii) to visualize dye
staining (Brilliant blue, BB) patterns within the vertical soil column
sections, (iii) to estimate and visualize pore architecture using X-ray
tomography, iv) to compare the BB stained soil images and X-ray
images to distinguish flow in macropores and matrix pores, and v) toapproximate dye transport using the two conceptual models (singleporosity and dual-permeability) and thus to identify the more suitablemodel.
2. Materials and methods
2.1. Study site and soil column preparation
The soil column study was performed in the laboratory of the
Department of Soil Amelioration (MELILAB), University of Zagreb
2Geoderma 380 (2020) 114699

V. Filipović, et al.

and are shown in the Results section (Table 2).
SHPs were described using the van Genuchten-Mualem model(VGM, van Genuchten, 1980):

Table 1
Average physicochemical properties of the vineyard Calcaric Leptosol soil
sampled at multiple depths at the Baštica experimental station located near
Zadar, Croatia.
Depth [cm]

0–25
25–50
50–75
75–125

Organic matter [%]

1.14
1.00
0.76
0.69

Soil texture [%]
Sand

Silt

Clay

pHH2O

EC* [dS m−1]

71
70
65
63

19
17
19
18

10
13
16
20

6.82
6.80
7.20
7.22

0.023
0.023
0.022
0.027

s

(h ) =r +(h ) =s forhSe =Faculty of Agriculture (UNIZG-AGR), using soil samples from theExperimental Station Baštica (44°09′19″N, 15°26′34″E), located near
Zadar, Croatia, and managed by the Department of Viticulture andEnology, UNIZG-AGR. The vineyard site occupies 6.5 ha and has a longterm average annual temperature of 15.5 °C and rainfall of 879 mm.
The vineyard was planted in 2008 with native Dalmatian varieties of V.
vinifera L. cv. Pošip and Maraština. For the determination of soil textural composition and physicochemical properties, disturbed soil sampling was performed (at 0–25, 25–50, 50–75, and 75–125 cm soil
depths) in triplicates (Table 1). Undisturbed soil columns were taken
within vineyard plant lines by gently pushing the PVC tubes (n = 7,
40 cm long with a 16 cm diameter) into the soil to minimize the disturbance of micro and macroporosity. Before the sampling, in order toprevent the edge effect of water leaching alongside column walls, a
non-reactive glue was applied to the column walls (e.g., Kodešová et al.,2010). One column was used for X-ray scanning, while the remaining
six columns were used for leaching experiments, three of which were
subsequently used for dye staining. The soil pH was measured using theMettler Toledo MPC 227 conductivity/pH-meter in water (pH H2O)
(HRN ISO 10390:2005). The organic matter content was determined bysulfochromic oxidation (HRN ISO 14235:1998). The soil particle size
distribution (fractions of sand, silt, and clay) was determined using thepipette method (Gee and Or, 2002). The soil was classified according to
the IUSS (2014) as Calcaric Leptosol. Soil pH, texture, bulk density, andorganic matter content are presented in Table 1.

for h < 0

(1)0

K (h) = Ks Sel (1* EC electrical conductivity.

r(1 + | h|n )m1

Sem ) m)2

(1

(2)

r
s

m=1

r

(3)

1
;n > 1n(4)where θr and θs denote residual and saturated volumetric watercontents [L3 L–3], respectively, h is the pressure head [L], Se is the effective saturation [–], α [L–1] and n [–] are shape parameters, and l [–]is a pore connectivity parameter. Pore connectivity parameter, l, is fixedin this study to a value of 0.5 (Mualem, 1976).
Undisturbed soil columns were saturated and left to drain in orderto reach the steady-state (matching) conditions before the intermittent
leaching experiment (n = 6). The bottom of the soil columns wascovered with a fine plastic mesh and tightly taped to prevent any soilmaterial disturbance. The experiment was conducted in seven cycles
carried out over 23 days. Approximately 350 mL of irrigation water wasapplied during each application (divided to 50 mL doses that wereapplied over a few hours to prevent excessive waterlogging) at the topof the columns in the form of mist using a hand-help sprayer to ensure
uniform water distribution. Water samples (leachate) were collected at
the bottom of the column after each application and volumes measured.
2.3. Numerical modeling
Numerical modeling was performed using the HYDRUS-1D program(version 4.0) (Šimůnek et al., 2016) that uses the Galerkin-type linear
finite element scheme to solve the governing flow equations for singleporosity (SP) and dual-permeability (DP) models. Water flow simulations (of intermittent leaching experiments) were based on a numericalsolution of the Richards equation for variably-saturated water flow in asingle-porosity (SP) media representing soil columns:

2.2. Soil hydraulic parameters and water flow estimation
Undisturbed soil cores (250 cm3, n = 5) were taken at multiple
depths of the topsoil horizon, up to 40 cm (5–10 cm, 10–15 cm,
20–25 cm, 30–35 cm, and 35–40 cm). The saturated hydraulic conductivity Ks [LT−1], was determined using the KSAT instrument(METER Group, USA) on these soil cores. The working principle of the
KSAT instrument is based on Darcy’s equation. After Ks estimation, theremaining soil hydraulic parameters (SHP) were estimated on the same
undisturbed soil cores (250 cm3) using the simplified evaporation
method (Schindler and Müller, 2017) and the HYPROP automatizedsystem applicable to the most soil types (Haghverdi et al., 2018). Thismethod takes into account the change of the sample weight and amatrix potential in the soil sample during the evaporation drying process (METER Group, USA). Soil hydraulic properties (SHPs) for Calcaric
Leptosol soil columns were determined using the HYPROP-FIT program

(h )
h
=
K (h )
+1
t
z
z

(5)
3 -3

where θ is the volumetric water content [L L ], h is the waterpressure head [L], K is the hydraulic conductivity [LT−1], t is time [T],and z is the vertical coordinate (positive upward) [L].
The advection–dispersion equation was used to simulate the solutetransport (Brilliant blue dye):

c
=
t
z

D

c
z

qc

(6)
-3

where c is the solute concentration in the liquid phase [M L ], q
represents the volumetric water flux density [L T−1] and D is the dispersion coefficient [L2T−1] defined as follows:

Table 2
VGM soil hydraulic parameters for Calcaric Leptosol, including porosity, determined using the HYPROP-FIT program, together with the goodness of fit (RMSE).
Sample No.

Depth [cm]Porosity [%]θr[cm3 cm−3]

θs[cm3 cm−3]α[cm−1]n [-]Ks[cm d−1]

RMSE_TH [cm3 cm−3]

RMSE_K [cm d−1]

1
2
3
4
5

5–10
10–15
20–25
30–35
35–40

43
42
40
40
37

0.031
0.087
0.070
0.088
0.057

0.395
0.371
0.361
0.380
0.358

0.019
0.020
0.029
0.018
0.016

1.56
1.50
1.53
2.00
1.61

3.00
4.20
12.80
2.90
3.86

0.017
0.002
0.006
0.009
0.01

0.82
0.16
0.25
0.59
0.39

3

Geoderma 380 (2020) 114699

V. Filipović, et al.

(7)

D = DL |q| + Dm

evaporation experiment. After a single-porosity model run, the dualpermeability model was applied to capture preferential flow leaching.
The fraction of the macropore domain (fw = 0.07) was determinedbased on the average pore distribution (macropores) derived from theCT images of the scanned soil column (see the Results and Discussion
Section below). Parameters characterizing the aggregate geometry wereset following the study of Kodešová et al. (2010). The following parameters were used to define the structure of both domains: β = 3 (theshape factor characterizing cubic aggregates), a = 2 cm (the characteristic length of an aggregate defining elongated pathways), and
Ksa = 1.0 × 10-6 cm h−1. The effective saturated hydraulic conductivity of the macropore-matrix interface was set relatively low to
avoid any numerical instabilities (Kodešová et al., 2010). The retention
curve parameters (θr, θs, α, and n) were used to characterize matrix(subscript m) and macropore (subscript f) domains where saturated
volumetric water contents were kept the same for both domains while
residual water contents were set to 0 (Gerke and Kohne, 2004;
Kodešová et al., 2008, 2010). The parameters αf and nf were increased
to 0.1 cm−1 and 3 for the entire column domain, respectively. It has
been found previously that such values increase macropore flow
(Kodešová et al., 2010). While the saturated hydraulic conductivity ofthe matrix domain, Ksm, was determined using the evaporation experiment (Table 2), the saturated hydraulic conductivity of the macropore domain, Ksf, was increased by 100 to account for macropore
flow. A related Ksf to Ksm ratio was estimated in Kodešová et al. (2008).
Solute mass transfer coefficient, s , was set to 0.02 d-1 (Gerke et al.,2013). It should be noted here that the objective of this work was not an
in-depth analysis of the dual-permeability model parametrization.
Initial conditions for modeling were set to the pressure head of 0 at
the bottom and −40 cm at the top (hydrostatic equilibrium) of the soilcolumn. A time-variable atmospheric flux boundary condition (applied
irrigation and evaporation) was used at the top, and a seepage faceboundary condition was applied at the bottom of the column (replicating experimental procedure in Section 2.2). Dye simulation followed the same initial and boundary condition with additionally adding
solute transport and adapting input dosage and amount as explained in
Section 2.4.
1 fw The performance of the calibrated water flow model was
evaluated using the root mean square error (RMSE, for evaporation andleaching experiments) and the coefficient of determination (R2, forleaching experiments). The model error for water retention (RMSE_TH)
was calculated separately from the model error for hydraulic conductivity (RMSE_K).

2 −1

where Dm is the coefficient of molecular diffusion [ML T ], τ is the
tortuosity factor [–], and DL represents the longitudinal dispersivity [L].
Additional modeling of intermittent leaching experiments was also
performed using a dual-permeability model, in which the Richards
equation is applied separately to two overlapping pore regions, i.e., themacropore and matrix domains (Gerke and van Genuchten, 1993). A setof two coupled Richards equations for flow in the fracture and matrix
domains, respectively, is defined as follows:f (h f )=tm (h m )t

K f (h f )

z

=

z

K m (h m )

hf
z

w

+1

(8)

fw

hm
+1
z

+

w

1

(9)

fw

where the subscripts f and m refer to fracture (macropore) andmatrix (micropore) domains, respectively, w is the water transfer term
for water exchange between the macropore and matrix domains [T−1],
and fw is the macropore domain fraction [-]. The matrix domain fraction is specified as 1 fw . The total composite soil water content andhydraulic conductivity are defined as the sums of the soil water contents and hydraulic conductivities of each domain multiplied by corresponding domain fractions.
The mass exchange between the matrix and macropore regions, w ,is expressed as:
w =

a2 w

K a (h a )

w

(hf

hm )

(10)

in which w (=0.4) is a dimensionless scaling factor (Gerke and van
Genuchten, 1996), and Ka is the effective hydraulic conductivity of theinterface between the two pore domains [L T−1], defined as:(11)

K a (ha) = Ksa Krm (ha )where Ksa is the saturated hydraulic conductivity of the interface
between the macropore and soil matrix domains, Krm(ha) is the relative
unsaturated hydraulic conductivity of the matrix domain calculatedusing the van Genuchten model for the pressure head of the interface ha[taken as (hf + hm)/2]. The dimensionless shape factor β (15 forspherical aggregates, 3 for cubic aggregates), and the characteristic
length of an aggregate a [L] (a sphere radius or a half size of the cubeedge) are the parameters describing the aggregate shapes (Gerke andvan Genuchten, 1996).
The dual-permeability formulation for solute transport is based on
two coupled advection–dispersion equations (Gerke and van
Genuchten, 1993):
f cf

=

t
m cm

t

=

f Df

z

z

cf

f cf

s

z

z

fw

m Dm

cm
z

m cm

z

+

(12)

1

(13)

fw
−1

where s is the solute mass transfer term [M L T
Genuchten, 1996):s (1fw ) m (cfcm ) +

R2 =

s

-3

s =

RMSE =

(14)

a2

Da

N
(Oi
i=1

O )(Pi0.5O )2] [

P)

N
(Pi
i=1

0.5

P )2]

(17)

2.4. Dye staining experiments and image analysis

where c* is either cf or cm depending on the exchange direction, and−1] of thes is the first-order solute mass transfer rate coefficient [T
form:
s =

(16)

N
N
(Oi
i=1

[

Pi ) 2where Oi are observation data points, Pi are model predictions, Ō isan average observation, P is an average prediction, and N is the sample
size.

] (Gerke and van

wc

N
(Oi
i=1

The dye experiments with soil columns were performed undercontrolled conditions in MELILAB at the UNIZG-AGR after the intermittent leaching experiment (Fig. 1). Brilliant blue dye (Brilliant Blue
FCF) in the form of solute (5 g L-1) was added to the soil columns(n = 3). The dye solution was dissolved in 500 mL of distilled waterand applied (divided into 250 mL dosages) with a handheld sprayer sothat approximately 1 cm of ponding occurred at the column top. After
24 h, columns were cut along the longer axis (vertically) in half. Fromeach column side, 2 cm slices (40 cm long) were prepared, resulting

(15)
2 −1

in which Da is an effective diffusion coefficient [L T ], which represents the diffusion properties of the fracture–matrix interface (e.g.,
accounting for organic coatings).
The parameters for the single-porosity model were derived from the4Geoderma 380 (2020) 114699

V. Filipović, et al.

Fig. 1. a) In-field sampling of undisturbed
soil columns at the Baštica experimental
vineyard, (b) the preparation for the dye
tracer experiment, and (c) the application of
the Brilliant Blue FCF dye at the top of thesoil column. (For interpretation of the references to color in this figure legend, thereader is referred to the web version of thisarticle.)altogether to eight slices per column. Soil slices were then photographed using a high-resolution camera (Samsung NX1000, 20.3 MP;
20–50 mm f/3.5–5.6 ED II Lens). Image analysis was performed using
the ImageJ software, commonly used for image processing and analysis(Ferreira and Rasband, 2012). The images were “thresholded” after
correction to create binary images with stained and unstained areas.
Dye coverage was calculated from the percentage of stained pixels vs.
depth, according to Flury et al. (1994) and Sander et al. (2007). ImageJ
was used for the determination of soil stained by Brilliant blue dye andfor the quantification of the pore volume and distribution inside the soilcolumn. In these analyses, only the upper 30 cm of soil samples weretaken into account to avoid disturbed soil at the bottom of the column.
Also, the applied Brilliant Blue FCF dye colored the soil down to the
30 cm depth.

estimation of the pore origin, e.g., bio pores or cracks.
3. Results and discussion
3.1. Soil hydraulic parameters and water flow
The SHPs of the VGM model are shown in Table 2. Errors in thewater retention and hydraulic conductivity model were below 0.02 cm3
cm−3 and 0.6 cm day−1, respectively, indicating that the unimodal van
Genuchten (VGM) model was appropriate for representing the soil hydraulic properties during the evaporation experiment. Table 2 showsthat the saturated hydraulic conductivity decreased in the top andbottom of the soil column, where values were 3.0 and 3.86 cm day−1,
respectively. The highest hydraulic conductivity was at a depth of
20–25 cm, where it was 12.8 cm day−1. The soil surface appears to beaffected by heavy trafficking in the vineyard, and the absence of roots
and macrofauna. The middle layer had higher porosity, resulting inhigher hydraulic conductivity. It is interesting that the bi-modal character of the capillary pore size distribution, usually present in structured soils (Kodešová et al., 2009, 2012), was not detected here using
the classical approach (Fig. 2). Also, one layer (i.e., the first 5–10 cm
layer) did not follow the VGM fit well. Indeed, when the bi-modal VGM
model was used in this soil layer, the improved fit was observed(RMSE_TH = 0.0035 cm3 cm−3 and RMSE_K = 0.3143 cm day−1).
However, such improvements were not found in other soil layers, andthus the single-porosity VGM model was applied.
As reported in Bezerra-Coelho et al. (2018), most soil types ranging
from fine to very coarse textures can be described using the singleporosity VGM model. The dual-porosity model (Durner 1994) is applicable only in cases when a good fit using a single-porosity VGM
model cannot be obtained. When a dual-porosity model of Durner(Durner 1994) is used, an apparent double S-shape curve in the soil
water retention curve is expected to fit the experimental data. It should
be noted here that the bi-modal character of the pore structure can be
hidden (Kodešová et al., 2009, 2008) even when the pore system isindeed bi-modal and needs to be divided into two regions (Gerke andKohne, 2004; Gerke and van Genuchten, 1996; Kodešová et al., 2010;
Vogel et al., 2000). In the Introduction, we cite Haghverdi et al. (2020),
who reported for a large dataset that the bi-modal and alternative expressions (e.g., including Peters–Durner–Iden PDI) could provide a
better fit than unimodal models. Still, the differences are sometimes
almost negligible (e.g., VG: RMSE = 0.015, r (correlation) = 0.995 vs.
VG-PDI RMSE = 0.009 r = 0.998; Haghverdi et al., 2020). Thus, insuch situations, both hydraulic property models would provide a very
similar description of flow behavior. This issue is further explored
below.
Additionally, the calculated porosity on the column scale was 44.2%(n = 6), showing comparable results to those obtained using HYPROP(40.4%, n = 5). A decrease in soil porosity and an increase in the clay
content of the last observed soil column part (Table 2) were in line withdecreasing hydraulic conductivity.

2.5. X-ray imaging of soil columns
The selected soil column was scanned using X-ray computed tomography (CT) Siemens Somatom Definition AS to visualize andquantify internal soil porosity. The column was scanned using a medical-grade scanner at a voltage of 120 kV and a current of 350 mAs.
Resolution was 0.48 mm × 0.48 mm × 0.48 mm. Only the upper 600
radiographs were selected for the analysis as the soil column showedsigns of disturbance in the bottom part (likely occurring during transport and handling). Pictures recorded for the bottom of the soil column
were therefore discarded as the soil started to dissipate (due to thehigher sand content in that layer, Table 1). By conducting the image
analysis, the pore volume and distribution were determined, as well as
the 3D representation of the porous soil system. The radiographs were
subsequently inverted into 3D images and exported as TIFF-stacks (atagged image file format). Image processing of the X-ray and dye experiments was accomplished using the open-source software ImageJand the bundle of plugins distributed by FIJI (Schindelin et al., 2012).
Since all X-ray parameters and settings were kept the same for all
images, image reconstruction of all layers could be accomplished on asemi-automated basis. Fixed threshold values, to distinguish between
pores and the solid phase, were selected based on the method presented
in the SoilJ plugin (Koestel, 2017). Two gray reference values wereselected that correspond to objects of known or at least constant density. Typically, one of the reference values corresponds to the grayvalue of the column walls. The second gray value was chosen as aquantile of the histogram of gray-values within individual horizontal
cross-sections. A very low quantile, e.g., 0.001, was chosen to represent
the least dense phase in the imaged object (i.e., air). X-ray scanning andimage analyses were performed using the methodology presented by
Larsbo et al. (2014). Image processing (thresholding and segmentation)
and analysis were carried out using the Fiji distribution (Schindelin
et al., 2012) of the ImageJ software (Abramoff et al., 2004), including
the BoneJ (Doube et al., 2010) and SoilJ (Koestel, 2017) plugins. The
main output of X-ray CT scanning was the visualization of the 3D poresystem architecture, the quantification of the pore system, and the5Geoderma 380 (2020) 114699

V. Filipović, et al.

Fig. 2. Soil water retention (a) and hydraulic conductivity curves (b, c) determined using the HYPROP-Fit software and fitted using the VGM model (for topsoil
samples taken at different depths down to the 40-cm depth (n = 5)).

Fig. 3. a) Measured and simulated (HYDRUS-1D, using the single-porosity (SP) and dual-permeability (DP) models) soil column cumulative water outflow (n = 6;
C1–C6); b) Simulated water outflow rates obtained using the SP and DP models.

Although the single-porosity VGM model provided an excellent
overall description of water flow (Fig. 3a, R2 of 0.984), the dye staining
and X-ray imaging presented evidence of extended macropore network(see Chapters 3.2.1 and 3.2.2), which should result in preferential flow
for conditions close to saturation. Thus, as reported previously(Kodešová et al., 2006; Kodešová et al., 2008; Filipović et al., 2019), a
different modeling approach, i.e., a dual-permeability model, can differentiate between flows in soil macropore and matrix regions. Additional dual-permeability modeling was thus performed using HYDRUS1D. Similar results for the intermittent leaching experiments (Fig. 3)
were obtained using the SP and DP models in terms of goodness of fit(R2 SP = 0.984; R2 DP = 0.987). Both modeling approaches, SP andDP, showed a similar water flow behavior and resulted in a good fit
when compared to the experimental data involving cumulative outflow.
A comparable observation was found by Kodešová et al. (2008), whoused the SP and DP models for multistep outflow data collected on

Greyic Phaeozem soil. In our leaching experiments, the use of the DPmodel, compared to the SP model, resulted in a faster response ofoutflow to irrigation and higher outflow rate peaks. Unfortunately, theexperimental data collected during the performed leaching experiments
did not provide sufficient evidence to distinguish between the twomodels and between equilibrium and nonequilibrium flow patterns.
While one would not expect to obtain such evidence from evaporation
experiments, which is run under unsaturated conditions, such evidence
could have been obtained from leaching experiments, had more detailed outflow dynamics been observed (more, rather than only two,
data points for each leaching cycle), or had environmental tracers beenused in infiltrating water (observing outflow breakthrough curves witharrival times, peaks, and possible tailing).
It is also well known that the soil hydraulic properties estimatedusing different methods can be quite different even in the case of sandy
soil, known for its uniformity (e.g., Gribb et al., 2004). Thus, indicating
6

Geoderma 380 (2020) 114699

V. Filipović, et al.

the necessity of implementing various methods and experiments duringtransient flow and transport studies in arable soils. Besides the effect ofmultimodality of the soil porous systems on soil water regimes, anotherreason may be the hysteresis in soil hydraulic functions. Other factors
affecting resulting soil hydraulic properties may include the characterof the experiment (i.e., steady-state or transient flow experiments), thesize of the soil sample, dimensions of the flow domain (1D, 2D, or 3D),and methods used for evaluation of soil hydraulic properties (analyticalequations or inverse modeling). Kutílek et al. (2009) tested empirical,
semi-empirical, and physically-based models of soil hydraulic functions
derived for bi-modal soils aimed at capturing preferential flow. They
concluded that, although many empirical models exist that could be
modified for multimodal functions, simple models often provide a good
agreement between modeled and measured flow data. This implies thatnumerical models simulating water flow using unimodal models canstill provide a good agreement with experimental data on water flow
collected in structured soils. However, such models then often fail toprovide a good description of the solute transport behavior.

Cambisol.
In our study, the plow pan effect was observed below the depth of
30 cm, which is probably due to the vineyard tillage management,
which is usually limited to the approximately upper 20–30 cm soil
depth between the rows (e.g., Comino et al., 2017). This is illustrated indye staining images and supported by lower porosities in deeper soil
layers. Sander and Gerke (2007) reported that dye penetration indrained paddy soils was directed horizontally when plow pan wasreached, while further penetration to deeper layers occurred throughbiopores present in the plow pan.
By additionally estimating dye staining soil patterns, information
regarding macropores origin could help in clarifying the preferential
flow behavior. This is further explored in the next section, where thedye staining and X-ray imaging results are compared. The aggregate
coatings also play an important role in dye transport and visualization.
Pore geometry in the soil changes with depth, thus influencing the
spatial orientation of the preferential flow (Guo and Lin, 2018).
Kodešová et al. (2015) found that dye distribution improved whenstraw residues, earthworms, and roots were present in the soil, but alsothat topology, macropores, and water content have an impact on dye
distribution as well. Presented dye staining images indicate that themain part of the macropore network is formed by root channels andearthworm burrows. Although macropore flow can be determined using
this method, limitations exist due to potential dye sorption onto soil
particles. Nevertheless, preferential flow and infiltration patterns in thevadose zone can be reliably determined by the application of Brilliant
Blue dye (Wu et al., 2015).

3.2. Identification of flow in preferential pathways and matrix pores
3.2.1. Soil column dye staining
Dye staining patterns in vertical cross-sections of the soil columns
indicated the highly heterogeneous nature of the flow patterns in thisarable soil. Percentages of stained areas for column 1 are in the range
from 17.2% to 33.4%, for column 2 from 18.5% to 31.9%, and forcolumn 3 from 18.4% to 25.3% (Table 3). Figs. 4 to 6 illustrate theuneven dye distribution caused by the enhanced dye penetration along
biopores and in regions with a higher fraction of larger interaggregate
pores formed in the root zone. A significant preferential flow could beobserved below the soil depth of approximately 10 cm, except in Fig. 4,
sections 1 − 4b, where preferential flow started to occur below 5 cm.
This was also confirmed when analyzing SHPs since the middle layer(20–25 cm) had higher hydraulic conductivity (Table 2). Similar observations were reported on a field scale by Kodešová et al. (2012).
They found the occurrence of extensive preferential flow in the subsurface horizons of arable soils, i.e., Haplic Luvisol and Haplic

3.2.2. X-ray Imaging, pore architecture, and bimodal system
conceptualization
From 600 images produced using the X-ray CT scanning and processed using the Image J software, image No. 8 had the largest fractionof the soil matrix within the sample column. In contrast, images No.
139 and 140 had the smallest fraction (Fig. 7). There were twenty
column sections with average values of 6.6% of macropores. Despite
equal percentages of the soil matrix in images No. 87 and 568, it isvisible that the pore arrangement in these column sections is entirely
different. The difference in the pore arrangement is caused by various
factors, mostly by the presence of flora and fauna in the soil (e.g.,
Menta, 2012). Images 1, 200, 400, and 600 (Fig. 8a) are chosen to illustrate changes in the appearance and distribution of pores at different
soil depths. The fractions of the macropores within the selected column
are in the range from 2.3 to 1%. Additionally, an unprocessed image
from X-ray scanning is shown in Fig. 8b to illustrate the scale of themacropore network. The largest pores had diameters from 5.2 up to
10.6 mm, suggesting that indeed the pore origin is from earthworm
burrows and roots. The preferential flow through the large diameter
macropores is well known to occur in arable soils. These macropores
are usually described as biopores (e.g., Beven and Germann, 2013; Guoand Lin, 2018; Sander and Gerke, 2007).
If we compare the results obtained by dye staining and CT, i.e., theaverage percentages of the pore systems identified by the two methods(23% compared to 7%, respectively), there is a large difference of 16%.
This difference can be related to dye coloring also the soil aroundmacropores, i.e., it can be related to dye penetration into the pores
smaller than the CT threshold. Again, this illustrates the difficulty of
selecting an appropriate method to quantify preferential flow in
structured soils. Soto-Gómez et al. (2019) presented a novel method of
combining fluorescence macrophotography and X-ray computed tomography (CT) in preferential flow estimation. They concluded that
preferential flow occurred in only a small fraction of the total pore
network and was controlled by pores connected to the soil surface andby matrix density. This can be further illustrated by overlapping the dye
staining and X-ray imaging. Fig. 9a shows the macroporosity detectedusing the X-ray scanning and the percentage of soil stained with Brilliant Blue FCF dye. The latter shows that dye was transported from theTable 3
Dye stained percentages of colored areas of the soil columns calculated afterimage processing using the Image J software.
Column No.

Image No.

Stained area [%]Column 1

1a
2a
3a
4a
1b
2b
3b
4b
1a
2a
3a
4a
1b
2b
3b
4b
1a
2a
3a
4a
1b
2b
3b
4b

26.7
31.1
33.4
17.2
23.4
26.4
19.4
18.2
23.9
23.2
31.9
26.3
24.3
30.8
29.2
18.5
21.8
21.2
20.4
21.0
20.3
25.3
18.4
21.3
23.9
4.62
0.19

Column 2

Column 3

AverageSD
CV

7

Geoderma 380 (2020) 114699

V. Filipović, et al.

Fig. 4. Images of soil after the application of the brilliant blue dye (Brilliant Blue FCF) (left image with brown background) and the ImageJ processing (right image)for soil Column 1 cut in half (marked as a and b) and then sliced to 4 column slices (1–4). (For interpretation of the references to color in this figure legend, the reader
is referred to the web version of this article.)column surface down to a depth of about 8–10 cm in the soil matrix
pores, while the deeper presence of dye can be explained only by thetransport in macropores. Numerical modeling of the Brilliant Bluetransport (Fig. 9b) with HYDRUS-1D indicate that a single-porosity
model cannot explain dye coloring in larger depths (i.e., below 15 cm)
while the dual-permeability model can capture the preferential flow inthe macropores as well (i.e., in entire soil column). It should be noted
here that Fig. 9a represents dye coloring which is considerable lighter atthe 10–20 cm depth (Figs. 4 to 6). Contrary Fig. 9b represents Brilliant
Blue mass which should explain the discrepancy between the two.

In some instances, according to Nielsen (2004), porosity determinedby X-ray CT differs from values obtained in the laboratory. However, Xray CT scans can be used to identify the largest pores present in thesystem and for macropore detection. By decreasing the image quality,
the accuracy of the porosity determination is reduced as well. This iswell known, especially when relatively large soil samples are used, as
was done in our study. A similar conclusion was made by Rab et al.
(2014). They concluded that sample scanning at a higher resolution
would result in more macropores detected, compared with sample
scanning at a lower resolution. One of the possible error origins is theFig. 5. Images of soil after the application of the brilliant blue dye (Brilliant Blue FCF) (left image with brown background) and the ImageJ processing (right image)for soil Column 2 cut in half (marked as a and b) and then sliced to 4 column slices (1–4). (For interpretation of the references to color in this figure legend, the reader
is referred to the web version of this article.)8Geoderma 380 (2020) 114699

V. Filipović, et al.

Fig. 6. Images of soil after the application of the brilliant blue dye (Brilliant Blue FCF) (left image with brown background) and the ImageJ processing (right image)for soil Column 3 cut in half (marked as a and b) and then sliced to 4 column slices (1–4). (For interpretation of the references to color in this figure legend, the reader
is referred to the web version of this article.)determination of thresholds between the solid soil phase and pores
during image processing, due to the subjective perception involved inthis step. However, in our study, the aim was to capture the pore architecture on a column scale, thus accepting the lower resolution of X-

ray images. Fig. 10 displays a three-dimensional (3D) pore system architecture obtained using X-ray on a column taken from the Calcaric
Leptosol soil. From the 3D visualization, it is clear that biopores originating from earthworm burrows or roots represent most of the poreFig. 7. Column X-ray images with the maximum (image No. 8 at 3.84 cm depth), minimum (image No. 139 at 6.67 cm depth), and average (images No. 87 at 4.18 cm
depth and 568 at 27.26 cm depth) solid phase (soil matrix) percentages. The soil pore system is represented in white color.
9Geoderma 380 (2020) 114699

V. Filipović, et al.

Fig. 8. a) Variations in the soil pore system (in white color) obtained using X-ray scanning at various depths of the soil column. (b) An unprocessed X-ray image withindicated macropore diameters.

system.
Leue et al. (2019) reported that the quantitative techniques, such as,e.g., X-ray CT combined with subsequent (3D) image analysis, areneeded to differentiate between biopores and cracks. Here, a more
holistic approach was used, in which a combination of methods resulted
not only in an overall description of the pore system but also of the flow
behavior. A relatively large column size (40 cm long × 16 cm in diameter) did not allow for much better image resolution using the medical scanner. However, by combining all presented methods, the preferential flow was quantified, and its origin determined. Gao and Lin(2018) reported that 25% of 190 analyzed case studies on preferential
flow combined multiple methods for its quantification. As also provenin this study, the integration of different methods overcomes some of

the limitations of each technique and helps to constrain the detected
preferential flow patterns and reveals some hidden information.
4. Conclusions
This study combined laboratory flow experiments and numericalmethods with dye staining and X-ray imaging to assess pore systemarchitecture and preferential flow in the vineyard soil. Soil hydraulicproperties measured in the laboratory and also determined using numerical modeling did not provide enough evidence supporting thenonequilibrium flow patterns. The numerical simulations performed
well when the single-porosity model was used to simulate leachingexperiments. However, because the extensive macropore network wasFig. 9. a) Percentages of flow pathways: Xray – potential preferential flow in macropores detected using the X-ray scanning, BB– flow in the macropores and soil matrix
derived from images of soil stained with
Brilliant Blue FCF dye. b) Brilliant Blue FCF
mass simulated using HYDRUS-1D and thedual-porosity (DP_BB; indicating matrix andmacropore transport, i.e., DP_m_BB and
DP_f_BB) or single-porosity (SP_BB) models.
(For interpretation of the references to color
in this figure legend, the reader is referred to
the web version of this article.)10

Geoderma 380 (2020) 114699

V. Filipović, et al.

Fig. 10. Three-dimensional (3D) visualization of the soil pore system (28.8 cm depth) of the vineyard soil obtained using computerized X-ray tomography (CT) andreproduced in the ImageJ software.

revealed using dye staining and X-ray imaging, the dual-permeability
model was used as well. The dual-permeability model also produced ahigh level of agreement (R2 = 0.987), while additionally providing afaster response to applied irrigation events. However, collected experimental data did not allow us to differentiate the performance of thetwo models. The study showed that the bi-modal character of pore
structure could often be hidden when a limited number of methods or
non-adequate methods (e.g., by relying only on laboratory-derived
SHPs and numerical modeling) are applied for its quantification.
The dye patterns, which represented about 17 to 33% (average23%) of the area within the vertical column sections, indicated thepresence of preferential (macropore) flow in the arable Calcaric
Leptosol vineyard soil. Additionally, a combination of dye staining and
X-ray imaging revealed a complex pore architecture network. Using
only dye staining can also be misleading due to dye sorption propertiesor relatively large mass exchange in preferential flow quantification
studies. The preferential flow was mostly caused by biopores, i.e.,
earthworm burrows and root channels (7% of the soil column). The
macropores with a diameter of up to 10.6 mm were found during theimage analysis, indicating their origin as stated above. A 3D visualization showed complex pore structure in the arable vineyard soil, withlarge vertical and horizontal biopores, confirming the importance ofcombining various techniques for preferential flow determination. By
combining the dye staining and X-ray imaging, different flow patterns
were revealed in the soil matrix and macropores. Additional modelingof the transport of Brilliant Blue dye revealed the necessity of using
dual-region transport models to capture non-equilibrium flow patterns
in agricultural soils with complex pore systems.

Croatia, which provided know-how and the equipment for X-ray computed tomography.
References
Abramoff, M.D., Magalhaes, P.J., Ram, S.J., 2004. Image processing with ImageJ.
Biophotonics Int. 11, 36–42.
Bebej, J., Homolák, M., Orfánus, T., 2017. Interaction of Brilliant Blue dye solution withsoil and its effect on mobility of compounds around the zones of preferential flows at
spruce stand. Cent. Eur. For. J. 63 (2–3), 79–90.
Beven, K., Germann, P., 2013. Macropores and water flow in soils revisited. Water Resour.
Res. 49, 3071–3092.
Bezerra-Coelho, C.R., Zhuang, L., Barbosa, M.C., Soto, M.A., van Genuchten, M.T., 2018.
Further tests of the HYPROP evaporation method for estimating the unsaturated soil
hydraulic properties. J. Hydrol. Hydromech. 66 (2), 161–169.
Bogunović, I., Telak, L.J., Pereira, P., 2020. Agriculture management impacts on soil
properties and hydrological response in Istria (Croatia). Agronomy 10 (2), 282.
Cercioglu, M., 2018. Imaging soil pore characteristics using computed tomography as
influenced by agroecosystems. Eurasian J. Soil Sci. 7 (3), 195–202.
Cerdan, O., Govers, G., Le Bissonnais, Y., Van Oost, K., Poesen, J., Saby, N., Gobin, A.,
Vacca, A., Quinton, J., Auerwald, K., Klik, A., Kwaad, F.J.P.M., Raclot, D., Ionita, I.,
Rejman, J., Rousseva, S., Muxart, T., Roxo, M.J., Dostal, T., 2010. Rates and spatial
variations of soil erosion in Europe: A study based on erosion plot data.
Geomorphology 122, 167–177.
Comino, J.R., Bogunović, I., Mohajerani, H., Pereira, P., Cerdà, A., Ruiz Sinoga, J.D., Ries,
J.B., 2017. The impact of vineyard abandonment on soil properties and hydrologicalprocesses. Vadose Zone J. 16 (12).
Doube, M., Klosowski, M.M., Arganda-Carreras, I., Cordelieres, F.P., Dougherty, R.P.,
Jackson, J.S., et al., 2010. BoneJ: free and extensible bone image analysis in ImageJ.
Bone 47, 1076–1079.
Durner, W., 1994. Hydraulic conductivity estimation for soils with heterogeneous pore
structure. Water Resour. Res. 32 (9), 211–223.
Fér, M., Kodešová, R., Hroníková, S., Nikodem, A., 2019. The effect of 12-year ecologicalfarming on the soil hydraulic properties and repellency index. Biologia. https://doi.
org/10.2478/s11756-019-00373-1.
Fér, M., Kodešová, R., Nikodem, A., Jirků, V., Jakšík, O., Němeček, K., 2016. The impactof the permanent grass cover or conventional tillage on hydraulic properties of Haplic
Cambisol developed on paragneiss substrate. Biologia 71 (10). https://doi.org/10.
1515/biolog-2016-0133.
Ferreira, T., Rasband, W., 2012. ImageJ User Guide, IJ, 1.46r, revised edition,. pp. 198.
Filipović, V., Coquet, Y., Gerke, H.H., 2019. Representation of Plot-Scale Soil
Heterogeneity in Dual-Domain Effective Flow and Transport Models with Mass
Exchange. Vadose Zone J. 18 (1), 1–14.
Flury, M., Flühler, H., 1994. Brilliant blue FCF as a dye tracer for solute transport studies –A toxicological overview. J. Environ. Qual. 23, 1108–1112.
Flury, M., Flühler, H., 1995. Tracer characteristics of brilliant blue FCF. Soil Sci. Soc. Am.
J. 59, 22–27.
Flury, M., Flühler, H., Jury, W.A., Leuenberger, J., 1994. Susceptibility of soils to preferential flow of water: A field study. Water Resour. Res. 30, 1945–1954.
Fuentes, I., Casanova, M., Seguel, O., Nájera, F., Salazar, O., 2014. Morphophysicalpedotransfer functions for groundwater pollution by nitrate leaching in Central Chile.
Chil. J. Agr. Res. 74 (3), 340–348.
Galati, A., Gristina, L., Crescimanno, M., Barone, E., Novara, A., 2015. Towards more
efficient incentives for agri-environment measures in degraded and eroded vineyards.

Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
The research was funded by the Croatian Waters project: Spatial andtemporal distribution of contaminants from agriculture (nitrates,
phosphates, pesticides, heavy metals) under different agroecologicalconditions (Grant No. 10-059/17). We thank Professor Boris Brkljačić
M.D. and medical staff at the Clinical Hospital Dubrava, Zagreb,
11Geoderma 380 (2020) 114699

V. Filipović, et al.

Marcelino, V., Mees, F. (Eds.), interpretation of micromorphological features of soils
and regoliths (Second Edition). Elsevier, pp. 447–469.
Kutílek, M., Jendele, L., Krejča, M., 2009. Comparison of empirical, semi-empirical andphysically based models of soil hydraulic functions derived for bi-modal soils. J.
Contam. Hydrol. 104 (1–4), 84–89. https://doi.org/10.1016/j.jconhyd.2008.10.007.
Larsbo, M., Koestel, J., Jarvis, N., 2014. Relations between macropore network characteristics and the degree of preferential solute transport. Hydrol. Earth Syst. Sc. 18,
5255–5269.
Leue, M., Uteau-Puschmann, D., Peth, S., Nellesen, J., Kodešová, R., Gerke, H.H., 2019.
Separation of soil macropore types in three-dimensional x-ray computed tomography
images based on pore geometry characteristics. Vadose Zone Journal 18 (1), 1–13.
https://doi.org/10.2136/vzj2018.09.0170.
Menta, C., 2012. Soil Fauna Diversity – Function, soil degradation, biological indices, soil
restoration. In Lameed G.A. (Ed.): Biodiversity conservation and utilization in a diverse World, IntechOpen. doi: 10.5772/51091.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated
porous media. Water Resour. Res. 12 (3), 513–522.
Nielsen, B.D., 2004. Non-destructive Soil Testing Using X-Ray Computed Tomography.
Master thesis. Montana State University, Bozeman, Montana, pp. 137.
Pagliai, M., Vignozzi, N., Pellegrini, S., 2004. Soil structure and the effect of management
practices. Soil. Till. Res. 79 (2), 131–143.
Périard, Y., José Gumiere, S., Long, B., Rousseau, A.N., Caron, J., 2016. Use of X-ray CT
scan to characterize the evolution of the hydraulic properties of a soil under drainage
conditions. Geoderma 279, 22–30.
Prosdocimi, M., Cerdà, A., Tarolli, P., 2016a. Soil water erosion on Mediterranean vineyards: a review. Catena 141, 1–21.
Prosdocimi, M., Jordán, A., Tarolli, P., Keesstra, S., Novara, A., Cerdà, A., 2016b. The
immediate effectiveness of barley straw mulch in reducing soil erodibility and surface
runoff generation in Mediterranean vineyards. Sci. Total Environ. 547, 323–330.
Rab, M.A., Haling, R.E., Aarons, S.R., Hannah, M., Young, I.M., Gibson, D., 2014.
Evaluation of X-ray computed tomography for quantifying macroporosity of loamy
pasture soils. Geoderma 213, 460–470.
Ranjit Kumar, M., Meenambal, T., Kumar, V., 2017. Macropore flow as a groundwatercomponent in hydrologic simulation: modelling, applications and results. Curr.
Science 112 (6), 1197–1207.
Rogasik, H., Onasch, I., Brunotte, J., Jegou, D., Wendroth, O., 2003. Assessment of soil
structure using X-ray computed tomography. Geological Society, London, Special
Publications 215 (1), 151–165.
Salvador, M.M.S., Köhne, S., Köhne, J.M., Lennartz, B., Libardi, P.L., 2011. Dye tracer andmorphophysical properties to observe water flow in a Gleyic Luvisol. Sci. Agr. 68 (2),160–166.
Sander, T., Gerke, H.H., 2007. Preferential flow patterns in paddy fields using a dye
tracer. Vadose Zone J. 6, 105–115.
Schindelin, J., Arganda-Carreras, I., Frise, E., Kaynig, V., Longair, M., Pietzsch, T.,
Preibisch, S., Rueden, C., Saalfeld, S., Schmid, B., Tinevez, J.-Y., White, D.J.,
Hartenstein, V., Eliceiri, K., Tomancak, P., Cardona, A., 2012. Fiji: an open-source
platform for biological-image analysis. Nat. Methods 9, 676–682.
Schindler, U., Müller, L., 2017. Soil hydraulic functions of international soils measured
with the extended evaporation method (EEM) and the HYPROP device. Open Data
Journal for Agricultural Research 3, 10–16.
Schneider, A., Hirsch, F., Raab, A., Raab, T., 2018. Dye tracer visualization of infiltration
patterns in soils on relict charcoal hearths. Front. Environ. Sci. 6. https://doi.org/10.
3389/fenvs.2018.00143.
Soto-Gómez, D., Pérez-Rodríguez, P., Vázquez Juíz, L., López-Periago, J.E., Paradelo
Pérez, M., 2019. A new method to trace colloid transport pathways in macroporous
soils using X-ray computed tomography and fluorescence macrophotography. Eur. J.
Soil Sci. 70 (3), 431–442. https://doi.org/10.1111/ejss.12783.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2016. Recent developments and applications of the HYDRUS computer software packages. Vadose Zone J. 15 (8), 25.
Strudley, M.W., Green, T.R., Ascough, J.C., 2008. Tillage effects on soil hydraulic properties in space and time: state of the science. Soil Tillage Res. 99 (1), 4–48.
Stumpp, C., Maloszewski, P., 2010. Quantification of preferential flow and flow heterogeneities in an unsaturated soil planted with different crops using the environmental
isotope delta O-18. J. Hydrol. 394, 407–415.
Tracy, S.R., Daly, K.R., Sturrock, C.J., Crout, N.M.J., Mooney, S.J., Roose, T., 2015. Threedimensional quantification of soil hydraulic properties using X-ray computed tomography and image-based modeling. Water Resour. Res. 51, 1006–1022.
van Genuchten, M.T., 1980. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–1037.
Vogel, T., Gerke, H.H., Zhang, R., van Genuchten, M.T., 2000. Modeling flow andtransport in a two-dimensional dual-permeability system with spatially variable hydraulic properties. J. Hydrol. 238 (1–2), 78–89. https://doi.org/10.1016/s00221694(00)00327-9.
Wu, Q., Liu, C., Lin, W., Zhang, M., Wang, G., Zhang, F., 2015. Quantifying the preferential flow by dye tracer in the North China Plain. J. Earth Sci. 26, 435–444.

Land Degrad. Dev. 26 (6), 557–564.
Gee, G.W., Or, D., 2002. Particle size analysis. In: Dane, J.H., Topp, G.C. (Eds.): methodsof soil analysis. Part 4. Physical Methods, Chapter 3.3. Soil Sci. Soc. Am., Madison,
WI, pp. 255–293.
Gerke, H.H., 2006. Preferential-flow descriptions for structured soils. J. Plant Nutr. Soil
Sci. 169, 1–19.
Gerke, H.H., Dusek, J., Vogel, T., 2013. Solute mass transfer effects in two-dimensional
dual-permeability modeling of bromide leaching from a tile-drained field. Vadose
Zone J. 12 (2). https://doi.org/10.2136/vzj2012.0091.
Gerke, H.H., Kohne, J.M., 2004. Dual-permeability modeling of preferential bromide
leaching from a tile-drained glacial till agricultural field. J. Hydrol. 289, 239–257.
Gerke, H.H., van Genuchten, M.T., 1993. A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media. Water Resour.
Res. 29, 305–319.
Gerke, H.H., van Genuchten, M.T., 1996. Macroscopic representation of structural geometry for simulating water and solute movement in dual-porosity media. Adv. Water
Resour. 19, 343–357.
Germán-Heins, J., Flury, M., 2000. Sorption of Brilliant Blue FCF in soils as affected by pH
and ionic strength. Geoderma 97, 87–101.
Gribb, M.M., Kodešová, R., Ordway, S.E., 2004. comparison of soil hydraulic property
measurement methods. Journal of Geotechnical and Geoenvironmental Engineering
130 (10), 1084–1095. https://doi.org/10.1061/(asce)1090-0241(2004)
130:10(1084).
Guo, L., Lin, H., 2018. Addressing two bottlenecks to advance the understanding ofpreferential flow in soils. Adv. Agron. 147, 61–117.
Guo, L., Lin, H., Fan, B., Nyquist, J., Toran, L., Mount, G.J., 2019. Preferential flow
through shallow fractured bedrock and a 3D fill-and-spill model of hillslope subsurface hydrology. J. Hydrol. 576, 430–442. https://doi.org/10.1016/j.jhydrol.2019.
06.070.
Guo, L., Mount, G.J., Hudson, S., Lin, H., Levia, D., 2020. Pairing geophysical techniques
improves understanding of the near-surface Critical Zone: Visualization of preferential routing of stemflow along coarse roots. Geoderma 357, 113953. https://doi.
org/10.1016/j.geoderma.2019.113953.
Haghverdi, A., Öztürk, H.S., Durner, W., 2018. Measurement and estimation of the soil
water retention curve using evaporation method and pseudo continuous pedotransfer
function. J. Hydrol. 563, 251–259.
Haghverdi, A., Najarchi, M., Öztürk, H.S., Durner, W., 2020. Studying unimodal, bimodal,
pdi and bimodal-pdi variants of multiple soil water retention models: i. direct model
fit using the extended evaporation and dewpoint methods. Water 12, 900. https://doi.org/10.3390/w12030900.
Horn, R., 2004. Time dependence of soil mechanical properties and pore functions forarable soils. Soil Sci. Soc. Am. J. 68 (4), 1131–1137.
IUSS Working Group WRB. World Reference Base for Soil Resources 2014. Update 2015:
International soil classification system for naming soils and creating legends for soil
maps; World Soil Resources Reports No. 106; FAO: Rome, Italy, 2015; pp. 192.
Jarvis, N.J., 2007. A review of non-equilibrium water flow and solute transport in soil
macropores: principles, controlling factors and consequences for water quality. Eur.
J. Soil Sci. 58, 523–546.
Jirků, V., Kodešová, R., Nikodem, A., Mühlhanselová, M., Žigová, A., 2013. Temporalvariability of structure and hydraulic properties of topsoil of three soil types.
Geoderma 204–205, 43–58.
Kodešová, R., Jirků, V., Kodeš, V., Mühlhanselová, M., Nikodem, A., Žigová, A., 2011. Soil
structure and soil hydraulic properties of Haplic Luvisol used as arable land andgrassland. Soil Till. Res. 111 (2), 154–161.
Kodešová, R., Kodeš, V., Žigová, A., Šimůnek, J., 2006. Impact of plant roots and soilorganisms on soil micromorphology and hydraulic properties. Biologia 61 (19).
https://doi.org/10.2478/s11756-006-0185-7.
Kodešová, R., Kočárek, M., Kodeš, V., Šimůnek, J., Kozák, J., 2008. Impact of soil micromorphological features on water flow and herbicide transport in soils. Vadose
Zone J. 7 (2), 798–809.
Kodešová, R., Vignozzi, N., Rohošková, M., Hájková, T., Kočárek, M., Pagliai, M., Kozák,J., Šimůnek, J., 2009. Impact of varying soil structure on transport processes indifferent diagnostic horizons of three soil types. J. Contam. Hydrol. 104, 107–125.
Kodešová, R., Němeček, K., Kodeš, V., Žigová, A., 2012. Using dye tracer for visualizationof preferential flow at macro- and microscales. Vadose Zone J. 11 (1). https://doi.
org/10.2136/vzj2011.0088.
Kodešová, R., Nĕmeček, K., Žigová, A., Nikodem, A., Fér, M., 2015. Using dye tracer forvisualizing roots impact on soil structure and soil porous system. Biologia 70 (11),1439–1443.
Kodešová, R., Šimůnek, J., Nikodem, A., Jirků, V., 2010. Estimation of the dual-permeability model parameters using tension disk infiltrometer and guelph permeameter.
Vadose Zone J. 9 (2), 213–225.
Koestel, J., 2017. SoilJ: An ImageJ plugin for the semiautomatic processing of three-dimensional X-ray images of soils. Vadose Zone J. 17 (1), 1–7.
Kooistra, M.J., Pulleman, M.M., 2018. Features Related to Faunal Activity. In: Soops, G.,

12